Merge branch 'dev'

.travis.yml

@@ -1,13 +1,23 @@
 sudo: false
 
+os:
+  - linux
+
+dist:
+  - xenial
+
 language: python
 
 python:
   - "3.6"
+  - "3.7"
+
+cache: pip
+
+services:
+  - xvfb
 
 before_install:
-  - export DISPLAY=:99.0
-  - sh -e /etc/init.d/xvfb start
   - wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh
   - bash miniconda.sh -b -p $HOME/miniconda
   - export PATH="$HOME/miniconda/bin:$PATH"

porespy/__init__.py

@@ -79,7 +79,7 @@
 
 '''
 
-__version__ = "1.1.2"
+__version__ = "1.2.0"
 
 from . import tools
 from . import filters

porespy/filters/__funcs__.py

@@ -1,6 +1,7 @@
 from collections import namedtuple
 import scipy as sp
 import numpy as np
+import operator as op
 import scipy.ndimage as spim
 import scipy.spatial as sptl
 import warnings
@@ -16,6 +17,66 @@
 from porespy.tools import _create_alias_map
 
 
+def trim_small_clusters(im, size=1):
+    r"""
+    Remove isolated voxels or clusters smaller than a given size
+
+    Parameters
+    ----------
+    im : ND-array
+        The binary image from which voxels are to be removed
+    size : scalar
+        The threshold size of clusters to trim.  As clusters with this many
+        voxels or fewer will be trimmed.  The default is 1 so only single
+        voxels are removed.
+
+    Returns
+    -------
+    im : ND-image
+        A copy of ``im`` with clusters of voxels smaller than the given
+        ``size`` removed.
+
+    """
+    if im.dims == 2:
+        strel = disk(1)
+    elif im.ndims == 3:
+        strel = ball(1)
+    else:
+        raise Exception('Only 2D or 3D images are accepted')
+    filtered_array = sp.copy(im)
+    labels, N = spim.label(filtered_array, structure=strel)
+    id_sizes = sp.array(spim.sum(im, labels, range(N + 1)))
+    area_mask = (id_sizes <= size)
+    filtered_array[area_mask[labels]] = 0
+    return filtered_array
+
+
+def hold_peaks(im, axis=-1):
+    r"""
+    Replaces each voxel with the highest value along the given axis
+
+    Parameters
+    ----------
+    im : ND-image
+        A greyscale image whose peaks are to
+    """
+    A = im
+    B = np.swapaxes(A, axis, -1)
+    updown = np.empty((*B.shape[:-1], B.shape[-1]+1), B.dtype)
+    updown[..., 0], updown[..., -1] = -1, -1
+    np.subtract(B[..., 1:], B[..., :-1], out=updown[..., 1:-1])
+    chnidx = np.where(updown)
+    chng = updown[chnidx]
+    pkidx, = np.where((chng[:-1] > 0) & (chng[1:] < 0) | (chnidx[-1][:-1] == 0))
+    pkidx = (*map(op.itemgetter(pkidx), chnidx),)
+    out = np.zeros_like(A)
+    aux = out.swapaxes(axis, -1)
+    aux[(*map(op.itemgetter(slice(1, None)), pkidx),)] = np.diff(B[pkidx])
+    aux[..., 0] = B[..., 0]
+    result = out.cumsum(axis=axis)
+    return result
+
+
 def distance_transform_lin(im, axis=0, mode='both'):
     r"""
     Replaces each void voxel with the linear distance to the nearest solid
@@ -985,10 +1046,11 @@ def apply_chords_3D(im, spacing=0, trim_edges=True):
 
 def local_thickness(im, sizes=25, mode='hybrid'):
     r"""
-    For each voxel, this functions calculates the radius of the largest sphere
-    that both engulfs the voxel and fits entirely within the foreground. This
-    is not the same as a simple distance transform, which finds the largest
-    sphere that could be *centered* on each voxel.
+    For each voxel, this function calculates the radius of the largest sphere
+    that both engulfs the voxel and fits entirely within the foreground.
+
+    This is not the same as a simple distance transform, which finds the
+    largest sphere that could be *centered* on each voxel.
 
     Parameters
     ----------
@@ -1111,7 +1173,7 @@ def porosimetry(im, sizes=25, inlets=None, access_limited=True,
 
     Notes
     -----
-    There are many ways to perform this filter, and PoreSpy offer 3, which
+    There are many ways to perform this filter, and PoreSpy offers 3, which
     users can choose between via the ``mode`` argument.  These methods all
     work in a similar way by finding which foreground voxels can accomodate
     a sphere of a given radius, then repeating for smaller radii.
@@ -1190,7 +1252,7 @@ def trim_disconnected_blobs(im, inlets):
     ----------
     im : ND-array
         The array to be trimmed
-    inlets : ND-array of tuple of indices
+    inlets : ND-array or tuple of indices
         The locations of the inlets.  Any voxels *not* connected directly to
         the inlets will be trimmed
 
@@ -1200,11 +1262,14 @@ def trim_disconnected_blobs(im, inlets):
         An array of the same shape as ``im``, but with all foreground
         voxels not connected to the ``inlets`` removed.
     """
-    temp = sp.zeros_like(im)
-    temp[inlets] = True
-    labels, N = spim.label(im + temp)
-    im = im ^ (clear_border(labels=labels) > 0)
-    return im
+    labels = spim.label(im)[0]
+    keep = sp.unique(labels[inlets])
+    keep = keep[keep > 0]
+    if len(keep) > 0:
+        im2 = sp.reshape(sp.in1d(labels, keep), newshape=im.shape)
+    else:
+        im2 = sp.zeros_like(im)
+    return im2
 
 
 def _get_axial_shifts(ndim=2, include_diagonals=False):
@@ -1319,3 +1384,66 @@ def nphase_border(im, include_diagonals=False):
         return out[1:-1, 1:-1].copy()
     else:
         return out[1:-1, 1:-1, 1:-1].copy()
+
+
+def prune_branches(skel, branch_points=None, iterations=1):
+    r"""
+    Removes all dangling ends or tails of a skeleton.
+
+    Parameters
+    ----------
+    skel : ND-image
+        A image of a full or partial skeleton from which the tails should be
+        trimmed.
+
+    branch_points : ND-image, optional
+        An image the same size ``skel`` with True values indicating the branch
+        points of the skeleton.  If this is not provided it is calculated
+        automatically.
+
+    Returns
+    -------
+    An ND-image containing the skeleton with tails removed.
+
+    """
+    skel = skel > 0
+    if skel.ndim == 2:
+        from skimage.morphology import square as cube
+    else:
+        from skimage.morphology import cube
+    # Create empty image to house results
+    im_result = sp.zeros_like(skel)
+    # If branch points are not supplied, attempt to find them
+    if branch_points is None:
+        branch_points = spim.convolve(skel*1.0, weights=cube(3)) > 3
+        branch_points = branch_points*skel
+    # Store original branch points before dilating
+    pts_orig = branch_points
+    # Find arcs of skeleton by deleting branch points
+    arcs = skel*(~branch_points)
+    # Label arcs
+    arc_labels = spim.label(arcs, structure=cube(3))[0]
+    # Dilate branch points so they overlap with the arcs
+    branch_points = spim.binary_dilation(branch_points, structure=cube(3))
+    pts_labels = spim.label(branch_points, structure=cube(3))[0]
+    # Now scan through each arc to see if it's connected to two branch points
+    slices = spim.find_objects(arc_labels)
+    label_num = 0
+    for s in slices:
+        label_num += 1
+        # Find branch point labels the overlap current arc
+        hits = pts_labels[s]*(arc_labels[s] == label_num)
+        # If image contains 2 branch points, then it's not a tail.
+        if len(sp.unique(hits)) == 3:
+            im_result[s] += arc_labels[s] == label_num
+    # Add missing branch points back to arc image to make complete skeleton
+    im_result += skel*pts_orig
+    if iterations > 1:
+        iterations -= 1
+        im_temp = sp.copy(im_result)
+        im_result = prune_branches(skel=im_result,
+                                   branch_points=None,
+                                   iterations=iterations)
+        if sp.all(im_temp == im_result):
+            iterations = 0
+    return im_result

porespy/filters/__init__.py

@@ -16,21 +16,26 @@
     porespy.filters.apply_chords
     porespy.filters.apply_chords_3D
     porespy.filters.distance_transform_lin
+    porespy.filters.fftmorphology
     porespy.filters.fill_blind_pores
     porespy.filters.find_disconnected_voxels
     porespy.filters.find_dt_artifacts
     porespy.filters.find_peaks
     porespy.filters.flood
     porespy.filters.local_thickness
     porespy.filters.porosimetry
+    porespy.filters.prune_branches
+    porespy.filters.reduce_peaks
     porespy.filters.region_size
     porespy.filters.snow_partitioning
+    porespy.filters.snow_partitioning_n
     porespy.filters.trim_disconnected_blobs
     porespy.filters.trim_extrema
     porespy.filters.trim_floating_solid
     porespy.filters.trim_nearby_peaks
     porespy.filters.trim_nonpercolating_paths
     porespy.filters.trim_saddle_points
+    porespy.filters.trim_small_cluster
 
 
 .. autofunction:: apply_chords
@@ -42,21 +47,25 @@
 .. autofunction:: find_dt_artifacts
 .. autofunction:: find_peaks
 .. autofunction:: flood
+.. autofunction:: hold_peaks
 .. autofunction:: local_thickness
 .. autofunction:: nphase_border
 .. autofunction:: porosimetry
+.. autofunction:: prune_branches
+.. autofunction:: reduce_peaks
 .. autofunction:: region_size
 .. autofunction:: snow_partitioning
 .. autofunction:: snow_partitioning_n
-.. autofunction:: trim_extrema
-.. autofunction:: reduce_peaks
 .. autofunction:: trim_disconnected_blobs
+.. autofunction:: trim_extrema
 .. autofunction:: trim_floating_solid
 .. autofunction:: trim_nearby_peaks
 .. autofunction:: trim_nonpercolating_paths
 .. autofunction:: trim_saddle_points
+.. autofunction:: trim_small_clusters
 
 """
+
 from .__funcs__ import apply_chords
 from .__funcs__ import apply_chords_3D
 from .__funcs__ import distance_transform_lin
@@ -66,9 +75,11 @@
 from .__funcs__ import find_dt_artifacts
 from .__funcs__ import find_peaks
 from .__funcs__ import flood
+from .__funcs__ import hold_peaks
 from .__funcs__ import local_thickness
 from .__funcs__ import nphase_border
 from .__funcs__ import porosimetry
+from .__funcs__ import prune_branches
 from .__funcs__ import reduce_peaks
 from .__funcs__ import region_size
 from .__funcs__ import snow_partitioning
@@ -79,3 +90,4 @@
 from .__funcs__ import trim_nonpercolating_paths
 from .__funcs__ import trim_nearby_peaks
 from .__funcs__ import trim_saddle_points
+from .__funcs__ import trim_small_clusters

porespy/generators/__imgen__.py

@@ -704,7 +704,7 @@ def blobs(shape: List[int], porosity: float = 0.5, blobiness: int = 1):
 
 
 def cylinders(shape: List[int], radius: int, ncylinders: int,
-              phi_max: float = 0, theta_max: float = 90):
+              phi_max: float = 0, theta_max: float = 90, length: float = None):
     r"""
     Generates a binary image of overlapping cylinders.  This is a good
     approximation of a fibrous mat.
@@ -722,14 +722,21 @@ def cylinders(shape: List[int], radius: int, ncylinders: int,
         cylinders overlap and intersect different fractions of the domain.
     theta_max : scalar
         A value between 0 and 90 that controls the amount of rotation *in the*
-        XY plane, with 0 meaning all fibers point in the X-direction, and
+        XY plane, with 0 meaning all cylinders point in the X-direction, and
         90 meaning they are randomly rotated about the Z axis by as much
         as +/- 90 degrees.
     phi_max : scalar
-        A value between 0 and 90 that controls the amount that the fibers
-        lie *out of* the XY plane, with 0 meaning all fibers lie in the XY
-        plane, and 90 meaning that fibers are randomly oriented out of the
+        A value between 0 and 90 that controls the amount that the cylinders
+        lie *out of* the XY plane, with 0 meaning all cylinders lie in the XY
+        plane, and 90 meaning that cylinders are randomly oriented out of the
         plane by as much as +/- 90 degrees.
+    length : scalar
+        The length of the cylinders to add.  If ``None`` (default) then the
+        cylinders will extend beyond the domain in both directions so no ends
+        will exist. If a scalar value is given it will be interpreted as the
+        Euclidean distance between the two ends of the cylinder.  Note that
+        one or both of the ends *may* still lie outside the domain, depending
+        on the randomly chosen center point of the cylinder.
 
     Returns
     -------
@@ -741,7 +748,10 @@ def cylinders(shape: List[int], radius: int, ncylinders: int,
         shape = sp.full((3, ), int(shape))
     elif sp.size(shape) == 2:
         raise Exception("2D cylinders don't make sense")
-    R = sp.sqrt(sp.sum(sp.square(shape))).astype(int)
+    if length is None:
+        R = sp.sqrt(sp.sum(sp.square(shape))).astype(int)
+    else:
+        R = length/2
     im = sp.zeros(shape)
     # Adjust max angles to be between 0 and 90
     if (phi_max > 90) or (phi_max < 0):

porespy/metrics/__funcs__.py

@@ -197,12 +197,12 @@ def porosity(im):
     solid phase.
 
     All other values are ignored, so this can also return the relative
-    fraction of a phase of interest.
+    fraction of a phase of interest in trinary or multiphase images.
 
     Parameters
     ----------
     im : ND-array
-        Image of the void space with 1's indicating void space (or True) and
+        Image of the void space with 1's indicating void phase (or True) and
         0's indicating the solid phase (or False).
 
     Returns
@@ -387,7 +387,10 @@ def pore_size_distribution(im, bins=10, log=True, voxel_size=1):
         be automatically generated that span the data range.
     log : boolean
         If ``True`` (default) the size data is converted to log (base-10)
-        values before processing.  This can help
+        values before processing.  This can help to plot wide size
+        distributions or to better visualize the in the small size region.
+        Note that you can anti-log the radii values in the retunred ``tuple``,
+        but the binning is performed on the logged radii values.
     voxel_size : scalar
         The size of a voxel side in preferred units.  The default is 1, so the
         user can apply the scaling to the returned results after the fact.
@@ -416,8 +419,8 @@ def pore_size_distribution(im, bins=10, log=True, voxel_size=1):
 
     Notes
     -----
-    (1) To ensure the returned values represent actual sizes be sure to scale
-    the distance transform by the voxel size first (``dt *= voxel_size``)
+    (1) To ensure the returned values represent actual sizes you can manually
+    scale the input image by the voxel size first (``im *= voxel_size``)
 
     plt.bar(psd.R, psd.satn, width=psd.bin_widths, edgecolor='k')
 

porespy/networks/__getnet__.py

@@ -81,7 +81,7 @@ def regions_to_network(im, dt=None, voxel_size=1):
         p_label[pore] = i
         p_coords[pore, :] = spim.center_of_mass(pore_im) + s_offset
         p_volume[pore] = sp.sum(pore_im)
-        p_dia_local[pore] = 2*sp.amax(pore_dt)
+        p_dia_local[pore] = (2*sp.amax(pore_dt)) - sp.sqrt(3)
         p_dia_global[pore] = 2*sp.amax(sub_dt)
         p_area_surf[pore] = sp.sum(pore_dt == 1)
         im_w_throats = spim.binary_dilation(input=pore_im, structure=struc_elem(1))