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/_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/_porosimetry.py b/porespy/filters/_porosimetry.py new file mode 100644 index 000000000..ea54817e0 --- /dev/null +++ b/porespy/filters/_porosimetry.py @@ -0,0 +1,139 @@ +import numpy as np +from edt import edt +import scipy.ndimage as spim +from porespy.tools import _insert_disk_at_points_parallel, get_tqdm + + +__all__ = [ + 'porosimetry_si', + 'porosimetry_dt', +] + + +tqdm = get_tqdm() + + +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 + + +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(tqdm(vals)): + seeds = dt >= r + 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=True, + overwrite=False, + ) + lt = lt*im + return lt + + +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 diff --git a/test/unit/test_filters.py b/test/unit/test_filters.py index bbc318458..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 @@ -563,6 +564,22 @@ 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()