Geometrical Tortousity

[aayushi363]

Hi all,
I am trying to find out Geometrical Tortousity of a2d Micro ct image of Carbonate. For that I am trying to find the shortest distance between the Surface pores on left side and right side (similarly for y direction) using find_path.

I have written the following code to find the left and right surface pores of network that was calculated using SNOW2. When I am trying to vizualize it in paraview, I am not able to see the pore on the coordinates that the program returns.

Am I missing something?

import numpy as np
import porespy as ps
import openpnm as op
from skimage import io
import matplotlib.pyplot as plt

ps.visualization.set_mpl_style()

image_path = "/Users/aayushigautam/Geometrical_Tortuosity/Sample/UofC0005.tif"
im = io.imread(image_path)

snow_output = ps.networks.snow2(im, voxel_size=1)
pn = op.io.network_from_porespy(snow_output.network)
print(pn)
path_to_file = "/Users/aayushigautam/Geometrical_Tortuosity/"
op.topotools.find_surface_pores(network=pn)  

Surface_Pores = pn.pores('surface')
print("The follwoing pores are surface pores: ", Surface_Pores)
Surface_Pores_coords = pn['pore.coords'][Surface_Pores]
print("The coordinates of the surface pores are: ", Surface_Pores_coords)
Surface_Pores_coords = np.array(Surface_Pores_coords)
pore_coords = pn['pore.coords']
max_x = np.max(pore_coords[:, 0])  # Max x-coordinate (column index 0)
max_y = np.max(pore_coords[:, 1])  # Max y-coordinate (column index 1)
min_x = np.min(pore_coords[:, 0])  # Min x-coordinate (column index 0)
min_y = np.min(pore_coords[:, 1])  # Min y-coordinate (column index 1)
tolerance = 1e-8

left_surface_pores = Surface_Pores_coords[np.isclose(Surface_Pores_coords[:, 0],min_x , atol=tolerance)]
pn.set_label(label='left', pores=left_surface_pores)
left_surface_pores = np.array(left_surface_pores)
print("The left surface pores are: ", left_surface_pores)

right_surface_pores = Surface_Pores_coords[np.isclose(Surface_Pores_coords[:, 0],max_x , atol=tolerance)]
pn.set_label(label='right', pores=right_surface_pores)
right_surface_pores = np.array(right_surface_pores)
print("The right surface pores are: ", right_surface_pores)
op.io._vtk.project_to_vtk(pn.project, filename=path_to_file+'/Paraview_network')

This is the paraview visualization:
The left surface pores are: [[ 1. 10. 0. ]
[ 1. 46. 0. ]
[ 1. 95. 0. ]
[ 1. 137. 0. ]
[ 1. 182. 0. ]
[ 1. 216. 0. ]
[ 1. 250.5 0. ]
[ 1. 275.5 0. ]
[ 1. 298. 0. ]
[ 1. 314. 0. ]
[ 1. 327.5 0. ]
[ 1. 337.5 0. ]
[ 1. 365. 0. ]
[ 1. 405.5 0. ]
[ 1. 424.5 0. ]
[ 1. 453.5 0. ]
[ 1. 475.5 0. ]
[ 1. 492.5 0. ]]
The right surface pores are: [[504. 4. 0. ]
[504. 18. 0. ]
[504. 40. 0. ]
[504. 75. 0. ]
[504. 115.5 0. ]
[504. 138. 0. ]
[504. 166.5 0. ]
[504. 191.5 0. ]
[504. 208. 0. ]
[504. 242. 0. ]
[504. 269. 0. ]
[504. 287. 0. ]
[504. 296.5 0. ]
[504. 320.5 0. ]
[504. 349.5 0. ]
[504. 389.5 0. ]
[504. 422.5 0. ]
[504. 451.5 0. ]
[504. 474.5 0. ]
[504. 485.5 0. ]
[504. 498. 0. ]]

[jgostick]

Geometric tortuosity is not a very robust metric and I don't recommend using it, but if you must...then you'll want to use Djikstra's Algorithm, which is available in scipy.sparse.csgraph. This tells you the distance between all pairs of pores, so you can ask what is the distance between a pore on the left edge and the nearest pore on the right edge. You can tell the algorithm to use the throat lengths as the weights.

[jgostick]

The reason that I put nearest in italics is because this is one of the reasons that the geometric tortuosity fails...real tortuosity is not between two specific points...it's and average between sets of points. So the question is if you find the distances from point A on the left and all points of the right, how do you average or weight these values?