From 103e7a04579ba3e818f90699eabe95c8bd2c7c5b Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Thu, 21 Dec 2023 22:23:33 -0500 Subject: [PATCH 1/5] doing lt with spheres, not actually faster --- porespy/filters/_funcs.py | 2 +- porespy/filters/_lt_mip.py | 88 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 89 insertions(+), 1 deletion(-) create mode 100644 porespy/filters/_lt_mip.py diff --git a/porespy/filters/_funcs.py b/porespy/filters/_funcs.py index ceefecd70..a909d2312 100644 --- a/porespy/filters/_funcs.py +++ b/porespy/filters/_funcs.py @@ -1161,7 +1161,7 @@ def porosimetry(im, sizes=25, inlets=None, access_limited=True, mode='hybrid', overlap=int(r) + 1, parallel=0, cores=settings.ncores, divs=divs) < r else: - imtemp = edt(~imtemp) < r + imtemp = edt(~imtemp, parallel=-1) < r imresults[(imresults == 0) * imtemp] = r elif mode == "hybrid": imresults = np.zeros(np.shape(im)) diff --git a/porespy/filters/_lt_mip.py b/porespy/filters/_lt_mip.py new file mode 100644 index 000000000..127e05f4f --- /dev/null +++ b/porespy/filters/_lt_mip.py @@ -0,0 +1,88 @@ +import numpy as np +from edt import edt +import scipy.ndimage as spim +import porespy as ps +from porespy.tools import _insert_disk_at_points_parallel +import matplotlib.pyplot as plt + + +im = ps.generators.blobs([300, 300, 300], seed=0) +dt = edt(im, parallel=-1) + + +def local_thickness_2(im, dt): + vals = np.arange(np.floor(dt.max()).astype(int), 0, -1) + lt = np.zeros_like(dt, dtype=int) + for i, r in enumerate(vals): + coords = np.vstack(np.where(dt == r)) + _insert_disk_at_points_parallel( + im=lt, + coords=coords, + r=r, + v=r, + smooth=False, + overwrite=False, + ) + lt = lt*im + return lt + + +# a = local_thickness_2(im=im, dt=dt) +# b = ps.filters.local_thickness(im=im, sizes=range(1, a.max())) +# c = a - b +# plt.imshow(c) + + +def trim_disconnected_seeds(im, inlets): + labels, N = spim.label(im+inlets) + labels = labels*im + keep = np.unique(labels[inlets]) + keep = keep[keep > 0] + seeds = np.isin(labels, keep) + return seeds + + +inlets = np.zeros_like(im, dtype=bool) +inlets[0, :] = True + + +def porosimetry_2(im, dt, inlets): + vals = np.arange(np.floor(dt.max()).astype(int), 0, -1) + lt = np.zeros_like(dt, dtype=int) + for i, r in enumerate(vals): + seeds = dt >= r + seeds = trim_disconnected_seeds(im=seeds, inlets=inlets) + coords = np.vstack(np.where(seeds)) + _insert_disk_at_points_parallel( + im=lt, + coords=coords, + r=r, + v=r, + smooth=False, + overwrite=False, + ) + lt = lt*im + return lt + +ps.tools.tic() +d = porosimetry_2(im=im, dt=dt, inlets=inlets) +ps.tools.toc() +ps.tools.tic() +vals = np.arange(np.floor(dt.max()).astype(int), 0, -1) +e = ps.filters.porosimetry(im=im, sizes=vals, inlets=inlets, mode='dt') +ps.tools.toc() +f = d - e +# plt.imshow(f) + + + + + + + + + + + + + From d686e03137085b98929ba18949f0774a14f44b1e Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Thu, 4 Jan 2024 21:39:27 -0500 Subject: [PATCH 2/5] after a lot of fiddling, the si method is only slighly faster than dt --- porespy/filters/_lt_mip.py | 145 +++++++++++++++++++++++++------------ 1 file changed, 98 insertions(+), 47 deletions(-) diff --git a/porespy/filters/_lt_mip.py b/porespy/filters/_lt_mip.py index 127e05f4f..ea54817e0 100644 --- a/porespy/filters/_lt_mip.py +++ b/porespy/filters/_lt_mip.py @@ -1,36 +1,16 @@ import numpy as np from edt import edt import scipy.ndimage as spim -import porespy as ps -from porespy.tools import _insert_disk_at_points_parallel -import matplotlib.pyplot as plt +from porespy.tools import _insert_disk_at_points_parallel, get_tqdm -im = ps.generators.blobs([300, 300, 300], seed=0) -dt = edt(im, parallel=-1) +__all__ = [ + 'porosimetry_si', + 'porosimetry_dt', +] -def local_thickness_2(im, dt): - vals = np.arange(np.floor(dt.max()).astype(int), 0, -1) - lt = np.zeros_like(dt, dtype=int) - for i, r in enumerate(vals): - coords = np.vstack(np.where(dt == r)) - _insert_disk_at_points_parallel( - im=lt, - coords=coords, - r=r, - v=r, - smooth=False, - overwrite=False, - ) - lt = lt*im - return lt - - -# a = local_thickness_2(im=im, dt=dt) -# b = ps.filters.local_thickness(im=im, sizes=range(1, a.max())) -# c = a - b -# plt.imshow(c) +tqdm = get_tqdm() def trim_disconnected_seeds(im, inlets): @@ -42,47 +22,118 @@ def trim_disconnected_seeds(im, inlets): return seeds -inlets = np.zeros_like(im, dtype=bool) -inlets[0, :] = True - - -def porosimetry_2(im, dt, inlets): +def porosimetry_si(im, dt, inlets=None): + r""" + Perform image-based porosimetry using sphere insertion + + Parameters + ---------- + im : ndarray + The boolean array with `True` indicating the phase of interest + dt : ndarray + The distance transform of `im`. If not provided it will be calculated + inlets : ndarray + A boolean array with `True` indicating the locations where non-wetting + fluid enters the domain. If `None` (default) then access limitations are + ignored and the result will correspond to the local thickness. + + Returns + ------- + results : ndarray + An array with the each voxel containing the radius of the largest sphere + that overlaps it. + + Notes + ----- + This function use direct sphere insertion to draw spheres at every location + where one can fit. + + """ + if dt is None: + dt = edt(im, parallel=-1) vals = np.arange(np.floor(dt.max()).astype(int), 0, -1) lt = np.zeros_like(dt, dtype=int) - for i, r in enumerate(vals): + for i, r in enumerate(tqdm(vals)): seeds = dt >= r - seeds = trim_disconnected_seeds(im=seeds, inlets=inlets) + if inlets is not None: + seeds = trim_disconnected_seeds(im=seeds, inlets=inlets) + lt[(lt == 0)*seeds] = r + seeds = seeds * (dt < (r + 1)) coords = np.vstack(np.where(seeds)) _insert_disk_at_points_parallel( im=lt, coords=coords, r=r, v=r, - smooth=False, + smooth=True, overwrite=False, ) lt = lt*im return lt -ps.tools.tic() -d = porosimetry_2(im=im, dt=dt, inlets=inlets) -ps.tools.toc() -ps.tools.tic() -vals = np.arange(np.floor(dt.max()).astype(int), 0, -1) -e = ps.filters.porosimetry(im=im, sizes=vals, inlets=inlets, mode='dt') -ps.tools.toc() -f = d - e -# plt.imshow(f) - - - - +def porosimetry_dt(im, dt, inlets=None): + r""" + Perform image-based porosimetry using distance transforms + + Parameters + ---------- + im : ndarray + The boolean array with `True` indicating the phase of interest + dt : ndarray + The distance transform of `im`. If not provided it will be calculated + inlets : ndarray + A boolean array with `True` indicating the locations where non-wetting + fluid enters the domain. If `None` (default) then access limitations are + ignored and the result will correspond to the local thickness. + + Returns + ------- + results : ndarray + An array with the each voxel containing the radius of the largest sphere + that overlaps it. + + Notes + ----- + This function use distance transforms to draw spheres + + """ + if dt is None: + dt = edt(im, parallel=-1) + vals = np.arange(np.floor(dt.max()).astype(int)) + lt = np.zeros_like(dt, dtype=int) + for i, r in enumerate(tqdm(vals)): + seeds = dt >= r + if inlets is not None: + seeds = trim_disconnected_seeds(im=seeds, inlets=inlets) + blobs = edt(~seeds, parallel=-1) < r + lt[blobs] = r + lt = lt*im + return lt +if __name__ == "__main__": + import porespy as ps + im = ps.generators.blobs([100, 100, 100], seed=0, porosity=0.7) + dt = edt(im, parallel=-1) + inlets = np.zeros_like(im, dtype=bool) + inlets[0, :] = True + ps.tools.tic() + d = porosimetry_si(im=im, dt=dt, inlets=inlets) + ps.tools.toc() + ps.tools.tic() + vals = np.arange(np.floor(dt.max()).astype(int), 0, -1) + e = ps.filters.porosimetry(im=im, sizes=vals, inlets=inlets, mode='dt') + ps.tools.toc() + ps.tools.tic() + f = porosimetry_dt(im=im, dt=dt, inlets=inlets) + ps.tools.toc() + # Make sure all three functions return exact same result + assert np.sum(d - e) == 0 + assert np.sum(e - f) == 0 From d9c8ff35d954f5b861b85604851d0ff4b2fc66a5 Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Thu, 4 Jan 2024 21:41:19 -0500 Subject: [PATCH 3/5] adding to init file --- porespy/filters/__init__.py | 3 ++- porespy/filters/{_lt_mip.py => _porosimetry.py} | 0 2 files changed, 2 insertions(+), 1 deletion(-) rename porespy/filters/{_lt_mip.py => _porosimetry.py} (100%) diff --git a/porespy/filters/__init__.py b/porespy/filters/__init__.py index 25eaf9411..ff865efcb 100644 --- a/porespy/filters/__init__.py +++ b/porespy/filters/__init__.py @@ -5,7 +5,7 @@ This module contains a variety of functions for altering images based on the structural characteristics, such as pore sizes. A definition of a -*filter* is a function that returns an image the shape as the original +*filter* is a function that returns an image the same shape as the original image, but with altered values. .. currentmodule:: porespy @@ -62,3 +62,4 @@ from ._nlmeans import * from ._fftmorphology import * from . import imagej +from ._porosimetry import * diff --git a/porespy/filters/_lt_mip.py b/porespy/filters/_porosimetry.py similarity index 100% rename from porespy/filters/_lt_mip.py rename to porespy/filters/_porosimetry.py From 4a45ebf067f9ba39e71c30a060f27775b39bb2f0 Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Thu, 4 Jan 2024 21:43:47 -0500 Subject: [PATCH 4/5] added test --- test/unit/test_filters.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/test/unit/test_filters.py b/test/unit/test_filters.py index bbc318458..76b3fe55f 100644 --- a/test/unit/test_filters.py +++ b/test/unit/test_filters.py @@ -563,6 +563,21 @@ def test_regions_size(self): hits = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 15, 16, 17, 19, 31, 32, 37] assert np.all(hits == np.unique(s)[1:]) + def test_porosimetry_dt_and_si(self): + im = ps.generators.blobs([100, 100, 100], seed=0, porosity=0.7) + dt = edt(im, parallel=-1) + + inlets = np.zeros_like(im, dtype=bool) + inlets[0, :] = True + + d = ps.filters.porosimetry_si(im=im, dt=dt, inlets=inlets) + vals = np.arange(np.floor(dt.max()).astype(int), 0, -1) + e = ps.filters.porosimetry(im=im, sizes=vals, inlets=inlets, mode='dt') + f = ps.filters.porosimetry_dt(im=im, dt=dt, inlets=inlets) + + # Make sure all three functions return exact same result + assert np.sum(d - e) == 0 + assert np.sum(e - f) == 0 if __name__ == '__main__': t = FilterTest() From 46df46d8ad6f9fc837cdd95ed703db55e937adea Mon Sep 17 00:00:00 2001 From: Jeff Gostick Date: Sat, 6 Jan 2024 23:01:29 -0500 Subject: [PATCH 5/5] fixing pep8 --- test/unit/test_filters.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/unit/test_filters.py b/test/unit/test_filters.py index 76b3fe55f..28881cf58 100644 --- a/test/unit/test_filters.py +++ b/test/unit/test_filters.py @@ -5,7 +5,8 @@ import scipy.ndimage as spim from skimage.morphology import disk, ball, skeletonize_3d from skimage.util import random_noise -from scipy.stats import norm + + ps.settings.tqdm['disable'] = True @@ -579,6 +580,7 @@ def test_porosimetry_dt_and_si(self): assert np.sum(d - e) == 0 assert np.sum(e - f) == 0 + if __name__ == '__main__': t = FilterTest() self = t