From e8faa2ede81eaf5c9245da001c101f7fa762d7be Mon Sep 17 00:00:00 2001 From: jgostick Date: Fri, 20 Jun 2025 17:12:49 -0400 Subject: [PATCH 1/9] removing white space --- src/porespy/filters/_funcs.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/porespy/filters/_funcs.py b/src/porespy/filters/_funcs.py index 9f9430cf4..f3d579246 100644 --- a/src/porespy/filters/_funcs.py +++ b/src/porespy/filters/_funcs.py @@ -688,21 +688,21 @@ def local_thickness( then parallel processing does not occur. `2` is equivalent to `[2, 2, 2]` for a 3D image. The number of cores used is specified in `porespy.settings.ncores` and defaults to all cores. - + parallel_kw : dict Dictionary containing the settings for parallelization by chunking. The optional settings include divs (scalar or list of scalars, default = 1), overlap (scalar or list of scalars, unused in this func), and cores (scalar, default is all available cores). - + Divs is the number of times to divide the image for parallel processing. If `1` then parallel processing does not occur. `2` is equivalent to `[2, 2, 2]` for a 3D image. - + Overlap is the amount of overlap to apply between chunks. In this function, the overlap is controlled by the distance transform and cannot be altered in any way. - + Cores is the number of cores that will be used to parallel process all domains. If ``None`` then all cores will be used but user can specify any integer values to control the memory usage. Setting value to 1 will @@ -865,21 +865,21 @@ def porosimetry( equivalent to `[2, 2, 2]` for a 3D image. The number of cores used is specified in `porespy.settings.ncores` and defaults to all cores. - + parallel_kw : dict Dictionary containing the settings for parallelization by chunking. The optional settings include divs (scalar or list of scalars, default = 1), overlap (scalar or list of scalars, unused in this func), and cores (scalar, default is all available cores). - + Divs is the number of times to divide the image for parallel processing. If `1` then parallel processing does not occur. `2` is equivalent to `[2, 2, 2]` for a 3D image. - + Overlap is the amount of overlap to apply between chunks. In this function, the overlap is controlled by the distance transform and cannot be altered in any way. - + Cores is the number of cores that will be used to parallel process all domains. If ``None`` then all cores will be used but user can specify any integer values to control the memory usage. Setting value to 1 will @@ -1221,23 +1221,23 @@ def chunked_func( func : function handle The function which should be applied to each chunk, such as `spipy.ndimage.binary_dilation`. - + parallel_kw : dict Dictionary containing the settings for parallelization by chunking. The optional settings include divs (scalar or list of scalars, default = [2, 2, 2]), overlap (scalar or list of scalars, optional), and cores (scalar, default is all available cores). - + Divs is the number of times to divide the image for parallel processing. If `1` then parallel processing does not occur. `2` is equivalent to `[2, 2, 2]` for a 3D image. - + Overlap is the amount of overlap to include when dividing up the image. This value will almost always be the size (i.e. raduis) of the structuring element. If not specified then the amount of overlap is inferred from the size of the structuring element, in which case the `strel_arg` must be specified. - + Cores is the number of cores that will be used to parallel process all domains. If ``None`` then all cores will be used but user can specify any integer values to control the memory usage. Setting value to 1 will From e4a0e8c2787ff1080f039b3398580a74d6a99f89 Mon Sep 17 00:00:00 2001 From: jgostick Date: Fri, 20 Jun 2025 17:13:32 -0400 Subject: [PATCH 2/9] moving magnet script to tests/integration --- .../integration/magnet_on_bentheimer.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename magnet_on_bentheimer.py => test/integration/magnet_on_bentheimer.py (100%) diff --git a/magnet_on_bentheimer.py b/test/integration/magnet_on_bentheimer.py similarity index 100% rename from magnet_on_bentheimer.py rename to test/integration/magnet_on_bentheimer.py From a9d0d144a59683a5e445226f699dbb04bd3a18af Mon Sep 17 00:00:00 2001 From: jgostick Date: Fri, 20 Jun 2025 17:14:08 -0400 Subject: [PATCH 3/9] added porosity_threshold function for finding percolation threshold of an image --- src/porespy/metrics/_funcs.py | 109 +++++++++++++++++++++++++++++++++- 1 file changed, 107 insertions(+), 2 deletions(-) diff --git a/src/porespy/metrics/_funcs.py b/src/porespy/metrics/_funcs.py index f4d29c7d6..7b5e0f645 100644 --- a/src/porespy/metrics/_funcs.py +++ b/src/porespy/metrics/_funcs.py @@ -5,15 +5,16 @@ import scipy.ndimage as spim import scipy.spatial as sptl import scipy.stats as spst -from skimage.morphology import skeletonize from numba import njit from scipy import fft as sp_ft from skimage.measure import regionprops +from skimage.morphology import skeletonize, ball, disk, square, cube from porespy import settings from porespy.filters import ( local_thickness, pc_to_seq, fill_closed_pores, + find_invalid_pores, ) from porespy.tools import ( Results, @@ -21,7 +22,8 @@ extend_slice, get_tqdm, get_slices_slabs, - get_edt + get_edt, + ps_round, ) @@ -34,6 +36,7 @@ "pore_size_distribution", "radial_density_distribution", "porosity", + "porosity_threshold", "porosity_profile", "satn_profile", "two_point_correlation", @@ -47,6 +50,108 @@ edt = get_edt() tqdm = get_tqdm() logger = logging.getLogger(__name__) +strel = {2: {'min': disk(1), 'max': square(3)}, 3: {'min': ball(1), 'max': cube(3)}} + + +def porosity_threshold(im, axis=0, method='repeated', conn='min'): + r""" + Find the porosity of the image at the percolation threshold + + Parameters + ---------- + im : ndarray + Image of the void space + axis : int + The axis along which percolation is checked + method : str + The method to use when eroding void space to find percolating threshold. + Options are: + + =========== ================================================================ + Method Description + =========== ================================================================ + repeated Erodes void space repeatedly using a round structuring element + of radius 1. + increasing Erodes void space using round structuring elements of increasing + radius, starting with 1. + =========== ================================================================ + conn : str + Can be either `'min'` or `'max'` and controls the shape of the structuring + element used to determine voxel connectivity. The default if `'min'` which + imposes the strictest criteria, so that voxels must share a face to be + considered connected. + + Returns + ------- + results + A results object with the following attributes: + + ================ =========================================================== + Attribute Description + ================ =========================================================== + eps_orig The total porosity of the original image, including closed + and surface pores + eps_orig_perc The percolating porosity of the original image (i.e. with + closed and surface pores filled) + eps_thresh The porosity of the image after eroding the void space + results in no percolating paths + eps_thresh_perc The percolating porosity of the eroded image (with closed + and surface pores filled). Note that technically NO voids + are percolating at the threshold, so this result uses + the percolating voids from the previous step as a mask + when counting the void volume of percolating voids at the + threshold. + ================ =========================================================== + """ + im = np.swapaxes(im, 0, axis) == 1 + im2 = np.swapaxes(im, 0, axis) == 1 + inlets = np.zeros_like(im) + inlets[0, ...] = True + inlets *= im + outlets = np.zeros_like(im) + outlets[-1, ...] = True + outlets *= im + + Vbulk = np.sum(im >= 0, dtype=np.int64) + invalid = find_invalid_pores(im, conn=conn) + invalid[~im] = -1 + valid = invalid == 0 + eps_orig = np.sum(im > 0, dtype=np.int64)/Vbulk + eps_orig_perc = np.sum(invalid[im] == 0, dtype=np.int64)/Vbulk + + R = 1 + while True: + # Erode void space + if method.startswith('repeat'): + se = ps_round(r=R, ndim=im2.ndim, smooth=False) + im2 = spim.binary_erosion(im2, structure=se, border_value=True) + elif method.startswith('increasing'): + se = ps_round(r=R, ndim=im2.ndim, smooth=False) + R = R + 1 + im2 = spim.binary_erosion(im, structure=se, border_value=True) + # Find labels present on inlet and outlet faces + labels, N = spim.label(im2, structure=strel[im.ndim][conn]) + a = np.unique(labels[inlets]) + a = a[a > 0] + b = np.unique(labels[outlets]) + b = b[b > 0] + if not np.any(np.isin(a, b)): + Vpore = im2 > 0 # Pore volume at current iteration + eps_thresh = np.sum(Vpore, dtype=np.int64)/Vbulk + eps_thresh_perc = np.sum(Vpore*valid, dtype=np.int64)/Vbulk + break + # Find currently percolating void space to use as a mask + invalid = find_invalid_pores(im2, conn=conn) + invalid[~im] = -1 + valid = invalid == 0 + + from porespy.tools import Results + r = Results() + r.eps_orig = eps_orig + r.eps_orig_perc = eps_orig_perc + r.eps_thresh = eps_thresh + r.eps_thresh_perc = eps_thresh_perc + return r def boxcount(im, bins=10): From 1c4c39d422584520b4cac6618ef31026b3a6d30a Mon Sep 17 00:00:00 2001 From: jgostick Date: Fri, 20 Jun 2025 23:14:08 -0400 Subject: [PATCH 4/9] adding is_percolating function, plus changes to porosity_threshold --- src/porespy/metrics/_funcs.py | 92 ++++++++++++++++++++++++----------- 1 file changed, 64 insertions(+), 28 deletions(-) diff --git a/src/porespy/metrics/_funcs.py b/src/porespy/metrics/_funcs.py index 7b5e0f645..8c925de15 100644 --- a/src/porespy/metrics/_funcs.py +++ b/src/porespy/metrics/_funcs.py @@ -32,6 +32,7 @@ "chord_counts", "chord_length_distribution", "find_h", + "is_percolating", "lineal_path_distribution", "pore_size_distribution", "radial_density_distribution", @@ -53,28 +54,65 @@ strel = {2: {'min': disk(1), 'max': square(3)}, 3: {'min': ball(1), 'max': cube(3)}} -def porosity_threshold(im, axis=0, method='repeated', conn='min'): +def is_percolating(im, axis=None, conn='min'): + R""" + + Parameters + ---------- + im : ndarray + Image of the void space with `True` indicating void space + axis : int + The axis along which percolation is checked. If `None` (default) then + percolation is checked in all dimensions + conn : str + Can be either `'min'` or `'max'` and controls the shape of the structuring + element used to determine voxel connectivity. The default if `'min'` which + imposes the strictest criteria, so that voxels must share a face to be + considered connected. + + Returns + ------- + percolating : bool or list of bools + A boolean value indicating if the domain percolates in the given direction. + If `axis=None` then all directions are checked and the result is returned + as a list like `[True, False, True]` indicating that the domain percolates + in the `x` and `z` directions, but not `y`. + + """ + if axis is not None: + im = np.swapaxes(im, 0, axis) == 1 + else: + ans = [] + for ax in range(im.ndim): + ans.append(is_percolating(im, axis=ax, conn=conn)) + return ans + + inlets = np.zeros_like(im) + inlets[0, ...] = True + inlets *= im + outlets = np.zeros_like(im) + outlets[-1, ...] = True + outlets *= im + + labels, N = spim.label(im, structure=strel[im.ndim][conn]) + a = np.unique(labels[inlets]) + a = a[a > 0] + b = np.unique(labels[outlets]) + b = b[b > 0] + hits = np.isin(a, b) + return np.any(hits) + + +def porosity_threshold(im, axis=0, conn='min'): r""" Find the porosity of the image at the percolation threshold Parameters ---------- im : ndarray - Image of the void space + Image of the void space with `True` indicating void space axis : int The axis along which percolation is checked - method : str - The method to use when eroding void space to find percolating threshold. - Options are: - - =========== ================================================================ - Method Description - =========== ================================================================ - repeated Erodes void space repeatedly using a round structuring element - of radius 1. - increasing Erodes void space using round structuring elements of increasing - radius, starting with 1. - =========== ================================================================ conn : str Can be either `'min'` or `'max'` and controls the shape of the structuring element used to determine voxel connectivity. The default if `'min'` which @@ -104,7 +142,7 @@ def porosity_threshold(im, axis=0, method='repeated', conn='min'): ================ =========================================================== """ im = np.swapaxes(im, 0, axis) == 1 - im2 = np.swapaxes(im, 0, axis) == 1 + dt = edt(im) inlets = np.zeros_like(im) inlets[0, ...] = True inlets *= im @@ -119,31 +157,29 @@ def porosity_threshold(im, axis=0, method='repeated', conn='min'): eps_orig = np.sum(im > 0, dtype=np.int64)/Vbulk eps_orig_perc = np.sum(invalid[im] == 0, dtype=np.int64)/Vbulk - R = 1 - while True: + Rmax = int(dt.max()) + desc = inspect.currentframe().f_code.co_name # Get current func name + pbar = tqdm(desc=desc, **settings.tqdm) + for R in range(1, Rmax, 1): # Erode void space - if method.startswith('repeat'): - se = ps_round(r=R, ndim=im2.ndim, smooth=False) - im2 = spim.binary_erosion(im2, structure=se, border_value=True) - elif method.startswith('increasing'): - se = ps_round(r=R, ndim=im2.ndim, smooth=False) - R = R + 1 - im2 = spim.binary_erosion(im, structure=se, border_value=True) + im2 = dt > R # Find labels present on inlet and outlet faces + # Note: This could use is_percolating but labels is needed below labels, N = spim.label(im2, structure=strel[im.ndim][conn]) a = np.unique(labels[inlets]) a = a[a > 0] b = np.unique(labels[outlets]) b = b[b > 0] - if not np.any(np.isin(a, b)): + hits = np.isin(a, b) + if not np.any(hits): Vpore = im2 > 0 # Pore volume at current iteration eps_thresh = np.sum(Vpore, dtype=np.int64)/Vbulk eps_thresh_perc = np.sum(Vpore*valid, dtype=np.int64)/Vbulk break # Find currently percolating void space to use as a mask - invalid = find_invalid_pores(im2, conn=conn) - invalid[~im] = -1 - valid = invalid == 0 + valid = np.isin(labels, a[hits]) + pbar.update(1) + pbar.close() from porespy.tools import Results r = Results() From 3bc08c3bbb0d7a384a42213aee752536933ea1f5 Mon Sep 17 00:00:00 2001 From: jgostick Date: Sat, 21 Jun 2025 16:28:55 -0400 Subject: [PATCH 5/9] removing magnet script from tests --- test/integration/magnet_on_bentheimer.py | 178 ----------------------- 1 file changed, 178 deletions(-) delete mode 100644 test/integration/magnet_on_bentheimer.py diff --git a/test/integration/magnet_on_bentheimer.py b/test/integration/magnet_on_bentheimer.py deleted file mode 100644 index fc132b56e..000000000 --- a/test/integration/magnet_on_bentheimer.py +++ /dev/null @@ -1,178 +0,0 @@ -import pandas as pd -import porespy as ps -import openpnm as op -import numpy as np -import scipy.ndimage as spim -from edt import edt -from skimage.morphology import square, cube - - -# %% Set simulation parameters -np.random.seed(10) -export = False -fill_blind_pores = True -padding = None -l_max=None -bw = None -endpoints=True -name = 'Bentheimer' -res = 5.82e-6 -Kx, Ky, Kz, Kavg = [], [], [], [] -porosity = [] - -# %% Import and prepare image -f = r"C:\Users\jeff\OneDrive - University of Waterloo\Research\Elidek - Bentheimer\bentheimer_binary.npz" -imb = np.load(f)['im'][:500, :500, :500] -imb = ps.filters.fill_blind_pores(imb, conn=6, surface=True) -imb = ps.filters.trim_floating_solid(imb, conn=6, surface=True) -im = imb -porosity.append(np.sum(im)/np.product(im.shape)*100) -print(f'porosity: {np.sum(im)/np.product(im.shape)*100}') - -# %% Get distance transform -dt = edt(im, parallel=8) - -# %% Get Skeleton -# sk = ps.networks.skeletonize_magnet2(im) -sk, im = ps.networks.skeleton(im, padding=None, surface=True, parallel=False) -labels, N = spim.label(sk, structure=ps.tools.ps_rect(3, 3)) -values, counts = np.unique(labels, return_counts=True) -if N > 1: - print(f"There are {N} separated segments in the skeleton") - print(f"Their sizes are {counts[1:]}") - -# %% Find junctions and end points on skeleton -juncs, endpts = ps.networks.analyze_skeleton(sk) -pores, throats = ps.networks.partition_skeleton(sk, juncs + endpts, dt) -new_juncs, pores, new_throats =\ - ps.networks.find_throat_junctions(im=im, pores=pores, throats=throats, dt=dt) -# pores, throats = ps.networks.partition_skeleton(sk, juncs + new_juncs + endpts, dt=dt) - -# net = ps.networks.sk_to_network(pores, throats, dt) - -# %% insert pores at junction points -pts = (pores + new_juncs) > 0 -fbd = ps.networks.insert_pore_bodies(sk, dt, pts) -# convert spheres to network dictionary -net = ps.networks.spheres_to_network( - sk=sk, dt=dt, fbd=fbd, voxel_size=res, boundary_width=0) - - -# %% import network to openpnm -net = op.io.network_from_porespy(net) -net['pore.diameter'] = net['pore.radius']*2 -net['throat.diameter'] = net['throat.radius']*2 - -# check network health -h = op.utils.check_network_health(net) -op.topotools.trim(net, pores=np.append(h['disconnected_pores'], h['isolated_pores'])) -h = op.utils.check_network_health(net) - -# add geometry models -geo = op.models.collections.geometry.cubes_and_cuboids.copy() -del geo['pore.diameter'], geo['throat.diameter'] -net['pore.diameter'] = net['pore.diameter'].copy() -net['throat.diameter'] = net['throat.diameter'].copy() -net.add_model_collection(geo) -net.regenerate_models() - -# add phase -phase = op.phase.Phase(network=net) -phase['pore.viscosity'] = 1e-3 - -# physics -phys = op.models.collections.physics.basic.copy() -phase.add_model_collection(phys) -phase.regenerate_models() - -# stokes flow simulation to estimate permeability -Pin = 5e-6 -Pout = 0 -A = (im.shape[0]*im.shape[1]) * res**2 -L = im.shape[1] * res -mu = phase['pore.viscosity'].max() - -# label boundary pores -xmin = net.coords[:, 0] <= net.coords[:, 0].max()*0.02 -xmax = net.coords[:, 0] >= net.coords[:, 0].max()*0.98 -flow_x = op.algorithms.StokesFlow(network=net, phase=phase) -flow_x.set_value_BC(pores=xmin, values=Pin) -flow_x.set_value_BC(pores=xmax, values=Pout) -flow_x.run() - -Q_x = flow_x.rate(pores=xmin, mode='group')[0] -K_x = Q_x * L * mu / (A * (Pin - Pout))/0.98e-12*1000 -print(f'K_x is: {K_x:.2f} mD') - -ymin = net.coords[:, 1] <= net.coords[:, 1].max()*0.02 -ymax = net.coords[:, 1] >= net.coords[:, 1].max()*0.98 -flow_y = op.algorithms.StokesFlow(network=net, phase=phase) -flow_y.set_value_BC(pores=ymin, values=Pin) -flow_y.set_value_BC(pores=ymax, values=Pout) -flow_y.run() - -Q_y = flow_y.rate(pores=ymin, mode='group')[0] -K_y = Q_y * L * mu / (A * (Pin - Pout))/0.98e-12*1000 -print(f'K_y is: {K_y:.2f} mD') - -zmin = net.coords[:, 2] <= net.coords[:, 2].max()*0.02 -zmax = net.coords[:, 2] >= net.coords[:, 2].max()*0.98 -flow_z = op.algorithms.StokesFlow(network=net, phase=phase) -flow_z.set_value_BC(pores=zmin, values=Pin) -flow_z.set_value_BC(pores=zmax, values=Pout) -flow_z.run() - -Q_z = flow_z.rate(pores=zmin, mode='group')[0] -K_z = Q_z * L * mu / (A * (Pin - Pout))/0.98e-12*1000 -print(f'K_z is: {K_z:.2f} mD') - -K = np.average([K_x, K_y, K_z]) -print(f'K is: {K:.2f} mD') - -# number of pore vs. skeleton clusters in network -from scipy.sparse import csgraph as csg -am = net.create_adjacency_matrix(fmt='coo', triu=True) -Np, cluster_num = csg.connected_components(am, directed=False) -print('Pore clusters:', Np) -# number of skeleton pieces -b = square(3) if im.ndim == 2 else cube(3) -_, Ns = spim.label(input=sk.astype('bool'), structure=b) -print('Skeleton clusters:', Ns) - -print('The expected permeability is 2490 mD') - -# append data -Kx.append(K_x) -Ky.append(K_y) -Kz.append(K_z) -Kavg.append(K) - - -# export -if export: - ps.io.to_stl(~sk, name + '-sk') - ps.io.to_stl(im, name + '-im') - net['pore.coords'] = net['pore.coords']/res - net['pore.diameter'] = net['pore.diameter']/res - net['throat.radius'] = net['throat.radius']/res - net['pore.coords'] += 10 # res - proj = net.project - op.io.project_to_xdmf(proj, filename=name + '-network') - -# create and export data frame -if 0: - data = { - "Name": name, - "Resolution": res, - "boundary_width": bw, - "l_max": l_max, - "padding": str(padding), - "endpoints": endpoints, - "porosity": porosity, - "K_x": Kx, - "K_y": Ky, - "K_z": Kz, - "K": Kavg - } - df = pd.DataFrame(data=data) - df.to_csv('MAGNET Permeability.csv', index=False) From 34f1570d4b208ee94eb5c286b013f5a9e986c250 Mon Sep 17 00:00:00 2001 From: jgostick Date: Mon, 23 Jun 2025 21:29:55 -0400 Subject: [PATCH 6/9] added percolating_porosity function, renamed porosity_threshold to find_porosity_threshold --- src/porespy/metrics/_funcs.py | 142 +++++++++++++++++++--------------- 1 file changed, 81 insertions(+), 61 deletions(-) diff --git a/src/porespy/metrics/_funcs.py b/src/porespy/metrics/_funcs.py index 8c925de15..82807cfdb 100644 --- a/src/porespy/metrics/_funcs.py +++ b/src/porespy/metrics/_funcs.py @@ -37,10 +37,11 @@ "pore_size_distribution", "radial_density_distribution", "porosity", - "porosity_threshold", + "find_porosity_threshold", "porosity_profile", "satn_profile", "two_point_correlation", + "percolating_porosity", "phase_fraction", "pc_curve", "pc_map_to_pc_curve", @@ -54,16 +55,18 @@ strel = {2: {'min': disk(1), 'max': square(3)}, 3: {'min': ball(1), 'max': cube(3)}} -def is_percolating(im, axis=None, conn='min'): +def is_percolating(im, axis=None, inlets=None, outlets=None, conn='min'): R""" + Determines if a percolating path exists across the domain (in the specified + direction) or between given inlets and outlets. Parameters ---------- im : ndarray - Image of the void space with `True` indicating void space + Image of the void space with `True` indicating void space. axis : int The axis along which percolation is checked. If `None` (default) then - percolation is checked in all dimensions + percolation is checked in all dimensions. conn : str Can be either `'min'` or `'max'` and controls the shape of the structuring element used to determine voxel connectivity. The default if `'min'` which @@ -79,21 +82,22 @@ def is_percolating(im, axis=None, conn='min'): in the `x` and `z` directions, but not `y`. """ - if axis is not None: + if (inlets is not None) and (outlets is not None): + pass + elif axis is not None: im = np.swapaxes(im, 0, axis) == 1 + inlets = np.zeros_like(im) + inlets[0, ...] = True + inlets *= im + outlets = np.zeros_like(im) + outlets[-1, ...] = True + outlets *= im else: ans = [] for ax in range(im.ndim): ans.append(is_percolating(im, axis=ax, conn=conn)) return ans - inlets = np.zeros_like(im) - inlets[0, ...] = True - inlets *= im - outlets = np.zeros_like(im) - outlets[-1, ...] = True - outlets *= im - labels, N = spim.label(im, structure=strel[im.ndim][conn]) a = np.unique(labels[inlets]) a = a[a > 0] @@ -103,10 +107,13 @@ def is_percolating(im, axis=None, conn='min'): return np.any(hits) -def porosity_threshold(im, axis=0, conn='min'): +def find_porosity_threshold(im, axis=0, conn='min'): r""" Find the porosity of the image at the percolation threshold + This function progressively dilates the solid and reports the porosity at the + step just before there are no percolating paths (in the specified direction) + Parameters ---------- im : ndarray @@ -131,65 +138,78 @@ def porosity_threshold(im, axis=0, conn='min'): and surface pores eps_orig_perc The percolating porosity of the original image (i.e. with closed and surface pores filled) - eps_thresh The porosity of the image after eroding the void space - results in no percolating paths + eps_thresh The total porosity of the image after eroding the void + space results in no percolating paths eps_thresh_perc The percolating porosity of the eroded image (with closed - and surface pores filled). Note that technically NO voids - are percolating at the threshold, so this result uses - the percolating voids from the previous step as a mask - when counting the void volume of percolating voids at the - threshold. + and surface pores filled) ================ =========================================================== """ - im = np.swapaxes(im, 0, axis) == 1 + def _check_percolation(dt, R, step, axis, conn): + while True: + im2 = dt >= R + if not is_percolating(im2, axis=axis, conn=conn): + break + R += step + return R + dt = edt(im) - inlets = np.zeros_like(im) - inlets[0, ...] = True - inlets *= im - outlets = np.zeros_like(im) - outlets[-1, ...] = True - outlets *= im - - Vbulk = np.sum(im >= 0, dtype=np.int64) - invalid = find_invalid_pores(im, conn=conn) - invalid[~im] = -1 - valid = invalid == 0 - eps_orig = np.sum(im > 0, dtype=np.int64)/Vbulk - eps_orig_perc = np.sum(invalid[im] == 0, dtype=np.int64)/Vbulk - - Rmax = int(dt.max()) - desc = inspect.currentframe().f_code.co_name # Get current func name - pbar = tqdm(desc=desc, **settings.tqdm) - for R in range(1, Rmax, 1): - # Erode void space - im2 = dt > R - # Find labels present on inlet and outlet faces - # Note: This could use is_percolating but labels is needed below - labels, N = spim.label(im2, structure=strel[im.ndim][conn]) - a = np.unique(labels[inlets]) - a = a[a > 0] - b = np.unique(labels[outlets]) - b = b[b > 0] - hits = np.isin(a, b) - if not np.any(hits): - Vpore = im2 > 0 # Pore volume at current iteration - eps_thresh = np.sum(Vpore, dtype=np.int64)/Vbulk - eps_thresh_perc = np.sum(Vpore*valid, dtype=np.int64)/Vbulk - break - # Find currently percolating void space to use as a mask - valid = np.isin(labels, a[hits]) - pbar.update(1) - pbar.close() + + R = _check_percolation(dt, R=1, step=10, axis=axis, conn=conn) + R = _check_percolation(dt, R=max(1, R-10), step=5, axis=axis, conn=conn) + R = _check_percolation(dt, R=max(1, R-5), step=1, axis=axis, conn=conn) + + im2 = dt >= (R - 1) + eps_thresh_total = porosity(im2) + eps_thresh_perc = percolating_porosity(im2, axis=axis) from porespy.tools import Results r = Results() - r.eps_orig = eps_orig - r.eps_orig_perc = eps_orig_perc - r.eps_thresh = eps_thresh + r.eps_orig = porosity(im) + r.eps_orig_perc = percolating_porosity(im, axis=axis) + r.eps_thresh = eps_thresh_total r.eps_thresh_perc = eps_thresh_perc return r +def percolating_porosity(im, axis=0, inlets=None, outlets=None, conn='min'): + r""" + Finds volume fraction of void space which belongs to percolating paths + across the domain in the direction specified. + + Parameters + ---------- + im : ndarray + Image of the void space with `True` indicating void space + axis : int + The axis along which percolation is checked + conn : str + Can be either `'min'` or `'max'` and controls the shape of the structuring + element used to determine voxel connectivity. The default if `'min'` which + imposes the strictest criteria, so that voxels must share a face to be + considered connected. + inlets, outlets : ndarray + Boolean arrays indicating the locations of the inlets and outlets. These + are useful if the domain is not cubic or if special inlet and outlet + locations are desired. + """ + se = strel[im.ndim][conn] + if (inlets is None) and (outlets is None): + im2 = np.swapaxes(im, 0, axis) + inlets = np.zeros_like(im2, dtype=bool) + inlets[0, ...] = True + outlets = np.zeros_like(im2, dtype=bool) + outlets[-1, ...] = True + labels, N = spim.label(im2, structure=se) + a = np.unique(labels*inlets) + a = a[a > 0] + b = np.unique(labels*outlets) + b = b[b > 0] + hits = np.intersect1d(a, b) + im3 = np.isin(labels, hits) + eps = porosity(im3) + return eps + + def boxcount(im, bins=10): r""" Calculates the fractal dimension of an image using the tiled box counting From 8f441bf02c3e52e8479c7d1c5ddb8ef5b0f55a99 Mon Sep 17 00:00:00 2001 From: jgostick Date: Tue, 24 Jun 2025 10:56:08 -0400 Subject: [PATCH 7/9] adding tests --- src/porespy/metrics/_funcs.py | 16 ++++++++++++---- test/unit/test_metrics.py | 28 ++++++++++++++++++++++++++++ 2 files changed, 40 insertions(+), 4 deletions(-) diff --git a/src/porespy/metrics/_funcs.py b/src/porespy/metrics/_funcs.py index 82807cfdb..e84f6c67a 100644 --- a/src/porespy/metrics/_funcs.py +++ b/src/porespy/metrics/_funcs.py @@ -16,6 +16,10 @@ fill_closed_pores, find_invalid_pores, ) +from porespy.generators import ( + faces, + borders, +) from porespy.tools import ( Results, _check_for_singleton_axes, @@ -109,7 +113,7 @@ def is_percolating(im, axis=None, inlets=None, outlets=None, conn='min'): def find_porosity_threshold(im, axis=0, conn='min'): r""" - Find the porosity of the image at the percolation threshold + Finds the porosity of the image at the percolation threshold This function progressively dilates the solid and reports the porosity at the step just before there are no percolating paths (in the specified direction) @@ -144,10 +148,14 @@ def find_porosity_threshold(im, axis=0, conn='min'): and surface pores filled) ================ =========================================================== """ + if axis is None: + raise Exception('axis must be specified') + def _check_percolation(dt, R, step, axis, conn): while True: im2 = dt >= R - if not is_percolating(im2, axis=axis, conn=conn): + check = np.array(is_percolating(im2, axis=axis, conn=conn)) + if not np.all(check): break R += step return R @@ -155,8 +163,8 @@ def _check_percolation(dt, R, step, axis, conn): dt = edt(im) R = _check_percolation(dt, R=1, step=10, axis=axis, conn=conn) - R = _check_percolation(dt, R=max(1, R-10), step=5, axis=axis, conn=conn) - R = _check_percolation(dt, R=max(1, R-5), step=1, axis=axis, conn=conn) + R = _check_percolation(dt, R=max(1, R-10), step=4, axis=axis, conn=conn) + R = _check_percolation(dt, R=max(1, R-4), step=1, axis=axis, conn=conn) im2 = dt >= (R - 1) eps_thresh_total = porosity(im2) diff --git a/test/unit/test_metrics.py b/test/unit/test_metrics.py index 9f16da3a3..686186273 100644 --- a/test/unit/test_metrics.py +++ b/test/unit/test_metrics.py @@ -514,6 +514,34 @@ def test_bond_number(self): source='dt', method='max') assert np.around(bo, decimals=5) == 0.79461 + def test_is_percolating(self): + im = np.ones([20, 20], dtype=bool) + im[:, 10:12] = False + im[10, 10] = True + im[11, 11] = True + assert ps.metrics.is_percolating(im=im, axis=0) + assert not ps.metrics.is_percolating(im=im, axis=1, conn='min') + assert ps.metrics.is_percolating(im=im, axis=1, conn='max') + assert ps.metrics.is_percolating(im=im, axis=None) == [True, False] + assert ps.metrics.is_percolating(im=im, axis=None, conn='max') == [True, True] + + def test_pecolating_porosity(self): + im = np.ones([20, 20], dtype=bool) + im[:10, 10] = False + im[10:, 11] = False + e = im.sum()/im.size + assert ps.metrics.percolating_porosity(im, axis=0, conn='min') == e + assert ps.metrics.percolating_porosity(im, axis=1, conn='min') == 0 + assert ps.metrics.percolating_porosity(im, axis=0, conn='max') == e + assert ps.metrics.percolating_porosity(im, axis=1, conn='max') == e + + def test_find_porosity_threshold(self): + im = ps.generators.lattice_spheres([31, 51], r=15, offset=5, spacing=30) + ep = ps.metrics.find_porosity_threshold(im, axis=0) + assert ep.eps_thresh == 0.4484503478810879 + ep = ps.metrics.find_porosity_threshold(im, axis=1) + assert ep.eps_thresh == 0.05439595192915876 + if __name__ == '__main__': t = MetricsTest() From 9c9d374a0aae9ec8b5adb8b0c337eadb79f7e049 Mon Sep 17 00:00:00 2001 From: jgostick Date: Tue, 24 Jun 2025 10:58:58 -0400 Subject: [PATCH 8/9] adding R to results object so the percolating image can be recovered --- src/porespy/metrics/_funcs.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/porespy/metrics/_funcs.py b/src/porespy/metrics/_funcs.py index e84f6c67a..4d4c981e9 100644 --- a/src/porespy/metrics/_funcs.py +++ b/src/porespy/metrics/_funcs.py @@ -133,7 +133,7 @@ def find_porosity_threshold(im, axis=0, conn='min'): Returns ------- results - A results object with the following attributes: + A Results object with the following attributes: ================ =========================================================== Attribute Description @@ -146,6 +146,8 @@ def find_porosity_threshold(im, axis=0, conn='min'): space results in no percolating paths eps_thresh_perc The percolating porosity of the eroded image (with closed and surface pores filled) + R The threshold to apply to the distance transform to + obtain the percolating image (i.e. im = dt >= R) ================ =========================================================== """ if axis is None: @@ -176,6 +178,7 @@ def _check_percolation(dt, R, step, axis, conn): r.eps_orig_perc = percolating_porosity(im, axis=axis) r.eps_thresh = eps_thresh_total r.eps_thresh_perc = eps_thresh_perc + r.R = R - 1 return r From 84fb9946571c304306e38f48e3807632a4a75e38 Mon Sep 17 00:00:00 2001 From: jgostick Date: Tue, 24 Jun 2025 16:24:20 -0400 Subject: [PATCH 9/9] removing unused imports --- src/porespy/metrics/_funcs.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/porespy/metrics/_funcs.py b/src/porespy/metrics/_funcs.py index 4d4c981e9..1094e1369 100644 --- a/src/porespy/metrics/_funcs.py +++ b/src/porespy/metrics/_funcs.py @@ -16,10 +16,6 @@ fill_closed_pores, find_invalid_pores, ) -from porespy.generators import ( - faces, - borders, -) from porespy.tools import ( Results, _check_for_singleton_axes,