Parallel extract

src/porespy/filters/_funcs.py

@@ -2,7 +2,6 @@
 import logging
 import dask
 import numpy as np
-from edt import edt
 import operator as op
 import scipy.ndimage as spim
 from deprecated import deprecated
@@ -16,6 +15,10 @@
 from porespy import settings
 from porespy.tools import get_tqdm
 from typing import Literal
+try:
+    from pyedt import edt
+except ModuleNotFoundError:
+    from edt import edt
 
 
 __all__ = [
@@ -1183,7 +1186,7 @@
                                                  strel=strel_2(1))
             if np.any(imtemp):
                 if parallel:
-                    imtemp = chunked_func(func=edt,
+                    imtemp = chunked_func(func=lambda x: edt(x),
                                           data=~imtemp, im_arg='data',
                                           overlap=int(r) + 1, parallel=0,
                                           cores=settings.ncores, divs=divs) < r

src/porespy/filters/_snows.py

@@ -3,7 +3,6 @@
 import logging
 import numpy as np
 from numba import njit, prange
-from edt import edt
 import scipy.ndimage as spim
 import scipy.spatial as sptl
 from skimage.segmentation import watershed
@@ -14,6 +13,10 @@
 from porespy.tools import get_tqdm
 from porespy.filters import chunked_func
 from porespy import settings
+try:
+    from pyedt import edt
+except ModuleNotFoundError:
+    from edt import edt
 
 
 __all__ = [
@@ -106,7 +109,7 @@
         logger.info("Converting supplied image to boolean")
         im = im > 0
     if dt is None:
-        logger.info("Peforming distance transform")
+        logger.info("Performing distance transform")
         if np.any(im_shape == 1):
             dt = edt(im.squeeze())
             dt = dt.reshape(im_shape)
@@ -566,7 +569,10 @@
     labels, N = spim.label(peaks > 0, structure=cube(3))
     crds = spim.measurements.center_of_mass(peaks > 0, labels=labels,
                                             index=np.arange(1, N + 1))
-    crds = np.vstack(crds).astype(int)  # Convert to numpy array of ints
+    try:
+        crds = np.vstack(crds).astype(int)  # Convert to numpy array of ints
+    except ValueError:
+        return peaks
     L = dt[tuple(crds.T)]  # Get distance to solid for each peak
     # Add tiny amount to joggle points to avoid equal distances to solid
     # arange was added instead of random values so the results are repeatable
@@ -597,14 +603,14 @@
     if mode == 'watershed':
         rev = spim.interpolation.zoom(im, zoom=zoom, order=0)
         rev = rev > 0
-        dt = edt(rev, parallel=0)
+        dt = edt(rev)
         rev_snow = snow_partitioning(rev, dt=dt)
         labels, counts = np.unique(rev_snow, return_counts=True)
         node = np.where(counts == counts[1:].max())[0][0]
         slices = spim.find_objects(rev_snow)
         overlap = max(rev_snow[slices[node - 1]].shape) / (zoom * 2.0)
     if mode == 'dt':
-        dt = edt((im > 0), parallel=0)
+        dt = edt((im > 0))
         overlap = dt.max()
     return overlap
 
@@ -679,7 +685,7 @@
     overlap = overlap / 2.0
     logger.debug(f'Overlap thickness: {int(2 * overlap)} voxels')
 
-    dt = edt((im > 0), parallel=0)
+    dt = edt((im > 0))
 
     # Get overlap and trim depth of all image dimension
     depth = {}
@@ -1089,7 +1095,10 @@
     dt2 = spim.gaussian_filter(input=dt, sigma=sigma)
     peaks = find_peaks(dt=dt2, r_max=r_max)
     peaks = trim_saddle_points(peaks=peaks, dt=dt)
-    peaks = trim_nearby_peaks(peaks=peaks, dt=dt)
-    peaks, N = spim.label(peaks > 0)
-    regions = watershed(image=-dt, markers=peaks)
+    if len(peaks) > 0:
+        peaks = trim_nearby_peaks(peaks=peaks, dt=dt)
+        peaks, N = spim.label(peaks > 0)
+        regions = watershed(image=-dt, markers=peaks)
+    else:
+        regions = np.ones_like(dt2)
     return regions * (dt > 0)

src/porespy/generators/_imgen.py

@@ -1,20 +1,30 @@
 import logging
 import numpy as np
 import inspect as insp
-from edt import edt
-import porespy as ps
 from numba import njit
 import scipy.spatial as sptl
 import scipy.ndimage as spim
 import scipy.stats as spst
-from deprecated import deprecated
-from porespy.tools import all_to_uniform, ps_ball, ps_disk, get_border, ps_round
-from porespy.tools import extract_subsection
-from porespy.tools import insert_sphere
-from porespy.tools import _insert_disk_at_points, _insert_disk_at_points_parallel
+from porespy import metrics
 from porespy import settings
 from typing import List, Literal
+from porespy.filters import chunked_func
+from porespy.tools import (
+    all_to_uniform,
+    ps_ball,
+    ps_disk,
+    get_border,
+    extract_subsection,
+    insert_sphere,
+    get_tqdm,
+    _insert_disk_at_points,
+    _insert_disk_at_points_parallel,
+)
 import numpy.typing as npt
+try:
+    from pyedt import edt
+except ModuleNotFoundError:
+    from edt import edt
 
 
 __all__ = [
@@ -34,7 +44,7 @@
 ]
 
 
-tqdm = ps.tools.get_tqdm()
+tqdm = get_tqdm()
 logger = logging.getLogger(__name__)
 
 
@@ -1040,9 +1050,9 @@ def blobs(
     im = np.random.random(shape)
     if parallel:
         overlap = max([int(s*4) for s in np.array(sigma, ndmin=1)])
-        im = ps.filters.chunked_func(func=spim.gaussian_filter,
-                                     input=im, sigma=sigma,
-                                     divs=divs, overlap=overlap)
+        im = chunked_func(func=spim.gaussian_filter,
+                          input=im, sigma=sigma,
+                          divs=divs, overlap=overlap)
     else:
         im = spim.gaussian_filter(im, sigma=sigma)
     im = all_to_uniform(im, scale=[0, 1])
@@ -1302,8 +1312,8 @@ def get_num_pixels(porosity):
             im = im * tmp
         n_fibers_added += n_fibers
         # Update parameters for next iteration
-        porosity = ps.metrics.porosity(im)
-        vol_added = get_num_pixels(porosity)
+        eps = metrics.porosity(im)
+        vol_added = get_num_pixels(eps)
         vol_fiber = vol_added / n_fibers_added
 
     logger.debug(f"{n_fibers_added} fibers added to reach target porosity.")

src/porespy/generators/_pseudo_packings.py

@@ -3,14 +3,22 @@
 import numpy as np
 import scipy.ndimage as spim
 import numpy.typing as npt
-from edt import edt
+from numba import njit
+from typing import Literal, List
 from skimage.morphology import disk, ball
 from porespy import settings
-from porespy.tools import get_tqdm, ps_round, get_border, unpad
-from porespy.tools import _insert_disk_at_point
+from porespy.tools import (
+    get_tqdm,
+    ps_round,
+    get_border,
+    unpad,
+    _insert_disk_at_point,
+)
 from porespy.filters import trim_disconnected_blobs
-from numba import njit
-from typing import Literal, List
+try:
+    from pyedt import edt
+except ModuleNotFoundError:
+    from edt import edt
 
 
 __all__ = [
@@ -369,7 +377,7 @@ def pseudo_electrostatic_packing(
         dt = spim.gaussian_filter(dt, sigma=0.5)
         strel = ps_round(r, ndim=im.ndim, smooth=True)
         sites = (spim.maximum_filter(dt, footprint=strel) == dt)*(mask > 0)
-        if np.any(dt == np.inf):  # In case above method failed.
+        if np.any(dt == np.inf) or np.all(dt == dt[0]):  # In case above method failed.
             sites = np.zeros_like(im)
             inds = tuple((np.array(im.shape)/2).astype(int))
             sites[inds] = True
@@ -461,13 +469,13 @@ def _do_packing(im, mask, q, r, value, clearance, smooth, maxiter):
         sites[150, 150] = True
         im = ps.generators.pseudo_electrostatic_packing(
             shape=[300, 300],
-            r=15,
-            sites=sites,
-            maxiter=30,
-            edges='contained',
-            clearance=0,
-            smooth=False,
-            compactness=0.1,
+            r=5,
+            # sites=sites,
+            # maxiter=1000,
+            edges='extended',
+            # clearance=0,
+            # smooth=True,
+            compactness=1.0,
         )
         ax[0][0].imshow(im)
     if 1:

src/porespy/io/_funcs.py

@@ -1,15 +1,16 @@
 import os
 import subprocess
-
 import numpy as np
 import scipy.ndimage as nd
 import skimage.measure as ms
-from edt import edt
 from skimage.morphology import ball
-
 from porespy.filters import reduce_peaks
 from porespy.networks import generate_voxel_image
 from porespy.tools import sanitize_filename
+try:
+    from pyedt import edt
+except ModuleNotFoundError:
+    from edt import edt
 
 
 def dict_to_vtk(data, filename, voxel_size=1, origin=(0, 0, 0)):
@@ -159,10 +160,10 @@ def to_palabos(im, filename, solid=0):
     bin_im = bin_im.astype(int)
     # Distance Transform computes Euclidean distance in lattice units to
     # Nearest fluid for every solid voxel
-    dt = nd.distance_transform_edt(bin_im)
-    dt[dt > np.sqrt(2)] = 2
-    dt[(dt > 0) * (dt <= np.sqrt(2))] = 1
-    dt = dt.astype(int)
+    dt = edt(bin_im)
+    dt[dt > 2] = 2
+    dt[(dt > 0) * (dt <= 2)] = 1
+    dt = np.sqrt(dt).astype(int)
     # Write out data
     with open(filename, "w") as f:
         out_data = dt.flatten().tolist()

src/porespy/metrics/_funcs.py

@@ -295,8 +295,7 @@ def radial_density_distribution(dt, bins=10, log=False, voxel_size=1):
     Parameters
     ----------
     dt : ndarray
-        A distance transform of the pore space (the ``edt`` package is
-        recommended).  Note that it is recommended to apply
+        A distance transform of the pore space.  Note that it is recommended to apply
         ``find_dt_artifacts`` to this image first, and set potentially
         erroneous values to 0 with ``dt[mask] = 0`` where
         ``mask = porespy.filters.find_dt_artifaces(dt)``.

src/porespy/metrics/_meshtools.py

@@ -22,7 +22,7 @@
 logger = logging.getLogger(__name__)
 
 
-def region_volumes(regions, mode='marching_cubes'):
+def region_volumes(regions, mode='marching_cubes', voxel_size=(1, 1, 1)):
     r"""
     Compute volume of each labelled region in an image
 
@@ -61,20 +61,24 @@ def region_volumes(regions, mode='marching_cubes'):
     for i, s in enumerate(tqdm(slices, desc=msg, **settings.tqdm)):
         region = regions[s] == (i + 1)
         if mode == 'marching_cubes':
-            vols[i] = mesh_volume(region)
+            vols[i] = mesh_volume(region, voxel_size=voxel_size)
         elif mode.startswith('voxel'):
             vols[i] = region.sum(dtype=np.int64)
     return vols
 
 
-def mesh_volume(region):
+def mesh_volume(region, voxel_size=(1., 1., 1.)):
     r"""
     Compute the volume of a single region by meshing it
 
     Parameters
     ----------
     region : ndarray
         An image with a single region labelled as ``True`` (or > 0)
+    voxel_size : tuple
+        The resolution of the image, expressed as the length of one side of a
+        voxel, so the volume of a voxel would be **voxel_size**-cubed.  The
+        default is 1.
 
     Returns
     -------
@@ -95,7 +99,8 @@ def mesh_volume(region):
     except ModuleNotFoundError:
         msg = 'The trimesh package can be installed with pip install trimesh'
         raise ModuleNotFoundError(msg)
-    mc = mesh_region(region > 0)
+
+    mc = mesh_region(region > 0, voxel_size=voxel_size)
     m = Trimesh(vertices=mc.verts, faces=mc.faces, vertex_normals=mc.norm)
     if m.is_watertight:
         vol = np.abs(m.volume)
@@ -104,7 +109,7 @@ def mesh_volume(region):
     return vol
 
 
-def region_surface_areas(regions, voxel_size=1, strel=None):
+def region_surface_areas(regions, voxel_size=(1, 1, 1), strel=None):
     r"""
     Extract the surface area of each region in a labeled image.
 
@@ -117,7 +122,7 @@ def region_surface_areas(regions, voxel_size=1, strel=None):
         An image of the pore space partitioned into individual pore regions.
         Note that zeros in the image will not be considered for area
         calculation.
-    voxel_size : scalar
+    voxel_size : tuple
         The resolution of the image, expressed as the length of one side of a
         voxel, so the volume of a voxel would be **voxel_size**-cubed.  The
         default is 1.
@@ -157,9 +162,11 @@ def region_surface_areas(regions, voxel_size=1, strel=None):
             s = extend_slice(slices[reg], im.shape)
             sub_im = im[s]
             mask_im = sub_im == i
-            mesh = mesh_region(region=mask_im, strel=strel)
+            mesh = mesh_region(region=mask_im,
+                               strel=strel,
+                               voxel_size=voxel_size)
             sa[reg] = mesh_surface_area(mesh)
-    result = sa * voxel_size**2
+    result = sa
     return result
 
 
@@ -291,13 +298,15 @@ def region_interface_areas(regions, areas, voxel_size=1, strel=None):
                                            slices[j][1].stop)]
                     merged_region = ((merged_region == reg + 1)
                                      + (merged_region == j + 1))
-                    mesh = mesh_region(region=merged_region, strel=strel)
+                    mesh = mesh_region(region=merged_region,
+                                       strel=strel,
+                                       voxel_size=voxel_size)
                     sa_combined.append(mesh_surface_area(mesh))
     # Interfacial area calculation
     cn = np.array(cn)
     ia = 0.5 * (sa[cn[:, 0]] + sa[cn[:, 1]] - sa_combined)
     ia[ia <= 0] = 1
     result = Results()
     result.conns = cn
-    result.area = ia * voxel_size**2
+    result.area = ia
     return result