Commit ⋅ PMEAL/PoreSpy

from collections import namedtuple
import sys
import dask
import dask.array as da
from dask.diagnostics import ProgressBar
import warnings
import numpy as np
from numba import njit, prange
from edt import edt
import operator as op
from tqdm import tqdm
import scipy.ndimage as spim
import scipy.spatial as sptl
import warnings
from collections import namedtuple
from scipy.signal import fftconvolve
from tqdm import tqdm
from numba import jit
from skimage.segmentation import clear_border
from skimage.segmentation import clear_border, watershed
from skimage.morphology import ball, disk, square, cube, diamond, octahedron
from skimage.morphology import reconstruction, watershed
from skimage.morphology import reconstruction
from porespy.tools import randomize_colors, fftmorphology
from porespy.tools import get_border, extend_slice, extract_subsection
from porespy.tools import ps_disk, ps_ball
from porespy.tools import _create_alias_map
from porespy.tools import ps_disk, ps_ball


def apply_padded(im, pad_width, func, pad_val=1, **kwargs):
    r"""
    Applies padding to an image before sending to ``func``, then extracts the
    result corresponding to the original image shape.

    Parameters
    ----------
    im : ND-image
        The image to which ``func`` should be applied
    pad_width : int or list of ints
        The amount of padding to apply to each axis.  Refer to ``numpy.pad``
        documentation for more details.
    pad_val : scalar
        The value to place into the padded voxels.  The default is 1 (or
        ``True``) which extends the pore space.
    func : function handle
        The function to apply to the padded image
    kwargs : additional keyword arguments
        All additional keyword arguments are collected and passed to ``func``.

    Notes
    -----
    A use case for this is when using ``skimage.morphology.skeletonize_3d``
    to ensure that the skeleton extends beyond the edges of the image.
    """
    padded = np.pad(im, pad_width=pad_width,
                    mode='constant', constant_values=pad_val)
    temp = func(padded, **kwargs)
    result = extract_subsection(im=temp, shape=im.shape)
    return result


def trim_small_clusters(im, size=1):
    r"""
    Remove isolated voxels or clusters smaller than a given size
    Remove isolated voxels or clusters of a given size or smaller

    Parameters
    ----------

@@ -36,16 +73,16 @@

        ``size`` removed.

    """
    if im.dims == 2:
    if im.ndim == 2:
        strel = disk(1)
    elif im.ndims == 3:
    elif im.ndim == 3:
        strel = ball(1)
    else:
        raise Exception('Only 2D or 3D images are accepted')
        raise Exception("Only 2D or 3D images are accepted")
    filtered_array = np.copy(im)
    labels, N = spim.label(filtered_array, structure=strel)
    id_sizes = np.array(spim.sum(im, labels, range(N + 1)))
    area_mask = (id_sizes <= size)
    area_mask = id_sizes <= size
    filtered_array[area_mask[labels]] = 0
    return filtered_array


@@ -61,12 +98,12 @@

61	98		"""
62	99		A = im
63	100		B = np.swapaxes(A, axis, -1)
64		-	updown = np.empty((*B.shape[:-1], B.shape[-1]+1), B.dtype)
	101	+	updown = np.empty((*B.shape[:-1], B.shape[-1] + 1), B.dtype)
65	102		updown[..., 0], updown[..., -1] = -1, -1
66	103		np.subtract(B[..., 1:], B[..., :-1], out=updown[..., 1:-1])
67	104		chnidx = np.where(updown)
68	105		chng = updown[chnidx]
69		-	pkidx, = np.where((chng[:-1] > 0) & (chng[1:] < 0) \| (chnidx[-1][:-1] == 0))
	106	+	(pkidx,) = np.where((chng[:-1] > 0) & (chng[1:] < 0) \| (chnidx[-1][:-1] == 0))
70	107		pkidx = (*map(op.itemgetter(pkidx), chnidx),)
71	108		out = np.zeros_like(A)
72	109		aux = out.swapaxes(axis, -1)

@@ -76,7 +113,7 @@

    return result


def distance_transform_lin(im, axis=0, mode='both'):
def distance_transform_lin(im, axis=0, mode="both"):
    r"""
    Replaces each void voxel with the linear distance to the nearest solid
    voxel along the specified axis.

@@ -110,33 +147,38 @@

        the nearest background along the specified axis.
    """
    if im.ndim != im.squeeze().ndim:
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
                      ' unexpected behavior.')
    if mode in ['backward', 'reverse']:
        warnings.warn(
            "Input image conains a singleton axis:"
            + str(im.shape)
            + " Reduce dimensionality with np.squeeze(im) to avoid"
            + " unexpected behavior."
        )
    if mode in ["backward", "reverse"]:
        im = np.flip(im, axis)
        im = distance_transform_lin(im=im, axis=axis, mode='forward')
        im = distance_transform_lin(im=im, axis=axis, mode="forward")
        im = np.flip(im, axis)
        return im
    elif mode in ['both']:
        im_f = distance_transform_lin(im=im, axis=axis, mode='forward')
        im_b = distance_transform_lin(im=im, axis=axis, mode='backward')
    elif mode in ["both"]:
        im_f = distance_transform_lin(im=im, axis=axis, mode="forward")
        im_b = distance_transform_lin(im=im, axis=axis, mode="backward")
        return np.minimum(im_f, im_b)
    else:
        b = np.cumsum(im > 0, axis=axis)
        c = np.diff(b*(im == 0), axis=axis)
        c = np.diff(b * (im == 0), axis=axis)
        d = np.minimum.accumulate(c, axis=axis)
        if im.ndim == 1:
            e = np.pad(d, pad_width=[1, 0], mode='constant', constant_values=0)
            e = np.pad(d, pad_width=[1, 0], mode="constant", constant_values=0)
        elif im.ndim == 2:
            ax = [[[1, 0], [0, 0]], [[0, 0], [1, 0]]]
            e = np.pad(d, pad_width=ax[axis], mode='constant', constant_values=0)
            e = np.pad(d, pad_width=ax[axis], mode="constant", constant_values=0)
        elif im.ndim == 3:
            ax = [[[1, 0], [0, 0], [0, 0]],
                  [[0, 0], [1, 0], [0, 0]],
                  [[0, 0], [0, 0], [1, 0]]]
            e = np.pad(d, pad_width=ax[axis], mode='constant', constant_values=0)
        f = im*(b + e)
            ax = [
                [[1, 0], [0, 0], [0, 0]],
                [[0, 0], [1, 0], [0, 0]],
                [[0, 0], [0, 0], [1, 0]],
            ]
            e = np.pad(d, pad_width=ax[axis], mode="constant", constant_values=0)
        f = im * (b + e)
        return f



@@ -205,36 +247,36 @@

    using marker-based watershed segmenation".  Physical Review E. (2017)

    """
    tup = namedtuple('results', field_names=['im', 'dt', 'peaks', 'regions'])
    print('_'*60)
    tup = namedtuple("results", field_names=["im", "dt", "peaks", "regions"])
    print("-" * 60)
    print("Beginning SNOW Algorithm")
    im_shape = np.array(im.shape)
    if im.dtype is not bool:
        print('Converting supplied image (im) to boolean')
        print("Converting supplied image (im) to boolean")
        im = im > 0
    if dt is None:
        print('Peforming Distance Transform')
        print("Peforming Distance Transform")
        if np.any(im_shape == 1):
            ax = np.where(im_shape == 1)[0][0]
            dt = spim.distance_transform_edt(input=im.squeeze())
            dt = edt(im.squeeze())
            dt = np.expand_dims(dt, ax)
        else:
            dt = spim.distance_transform_edt(input=im)
            dt = edt(im)

    tup.im = im
    tup.dt = dt

    if sigma > 0:
        print('Applying Gaussian blur with sigma =', str(sigma))
        print("Applying Gaussian blur with sigma =", str(sigma))
        dt = spim.gaussian_filter(input=dt, sigma=sigma)

    peaks = find_peaks(dt=dt, r_max=r_max)
    print('Initial number of peaks: ', spim.label(peaks)[1])
    print("Initial number of peaks: ", spim.label(peaks)[1])
    peaks = trim_saddle_points(peaks=peaks, dt=dt, max_iters=500)
    print('Peaks after trimming saddle points: ', spim.label(peaks)[1])
    print("Peaks after trimming saddle points: ", spim.label(peaks)[1])
    peaks = trim_nearby_peaks(peaks=peaks, dt=dt)
    peaks, N = spim.label(peaks)
    print('Peaks after trimming nearby peaks: ', N)
    print("Peaks after trimming nearby peaks: ", N)
    tup.peaks = peaks
    if mask:
        mask_solid = im > 0

@@ -333,30 +375,32 @@

    combined_dt = 0
    combined_region = 0
    num = [0]
    for i in phases_num:
        print('_' * 60)
    for i, j in enumerate(phases_num):
        print("_" * 60)
        if alias is None:
            print('Processing Phase {}'.format(i))
        else:
            print('Processing Phase {}'.format(al[i]))
        phase_snow = snow_partitioning(im == i,
                                       dt=None, r_max=r_max, sigma=sigma,
                                       return_all=return_all, mask=mask,
                                       randomize=randomize)
        if len(phases_num) == 1 and phases_num == 1:
            combined_dt = phase_snow.dt
            combined_region = phase_snow.regions
            print("Processing Phase {}".format(j))
        else:
            combined_dt += phase_snow.dt
            phase_snow.regions *= phase_snow.im
            phase_snow.regions += num[i - 1]
            phase_ws = phase_snow.regions * phase_snow.im
            phase_ws[phase_ws == num[i - 1]] = 0
            combined_region += phase_ws
            print("Processing Phase {}".format(al[j]))
        phase_snow = snow_partitioning(
            im == j,
            dt=None,
            r_max=r_max,
            sigma=sigma,
            return_all=return_all,
            mask=mask,
            randomize=randomize,
        )
        combined_dt += phase_snow.dt
        phase_snow.regions *= phase_snow.im
        phase_snow.regions += num[i]
        phase_ws = phase_snow.regions * phase_snow.im
        phase_ws[phase_ws == num[i]] = 0
        combined_region += phase_ws
        num.append(np.amax(combined_region))
    if return_all:
        tup = namedtuple('results', field_names=['im', 'dt', 'phase_max_label',
                                                 'regions'])
        tup = namedtuple(
            "results", field_names=["im", "dt", "phase_max_label", "regions"]
        )
        tup.im = im
        tup.dt = combined_dt
        tup.phase_max_label = num[1:]

@@ -366,24 +410,22 @@

        return combined_region


def find_peaks(dt, r_max=4, footprint=None):
def find_peaks(dt, r_max=4, footprint=None, **kwargs):
    r"""
    Returns all local maxima in the distance transform
    Finds local maxima in the distance transform

    Parameters
    ----------
    dt : ND-array
        The distance transform of the pore space.  This may be calculated and
        filtered using any means desired.

    r_max : scalar
        The size of the structuring element used in the maximum filter.  This
        controls the localness of any maxima. The default is 4 voxels.

    footprint : ND-array
        Specifies the shape of the structuring element used to define the
        neighborhood when looking for peaks.  If none is specified then a
        spherical shape is used (or circular in 2D).
        neighborhood when looking for peaks.  If ``None`` (the default) is
        specified then a spherical shape is used (or circular in 2D).

    Returns
    -------

@@ -399,14 +441,15 @@

    ``peaks = peak_local_max(image=dt, min_distance=r, exclude_border=0,
    indices=False)``

    This automatically uses a square structuring element which is significantly
    faster than using a circular or spherical element.
    The *skimage* function automatically uses a square structuring element
    which is significantly faster than using a circular or spherical element.
    """
    im = dt > 0
    if im.ndim != im.squeeze().ndim:
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
                      ' unexpected behavior.')
        warnings.warn("Input image conains a singleton axis:"
                      + str(im.shape)
                      + " Reduce dimensionality with np.squeeze(im) to avoid"
                      + " unexpected behavior.")
    if footprint is None:
        if im.ndim == 2:
            footprint = disk

@@ -414,8 +457,17 @@

            footprint = ball
        else:
            raise Exception("only 2-d and 3-d images are supported")
    mx = spim.maximum_filter(dt + 2*(~im), footprint=footprint(r_max))
    peaks = (dt == mx)*im
    parallel = kwargs.pop('parallel', False)
    cores = kwargs.pop('cores', None)
    divs = kwargs.pop('cores', 2)
    if parallel:
        overlap = max(footprint(r_max).shape)
        peaks = chunked_func(func=find_peaks, overlap=overlap,
                             im_arg='dt', dt=dt, footprint=footprint,
                             cores=cores, divs=divs)
    else:
        mx = spim.maximum_filter(dt + 2 * (~im), footprint=footprint(r_max))
        peaks = (dt == mx) * im
    return peaks



@@ -427,14 +479,14 @@

    Parameters
    ----------
    peaks : ND-image
        An image containing True values indicating peaks in the distance
        An image containing ``True`` values indicating peaks in the distance
        transform

    Returns
    -------
    image : ND-array
        An array with the same number of isolated peaks as the original image,
        but fewer total voxels.
        but fewer total ``True`` voxels.

    Notes
    -----

@@ -447,9 +499,9 @@

    else:
        strel = cube
    markers, N = spim.label(input=peaks, structure=strel(3))
    inds = spim.measurements.center_of_mass(input=peaks,
                                            labels=markers,
                                            index=np.arange(1, N+1))
    inds = spim.measurements.center_of_mass(
        input=peaks, labels=markers, index=np.arange(1, N + 1)
    )
    inds = np.floor(inds).astype(int)
    # Centroid may not be on old pixel, so create a new peaks image
    peaks_new = np.zeros_like(peaks, dtype=bool)

@@ -457,7 +509,7 @@

    return peaks_new


def trim_saddle_points(peaks, dt, max_iters=10):
def trim_saddle_points(peaks, dt, max_iters=10, verbose=1):
    r"""
    Removes peaks that were mistakenly identified because they lied on a
    saddle or ridge in the distance transform that was not actually a true

@@ -499,26 +551,27 @@

    slices = spim.find_objects(labels)
    for i in range(N):
        s = extend_slice(s=slices[i], shape=peaks.shape, pad=10)
        peaks_i = labels[s] == i+1
        peaks_i = labels[s] == i + 1
        dt_i = dt[s]
        im_i = dt_i > 0
        iters = 0
        peaks_dil = np.copy(peaks_i)
        while iters < max_iters:
            iters += 1
            peaks_dil = spim.binary_dilation(input=peaks_dil,
                                             structure=cube(3))
            peaks_max = peaks_dil*np.amax(dt_i*peaks_dil)
            peaks_extended = (peaks_max == dt_i)*im_i
            peaks_dil = spim.binary_dilation(input=peaks_dil, structure=cube(3))
            peaks_max = peaks_dil * np.amax(dt_i * peaks_dil)
            peaks_extended = (peaks_max == dt_i) * im_i
            if np.all(peaks_extended == peaks_i):
                break  # Found a true peak
            elif np.sum(peaks_extended*peaks_i) == 0:
            elif np.sum(peaks_extended * peaks_i) == 0:
                peaks_i = False
                break  # Found a saddle point
        peaks[s] = peaks_i
        if iters >= max_iters:
            print('Maximum number of iterations reached, consider'
                  + 'running again with a larger value of max_iters')
        if iters >= max_iters and verbose:
            print(
                "Maximum number of iterations reached, consider "
                + "running again with a larger value of max_iters"
            )
    return peaks



@@ -561,7 +614,7 @@

        from skimage.morphology import cube
    peaks, N = spim.label(peaks, structure=cube(3))
    crds = spim.measurements.center_of_mass(peaks, labels=peaks,
                                            index=np.arange(1, N+1))
                                            index=np.arange(1, N + 1))
    crds = np.vstack(crds).astype(int)  # Convert to numpy array of ints
    # Get distance between each peak as a distance map
    tree = sptl.cKDTree(data=crds)

@@ -583,7 +636,7 @@

    slices = spim.find_objects(input=peaks)
    for s in drop_peaks:
        peaks[slices[s]] = 0
    return (peaks > 0)
    return peaks > 0


def find_disconnected_voxels(im, conn=None):

@@ -597,11 +650,10 @@

    im : ND-image
        A Boolean image, with True values indicating the phase for which
        disconnected voxels are sought.

    conn : int
        For 2D the options are 4 and 8 for square and diagonal neighbors, while
        for the 3D the options are 6 and 26, similarily for square and diagonal
        neighbors.  The default is max
        neighbors.  The default is the maximum option.

    Returns
    -------

@@ -618,25 +670,32 @@


    """
    if im.ndim != im.squeeze().ndim:
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
                      ' unexpected behavior.')
        warnings.warn(
            "Input image conains a singleton axis:"
            + str(im.shape)
            + " Reduce dimensionality with np.squeeze(im) to avoid"
            + " unexpected behavior."
        )
    if im.ndim == 2:
        if conn == 4:
            strel = disk(1)
        elif conn in [None, 8]:
            strel = square(3)
        else:
            raise Exception("Received conn is not valid")
    elif im.ndim == 3:
        if conn == 6:
            strel = ball(1)
        elif conn in [None, 26]:
            strel = cube(3)
        else:
            raise Exception("Received conn is not valid")
    labels, N = spim.label(input=im, structure=strel)
    holes = clear_border(labels=labels) > 0
    return holes


def fill_blind_pores(im):
def fill_blind_pores(im, conn=None):
    r"""
    Fills all pores that are not connected to the edges of the image.


@@ -649,6 +708,10 @@

    -------
    image : ND-array
        A version of ``im`` but with all the disconnected pores removed.
    conn : int
        For 2D the options are 4 and 8 for square and diagonal neighbors, while
        for the 3D the options are 6 and 26, similarily for square and diagonal
        neighbors.  The default is the maximum option.

    See Also
    --------

@@ -656,12 +719,12 @@


    """
    im = np.copy(im)
    holes = find_disconnected_voxels(im)
    holes = find_disconnected_voxels(im, conn=conn)
    im[holes] = False
    return im


def trim_floating_solid(im):
def trim_floating_solid(im, conn=None):
    r"""
    Removes all solid that that is not attached to the edges of the image.


@@ -669,6 +732,10 @@

    ----------
    im : ND-array
        The image of the porous material
    conn : int
        For 2D the options are 4 and 8 for square and diagonal neighbors, while
        for the 3D the options are 6 and 26, similarily for square and diagonal
        neighbors.  The default is the maximum option.

    Returns
    -------

@@ -681,12 +748,13 @@


    """
    im = np.copy(im)
    holes = find_disconnected_voxels(~im)
    holes = find_disconnected_voxels(~im, conn=conn)
    im[holes] = True
    return im


def trim_nonpercolating_paths(im, inlet_axis=0, outlet_axis=0):
def trim_nonpercolating_paths(im, inlet_axis=0, outlet_axis=0,
                              inlets=None, outlets=None):
    r"""
    Removes all nonpercolating paths between specified edges


@@ -699,16 +767,20 @@

    im : ND-array
        The image of the porous material with ```True`` values indicating the
        phase of interest

    inlet_axis : int
        Inlet axis of boundary condition. For three dimensional image the
        number ranges from 0 to 2. For two dimensional image the range is
        between 0 to 1.

        between 0 to 1. If ``inlets`` is given then this argument is ignored.
    outlet_axis : int
        Outlet axis of boundary condition. For three dimensional image the
        number ranges from 0 to 2. For two dimensional image the range is
        between 0 to 1.
        between 0 to 1. If ``outlets`` is given then this argument is ignored.
    inlets : ND-image (optional)
        A boolean mask indicating locations of inlets.  If this argument is
        supplied then ``inlet_axis`` is ignored.
    outlets : ND-image (optional)
        A boolean mask indicating locations of outlets. If this argument is
        supplied then ``outlet_axis`` is ignored.

    Returns
    -------

@@ -723,46 +795,50 @@


    """
    if im.ndim != im.squeeze().ndim:
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
                      ' unexpected behavior.')
        warnings.warn(
            "Input image conains a singleton axis:"
            + str(im.shape)
            + " Reduce dimensionality with np.squeeze(im) to avoid"
            + " unexpected behavior."
        )
    im = trim_floating_solid(~im)
    labels = spim.label(~im)[0]
    inlet = np.zeros_like(im, dtype=int)
    outlet = np.zeros_like(im, dtype=int)
    if im.ndim == 3:
        if inlet_axis == 0:
            inlet[0, :, :] = 1
        elif inlet_axis == 1:
            inlet[:, 0, :] = 1
        elif inlet_axis == 2:
            inlet[:, :, 0] = 1

        if outlet_axis == 0:
            outlet[-1, :, :] = 1
        elif outlet_axis == 1:
            outlet[:, -1, :] = 1
        elif outlet_axis == 2:
            outlet[:, :, -1] = 1

    if im.ndim == 2:
        if inlet_axis == 0:
            inlet[0, :] = 1
        elif inlet_axis == 1:
            inlet[:, 0] = 1

        if outlet_axis == 0:
            outlet[-1, :] = 1
        elif outlet_axis == 1:
            outlet[:, -1] = 1
    IN = np.unique(labels*inlet)
    OUT = np.unique(labels*outlet)
    if inlets is None:
        inlets = np.zeros_like(im, dtype=bool)
        if im.ndim == 3:
            if inlet_axis == 0:
                inlets[0, :, :] = True
            elif inlet_axis == 1:
                inlets[:, 0, :] = True
            elif inlet_axis == 2:
                inlets[:, :, 0] = True
        if im.ndim == 2:
            if inlet_axis == 0:
                inlets[0, :] = True
            elif inlet_axis == 1:
                inlets[:, 0] = True
    if outlets is None:
        outlets = np.zeros_like(im, dtype=bool)
        if im.ndim == 3:
            if outlet_axis == 0:
                outlets[-1, :, :] = True
            elif outlet_axis == 1:
                outlets[:, -1, :] = True
            elif outlet_axis == 2:
                outlets[:, :, -1] = True
        if im.ndim == 2:
            if outlet_axis == 0:
                outlets[-1, :] = True
            elif outlet_axis == 1:
                outlets[:, -1] = True
    IN = np.unique(labels * inlets)
    OUT = np.unique(labels * outlets)
    new_im = np.isin(labels, list(set(IN) ^ set(OUT)), invert=True)
    im[new_im == 0] = True
    return ~im


def trim_extrema(im, h, mode='maxima'):
def trim_extrema(im, h, mode="maxima"):
    r"""
    Trims local extrema in greyscale values by a specified amount.


@@ -790,37 +866,40 @@


    """
    result = im
    if mode in ['maxima', 'extrema']:
        result = reconstruction(seed=im - h, mask=im, method='dilation')
    elif mode in ['minima', 'extrema']:
        result = reconstruction(seed=im + h, mask=im, method='erosion')
    if mode in ["maxima", "extrema"]:
        result = reconstruction(seed=im - h, mask=im, method="dilation")
    elif mode in ["minima", "extrema"]:
        result = reconstruction(seed=im + h, mask=im, method="erosion")
    return result


@jit(forceobj=True)
def flood(im, regions=None, mode='max'):
def flood(im, regions=None, mode="max"):
    r"""
    Floods/fills each region in an image with a single value based on the
    specific values in that region.  The ``mode`` argument is used to
    determine how the value is calculated.
    specific values in that region.

    The ``mode`` argument is used to determine how the value is calculated.

    Parameters
    ----------
    im : array_like
        An ND image with isolated regions containing 0's elsewhere.

        An image with isolated regions with numerical values in each voxel,
        and 0's elsewhere.
    regions : array_like
        An array the same shape as ``im`` with each region labeled.  If None is
        supplied (default) then ``scipy.ndimage.label`` is used with its
        default arguments.

        An array the same shape as ``im`` with each region labeled.  If
        ``None`` is supplied (default) then ``scipy.ndimage.label`` is used
        with its default arguments.
    mode : string
        Specifies how to determine which value should be used to flood each
        region.  Options are:

        'max' - Floods each region with the local maximum in that region
        'maximum' - Floods each region with the local maximum in that region

        'minimum' - Floods each region the local minimum in that region

        'median' - Floods each region the local median in that region

        'min' - Floods each region the local minimum in that region
        'mean' - Floods each region the local mean in that region

        'size' - Floods each region with the size of that region


@@ -841,24 +920,16 @@

    else:
        labels = np.copy(regions)
        N = labels.max()
    I = im.flatten()
    L = labels.flatten()
    if mode.startswith('max'):
        V = np.zeros(shape=N+1, dtype=float)
        for i in range(len(L)):
            if V[L[i]] < I[i]:
                V[L[i]] = I[i]
    elif mode.startswith('min'):
        V = np.ones(shape=N+1, dtype=float)*np.inf
        for i in range(len(L)):
            if V[L[i]] > I[i]:
                V[L[i]] = I[i]
    elif mode.startswith('size'):
        V = np.zeros(shape=N+1, dtype=int)
        for i in range(len(L)):
            V[L[i]] += 1
    im_flooded = np.reshape(V[labels], newshape=im.shape)
    im_flooded = im_flooded*mask
    mode = "sum" if mode == "size" else mode
    mode = "maximum" if mode == "max" else mode
    mode = "minimum" if mode == "min" else mode
    if mode in ["mean", "median", "maximum", "minimum", "sum"]:
        f = getattr(spim, mode)
        vals = f(input=im, labels=labels, index=range(0, N + 1))
        im_flooded = vals[labels]
        im_flooded = im_flooded * mask
    else:
        raise Exception(mode + " is not a recognized mode")
    return im_flooded



@@ -886,10 +957,10 @@

        the image.  Obviously, voxels with a value of zero have no error.

    """
    temp = np.ones(shape=dt.shape)*np.inf
    temp = np.ones(shape=dt.shape) * np.inf
    for ax in range(dt.ndim):
        dt_lin = distance_transform_lin(np.ones_like(temp, dtype=bool),
                                        axis=ax, mode='both')
                                        axis=ax, mode="both")
        temp = np.minimum(temp, dt_lin)
    result = np.clip(dt - temp, a_min=0, a_max=np.inf)
    return result

@@ -913,6 +984,10 @@

        A copy of ``im`` with each voxel value indicating the size of the
        region to which it belongs.  This is particularly useful for finding
        chord sizes on the image produced by ``apply_chords``.

    See Also
    --------
    flood
    """
    if im.dtype == bool:
        im = spim.label(im)[0]

@@ -962,20 +1037,23 @@


    """
    if im.ndim != im.squeeze().ndim:
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
                      ' unexpected behavior.')
        warnings.warn(
            "Input image conains a singleton axis:"
            + str(im.shape)
            + " Reduce dimensionality with np.squeeze(im) to avoid"
            + " unexpected behavior."
        )
    if spacing < 0:
        raise Exception('Spacing cannot be less than 0')
        raise Exception("Spacing cannot be less than 0")
    if spacing == 0:
        label = True
    result = np.zeros(im.shape, dtype=int)  # Will receive chords at end
    slxyz = [slice(None, None, spacing*(axis != i) + 1) for i in [0, 1, 2]]
    slices = tuple(slxyz[:im.ndim])
    slxyz = [slice(None, None, spacing * (axis != i) + 1) for i in [0, 1, 2]]
    slices = tuple(slxyz[: im.ndim])
    s = [[0, 1, 0], [0, 1, 0], [0, 1, 0]]  # Straight-line structuring element
    if im.ndim == 3:  # Make structuring element 3D if necessary
        s = np.pad(np.atleast_3d(s), pad_width=((0, 0), (0, 0), (1, 1)),
                   mode='constant', constant_values=0)
                   mode="constant", constant_values=0)
    im = im[slices]
    s = np.swapaxes(s, 0, axis)
    chords = spim.label(im, structure=s)[0]

@@ -1025,25 +1103,28 @@


    """
    if im.ndim != im.squeeze().ndim:
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
                      ' unexpected behavior.')
        warnings.warn(
            "Input image conains a singleton axis:"
            + str(im.shape)
            + " Reduce dimensionality with np.squeeze(im) to avoid"
            + " unexpected behavior."
        )
    if im.ndim < 3:
        raise Exception('Must be a 3D image to use this function')
        raise Exception("Must be a 3D image to use this function")
    if spacing < 0:
        raise Exception('Spacing cannot be less than 0')
        raise Exception("Spacing cannot be less than 0")
    ch = np.zeros_like(im, dtype=int)
    ch[:, ::4+2*spacing, ::4+2*spacing] = 1  # X-direction
    ch[::4+2*spacing, :, 2::4+2*spacing] = 2  # Y-direction
    ch[2::4+2*spacing, 2::4+2*spacing, :] = 3  # Z-direction
    chords = ch*im
    ch[:, :: 4 + 2 * spacing, :: 4 + 2 * spacing] = 1       # X-direction
    ch[:: 4 + 2 * spacing, :, 2::4 + 2 * spacing] = 2     # Y-direction
    ch[2::4 + 2 * spacing, 2::4 + 2 * spacing, :] = 3   # Z-direction
    chords = ch * im
    if trim_edges:
        temp = clear_border(spim.label(chords > 0)[0]) > 0
        chords = temp*chords
        chords = temp * chords
    return chords


def local_thickness(im, sizes=25, mode='hybrid'):
def local_thickness(im, sizes=25, mode="hybrid", **kwargs):
    r"""
    For each voxel, this function calculates the radius of the largest sphere
    that both engulfs the voxel and fits entirely within the foreground.

@@ -1104,19 +1185,20 @@

    function.  This is not needed in ``local_thickness`` however.

    """
    im_new = porosimetry(im=im, sizes=sizes, access_limited=False, mode=mode)
    im_new = porosimetry(im=im, sizes=sizes, access_limited=False, mode=mode,
                         **kwargs)
    return im_new


def porosimetry(im, sizes=25, inlets=None, access_limited=True,
                mode='hybrid'):
def porosimetry(im, sizes=25, inlets=None, access_limited=True, mode='hybrid',
                fft=True, **kwargs):
    r"""
    Performs a porosimetry simulution on the image
    Performs a porosimetry simulution on an image

    Parameters
    ----------
    im : ND-array
        An ND image of the porous material containing True values in the
        An ND image of the porous material containing ``True`` values in the
        pore space.

    sizes : array_like or scalar

@@ -1125,7 +1207,7 @@

        the min and max of the distance transform are used.

    inlets : ND-array, boolean
        A boolean mask with True values indicating where the invasion
        A boolean mask with ``True`` values indicating where the invasion
        enters the image.  By default all faces are considered inlets,
        akin to a mercury porosimetry experiment.  Users can also apply
        solid boundaries to their image externally before passing it in,

@@ -1134,7 +1216,7 @@


    access_limited : Boolean
        This flag indicates if the intrusion should only occur from the
        surfaces (``access_limited`` is True, which is the default), or
        surfaces (``access_limited`` is ``True``, which is the default), or
        if the invading phase should be allowed to appear in the core of
        the image.  The former simulates experimental tools like mercury
        intrusion porosimetry, while the latter is useful for comparison

@@ -1159,6 +1241,12 @@

        included mostly for comparison purposes.  The morphological operations
        are done using fft-based method implementations.

    fft : boolean (default is ``True``)
        Indicates whether to use the ``fftmorphology`` function in
        ``porespy.filters`` or to use the standard morphology functions in
        ``scipy.ndimage``.  Always use ``fft=True`` unless you have a good
        reason not to.

    Returns
    -------
    image : ND-array

@@ -1168,7 +1256,7 @@

        capillary pressure by applying a boolean comparison:
        ``inv_phase = im > r`` where ``r`` is the radius (in voxels) of the
        invading sphere.  Of course, ``r`` can be converted to capillary
        pressure using your favorite model.
        pressure using a preferred model.

    Notes
    -----

@@ -1179,19 +1267,19 @@


    See Also
    --------
    fftmorphology
    local_thickness

    """
    if im.ndim != im.squeeze().ndim:
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
                      ' unexpected behavior.')
        warnings.warn("Input image contains a singleton axis:"
                      + str(im.shape)
                      + " Reduce dimensionality with np.squeeze(im) to avoid"
                      + " unexpected behavior.")

    dt = spim.distance_transform_edt(im > 0)
    dt = edt(im > 0)

    if inlets is None:
        inlets = get_border(im.shape, mode='faces')
        inlets = get_border(im.shape, mode="faces")

    if isinstance(sizes, int):
        sizes = np.logspace(start=np.log10(np.amax(dt)), stop=0, num=sizes)

@@ -1203,59 +1291,100 @@

    else:
        strel = ps_ball

    if mode == 'mio':
    # Parse kwargs for any parallelization arguments
    parallel = kwargs.pop('parallel', False)
    cores = kwargs.pop('cores', None)
    divs = kwargs.pop('divs', 2)

    if mode == "mio":
        pw = int(np.floor(dt.max()))
        impad = np.pad(im, mode='symmetric', pad_width=pw)
        inletspad = np.pad(inlets, mode='symmetric', pad_width=pw)
        inlets = np.where(inletspad)
#        sizes = np.unique(np.around(sizes, decimals=0).astype(int))[-1::-1]
        impad = np.pad(im, mode="symmetric", pad_width=pw)
        inlets = np.pad(inlets, mode="symmetric", pad_width=pw)
        # sizes = np.unique(np.around(sizes, decimals=0).astype(int))[-1::-1]
        imresults = np.zeros(np.shape(impad))
        for r in tqdm(sizes):
            imtemp = fftmorphology(impad, strel(r), mode='erosion')
            if access_limited:
                imtemp = trim_disconnected_blobs(imtemp, inlets)
            imtemp = fftmorphology(imtemp, strel(r), mode='dilation')
            if np.any(imtemp):
                imresults[(imresults == 0)*imtemp] = r
        with tqdm(sizes) as pbar:
            for r in sizes:
                pbar.update()
                if parallel:
                    imtemp = chunked_func(func=spim.binary_erosion,
                                          input=impad, structure=strel(r),
                                          overlap=int(2*r) + 1,
                                          cores=cores, divs=divs)
                elif fft:
                    imtemp = fftmorphology(impad, strel(r), mode="erosion")
                else:
                    imtemp = spim.binary_erosion(input=impad,
                                                 structure=strel(r))
                if access_limited:
                    imtemp = trim_disconnected_blobs(imtemp, inlets)
                if parallel:
                    imtemp = chunked_func(func=spim.binary_dilation,
                                          input=imtemp, structure=strel(r),
                                          overlap=int(2*r) + 1,
                                          cores=cores, divs=divs)
                elif fft:
                    imtemp = fftmorphology(imtemp, strel(r), mode="dilation")
                else:
                    imtemp = spim.binary_dilation(input=imtemp,
                                                  structure=strel(r))
                if np.any(imtemp):
                    imresults[(imresults == 0) * imtemp] = r
        imresults = extract_subsection(imresults, shape=im.shape)
    elif mode == 'dt':
        inlets = np.where(inlets)
    elif mode == "dt":
        imresults = np.zeros(np.shape(im))
        for r in tqdm(sizes):
            imtemp = dt >= r
            if access_limited:
                imtemp = trim_disconnected_blobs(imtemp, inlets)
            if np.any(imtemp):
                imtemp = spim.distance_transform_edt(~imtemp) < r
                imresults[(imresults == 0)*imtemp] = r
    elif mode == 'hybrid':
        inlets = np.where(inlets)
        with tqdm(sizes) as pbar:
            for r in sizes:
                pbar.update()
                imtemp = dt >= r
                if access_limited:
                    imtemp = trim_disconnected_blobs(imtemp, inlets)
                if np.any(imtemp):
                    imtemp = edt(~imtemp) < r
                    imresults[(imresults == 0) * imtemp] = r
    elif mode == "hybrid":
        imresults = np.zeros(np.shape(im))
        for r in tqdm(sizes):
            imtemp = dt >= r
            if access_limited:
                imtemp = trim_disconnected_blobs(imtemp, inlets)
            if np.any(imtemp):
                imtemp = fftconvolve(imtemp, strel(r), mode='same') > 0.0001
                imresults[(imresults == 0)*imtemp] = r
        with tqdm(sizes) as pbar:
            for r in sizes:
                pbar.update()
                imtemp = dt >= r
                if access_limited:
                    imtemp = trim_disconnected_blobs(imtemp, inlets)
                if np.any(imtemp):
                    if parallel:
                        imtemp = chunked_func(func=spim.binary_dilation,
                                              input=imtemp, structure=strel(r),
                                              overlap=int(2*r) + 1,
                                              cores=cores, divs=divs)
                    elif fft:
                        imtemp = fftmorphology(imtemp, strel(r),
                                               mode="dilation")
                    else:
                        imtemp = spim.binary_dilation(input=imtemp,
                                                      structure=strel(r))
                    imresults[(imresults == 0) * imtemp] = r
    else:
        raise Exception('Unreckognized mode ' + mode)
        raise Exception("Unrecognized mode " + mode)
    return imresults


def trim_disconnected_blobs(im, inlets):
def trim_disconnected_blobs(im, inlets, strel=None):
    r"""
    Removes foreground voxels not connected to specified inlets

    Parameters
    ----------
    im : ND-array
        The array to be trimmed
        The image containing the blobs to be trimmed
    inlets : ND-array or tuple of indices
        The locations of the inlets.  Can either be a boolean mask the same
        shape as ``im``, or a tuple of indices such as that returned by the
        ``where`` function.  Any voxels *not* connected directly to
        the inlets will be trimmed.
    strel : ND-array
        The neighborhood over which connectivity should be checked. It must
        be symmetric and the same dimensionality as the image.  It is passed
        directly to the ``scipy.ndimage.label`` function as the ``structure``
        argument so refer to that docstring for additional info.

    Returns
    -------

@@ -1270,23 +1399,27 @@

    elif (inlets.shape == im.shape) and (inlets.max() == 1):
        inlets = inlets.astype(bool)
    else:
        raise Exception('inlets not valid, refer to docstring for info')
    labels = spim.label(inlets + (im > 0))[0]
        raise Exception("inlets not valid, refer to docstring for info")
    if im.ndim == 3:
        strel = cube
    else:
        strel = square
    labels = spim.label(inlets + (im > 0), structure=strel(3))[0]
    keep = np.unique(labels[inlets])
    keep = keep[keep > 0]
    if len(keep) > 0:
        im2 = np.reshape(np.in1d(labels, keep), newshape=im.shape)
    else:
        im2 = np.zeros_like(im)
    im2 = im2*im
    im2 = im2 * im
    return im2


def _get_axial_shifts(ndim=2, include_diagonals=False):
    r'''
    r"""
    Helper function to generate the axial shifts that will be performed on
    the image to identify bordering pixels/voxels
    '''
    """
    if ndim == 2:
        if include_diagonals:
            neighbors = square(3)

@@ -1311,40 +1444,37 @@



def _make_stack(im, include_diagonals=False):
    r'''
    r"""
    Creates a stack of images with one extra dimension to the input image
    with length equal to the number of borders to search + 1.
    Image is rolled along the axial shifts so that the border pixel is
    overlapping the original pixel. First image in stack is the original.
    Stacking makes direct vectorized array comparisons possible.
    '''
    """
    ndim = len(np.shape(im))
    axial_shift = _get_axial_shifts(ndim, include_diagonals)
    if ndim == 2:
        stack = np.zeros([np.shape(im)[0],
                          np.shape(im)[1],
                          len(axial_shift)+1])
        stack = np.zeros([np.shape(im)[0], np.shape(im)[1], len(axial_shift) + 1])
        stack[:, :, 0] = im
        for i in range(len(axial_shift)):
            ax0, ax1 = axial_shift[i]
            temp = np.roll(np.roll(im, ax0, 0), ax1, 1)
            stack[:, :, i+1] = temp
            stack[:, :, i + 1] = temp
        return stack
    elif ndim == 3:
        stack = np.zeros([np.shape(im)[0],
                          np.shape(im)[1],
                          np.shape(im)[2],
                          len(axial_shift)+1])
        stack = np.zeros(
            [np.shape(im)[0], np.shape(im)[1], np.shape(im)[2], len(axial_shift) + 1]
        )
        stack[:, :, :, 0] = im
        for i in range(len(axial_shift)):
            ax0, ax1, ax2 = axial_shift[i]
            temp = np.roll(np.roll(np.roll(im, ax0, 0), ax1, 1), ax2, 2)
            stack[:, :, :, i+1] = temp
            stack[:, :, :, i + 1] = temp
        return stack


def nphase_border(im, include_diagonals=False):
    r'''
    r"""
    Identifies the voxels in regions that border *N* other regions.

    Useful for finding triple-phase boundaries.

@@ -1365,17 +1495,20 @@

    image : ND-array
        A copy of ``im`` with voxel values equal to the number of uniquely
        different bordering values
    '''
    """
    if im.ndim != im.squeeze().ndim:
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
                      ' unexpected behavior.')
        warnings.warn(
            "Input image conains a singleton axis:"
            + str(im.shape)
            + " Reduce dimensionality with np.squeeze(im) to avoid"
            + " unexpected behavior."
        )
    # Get dimension of image
    ndim = len(np.shape(im))
    if ndim not in [2, 3]:
        raise NotImplementedError("Function only works for 2d and 3d images")
    # Pad image to handle edges
    im = np.pad(im, pad_width=1, mode='edge')
    im = np.pad(im, pad_width=1, mode="edge")
    # Stack rolled images for each neighbor to be inspected
    stack = _make_stack(im, include_diagonals)
    # Sort the stack along the last axis

@@ -1385,9 +1518,9 @@

    # Number of changes is number of unique bordering regions
    for k in range(np.shape(stack)[ndim])[1:]:
        if ndim == 2:
            mask = stack[:, :, k] != stack[:, :, k-1]
            mask = stack[:, :, k] != stack[:, :, k - 1]
        elif ndim == 3:
            mask = stack[:, :, :, k] != stack[:, :, :, k-1]
            mask = stack[:, :, :, k] != stack[:, :, :, k - 1]
        out += mask
    # Un-pad
    if ndim == 2:

@@ -1396,7 +1529,7 @@

        return out[1:-1, 1:-1, 1:-1].copy()


def prune_branches(skel, branch_points=None, iterations=1):
def prune_branches(skel, branch_points=None, iterations=1, **kwargs):
    r"""
    Removes all dangling ends or tails of a skeleton.


@@ -1421,20 +1554,28 @@

        from skimage.morphology import square as cube
    else:
        from skimage.morphology import cube
    parallel = kwargs.pop('parallel', False)
    divs = kwargs.pop('divs', 2)
    cores = kwargs.pop('cores', None)
    # Create empty image to house results
    im_result = np.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
        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)
    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))
    if parallel:
        branch_points = chunked_func(func=spim.binary_dilation,
                                     input=branch_points, structure=cube(3),
                                     overlap=3, divs=divs, cores=cores)
    else:
        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)

@@ -1442,18 +1583,750 @@

    for s in slices:
        label_num += 1
        # Find branch point labels the overlap current arc
        hits = pts_labels[s]*(arc_labels[s] == label_num)
        hits = pts_labels[s] * (arc_labels[s] == label_num)
        # If image contains 2 branch points, then it's not a tail.
        if len(np.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
    im_result += skel * pts_orig
    if iterations > 1:
        iterations -= 1
        im_temp = np.copy(im_result)
        im_result = prune_branches(skel=im_result,
                                   branch_points=None,
                                   iterations=iterations)
                                   iterations=iterations,
                                   parallel=parallel,
                                   divs=divs, cores=cores)
        if np.all(im_temp == im_result):
            iterations = 0
    return im_result


def chunked_func(func,
                 overlap=None,
                 divs=2,
                 cores=None,
                 im_arg=["input", "image", "im"],
                 strel_arg=["strel", "structure", "footprint"],
                 **kwargs):
    r"""
    Performs the specfied operation "chunk-wise" in parallel using dask

    This can be used to save memory by doing one chunk at a time (``cores=1``)
    or to increase computation speed by spreading the work across multiple
    cores (e.g. ``cores=8``)

    This function can be used with any operation that applies a structuring
    element of some sort, since this implies that the operation is local
    and can be chunked.

    Parameters
    ----------
    func : function handle
        The function which should be applied to each chunk, such as
        ``spipy.ndimage.binary_dilation``.
    overlap : scalar or list of scalars, optional
        The amount of overlap to include when dividing up the image.  This
        value will almost always be the size (i.e. diameter) 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.
    divs : scalar or list of scalars (default = [2, 2, 2])
        The number of chunks to divide the image into in each direction.  The
        default is 2 chunks in each direction, resulting in a quartering of
        the image and 8 total chunks (in 3D).  A scalar is interpreted as
        applying to all directions, while a list of scalars is interpreted
        as applying to each individual direction.
    cores : scalar
        The number of cores which should be used.  By default, all cores will
        be used, or as many are needed for the given number of chunks, which
        ever is smaller.
    im_arg : string
        The keyword used by ``func`` for the image to be operated on.  By
        default this function will look for ``image``, ``input``, and ``im``
        which are commonly used by *scipy.ndimage* and *skimage*.
    strel_arg : string
        The keyword used by ``func`` for the structuring element to apply.
        This is only needed if ``overlap`` is not specified. By default this
        function will look for ``strel``, ``structure``, and ``footprint``
        which are commonly used by *scipy.ndimage* and *skimage*.
    kwargs : additional keyword arguments
        All other arguments are passed to ``func`` as keyword arguments. Note
        that PoreSpy will fetch the image from this list of keywords using the
        value provided to ``im_arg``.

    Returns
    -------
    result : ND-image
        An image the same size as the input image, with the specified filter
        applied as though done on a single large image.  There should be *no*
        difference.

    Notes
    -----
    This function divides the image into the specified number of chunks, but
    also applies a padding to each chunk to create an overlap with neighboring
    chunks.  This way the operation does not have any edge artifacts. The
    amount of padding is usually equal to the radius of the structuring
    element but some functions do not use one, such as the distance transform
    and Gaussian blur.  In these cases the user can specify ``overlap``.

    See Also
    --------
    skikit-image.util.apply_parallel

    Examples
    --------
    >>> import scipy.ndimage as spim
    >>> import porespy as ps
    >>> from skimage.morphology import ball
    >>> im = ps.generators.blobs(shape=[100, 100, 100])
    >>> f = spim.binary_dilation
    >>> im2 = ps.filters.chunked_func(func=f, overlap=7, im_arg='input',
    ...                               input=im, structure=ball(3), cores=1)
    >>> im3 = spim.binary_dilation(input=im, structure=ball(3))
    >>> np.all(im2 == im3)
    True

    """

    @dask.delayed
    def apply_func(func, **kwargs):
        # Apply function on sub-slice of overall image
        return func(**kwargs)

    # Import the array_split methods
    from array_split import shape_split, ARRAY_BOUNDS

    # Determine the value for im_arg
    if type(im_arg) == str:
        im_arg = [im_arg]
    for item in im_arg:
        if item in kwargs.keys():
            im = kwargs[item]
            im_arg = item
            break
    # Fetch image from the kwargs dict
    im = kwargs[im_arg]
    # Determine the number of divisions to create
    divs = np.ones((im.ndim,), dtype=int) * np.array(divs)
    # If overlap given then use it, otherwise search for strel in kwargs
    if overlap is not None:
        halo = overlap * (divs > 1)
    else:
        if type(strel_arg) == str:
            strel_arg = [strel_arg]
        for item in strel_arg:
            if item in kwargs.keys():
                strel = kwargs[item]
                break
        halo = np.array(strel.shape) * (divs > 1)
    slices = np.ravel(
        shape_split(
            im.shape, axis=divs, halo=halo.tolist(), tile_bounds_policy=ARRAY_BOUNDS
        )
    )
    # Apply func to each subsection of the image
    res = []
    # print('Image will be broken into the following chunks:')
    for s in slices:
        # Extract subsection from image and input into kwargs
        kwargs[im_arg] = im[tuple(s)]
        # print(kwargs[im_arg].shape)
        res.append(apply_func(func=func, **kwargs))
    # Have dask actually compute the function on each subsection in parallel
    # with ProgressBar():
    #    ims = dask.compute(res, num_workers=cores)[0]
    ims = dask.compute(res, num_workers=cores)[0]
    # Finally, put the pieces back together into a single master image, im2
    im2 = np.zeros_like(im, dtype=im.dtype)
    for i, s in enumerate(slices):
        # Prepare new slice objects into main and sub-sliced image
        a = []  # Slices into main image
        b = []  # Slices into chunked image
        for dim in range(im.ndim):
            if s[dim].start == 0:
                ax = bx = 0
            else:
                ax = s[dim].start + halo[dim]
                bx = halo[dim]
            if s[dim].stop == im.shape[dim]:
                ay = by = im.shape[dim]
            else:
                ay = s[dim].stop - halo[dim]
                by = s[dim].stop - s[dim].start - halo[dim]
            a.append(slice(ax, ay, None))
            b.append(slice(bx, by, None))
        # Convert lists of slices to tuples
        a = tuple(a)
        b = tuple(b)
        # Insert image chunk into main image
        im2[a] = ims[i][b]
    return im2


def snow_partitioning_parallel(im,
                               overlap='dt',
                               divs=2,
                               mode='parallel',
                               num_workers=None,
                               crop=True,
                               zoom_factor=0.5,
                               r_max=5,
                               sigma=0.4,
                               return_all=False):
    r"""
    Perform SNOW algorithm in parallel and serial mode to reduce time and
    memory usage repectively by geomertirc domain decomposition of large size
    image.

    Parameters
    ----------
    im: ND_array
        A binary image of porous media with 'True' values indicating phase of
        interest

    overlap: float or int or str
        Overlapping thickness between two subdomains that is used to merge
        watershed segmented regions at intersection of two or more subdomains.
        If 'dt' the overlap will be calculated based on maximum
        distance transform in the whole image.
        If 'ws' the overlap will be calculated by finding the maximum dimension
        of the bounding box of largest segmented region. The image is scale down
        by 'zoom_factor' provided by user.
        If any real number of overlap is provided then this value will be
        considered as overlapping thickness.

    divs: list or int
        Number of domains each axis will be divided.
        If a scalar is provided then it will be assigned to all axis.
        If list is provided then each respective axis will be divided by its
        corresponding number in the list. For example [2, 3, 4] will divide
        z, y and x axis to 2, 3, and 4 respectively.

    mode: str
        if 'parallel' then all subdomains will be processed in number of cores
        provided as num_workers
        if 'serial' then all subdomains will be processed one by one in one core
        of CPU.

    num_workers: int or None
        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.

    crop: bool
        If True the image shape is cropped to fit specified division.

    zoom_factor: float or int
        The amount of zoom appiled to image to find overlap thickness using "ws"
        overlap mode.

    return_all : boolean
        If set to ``True`` a named tuple is returned containing the original
        image, the distance transform, and the final
        pore regions.  The default is ``False``

    Returns
    ----------
    regions: ND_array
        Partitioned image of segmentated regions with unique labels. Each
        region correspond to pore body while intersection with other region
        correspond throat area.
    """
    # --------------------------------------------------------------------------
    # Adjust image shape according to specified dimension
    tup = namedtuple("results", field_names=["im", "dt", "regions"])
    if isinstance(divs, int):
        divs = [divs for i in range(im.ndim)]
    shape = []
    for i in range(im.ndim):
        shape.append(divs[i] * (im.shape[i] // divs[i]))

    if shape != list(im.shape):
        if crop:
            for i in range(im.ndim):
                im = im.swapaxes(0, i)
                im = im[:shape[i], ...]
                im = im.swapaxes(i, 0)
            print(f'Image is cropped to shape {shape}.')
        else:
            print('-' * 80)
            print(f"Possible image shape for specified divisions is {shape}.")
            print("To crop the image please set crop argument to 'True'.")
            return
    # --------------------------------------------------------------------------
    # Get overlap thickness from distance transform
    chunk_shape = (np.array(shape) / np.array(divs)).astype(int)
    print('# Beginning parallel SNOW algorithm...')
    print('=' * 80)
    print('Calculating overlap thickness')
    if overlap == 'dt':
        dt = edt((im > 0), parallel=0)
        overlap = dt.max()
    elif overlap == 'ws':
        rev = spim.interpolation.zoom(im, zoom=zoom_factor, order=0)
        rev = rev > 0
        dt = edt(rev, parallel=0)
        rev_snow = snow_partitioning(rev, dt=dt, r_max=r_max, sigma=sigma)
        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_factor * 2.0)
        dt = edt((im > 0), parallel=0)
    else:
        overlap = overlap / 2.0
        dt = edt((im > 0), parallel=0)
    print('Overlap Thickness: ' + str(int(2.0 * overlap)) + ' voxels')
    # --------------------------------------------------------------------------
    # Get overlap and trim depth of all image dimension
    depth = {}
    trim_depth = {}
    for i in range(im.ndim):
        depth[i] = int(2.0 * overlap)
        trim_depth[i] = int(2.0 * overlap) - 1

    tup.im = im
    tup.dt = dt
    # --------------------------------------------------------------------------
    # Applying snow to image chunks
    im = da.from_array(dt, chunks=chunk_shape)
    im = da.overlap.overlap(im, depth=depth, boundary='none')
    im = im.map_blocks(chunked_snow, r_max=r_max, sigma=sigma)
    im = da.overlap.trim_internal(im, trim_depth, boundary='none')
    if mode == 'serial':
        num_workers = 1
    elif mode == 'parallel':
        num_workers = num_workers
    else:
        raise Exception('Mode of operation can either be parallel or serial')
    with ProgressBar():
        # print('-' * 80)
        print('Applying snow to image chunks')
        regions = im.compute(num_workers=num_workers)
    # --------------------------------------------------------------------------
    # Relabelling watershed chunks
    # print('-' * 80)
    print('Relabelling watershed chunks')
    regions = relabel_chunks(im=regions, chunk_shape=chunk_shape)
    # --------------------------------------------------------------------------
    # Stitching watershed chunks
    # print('-' * 80)
    print('Stitching watershed chunks')
    regions = watershed_stitching(im=regions, chunk_shape=chunk_shape)
    print('=' * 80)
    if return_all:
        tup.regions = regions
        return tup
    else:
        return regions

    return regions


def chunked_snow(im, r_max=5, sigma=0.4):
    r"""
    Partitions the void space into pore regions using a marker-based watershed
    algorithm, with specially filtered peaks as markers.

    The SNOW network extraction algorithm (Sub-Network of an Over-segmented
    Watershed) was designed to handle to perculiarities of high porosity
    materials, but it applies well to other materials as well.

    Parameters
    ----------
    im : array_like
        Distance transform of phase of interest in a binary image
    r_max : int
        The radius of the spherical structuring element to use in the Maximum
        filter stage that is used to find peaks.  The default is 5
    sigma : float
        The standard deviation of the Gaussian filter used in step 1.  The
        default is 0.4.  If 0 is given then the filter is not applied, which is
        useful if a distance transform is supplied as the ``im`` argument that
        has already been processed.

    Returns
    -------
    image : ND-array
        An image the same shape as ``im`` with the void space partitioned into
        pores using a marker based watershed with the peaks found by the
        SNOW algorithm [1].

    References
    ----------
    [1] Gostick, J. "A versatile and efficient network extraction algorithm
    using marker-based watershed segmenation".  Physical Review E. (2017)
    """

    dt = spim.gaussian_filter(input=im, sigma=sigma)
    peaks = find_peaks(dt=dt, r_max=r_max)
    peaks = trim_saddle_points(peaks=peaks, dt=dt, max_iters=99, verbose=0)
    peaks = trim_nearby_peaks(peaks=peaks, dt=dt)
    peaks, N = spim.label(peaks)
    regions = watershed(image=-dt, markers=peaks, mask=im > 0)

    return regions * (im > 0)


def pad(im, pad_width=1, constant_value=0):
    r"""
    Pad the image with a constant values and width.

    Parameters:
    ----------
    im : ND-array
        The image that requires padding
    pad_width : int
        The number of values that will be padded from the edges. Default values
        is 1.
    contant_value : int
        Pads with the specified constant value

    return:
    -------
    output: ND-array
        Padded image with same dimnesions as provided image
    """
    shape = np.array(im.shape)
    pad_shape = shape + (2 * pad_width)
    temp = np.zeros(pad_shape, dtype=np.uint32)
    if constant_value != 0:
        temp = temp + constant_value
    if im.ndim == 3:
        temp[pad_width: -pad_width,
             pad_width: -pad_width,
             pad_width: -pad_width] = im
    elif im.ndim == 2:
        temp[pad_width: -pad_width,
             pad_width: -pad_width] = im
    else:
        temp[pad_width: -pad_width] = im

    return temp


def relabel_chunks(im, chunk_shape):  # pragma: no cover
    r"""
    Assign new labels to each chunk or sub-domain of actual image. This
    prevents from two or more regions to have same label.

    Parameters:
    -----------

    im: ND-array
        Actual image that contains repeating labels in chunks or sub-domains

    chunk_shape: tuple
        The shape of chunk that will be relabeled in actual image. Note the
        chunk shape should be a multiple of actual image shape otherwise some
        labels will not be relabeled.

    return:
    -------
    output : ND-array
        Relabeled image with unique label assigned to each region.
    """
    im = pad(im, pad_width=1)
    im_shape = np.array(im.shape, dtype=np.uint32)
    max_num = 0
    c = np.array(chunk_shape, dtype=np.uint32) + 2
    num = (im_shape / c).astype(int)

    if im.ndim == 3:
        for z in range(num[0]):
            for y in range(num[1]):
                for x in range(num[2]):
                    chunk = im[z * c[0]: (z + 1) * c[0],
                               y * c[1]: (y + 1) * c[1],
                               x * c[2]: (x + 1) * c[2]]
                    chunk += max_num
                    chunk[chunk == max_num] = 0
                    max_num = chunk.max()
                    im[z * c[0]: (z + 1) * c[0],
                       y * c[1]: (y + 1) * c[1],
                       x * c[2]: (x + 1) * c[2]] = chunk
    else:
        for y in range(num[0]):
            for x in range(num[1]):
                chunk = im[y * c[0]: (y + 1) * c[0],
                           x * c[1]: (x + 1) * c[1]]
                chunk += max_num
                chunk[chunk == max_num] = 0
                max_num = chunk.max()
                im[y * c[0]: (y + 1) * c[0],
                   x * c[1]: (x + 1) * c[1]] = chunk

    return im


def trim_internal_slice(im, chunk_shape):  # pragma: no cover
    r"""
    Delete extra slices from image that were used to stitch two or more chunks
    together.

    Parameters:
    -----------

    im :  ND-array
        image that contains extra slices in x, y, z direction.

    chunk_shape : tuple
        The shape of the chunk from which image is subdivided.

    Return:
    -------
    output :  ND-array
        Image without extra internal slices. The shape of the image will be
        same as input image provided for  waterhsed segmentation.
    """
    im_shape = np.array(im.shape, dtype=np.uint32)
    c1 = np.array(chunk_shape, dtype=np.uint32) + 2
    c2 = np.array(chunk_shape, dtype=np.uint32)
    num = (im_shape / c1).astype(int)
    out_shape = num * c2
    out = np.empty((out_shape), dtype=np.uint32)

    if im.ndim == 3:
        for z in range(num[0]):
            for y in range(num[1]):
                for x in range(num[2]):
                    chunk = im[z * c1[0]: (z + 1) * c1[0],
                               y * c1[1]: (y + 1) * c1[1],
                               x * c1[2]: (x + 1) * c1[2]]

                    out[z * c2[0]: (z + 1) * c2[0],
                        y * c2[1]: (y + 1) * c2[1],
                        x * c2[2]: (x + 1) * c2[2]] = chunk[1:-1, 1:-1, 1:-1]
    else:
        for y in range(num[0]):
            for x in range(num[1]):
                chunk = im[y * c1[0]: (y + 1) * c1[0],
                           x * c1[1]: (x + 1) * c1[1]]

                out[y * c2[0]: (y + 1) * c2[0],
                    x * c2[1]: (x + 1) * c2[1]] = chunk[1:-1, 1:-1]

    return out


def watershed_stitching(im, chunk_shape):
    r"""
    Stitch individual sub-domains of watershed segmentation into one big
    segmentation with all boundary labels of each sub-domain relabeled to merge
    boundary regions.

    Parameters:
    -----------
    im : ND-array
        A worked image with watershed segmentation performed on all sub-domains
        individually.

    chunk_shape: tuple
        The shape of the sub-domain in which image segmentation is performed.

    return:
    -------
    output : ND-array
        Stitched watershed segmentation with all sub-domains merged to form a
        single watershed segmentation.
    """

    c_shape = np.array(chunk_shape)
    cuts_num = (np.array(im.shape) / c_shape).astype(np.uint32)

    for axis, num in enumerate(cuts_num):
        keys = []
        values = []
        if num > 1:
            im = im.swapaxes(0, axis)
            for i in range(1, num):
                sl = i * (chunk_shape[axis] + 3) - (i - 1)
                sl1 = im[sl - 3, ...]
                sl1_mask = sl1 > 0
                sl2 = im[sl - 1, ...] * sl1_mask
                sl1_labels = sl1.flatten()[sl1.flatten() > 0]
                sl2_labels = sl2.flatten()[sl2.flatten() > 0]
                if sl1_labels.size != sl2_labels.size:
                    raise Exception('The selected overlapping thickness is not '
                                    'suitable for input image. Change '
                                    'overlapping criteria '
                                    'or manually input value.')
                keys.append(sl1_labels)
                values.append(sl2_labels)
            im = replace_labels(array=im, keys=keys, values=values)
            im = im.swapaxes(axis, 0)
    im = trim_internal_slice(im=im, chunk_shape=chunk_shape)
    im = resequence_labels(array=im)

    return im


@njit(parallel=True)
def copy(im, output):  # pragma: no cover
    r"""
    The function copy the input array and make output array that is allocated
    in different memory space. This a numba version of copy function of numpy.
    Because each element is copied using parallel approach this implementation
    is facter than numpy version of copy.

    parameter:
    ----------
    array: ND-array
        Array that needs to be copied

    Return:
    -------
    output: ND-array
        Copied array
    """

    if im.ndim == 3:
        for i in prange(im.shape[0]):
            for j in prange(im.shape[1]):
                for k in prange(im.shape[2]):
                    output[i, j, k] = im[i, j, k]
    elif im.ndim == 2:
        for i in prange(im.shape[0]):
            for j in prange(im.shape[1]):
                output[i, j] = im[i, j]
    else:
        for i in prange(im.shape[0]):
            output[i] = im[i]

    return output


@njit(parallel=True)
def _replace(array, keys, values, ind_sort):  # pragma: no cover
    r"""
    This function replace keys elements in input array with new value elements.
    This function is used as internal function of replace_relabels.

    Parameter:
    ----------
    array : ND-array
        Array which requires replacing labels
    keys :  1D-array
        The unique labels that need to be replaced
    values : 1D-array
        The unique values that will be assigned to labels

    return:
    -------
    array : ND-array
        Array with replaced labels.
    """
    # ind_sort = np.argsort(keys)
    keys_sorted = keys[ind_sort]
    values_sorted = values[ind_sort]
    s_keys = set(keys)

    for i in prange(array.shape[0]):
        if array[i] in s_keys:
            ind = np.searchsorted(keys_sorted, array[i])
            array[i] = values_sorted[ind]


def replace_labels(array, keys, values):
    r"""
    Replace labels in array provided as keys to values.

    Parameter:
    ----------
    array : ND-array
        Array which requires replacing labels
    keys :  1D-array
        The unique labels that need to be replaced
    values : 1D-array
        The unique values that will be assigned to labels

    return:
    -------
    array : ND-array
        Array with replaced labels.
    """
    a_shape = array.shape
    array = array.flatten()
    keys = np.concatenate(keys, axis=0)
    values = np.concatenate(values, axis=0)
    ind_sort = np.argsort(keys)
    _replace(array, keys, values, ind_sort)

    return array.reshape(a_shape)


@njit()
def _sequence(array, count):  # pragma: no cover
    r"""
    Internal function of resequnce_labels method. This function resquence array
    elements in an ascending order using numba technique which is many folds
    faster than make contigious funcition.

    parameter:
    ----------
    array: 1d-array
        1d-array that needs resquencing
    count: 1d-array
        1d-array of zeros having same size as array

    return:
    -------
    array: 1d-array
        The input array with elements resequenced in ascending order
    Note: The output of this function is not same as make_contigous or
    relabel_sequential function of scikit-image. This function resequence and
    randomize the regions while other methods only do resequencing and output
    sorted array.
    """
    a = 1
    i = 0
    while i < (len(array)):
        data = array[i]
        if data != 0:
            if count[data] == 0:
                count[data] = a
                a += 1
        array[i] = count[data]
        i += 1


@njit(parallel=True)
def amax(array):  # pragma: no cover
    r"""
    Find largest element in an array using fast parallel numba technique

    Parameter:
    ----------
    array: ND-array
        array in which largest elements needs to be calcuted

    return:
    scalar: float or int
        The largest element value in the input array
    """

    return np.max(array)


def resequence_labels(array):
    r"""
    Resequence the lablels to make them contigious.

    Parameter:
    ----------
    array: ND-array
        Array that requires resequencing

    return:
    -------
    array : ND-array
        Resequenced array with same shape as input array
    """
    a_shape = array.shape
    array = array.ravel()
    max_num = amax(array) + 1
    count = np.zeros(max_num, dtype=np.uint32)
    _sequence(array, count)

    return array.reshape(a_shape)

@@ -12,7 +12,6 @@

              voxel_size=1,
              boundary_faces=['top', 'bottom', 'left', 'right', 'front', 'back'],
              marching_cubes_area=False):

    r"""
    Analyzes an image that has been partitioned into void and solid regions
    and extracts the void and solid phase geometry as well as network

@@ -1,11 +1,12 @@

import numpy as np
import scipy as sp
import numpy as np
import scipy.ndimage as spim
import warnings
from edt import edt
from collections import namedtuple
from skimage.morphology import ball, disk
from skimage.measure import marching_cubes_lewiner
from array_split import shape_split
from array_split import shape_split, ARRAY_BOUNDS
from scipy.signal import fftconvolve



@@ -27,9 +28,9 @@

        Returns a copy of ``im`` rotated accordingly.
    """
    if im.ndim != im.squeeze().ndim:
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
                      ' unexpected behavior.')
        warnings.warn(f'Input image conains a singleton axis: {im.shape}'
                      + ' Reduce dimensionality with np.squeeze(im) to avoid'
                      + ' unexpected behavior.')
    im = np.copy(im)
    if im.ndim == 2:
        im = (np.swapaxes(im, 1, 0))

@@ -112,9 +113,9 @@

        return t

    if im.ndim != im.squeeze().ndim:
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
                      ' unexpected behavior.')
        warnings.warn(f'Input image conains a singleton axis: {im.shape}'
                      + ' Reduce dimensionality with np.squeeze(im) to avoid'
                      + ' unexpected behavior.')

    # Perform erosion and dilation
    # The array must be padded with 0's so it works correctly at edges

@@ -141,7 +142,7 @@

    return result


def subdivide(im, divs=2):
def subdivide(im, divs=2, overlap=0, flatten=False):
    r"""
    Returns slices into an image describing the specified number of sub-arrays.


@@ -158,48 +159,53 @@

        The number of sub-divisions to create in each axis of the image.  If a
        scalar is given it is assumed this value applies in all dimensions.

    overlap : scalar or array_like
        The amount of overlap to use when dividing along each axis.  If a
        scalar is given it is assumed this value applies in all dimensions.

    flatten : boolean
        If set to ``True`` then the slice objects are returned as a flat
        list, while if ``False`` they are returned in a ND-array where each
        subdivision is accessed using row-col or row-col-layer indexing.

    Returns
    -------
    slices : 1D-array
        A 1-D array containing slice objects for indexing into ``im`` that
        extract the sub-divided arrays.
    slices : ND-array
        An ND-array containing sets of slice objects for indexing into ``im``
        that extract subdivisions of an image.  If ``flatten`` was ``True``,
        then this array is flat, suitable  for iterating.  If ``flatten`` was
        ``False`` then the slice objects must be accessed by row, col, layer
        indices.  An ND-array is the preferred container since it's shape can
        be easily queried.

    Notes
    -----
    This method uses the
    `array_split package <https://github.com/array-split/array_split>`_ which
    offers the same functionality as the ``split`` method of Numpy's ND-array,
    but supports the splitting multidimensional arrays in all dimensions.
    but supports the splitting of multidimensional arrays in all dimensions.

    See Also
    --------
    chunked_func

    Examples
    --------
    >>> import porespy as ps
    >>> import matplotlib.pyplot as plt
    >>> im = ps.generators.blobs(shape=[200, 200])
    >>> s = ps.tools.subdivide(im, divs=[2, 2])

    ``s`` contains an array with the shape given by ``divs``.  To access the
    first and last quadrants of ``im`` use:
    >>> print(im[tuple(s[0, 0])].shape)
    (100, 100)
    >>> print(im[tuple(s[1, 1])].shape)
    (100, 100)

    It can be easier to index the array with the slices by applying ``flatten``
    first:
    >>> s_flat = s.flatten()
    >>> for i in s_flat:
    ...     print(im[i].shape)
    (100, 100)
    (100, 100)
    (100, 100)
    (100, 100)
    >>> s = ps.tools.subdivide(im, divs=[2, 2], flatten=True)
    >>> print(len(s))
    4

    """
    # Expand scalar divs
    if isinstance(divs, int):
        divs = [divs for i in range(im.ndim)]
    s = shape_split(im.shape, axis=divs)
    return s
    divs = np.ones((im.ndim,), dtype=int) * np.array(divs)
    halo = overlap * (divs > 1)
    slices = shape_split(im.shape, axis=divs, halo=halo.tolist(),
                         tile_bounds_policy=ARRAY_BOUNDS).astype(object)
    if flatten is True:
        slices = np.ravel(slices)
    return slices


def bbox_to_slices(bbox):

@@ -265,16 +271,16 @@


    """
    if r == 0:
        dt = spim.distance_transform_edt(input=im)
        dt = edt(input=im)
        r = int(np.amax(dt)) * 2
    im_padded = np.pad(array=im, pad_width=r, mode='constant',
                       constant_values=True)
    dt = spim.distance_transform_edt(input=im_padded)
    dt = edt(input=im_padded)
    seeds = (dt >= r) + get_border(shape=im_padded.shape)
    # Remove seeds not connected to edges
    labels = spim.label(seeds)[0]
    mask = labels == 1  # Assume label of 1 on edges, assured by adding border
    dt = spim.distance_transform_edt(~mask)
    dt = edt(~mask)
    outer_region = dt < r
    outer_region = extract_subsection(im=outer_region, shape=im.shape)
    return outer_region

@@ -314,7 +320,7 @@

    dim = [range(int(-s / 2), int(s / 2) + s % 2) for s in im.shape]
    inds = np.meshgrid(*dim, indexing='ij')
    inds[axis] = inds[axis] * 0
    d = np.sqrt(np.sum(np.square(inds), axis=0))
    d = np.sqrt(np.sum(sp.square(inds), axis=0))
    mask = d < r
    im_temp = im*mask
    return im_temp

@@ -341,7 +347,7 @@

341	347
342	348		Examples
343	349		--------
344		-	>>> import numpy as np
	350	+	>>> import scipy as sp
345	351		>>> from porespy.tools import extract_subsection
346	352		>>> im = np.array([[1, 1, 1, 1], [1, 2, 2, 2], [1, 2, 3, 3], [1, 2, 3, 4]])
347	353		>>> print(im)

@@ -505,9 +511,9 @@

505	511		Examples
506	512		--------
507	513		>>> import porespy as ps
508		-	>>> import numpy as np
509		-	>>> np.random.seed(0)
510		-	>>> im = np.random.randint(low=0, high=5, size=[4, 4])
	514	+	>>> import scipy as sp
	515	+	>>> sp.random.seed(0)
	516	+	>>> im = sp.random.randint(low=0, high=5, size=[4, 4])
511	517		>>> print(im)
512	518		[[4 0 3 3]
513	519		[3 1 3 2]

@@ -529,7 +535,7 @@

    keep_vals = np.array(keep_vals)
    swap_vals = ~np.in1d(im_flat, keep_vals)
    im_vals = np.unique(im_flat[swap_vals])
    new_vals = np.random.permutation(im_vals)
    new_vals = sp.random.permutation(im_vals)
    im_map = np.zeros(shape=[np.amax(im_vals) + 1, ], dtype=int)
    im_map[im_vals] = new_vals
    im_new = im_map[im_flat]

@@ -566,7 +572,7 @@

566	572		Example
567	573		-------
568	574		>>> import porespy as ps
569		-	>>> import numpy as np
	575	+	>>> import scipy as sp
570	576		>>> im = np.array([[0, 2, 9], [6, 8, 3]])
571	577		>>> im = ps.tools.make_contiguous(im)
572	578		>>> print(im)

@@ -632,7 +638,7 @@

    Examples
    --------
    >>> import porespy as ps
    >>> import numpy as np
    >>> import scipy as sp
    >>> mask = ps.tools.get_border(shape=[3, 3], mode='corners')
    >>> print(mask)
    [[ True False  True]

@@ -669,7 +675,7 @@

            border[0::, t:-t, 0::] = False
            border[0::, 0::, t:-t] = False
    if return_indices:
        border = np.where(border)
        border = sp.where(border)
    return border



@@ -732,7 +738,7 @@

    return im


def functions_to_table(mod, colwidth=[27, 48]):
def _functions_to_table(mod, colwidth=[27, 48]):
    r"""
    Given a module of functions, returns a ReST formatted text string that
    outputs a table when printed.

@@ -800,9 +806,9 @@

    """
    im = region
    if im.ndim != im.squeeze().ndim:
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
                      ' unexpected behavior.')
        warnings.warn(f'Input image conains a singleton axis: {im.shape}'
                      + ' Reduce dimensionality with np.squeeze(im) to avoid'
                      + ' unexpected behavior.')
    if strel is None:
        if region.ndim == 3:
            strel = ball(1)

@@ -842,7 +848,7 @@

    rad = int(np.ceil(radius))
    other = np.ones((2 * rad + 1, 2 * rad + 1), dtype=bool)
    other[rad, rad] = False
    disk = spim.distance_transform_edt(other) < radius
    disk = edt(other) < radius
    return disk



@@ -863,7 +869,7 @@

    rad = int(np.ceil(radius))
    other = np.ones((2 * rad + 1, 2 * rad + 1, 2 * rad + 1), dtype=bool)
    other[rad, rad, rad] = False
    ball = spim.distance_transform_edt(other) < radius
    ball = edt(other) < radius
    return ball



@@ -901,7 +907,7 @@

    return im1


def insert_sphere(im, c, r):
def insert_sphere(im, c, r, v=True, overwrite=True):
    r"""
    Inserts a sphere of a specified radius into a given image


@@ -913,26 +919,48 @@

        The [x, y, z] coordinate indicating the center of the sphere
    r : int
        The radius of sphere to insert
    v : int
        The value to put into the sphere voxels.  The default is ``True``
        which corresponds to inserting spheres into a Boolean image.  If
        a numerical value is given, ``im`` is converted to the same type as
        ``v``.
    overwrite : boolean
        If ``True`` (default) then the sphere overwrites whatever values are
        present in ``im``.  If ``False`` then the sphere values are only
        inserted into locations that are 0 or ``False``.

    Returns
    -------
    image : ND-array
        The original image with a sphere inerted at the specified location
    """
    # Convert image to same type os v for eventual insertion
    if im.dtype != type(v):
        im = im.astype(type(v))
    # Parse the arugments
    r = int(sp.around(r, decimals=0))
    if r == 0:
        return im
    c = np.array(c, dtype=int)
    if c.size != im.ndim:
        raise Exception('Coordinates do not match dimensionality of image')

    # Define a bounding box around inserted sphere, minding imaage boundaries
    bbox = []
    [bbox.append(np.clip(c[i] - r, 0, im.shape[i])) for i in range(im.ndim)]
    [bbox.append(np.clip(c[i] + r, 0, im.shape[i])) for i in range(im.ndim)]
    bbox = np.ravel(bbox)
    # Obtain slices into image
    s = bbox_to_slices(bbox)
    temp = im[s]
    blank = np.ones_like(temp)
    # Generate sphere template within image boundaries
    blank = np.ones_like(im[s], dtype=int)
    blank[tuple(c - bbox[0:im.ndim])] = 0
    blank = spim.distance_transform_edt(blank) < r
    im[s] = blank
    sph = edt(blank) < r
    if overwrite:  # Clear voxles under sphere to be zero
        temp = im[s]*sph > 0
        im[s][temp] = 0
    else:  # Clear portions of sphere to prevent overwriting
        sph *= im[s] == 0
    im[s] = im[s] + sph*v
    return im



@@ -983,7 +1011,7 @@

    else:
        xyz_line_in_template_coords = [xyz_line[i] - xyz_min[i] for i in range(3)]
        template[tuple(xyz_line_in_template_coords)] = 1
        template = spim.distance_transform_edt(template == 0) <= r
        template = edt(template == 0) <= r

    im[xyz_min[0]:xyz_max[0]+1,
       xyz_min[1]:xyz_max[1]+1,

@@ -1023,9 +1051,9 @@

    add_boundary_regions
    """
    if im.ndim != im.squeeze().ndim:
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
                      ' unexpected behavior.')
        warnings.warn(f'Input image conains a singleton axis: {im.shape}'
                      + ' Reduce dimensionality with np.squeeze(im) to avoid'
                      + ' unexpected behavior.')
    f = faces
    if f is not None:
        if im.ndim == 2:

@@ -1070,18 +1098,37 @@

    """
    # -------------------------------------------------------------------------
    # Get alias if provided by user
    phases_num = np.unique(im * 1)
    phases_num = np.unique(im).astype(int)
    phases_num = np.trim_zeros(phases_num)
    al = {}
    wrong_labels = []
    for values in phases_num:
        al[values] = 'phase{}'.format(values)
    if alias is not None:
        alias_sort = dict(sorted(alias.items()))
        phase_labels = np.array([*alias_sort])
        al = alias
        if set(phase_labels) != set(phases_num):
            raise Exception('Alias labels does not match with image labels '
                            'please provide correct image labels')
        for i in phase_labels:
            if i == 0:
                raise Exception("Label 0 is not allowed in alias. "
                                + "Please specify alias with a positive "
                                  "integer")
            elif i not in phases_num:
                wrong_labels.append(i)
        if wrong_labels:
            raise Exception("Alias label(s) {} does not "
                            "match with image "
                            "label(s).".format(wrong_labels)
                            + "Please provide correct labels from image.")
        if phase_labels.size < phases_num.size:
            missed_labels = np.setdiff1d(phases_num, phase_labels)
            for i in missed_labels:
                warnings.warn(
                    "label_{} alias is not provided although it "
                    "exists in the input image.".format(i)
                    + "The default label alias phase{} is assigned to "
                      "label_{}".format(i, i))
                al[i] = 'phase{}'.format(i)
    return al



@@ -1110,7 +1157,7 @@

        labels = [labels]
    s = spim.find_objects(regions)
    im_new = np.zeros_like(regions)
    x_min, y_min, z_min = np.inf, np.inf, np.inf
    x_min, y_min, z_min = sp.inf, sp.inf, sp.inf
    x_max, y_max, z_max = 0, 0, 0
    for i in labels:
        im_new[s[i-1]] = regions[s[i-1]] == i

@@ -1125,3 +1172,81 @@

1125	1172		bbox = bbox_to_slices([x_min, y_min, x_max, y_max])
1126	1173		im_new = im_new[bbox]
1127	1174		return im_new
	1175	+
	1176	+
	1177	+	def size_to_seq(size, bins=None):
	1178	+	r"""
	1179	+	Converts an image of invasion size values into sequence values.
	1180	+
	1181	+	This is meant to accept the output of the ``porosimetry`` function.
	1182	+
	1183	+	Parameters
	1184	+	----------
	1185	+	size : ND-image
	1186	+	The image containing invasion size values in each voxel.
	1187	+	bins : array_like or int (optional)
	1188	+	The bins to use when converting sizes to sequence. The default is
	1189	+	to create 1 bin for each unique value in ``size``. If an int
	1190	+	is supplied it is interpreted as the number of bins between 0 and the
	1191	+	maximum value in ``size``. If an array is supplied it is used as
	1192	+	the bins directly.
	1193	+
	1194	+	Returns
	1195	+	-------
	1196	+	seq : ND-image
	1197	+	An ND-image the same shape as ``size`` with invasion size values
	1198	+	replaced by the invasion sequence. This assumes that the invasion
	1199	+	process occurs via increasing pressure steps, such as produced by
	1200	+	the ``porosimetry`` function.
	1201	+
	1202	+	"""
	1203	+	solid = size == 0
	1204	+	if bins is None:
	1205	+	bins = np.unique(size)
	1206	+	elif isinstance(bins, int):
	1207	+	bins = np.linspace(0, size.max(), bins)
	1208	+	vals = np.digitize(size, bins=bins, right=True)
	1209	+	# Invert the vals so smallest size has largest sequence
	1210	+	vals = -(vals - vals.max() - 1)*~solid
	1211	+	# In case too many bins are given, remove empty ones
	1212	+	vals = make_contiguous(vals)
	1213	+
	1214	+	# Possibly simpler way?
	1215	+	# vals = (-(size - size.max())).astype(int) + 1
	1216	+	# vals[vals > size.max()] = 0
	1217	+
	1218	+	return vals
	1219	+
	1220	+
	1221	+	def seq_to_satn(seq):
	1222	+	r"""
	1223	+	Converts an image of invasion sequence values to saturation values.
	1224	+
	1225	+	Parameters
	1226	+	----------
	1227	+	seq : ND-image
	1228	+	The image containing invasion sequence values in each voxel.
	1229	+	Note that the invasion steps must be positive integers, solid voxels
	1230	+	indicated by 0, and uninvaded voxels indicated by -1.
	1231	+
	1232	+	Returns
	1233	+	-------
	1234	+	satn : ND-image
	1235	+	An ND-iamge the same size as ``seq`` but with sequnece values replaced
	1236	+	by the fraction of pores invaded at or below the sequence number.
	1237	+	Solid voxels and uninvaded voxels are represented by 0 and -1
	1238	+	respectively.
	1239	+
	1240	+	"""
	1241	+	seq = np.copy(seq).astype(int)
	1242	+	solid = seq == 0
	1243	+	uninvaded = seq == -1
	1244	+	seq = np.clip(seq, a_min=0, a_max=None)
	1245	+	seq = make_contiguous(seq)
	1246	+	b = np.bincount(seq.flatten())
	1247	+	b[0] = 0
	1248	+	c = np.cumsum(b)
	1249	+	satn = c[seq]/((seq > 0).sum() + uninvaded.sum())
	1250	+	satn[solid] = 0.0
	1251	+	satn[uninvaded] = -1.0
	1252	+	return satn

@@ -3,6 +3,81 @@

3	3		from mpl_toolkits.mplot3d.art3d import Poly3DCollection
4	4
5	5
	6	+	def bar(tup, h='pdf', **kwargs):
	7	+	r"""
	8	+	Convenience wrapper for matplotlib's ``bar``.
	9	+
	10	+	This automatically:
	11	+	* fetches the ``bin_centers``
	12	+	* fetches the bin heights from the specified ``h``
	13	+	* sets the bin widths
	14	+	* sets the edges to black
	15	+
	16	+	Parameters
	17	+	----------
	18	+	tup : named-tuple
	19	+	The named tuple returned by various functions in the
	20	+	``porespy.metrics`` submodule, such as ``chord_length_distribution``.
	21	+	h : str
	22	+	The value to use for bin heights. The default is ``pdf``, but
	23	+	``cdf`` is another option. Depending on the function the named-tuple
	24	+	may have different options.
	25	+	kwargs : keyword arguments
	26	+	All other keyword arguments are passed to ``bar``, including
	27	+	``edgecolor`` if you wish to overwrite the default black.
	28	+
	29	+	Returns
	30	+	-------
	31	+	fig: Matplotlib figure handle
	32	+	"""
	33	+	if 'edgecolor' not in kwargs:
	34	+	kwargs['edgecolor'] = 'k'
	35	+	fig = plt.bar(x=tup.bin_centers, height=getattr(tup, h),
	36	+	width=tup.bin_widths, **kwargs)
	37	+	return fig
	38	+
	39	+
	40	+	def imshow(*im, ind=None, axis=None):
	41	+	r"""
	42	+	Convenience wrapper for matplotlib's ``imshow``.
	43	+
	44	+	This automatically:
	45	+	* slices a 3D image in the middle of the last axis
	46	+	* uses a masked array to make 0's white
	47	+	* sets the origin to 'lower' so bottom-left corner is [0, 0]
	48	+
	49	+	Parameters
	50	+	----------
	51	+	im : ND-array
	52	+	The 2D or 3D image (or images) to show. If 2D then all other
	53	+	arguments are ignored.
	54	+	ind : int
	55	+	The slice to show if ``im`` is 3D. If not given then the middle of
	56	+	the image is used.
	57	+	axis : int
	58	+	The axis to show if ``im`` is 3D. If not given, then the last
	59	+	axis of the image is used, so an 'lower' slice is shown.
	60	+
	61	+	Note
	62	+	----
	63	+	``im`` can also be a series of unnamed arguments, in which case all
	64	+	received images will be shown using ``subplot``.
	65	+	"""
	66	+	if not isinstance(im, tuple):
	67	+	im = tuple([im])
	68	+	for i, image in enumerate(im):
	69	+	if image.ndim == 3:
	70	+	if axis is None:
	71	+	axis = 2
	72	+	if ind is None:
	73	+	ind = int(image.shape[axis]/2)
	74	+	image = image.take(indices=ind, axis=axis)
	75	+	image = np.ma.array(image, mask=image == 0)
	76	+	fig = plt.subplot(1, len(im), i+1)
	77	+	plt.imshow(image, origin='lower')
	78	+	return fig
	79	+
	80	+
6	81		def show_mesh(mesh):
7	82		r"""
8	83		Visualizes the mesh of a region as obtained by ``get_mesh`` function in

@@ -1,3 +1,4 @@

	1	+	import sys
1	2		import numpy as np
2	3		import openpnm as op
3	4		from tqdm import tqdm

@@ -35,7 +36,7 @@

    directly to an OpenPNM network object using the ``update`` command.

    """
    print('_'*60)
    print('-'*60)
    print('Extracting pore and throat information from image')
    from skimage.morphology import disk, ball
    struc_elem = disk if im.ndim == 2 else ball

@@ -67,7 +68,7 @@

    # dt_shape = np.array(dt.shape)

    # Start extracting size information for pores and throats
    for i in tqdm(Ps):
    for i in tqdm(Ps, file=sys.stdout):
        pore = i - 1
        if slices[pore] is None:
            continue

@@ -0,0 +1,103 @@

1	+	import numpy as np
2	+	import openpnm as op
3	+	from porespy.filters import trim_nonpercolating_paths
4	+	import collections
5	+
6	+
7	+	def tortuosity(im, axis, return_im=False, **kwargs):
8	+	r"""
9	+	Calculates tortuosity of given image in specified direction
10	+
11	+	Parameters
12	+	----------
13	+	im : ND-image
14	+	The binary image to analyze with ``True`` indicating phase of interest
15	+	axis : int
16	+	The axis along which to apply boundary conditions
17	+	return_im : boolean
18	+	If ``True`` then the resulting tuple contains a copy of the input
19	+	image with the concentration profile.
20	+
21	+	Returns
22	+	-------
23	+	results : tuple
24	+	A named-tuple containing:
25	+	* ``tortuosity`` calculated using the ``effective_porosity`` as
26	+	* ``effective_porosity`` of the image after applying
27	+	``trim_nonpercolating_paths``. This removes disconnected
28	+	voxels which cause singular matrices.
29	+	:math:`(D_{AB}/D_{eff}) \cdot \varepsilon`.
30	+	* ``original_porosity`` of the image as given
31	+	* ``formation_factor`` found as :math:`D_{AB}/D_{eff}`.
32	+	* ``image`` containing the concentration values from the simulation.
33	+	This is only returned if ``return_im`` is ``True``.
34	+
35	+	"""
36	+	if axis > (im.ndim - 1):
37	+	raise Exception("Axis argument is too high")
38	+	# Obtain original porosity
39	+	porosity_orig = im.sum()/im.size
40	+	# removing floating pores
41	+	im = trim_nonpercolating_paths(im, inlet_axis=axis, outlet_axis=axis)
42	+	# porosity is changed because of trimmimg floating pores
43	+	porosity_true = im.sum()/im.size
44	+	if porosity_true < porosity_orig:
45	+	print('Caution, True porosity is:', porosity_true,
46	+	'and volume fraction filled:',
47	+	abs(porosity_orig-porosity_true)*100, '%')
48	+	# cubic network generation
49	+	net = op.network.CubicTemplate(template=im, spacing=1)
50	+	# adding phase
51	+	water = op.phases.Water(network=net)
52	+	water['throat.diffusive_conductance'] = 1 # dummy value
53	+	# running Fickian Diffusion
54	+	fd = op.algorithms.FickianDiffusion(network=net, phase=water)
55	+	# choosing axis of concentration gradient
56	+	inlets = net['pore.coords'][:, axis] <= 1
57	+	outlets = net['pore.coords'][:, axis] >= im.shape[axis]-1
58	+	# boundary conditions on concentration
59	+	C_in = 1.0
60	+	C_out = 0.0
61	+	fd.set_value_BC(pores=inlets, values=C_in)
62	+	fd.set_value_BC(pores=outlets, values=C_out)
63	+	# Use specified solver if given
64	+	if 'solver_family' in kwargs.keys():
65	+	fd.settings.update(kwargs)
66	+	else: # Use pyamg otherwise, if presnet
67	+	try:
68	+	import pyamg
69	+	fd.settings['solver_family'] = 'pyamg'
70	+	except ModuleNotFoundError: # Use scipy cg as last resort
71	+	fd.settings['solver_family'] = 'scipy'
72	+	fd.settings['solver_type'] = 'cg'
73	+	op.utils.tic()
74	+	fd.run()
75	+	op.utils.toc()
76	+	# calculating molar flow rate, effective diffusivity and tortuosity
77	+	rate_out = fd.rate(pores=outlets)[0]
78	+	rate_in = fd.rate(pores=inlets)[0]
79	+	if not np.allclose(-rate_out, rate_in):
80	+	raise Exception('Something went wrong, inlet and outlet rate do not match')
81	+	delta_C = C_in - C_out
82	+	L = im.shape[axis]
83	+	A = np.prod(im.shape)/L
84	+	N_A = A/L*delta_C
85	+	Deff = rate_in/N_A
86	+	tau = porosity_true/(Deff)
87	+	result = collections.namedtuple('tortuosity_result', ['tortuosity',
88	+	'effective_porosity',
89	+	'original_porosity',
90	+	'formation_factor',
91	+	'image'])
92	+	result.tortuosity = tau
93	+	result.formation_factor = 1/Deff
94	+	result.original_porosity = porosity_orig
95	+	result.effective_porosity = porosity_true
96	+	if return_im:
97	+	conc = np.zeros([im.size, ], dtype=float)
98	+	conc[net['pore.template_indices']] = fd['pore.concentration']
99	+	conc = np.reshape(conc, newshape=im.shape)
100	+	result.image = conc
101	+	else:
102	+	result.image = None
103	+	return result

@@ -3,9 +3,9 @@


def set_mpl_style():

    sfont = 25
    mfont = 35
    lfont = 45
    sfont = 15
    mfont = 15
    lfont = 15

    line_props = {'linewidth': 4,
                  'markersize': 10}

@@ -85,32 +85,32 @@

    return new_im


def sem(im, direction='X'):
def sem(im, axis=0):
    r"""
    Simulates an SEM photograph looking into the porous material in the
    specified direction.  Features are colored according to their depth into
    the image, so darker features are further away.
    Simulates an SEM photograph looking into the porous material.

    Features are colored according to their depth into the image, so
    darker features are further away.

    Parameters
    ----------
    im : array_like
        ND-image of the porous material with the solid phase marked as 1 or
        True

    direction : string
        Specify the axis along which the camera will point.  Options are
        'X', 'Y', and 'Z'.
    axis : int
        Specifes the axis along which the camera will point.

    Returns
    -------
    image : 2D-array
        A 2D greyscale image suitable for use in matplotlib\'s ```imshow```
        A 2D greyscale image suitable for use in matplotlib's ``imshow``
        function.

    """
    im = np.array(~im, dtype=int)
    if direction in ['Y', 'y']:
    if axis == 1:
        im = np.transpose(im, axes=[1, 0, 2])
    if direction in ['Z', 'z']:
    if axis == 2:
        im = np.transpose(im, axes=[2, 1, 0])
    t = im.shape[0]
    depth = np.reshape(np.arange(0, t), [t, 1, 1])

@@ -119,12 +119,12 @@

    return im


def xray(im, direction='X'):
def xray(im, axis=0):
    r"""
    Simulates an X-ray radiograph looking through the porouls material in the
    specfied direction.  The resulting image is colored according to the amount
    of attenuation an X-ray would experience, so regions with more solid will
    appear darker.
    Simulates an X-ray radiograph looking through the porous material.

    The resulting image is colored according to the amount of attenuation an
    X-ray would experience, so regions with more solid will appear darker.

    Parameters
    ----------

@@ -132,9 +132,8 @@

        ND-image of the porous material with the solid phase marked as 1 or
        True

    direction : string
        Specify the axis along which the camera will point.  Options are
        'X', 'Y', and 'Z'.
    axis : int
        Specifes the axis along which the camera will point.

    Returns
    -------

@@ -143,9 +142,9 @@

143	142		function.
144	143		"""
145	144		im = np.array(~im, dtype=int)
146		-	if direction in ['Y', 'y']:
	145	+	if axis == 1:
147	146		im = np.transpose(im, axes=[1, 0, 2])
148		-	if direction in ['Z', 'z']:
	147	+	if axis == 2:
149	148		im = np.transpose(im, axes=[2, 1, 0])
150	149		im = np.sum(im, axis=0)
151	150		return im

@@ -391,23 +391,23 @@

    num = [0, *num]
    phases_num = np.unique(im * 1)
    phases_num = np.trim_zeros(phases_num)
    for i in phases_num:
        loc1 = np.logical_and(conns1 >= num[i - 1], conns1 < num[i])
        loc2 = np.logical_and(conns2 >= num[i - 1], conns2 < num[i])
        loc3 = np.logical_and(label >= num[i - 1], label < num[i])
        net['throat.{}'.format(al[i])] = loc1 * loc2
        net['pore.{}'.format(al[i])] = loc3
        if i == phases_num[-1]:
    for i0, i1 in enumerate(phases_num):
        loc1 = np.logical_and(conns1 >= num[i0], conns1 < num[i0 + 1])
        loc2 = np.logical_and(conns2 >= num[i0], conns2 < num[i0 + 1])
        loc3 = np.logical_and(label >= num[i0], label < num[i0 + 1])
        net['throat.{}'.format(al[i1])] = loc1 * loc2
        net['pore.{}'.format(al[i1])] = loc3
        if i1 == phases_num[-1]:
            loc4 = np.logical_and(conns1 < num[-1], conns2 >= num[-1])
            loc5 = label >= num[-1]
            net['throat.boundary'] = loc4
            net['pore.boundary'] = loc5
        for j in phases_num:
            if j > i:
                pi_pj_sa = np.zeros_like(label)
                loc6 = np.logical_and(conns2 >= num[j - 1], conns2 < num[j])
        for j0, j1 in enumerate(phases_num):
            if j0 > i0:
                pi_pj_sa = np.zeros_like(label, dtype=float)
                loc6 = np.logical_and(conns2 >= num[j0], conns2 < num[j0 + 1])
                pi_pj_conns = loc1 * loc6
                net['throat.{}_{}'.format(al[i], al[j])] = pi_pj_conns
                net['throat.{}_{}'.format(al[i1], al[j1])] = pi_pj_conns
                if any(pi_pj_conns):
                    # ---------------------------------------------------------
                    # Calculates phase[i] interfacial area that connects with

@@ -434,9 +434,7 @@

                        s_pa_c = np.trim_zeros(s_pa_c)
                        pi_pj_sa[i_index] = p_sa_c
                        pi_pj_sa[j_index] = s_pa_c
                    net['pore.{}_{}_area'.format(al[i],
                                                 al[j])] = (pi_pj_sa *
                                                            voxel_size ** 2)
                    net[f'pore.{al[i1]}_{al[j1]}_area'] = pi_pj_sa * voxel_size ** 2
    return net



@@ -479,10 +477,12 @@

            dic['bottom'] = 1
        for i in f:
            if i in ['left', 'front', 'bottom']:
                network['pore.{}'.format(i)] = (coords[:, dic[i]] <
                                                min(condition[:, dic[i]]))
                network['pore.{}'.format(i)] = (
                    coords[:, dic[i]] < min(condition[:, dic[i]])
                )
            elif i in ['right', 'back', 'top']:
                network['pore.{}'.format(i)] = (coords[:, dic[i]] >
                                                max(condition[:, dic[i]]))
                network['pore.{}'.format(i)] = (
                    coords[:, dic[i]] > max(condition[:, dic[i]])
                )

    return network

@@ -1,10 +1,11 @@

import sys
import numpy as np
import scipy.ndimage as spim
from tqdm import tqdm
import scipy.ndimage as spim
from porespy.tools import extract_subsection, bbox_to_slices
from skimage.measure import regionprops
from skimage.measure import mesh_surface_area, marching_cubes_lewiner
from skimage.morphology import skeletonize_3d, ball
from skimage.measure import regionprops
from pandas import DataFrame



@@ -176,53 +177,54 @@

    which may be helpful.

    """
    print('_'*60)
    print('-'*60)
    print('Calculating regionprops')

    results = regionprops(im, coordinates='xy')
    for i in tqdm(range(len(results))):
        mask = results[i].image
        mask_padded = np.pad(mask, pad_width=1, mode='constant')
        temp = spim.distance_transform_edt(mask_padded)
        dt = extract_subsection(temp, shape=mask.shape)
        # ---------------------------------------------------------------------
        # Slice indices
        results[i].slice = results[i]._slice
        # ---------------------------------------------------------------------
        # Volume of regions in voxels
        results[i].volume = results[i].area
        # ---------------------------------------------------------------------
        # Volume of bounding box, in voxels
        results[i].bbox_volume = np.prod(mask.shape)
        # ---------------------------------------------------------------------
        # Create an image of the border
        results[i].border = dt == 1
        # ---------------------------------------------------------------------
        # Create an image of the maximal inscribed sphere
        r = dt.max()
        inv_dt = spim.distance_transform_edt(dt < r)
        results[i].inscribed_sphere = inv_dt < r
        # ---------------------------------------------------------------------
        # Find surface area using marching cubes and analyze the mesh
        tmp = np.pad(np.atleast_3d(mask), pad_width=1, mode='constant')
        tmp = spim.convolve(tmp, weights=ball(1))/5
        verts, faces, norms, vals = marching_cubes_lewiner(volume=tmp, level=0)
        results[i].surface_mesh_vertices = verts
        results[i].surface_mesh_simplices = faces
        area = mesh_surface_area(verts, faces)
        results[i].surface_area = area
        # ---------------------------------------------------------------------
        # Find sphericity
        vol = results[i].volume
        r = (3/4/np.pi*vol)**(1/3)
        a_equiv = 4*np.pi*(r)**2
        a_region = results[i].surface_area
        results[i].sphericity = a_equiv/a_region
        # ---------------------------------------------------------------------
        # Find skeleton of region
        results[i].skeleton = skeletonize_3d(mask)
        # ---------------------------------------------------------------------
        # Volume of convex image, equal to area in 2D, so just translating
        results[i].convex_volume = results[i].convex_area
    results = regionprops(im)
    with tqdm(range(len(results))) as pbar:
        for i in range(len(results)):
            pbar.update()
            mask = results[i].image
            mask_padded = np.pad(mask, pad_width=1, mode='constant')
            temp = spim.distance_transform_edt(mask_padded)
            dt = extract_subsection(temp, shape=mask.shape)
            # Slice indices
            results[i].slice = results[i]._slice
            # ----------------------------------------------------------------
            # Volume of regions in voxels
            results[i].volume = results[i].area
            # ----------------------------------------------------------------
            # Volume of bounding box, in voxels
            results[i].bbox_volume = np.prod(mask.shape)
            # ----------------------------------------------------------------
            # Create an image of the border
            results[i].border = dt == 1
            # ----------------------------------------------------------------
            # Create an image of the maximal inscribed sphere
            r = dt.max()
            inv_dt = spim.distance_transform_edt(dt < r)
            results[i].inscribed_sphere = inv_dt < r
            # ----------------------------------------------------------------
            # Find surface area using marching cubes and analyze the mesh
            tmp = np.pad(np.atleast_3d(mask), pad_width=1, mode='constant')
            tmp = spim.convolve(tmp, weights=ball(1))/5
            verts, faces, norms, vals = marching_cubes_lewiner(volume=tmp, level=0)
            results[i].surface_mesh_vertices = verts
            results[i].surface_mesh_simplices = faces
            area = mesh_surface_area(verts, faces)
            results[i].surface_area = area
            # ----------------------------------------------------------------
            # Find sphericity
            vol = results[i].volume
            r = (3/4/np.pi*vol)**(1/3)
            a_equiv = 4*np.pi*(r)**2
            a_region = results[i].surface_area
            results[i].sphericity = a_equiv/a_region
            # ----------------------------------------------------------------
            # Find skeleton of region
            results[i].skeleton = skeletonize_3d(mask)
            # ----------------------------------------------------------------
            # Volume of convex image, equal to area in 2D, so just translating
            results[i].convex_volume = results[i].convex_area

    return results

@@ -1,10 +1,15 @@

import porespy as ps
import numpy as np
from edt import edt
import porespy as ps
import scipy.spatial as sptl
import scipy.ndimage as spim
from porespy.tools import norm_to_uniform, ps_ball, ps_disk
from numba import njit, jit
from porespy.tools import norm_to_uniform, ps_ball, ps_disk, get_border
from porespy.tools import insert_sphere, fftmorphology
from typing import List

from numpy import array
from tqdm import tqdm
from skimage.morphology import ball, disk


def insert_shape(im, element, center=None, corner=None, value=1,

@@ -56,8 +61,8 @@

        for dim in range(im.ndim):
            r, d = np.divmod(element.shape[dim], 2)
            if d == 0:
                raise Exception('Cannot specify center point when element ' +
                                'has one or more even dimension')
                raise Exception('Cannot specify center point when element '
                                + 'has one or more even dimension')
            lower_im = np.amax((center[dim] - r, 0))
            upper_im = np.amin((center[dim] + r + 1, im.shape[dim]))
            s_im.append(slice(lower_im, upper_im))

@@ -87,13 +92,17 @@

    return im


def RSA(im: np.array, radius: int, volume_fraction: int = 1,
        mode: str = 'extended'):
def RSA(im: array, radius: int, volume_fraction: int = 1, n_max: int = None,
        mode: str = 'contained'):
    r"""
    Generates a sphere or disk packing using Random Sequential Addition

    This which ensures that spheres do not overlap but does not guarantee they
    are tightly packed.
    This algorithm ensures that spheres do not overlap but does not
    guarantee they are tightly packed.

    This function adds spheres to the background of the received ``im``, which
    allows iteratively adding spheres of different radii to the unfilled space,
    be repeatedly passing in the result of previous calls to RSA.

    Parameters
    ----------

@@ -101,97 +110,176 @@

        The image into which the spheres should be inserted.  By accepting an
        image rather than a shape, it allows users to insert spheres into an
        already existing image.  To begin the process, start with an array of
        zero such as ``im = np.zeros([200, 200], dtype=bool)``.
        zeros such as ``im = np.zeros([200, 200, 200], dtype=bool)``.
    radius : int
        The radius of the disk or sphere to insert.
    volume_fraction : scalar
    volume_fraction : scalar (default is 1.0)
        The fraction of the image that should be filled with spheres.  The
        spheres are addeds 1's, so each sphere addition increases the
        ``volume_fraction`` until the specified limit is reach.
    mode : string
        spheres are added as 1's, so each sphere addition increases the
        ``volume_fraction`` until the specified limit is reach.  Note that if
        ``n_max`` is reached first, then ``volume_fraction`` will not be
        acheived.
    n_max : int (default is 10,000)
        The maximum number of spheres to add.  By default the value of
        ``n_max`` is high so that the addition of spheres will go indefinately
        until ``volume_fraction`` is met, but specifying a smaller value
        will halt addition after the given number of spheres are added.
    mode : string (default is 'contained')
        Controls how the edges of the image are handled.  Options are:

        'extended' - Spheres are allowed to extend beyond the edge of the image

        'contained' - Spheres are all completely within the image

        'periodic' - The portion of a sphere that extends beyond the image is
        inserted into the opposite edge of the image (Not Implemented Yet!)
        'extended' - Spheres are allowed to extend beyond the edge of the
        image.  In this mode the volume fraction will be less that requested
        since some spheres extend beyond the image, but their entire volume
        is counted as added for computational efficiency.

    Returns
    -------
    image : ND-array
        A copy of ``im`` with spheres of specified radius *added* to the
        background.
        A handle to the input ``im`` with spheres of specified radius
        *added* to the background.

    Notes
    -----
    Each sphere is filled with 1's, but the center is marked with a 2.  This
    allows easy boolean masking to extract only the centers, which can be
    converted to coordinates using ``scipy.where`` and used for other purposes.
    The obtain only the spheres, use``im = im == 1``.

    This function adds spheres to the background of the received ``im``, which
    allows iteratively adding spheres of different radii to the unfilled space.
    This function uses Numba to speed up the search for valid sphere insertion
    points.  It seems that Numba does not look at the state of the scipy
    random number generator, so setting the seed to a known value has no
    effect on the output of this function. Each call to this function will
    produce a unique image.  If you wish to use the same realization multiple
    times you must save the array (e.g. ``numpy.save``).

    References
    ----------
    [1] Random Heterogeneous Materials, S. Torquato (2001)

    """
    # Note: The 2D vs 3D splitting of this just me being lazy...I can't be
    # bothered to figure it out programmatically right now
    # TODO: Ideally the spheres should be added periodically
    print(78*'―')
    print('RSA: Adding spheres of size ' + str(radius))
    d2 = len(im.shape) == 2
    mrad = 2*radius
    if d2:
        im_strel = ps_disk(radius)
        mask_strel = ps_disk(mrad)
    else:
        im_strel = ps_ball(radius)
        mask_strel = ps_ball(mrad)
    if np.any(im > 0):
        # Dilate existing objects by im_strel to remove pixels near them
        # from consideration for sphere placement
        mask = ps.tools.fftmorphology(im > 0, im_strel > 0, mode='dilate')
        mask = mask.astype(int)
    print(80*'-')
    print(f'RSA: Adding spheres of size {radius}')
    im = im.astype(bool)
    if n_max is None:
        n_max = 10000
    vf_final = volume_fraction
    vf_start = im.sum()/im.size
    print('Initial volume fraction:', vf_start)
    if im.ndim == 2:
        template_lg = ps_disk(radius*2)
        template_sm = ps_disk(radius)
    else:
        mask = np.zeros_like(im)
        template_lg = ps_ball(radius*2)
        template_sm = ps_ball(radius)
    vf_template = template_sm.sum()/im.size
    # Pad image by the radius of large template to enable insertion near edges
    im = np.pad(im, pad_width=2*radius, mode='edge')
    # Depending on mode, adjust mask to remove options around edge
    if mode == 'contained':
        mask = _remove_edge(mask, radius)
        border = get_border(im.shape, thickness=2*radius, mode='faces')
    elif mode == 'extended':
        pass
    elif mode == 'periodic':
        raise Exception('Periodic edges are not implemented yet')
        border = get_border(im.shape, thickness=radius+1, mode='faces')
    else:
        raise Exception('Unrecognized mode: ' + mode)
    vf = im.sum()/im.size
    free_spots = np.argwhere(mask == 0)
        raise Exception('Unrecognized mode: ', mode)
    # Remove border pixels
    im[border] = True
    # Dilate existing objects by strel to remove pixels near them
    # from consideration for sphere placement
    print('Dilating foreground features by sphere radius')
    dt = edt(im == 0)
    options_im = (dt >= radius)
    # ------------------------------------------------------------------------
    # Begin inserting the spheres
    vf = vf_start
    free_sites = np.flatnonzero(options_im)
    i = 0
    while vf <= volume_fraction and len(free_spots) > 0:
        choice = np.random.randint(0, len(free_spots), size=1)
        if d2:
            [x, y] = free_spots[choice].flatten()
            im = _fit_strel_to_im_2d(im, im_strel, radius, x, y)
            mask = _fit_strel_to_im_2d(mask, mask_strel, mrad, x, y)
            im[x, y] = 2
        else:
            [x, y, z] = free_spots[choice].flatten()
            im = _fit_strel_to_im_3d(im, im_strel, radius, x, y, z)
            mask = _fit_strel_to_im_3d(mask, mask_strel, mrad, x, y, z)
            im[x, y, z] = 2
        free_spots = np.argwhere(mask == 0)
        vf = im.sum()/im.size
    while (vf <= vf_final) and (i < n_max) and (len(free_sites) > 0):
        c, count = _make_choice(options_im, free_sites=free_sites)
        # The 100 below is arbitrary and may change performance
        if count > 100:
            # Regenerate list of free_sites
            print('Regenerating free_sites after', i, 'iterations')
            free_sites = np.flatnonzero(options_im)
        if all(np.array(c) == -1):
            break
        s_sm = tuple([slice(x - radius, x + radius + 1, None) for x in c])
        s_lg = tuple([slice(x - 2*radius, x + 2*radius + 1, None) for x in c])
        im[s_sm] += template_sm  # Add ball to image
        options_im[s_lg][template_lg] = False  # Update extended region
        vf += vf_template
        i += 1
    if vf > volume_fraction:
        print('Volume Fraction', volume_fraction, 'reached')
    if len(free_spots) == 0:
        print('No more free spots', 'Volume Fraction', vf)
    print('Number of spheres inserted:', i)
    # ------------------------------------------------------------------------
    # Get slice into returned image to retain original size
    s = tuple([slice(2*radius, d-2*radius, None) for d in im.shape])
    im = im[s]
    vf = im.sum()/im.size
    print('Final volume fraction:', vf)
    return im


@njit
def _make_choice(options_im, free_sites):
    r"""
    This function is called by _begin_inserting to find valid insertion points

    Parameters
    ----------
    options_im : ND-array
        An array with ``True`` at all valid locations and ``False`` at all
        locations where a sphere already exists PLUS a region of radius R
        around each sphere since these points are also invalid insertion
        points.
    free_sites : array_like
        A 1D array containing valid insertion indices.  This list is used to
        select insertion points from a limited which occasionally gets
        smaller.

    Returns
    -------
    coords : list
        The XY or XYZ coordinates of the next insertion point
    count : int
        The number of attempts that were needed to find valid point.  If
        this value gets too high, a short list of ``free_sites`` should be
        generated in the calling function.

    """
    choice = False
    count = 0
    upper_limit = len(free_sites)
    max_iters = upper_limit*20
    if options_im.ndim == 2:
        coords = [-1, -1]
        Nx, Ny = options_im.shape
        while not choice:
            if count >= max_iters:
                coords = [-1, -1]
                break
            ind = np.random.randint(0, upper_limit)
            # This numpy function is not supported by numba yet
            # c1, c2 = np.unravel_index(free_sites[ind], options_im.shape)
            # So using manual unraveling
            coords[1] = free_sites[ind] % Ny
            coords[0] = (free_sites[ind] // Ny) % Nx
            choice = options_im[coords[0], coords[1]]
            count += 1
    if options_im.ndim == 3:
        coords = [-1, -1, -1]
        Nx, Ny, Nz = options_im.shape
        while not choice:
            if count >= max_iters:
                coords = [-1, -1, -1]
                break
            ind = np.random.randint(0, upper_limit)
            # This numpy function is not supported by numba yet
            # c1, c2, c3 = np.unravel_index(free_sites[ind], options_im.shape)
            # So using manual unraveling
            coords[2] = free_sites[ind] % Nz
            coords[1] = (free_sites[ind] // Nz) % Ny
            coords[0] = (free_sites[ind] // (Nz * Ny)) % Nx
            choice = options_im[coords[0], coords[1], coords[2]]
            count += 1
    return coords, count


def bundle_of_tubes(shape: List[int], spacing: int):
    r"""
    Create a 3D image of a bundle of tubes, in the form of a rectangular

@@ -351,7 +439,7 @@

        if np.all(pts >= 0) and np.all(pts < im.shape):
            line_pts = line_segment(pts[0], pts[1])
            im[tuple(line_pts)] = True
    im = spim.distance_transform_edt(~im) > radius
    im = edt(~im) > radius
    return im



@@ -492,7 +580,7 @@

492	580		s+r:im.shape[1]-r:2*s,
493	581		s:im.shape[2]-r:2*s]
494	582		im[coords[0], coords[1], coords[2]] = 1
495		-	im = ~(spim.distance_transform_edt(~im) < r)
	583	+	im = ~(edt(~im) < r)
496	584		return im
497	585
498	586

@@ -543,8 +631,12 @@

    im = np.random.random(size=shape)

    # Helper functions for calculating porosity: phi = g(f(N))
    f = lambda N: spim.distance_transform_edt(im > N/bulk_vol) < radius
    g = lambda im: 1 - im.sum() / np.prod(shape)
    def f(N):
        return edt(im > N/bulk_vol) < radius

    def g(im):
        r"""Returns fraction of 0s, given a binary image"""
        return 1 - im.sum() / np.prod(shape)

    # # Newton's method for getting image porosity match the given
    # w = 1.0                         # Damping factor

@@ -574,92 +666,187 @@

    return ~f(N)


def generate_noise(shape: List[int], porosity=None, octaves: int = 3,
                   frequency: int = 32, mode: str = 'simplex'):
def perlin_noise(shape: List[int], porosity=None, octaves: int = 3,
                 frequency: List[int] = 2, persistence: float = 0.5):
    r"""
    Generate a field of spatially correlated random noise using the Perlin
    noise algorithm, or the updated Simplex noise algorithm.
    Generate a Perlin noise field

    Parameters
    ----------
    shape : array_like
        The size of the image to generate in [Nx, Ny, Nz] where N is the
        number of voxels.

        The shape of the desired image
    frequncy : array_like
        Controls the frequency of the noise, with higher values leading to
        smaller features or more tightly spaced undulations in the brightness.
    porosity : float
        If specified, this will threshold the image to the specified value
        prior to returning.  If no value is given (the default), then the
        scalar noise field is returned.

        If specified, the returned image will be thresholded to the specified
        porosity.  If not provided, the greyscale noise is returned (default).
    octaves : int
        Controls the *texture* of the noise, with higher octaves giving more
        complex features over larger length scales.

    frequency : array_like
        Controls the relative sizes of the features, with higher frequencies
        giving larger features.  A scalar value will apply the same frequency
        in all directions, given an isotropic field; a vector value will
        apply the specified values along each axis to create anisotropy.

    mode : string
        Which noise algorithm to use, either ``'simplex'`` (default) or
        ``'perlin'``.
        Controls the texture of the noise, with higher values giving more
        comlex features of larger length scales.
    persistence : float
        Controls how prominent each successive octave is.  Shoul be a number
        less than 1.

    Returns
    -------
    image : ND-array
        If porosity is given, then a boolean array with ``True`` values
        denoting the pore space is returned.  If not, then normally
        distributed and spatially correlated randomly noise is returned.
    An ND-array of the specified ``shape``.  If ``porosity`` is not given
    then the array contains greyscale values distributed normally about 0.
    Use ``porespy.tools.norm_to_uniform`` to create an well-scale image for
    thresholding.  If ``porosity`` is given then these steps are done
    internally and a boolean image is returned.

    Notes
    -----
    This method depends the a package called 'noise' which must be
    compiled. It is included in the Anaconda distribution, or a platform
    specific binary can be downloaded.
    The implementation used here is a bit fussy about the values of
    ``frequency`` and ``octaves``.  (1) the image ``shape`` must an integer
    multiple of ``frequency`` in each direction, and (2) ``frequency`` to the
    power of ``octaves`` must be less than or equal the``shape`` in each
    direction.  Exceptions are thrown if these conditions are not met.

    See Also
    --------
    porespy.tools.norm_to_uniform
    References
    ----------
    This implementation is taken from Pierre Vigier's
    `Github repo <https://github.com/pvigier/perlin-numpy>`_

    """
    try:
        import noise
    except ModuleNotFoundError:
        raise Exception("The noise package must be installed")
    # Parse args
    shape = np.array(shape)
    if np.size(shape) == 1:
        Lx, Ly, Lz = np.full((3, ), int(shape))
    elif len(shape) == 2:
        Lx, Ly = shape
        Lz = 1
    elif len(shape) == 3:
        Lx, Ly, Lz = shape
    if mode == 'simplex':
        f = noise.snoise3
    else:
        f = noise.pnoise3
    frequency = np.atleast_1d(frequency)
    if frequency.size == 1:
        freq = np.full(shape=[3, ], fill_value=frequency[0])
    elif frequency.size == 2:
        freq = np.concatenate((frequency, [1]))
    else:
        freq = np.array(frequency)
    im = np.zeros(shape=[Lx, Ly, Lz], dtype=float)
    for x in range(Lx):
        for y in range(Ly):
            for z in range(Lz):
                im[x, y, z] = f(x=x/freq[0], y=y/freq[1], z=z/freq[2],
                                octaves=octaves)
    im = im.squeeze()
    if porosity:
        im = norm_to_uniform(im, scale=[0, 1])
        im = im < porosity
    return im


def blobs(shape: List[int], porosity: float = 0.5, blobiness: int = 1):
    if shape.size == 1:  # Assume 3D
        shape = np.ones(3, dtype=int)*shape
    res = np.array(frequency)
    if res.size == 1:  # Assume shape as shape
        res = np.ones(shape.size, dtype=int)*res

    # Check inputs for various sins
    if res.size != shape.size:
        raise Exception('shape and res must have same dimensions')
    if np.any(np.mod(shape, res) > 0):
        raise Exception('res must be a multiple of shape along each axis')
    if np.any(shape/res**octaves < 1):
        raise Exception('(res[i])**octaves must be <= shape[i]')
    check = shape/(res**octaves)
    if np.any(check % 1):
        raise Exception("Image size must be factor of res**octaves")

    # Generate noise
    noise = np.zeros(shape)
    frequency = 1
    amplitude = 1
    for _ in tqdm(range(octaves)):
        if noise.ndim == 2:
            noise += amplitude * _perlin_noise_2D(shape, frequency*res)
        elif noise.ndim == 3:
            noise += amplitude * _perlin_noise_3D(shape, frequency*res)
        frequency *= 2
        amplitude *= persistence

    if porosity is not None:
        noise = norm_to_uniform(noise, scale=[0, 1])
        noise = noise > porosity

    return noise




366	410		return combined_region
367	411
368	412
369		-	def find_peaks(dt, r_max=4, footprint=None):
	413	+	def find_peaks(dt, r_max=4, footprint=None, **kwargs):
370	414		r"""
371		-	Returns all local maxima in the distance transform
	415	+	Finds local maxima in the distance transform
372	416
373	417		Parameters
374	418		----------
375	419		dt : ND-array
376	420		The distance transform of the pore space. This may be calculated and
377	421		filtered using any means desired.
378		-
379	422		r_max : scalar
380	423		The size of the structuring element used in the maximum filter. This
381	424		controls the localness of any maxima. The default is 4 voxels.
382		-
383	425		footprint : ND-array
384	426		Specifies the shape of the structuring element used to define the
385		-	neighborhood when looking for peaks. If none is specified then a
386		-	spherical shape is used (or circular in 2D).
	427	+	neighborhood when looking for peaks. If ``None`` (the default) is
	428	+	specified then a spherical shape is used (or circular in 2D).
387	429
388	430		Returns
389	431		-------

414	457		footprint = ball
415	458		else:
416	459		raise Exception("only 2-d and 3-d images are supported")
417		-	mx = spim.maximum_filter(dt + 2*(~im), footprint=footprint(r_max))
418		-	peaks = (dt == mx)*im
	460	+	parallel = kwargs.pop('parallel', False)
	461	+	cores = kwargs.pop('cores', None)
	462	+	divs = kwargs.pop('cores', 2)
	463	+	if parallel:
	464	+	overlap = max(footprint(r_max).shape)
	465	+	peaks = chunked_func(func=find_peaks, overlap=overlap,
	466	+	im_arg='dt', dt=dt, footprint=footprint,
	467	+	cores=cores, divs=divs)
	468	+	else:
	469	+	mx = spim.maximum_filter(dt + 2 * (~im), footprint=footprint(r_max))
	470	+	peaks = (dt == mx) * im
419	471		return peaks
420	472
421	473

427	479		Parameters
428	480		----------
429	481		peaks : ND-image
430		-	An image containing True values indicating peaks in the distance
	482	+	An image containing ``True`` values indicating peaks in the distance
431	483		transform
432	484
433	485		Returns
434	486		-------
435	487		image : ND-array
436	488		An array with the same number of isolated peaks as the original image,
437		-	but fewer total voxels.
	489	+	but fewer total ``True`` voxels.
438	490
439	491		Notes
440	492		-----

597	650		im : ND-image
598	651		A Boolean image, with True values indicating the phase for which
599	652		disconnected voxels are sought.
600		-
601	653		conn : int
602	654		For 2D the options are 4 and 8 for square and diagonal neighbors, while
603	655		for the 3D the options are 6 and 26, similarily for square and diagonal
604		-	neighbors. The default is max
	656	+	neighbors. The default is the maximum option.
605	657
606	658		Returns
607	659		-------

618	670
619	671		"""
620	672		if im.ndim != im.squeeze().ndim:
621		-	warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
622		-	' Reduce dimensionality with np.squeeze(im) to avoid' +
623		-	' unexpected behavior.')
	673	+	warnings.warn(
	674	+	"Input image conains a singleton axis:"
	675	+	+ str(im.shape)
	676	+	+ " Reduce dimensionality with np.squeeze(im) to avoid"
	677	+	+ " unexpected behavior."
	678	+	)
624	679		if im.ndim == 2:
625	680		if conn == 4:
626	681		strel = disk(1)
627	682		elif conn in [None, 8]:
628	683		strel = square(3)
	684	+	else:
	685	+	raise Exception("Received conn is not valid")
629	686		elif im.ndim == 3:
630	687		if conn == 6:
631	688		strel = ball(1)
632	689		elif conn in [None, 26]:
633	690		strel = cube(3)
	691	+	else:
	692	+	raise Exception("Received conn is not valid")
634	693		labels, N = spim.label(input=im, structure=strel)
635	694		holes = clear_border(labels=labels) > 0
636	695		return holes
637	696
638	697
639		-	def fill_blind_pores(im):
	698	+	def fill_blind_pores(im, conn=None):
640	699		r"""
641	700		Fills all pores that are not connected to the edges of the image.
642	701

1		-	from collections import namedtuple
	1	+	import sys
	2	+	import dask
	3	+	import dask.array as da
	4	+	from dask.diagnostics import ProgressBar
	5	+	import warnings
2	6		import numpy as np
	7	+	from numba import njit, prange
	8	+	from edt import edt
3	9		import operator as op
	10	+	from tqdm import tqdm
4	11		import scipy.ndimage as spim
5	12		import scipy.spatial as sptl
6		-	import warnings
	13	+	from collections import namedtuple
7	14		from scipy.signal import fftconvolve
8		-	from tqdm import tqdm
9		-	from numba import jit
10		-	from skimage.segmentation import clear_border
	15	+	from skimage.segmentation import clear_border, watershed
11	16		from skimage.morphology import ball, disk, square, cube, diamond, octahedron
12		-	from skimage.morphology import reconstruction, watershed
	17	+	from skimage.morphology import reconstruction
13	18		from porespy.tools import randomize_colors, fftmorphology
14	19		from porespy.tools import get_border, extend_slice, extract_subsection
15		-	from porespy.tools import ps_disk, ps_ball
16	20		from porespy.tools import _create_alias_map
	21	+	from porespy.tools import ps_disk, ps_ball
	22	+
	23	+
	24	+	def apply_padded(im, pad_width, func, pad_val=1, **kwargs):
	25	+	r"""
	26	+	Applies padding to an image before sending to ``func``, then extracts the
	27	+	result corresponding to the original image shape.
	28	+
	29	+	Parameters
	30	+	----------
	31	+	im : ND-image
	32	+	The image to which ``func`` should be applied
	33	+	pad_width : int or list of ints
	34	+	The amount of padding to apply to each axis. Refer to ``numpy.pad``
	35	+	documentation for more details.
	36	+	pad_val : scalar
	37	+	The value to place into the padded voxels. The default is 1 (or
	38	+	``True``) which extends the pore space.
	39	+	func : function handle
	40	+	The function to apply to the padded image
	41	+	kwargs : additional keyword arguments
	42	+	All additional keyword arguments are collected and passed to ``func``.
	43	+
	44	+	Notes
	45	+	-----
	46	+	A use case for this is when using ``skimage.morphology.skeletonize_3d``
	47	+	to ensure that the skeleton extends beyond the edges of the image.
	48	+	"""
	49	+	padded = np.pad(im, pad_width=pad_width,
	50	+	mode='constant', constant_values=pad_val)
	51	+	temp = func(padded, **kwargs)
	52	+	result = extract_subsection(im=temp, shape=im.shape)
	53	+	return result
17	54
18	55
19	56		def trim_small_clusters(im, size=1):
20	57		r"""
21		-	Remove isolated voxels or clusters smaller than a given size
	58	+	Remove isolated voxels or clusters of a given size or smaller
22	59
23	60		Parameters
24	61		----------

36	73		``size`` removed.
37	74
38	75		"""
39		-	if im.dims == 2:
	76	+	if im.ndim == 2:
40	77		strel = disk(1)
41		-	elif im.ndims == 3:
	78	+	elif im.ndim == 3:
42	79		strel = ball(1)
43	80		else:
44		-	raise Exception('Only 2D or 3D images are accepted')
	81	+	raise Exception("Only 2D or 3D images are accepted")
45	82		filtered_array = np.copy(im)
46	83		labels, N = spim.label(filtered_array, structure=strel)
47	84		id_sizes = np.array(spim.sum(im, labels, range(N + 1)))
48		-	area_mask = (id_sizes <= size)
	85	+	area_mask = id_sizes <= size
49	86		filtered_array[area_mask[labels]] = 0
50	87		return filtered_array
51	88

76	113		return result
77	114
78	115
79		-	def distance_transform_lin(im, axis=0, mode='both'):
	116	+	def distance_transform_lin(im, axis=0, mode="both"):
80	117		r"""
81	118		Replaces each void voxel with the linear distance to the nearest solid
82	119		voxel along the specified axis.

110	147		the nearest background along the specified axis.
111	148		"""
112	149		if im.ndim != im.squeeze().ndim:
113		-	warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
114		-	' Reduce dimensionality with np.squeeze(im) to avoid' +
115		-	' unexpected behavior.')
116		-	if mode in ['backward', 'reverse']:
	150	+	warnings.warn(
	151	+	"Input image conains a singleton axis:"
	152	+	+ str(im.shape)
	153	+	+ " Reduce dimensionality with np.squeeze(im) to avoid"
	154	+	+ " unexpected behavior."
	155	+	)
	156	+	if mode in ["backward", "reverse"]:
117	157		im = np.flip(im, axis)
118		-	im = distance_transform_lin(im=im, axis=axis, mode='forward')
	158	+	im = distance_transform_lin(im=im, axis=axis, mode="forward")
119	159		im = np.flip(im, axis)
120	160		return im
121		-	elif mode in ['both']:
122		-	im_f = distance_transform_lin(im=im, axis=axis, mode='forward')
123		-	im_b = distance_transform_lin(im=im, axis=axis, mode='backward')
	161	+	elif mode in ["both"]:
	162	+	im_f = distance_transform_lin(im=im, axis=axis, mode="forward")
	163	+	im_b = distance_transform_lin(im=im, axis=axis, mode="backward")
124	164		return np.minimum(im_f, im_b)
125	165		else:
126	166		b = np.cumsum(im > 0, axis=axis)
127		-	c = np.diff(b*(im == 0), axis=axis)
	167	+	c = np.diff(b * (im == 0), axis=axis)
128	168		d = np.minimum.accumulate(c, axis=axis)
129	169		if im.ndim == 1:
130		-	e = np.pad(d, pad_width=[1, 0], mode='constant', constant_values=0)
	170	+	e = np.pad(d, pad_width=[1, 0], mode="constant", constant_values=0)
131	171		elif im.ndim == 2:
132	172		ax = [[[1, 0], [0, 0]], [[0, 0], [1, 0]]]
133		-	e = np.pad(d, pad_width=ax[axis], mode='constant', constant_values=0)
	173	+	e = np.pad(d, pad_width=ax[axis], mode="constant", constant_values=0)
134	174		elif im.ndim == 3:
135		-	ax = [[[1, 0], [0, 0], [0, 0]],
136		-	[[0, 0], [1, 0], [0, 0]],
137		-	[[0, 0], [0, 0], [1, 0]]]
138		-	e = np.pad(d, pad_width=ax[axis], mode='constant', constant_values=0)
139		-	f = im*(b + e)
	175	+	ax = [
	176	+	[[1, 0], [0, 0], [0, 0]],
	177	+	[[0, 0], [1, 0], [0, 0]],
	178	+	[[0, 0], [0, 0], [1, 0]],
	179	+	]
	180	+	e = np.pad(d, pad_width=ax[axis], mode="constant", constant_values=0)
	181	+	f = im * (b + e)
140	182		return f
141	183
142	184

205	247		using marker-based watershed segmenation". Physical Review E. (2017)
206	248
207	249		"""
208		-	tup = namedtuple('results', field_names=['im', 'dt', 'peaks', 'regions'])
209		-	print('_'*60)
	250	+	tup = namedtuple("results", field_names=["im", "dt", "peaks", "regions"])
	251	+	print("-" * 60)
210	252		print("Beginning SNOW Algorithm")
211	253		im_shape = np.array(im.shape)
212	254		if im.dtype is not bool:
213		-	print('Converting supplied image (im) to boolean')
	255	+	print("Converting supplied image (im) to boolean")
214	256		im = im > 0
215	257		if dt is None:
216		-	print('Peforming Distance Transform')
	258	+	print("Peforming Distance Transform")
217	259		if np.any(im_shape == 1):
218	260		ax = np.where(im_shape == 1)[0][0]
219		-	dt = spim.distance_transform_edt(input=im.squeeze())
	261	+	dt = edt(im.squeeze())
220	262		dt = np.expand_dims(dt, ax)
221	263		else:
222		-	dt = spim.distance_transform_edt(input=im)
	264	+	dt = edt(im)
223	265
224	266		tup.im = im
225	267		tup.dt = dt
226	268
227	269		if sigma > 0:
228		-	print('Applying Gaussian blur with sigma =', str(sigma))
	270	+	print("Applying Gaussian blur with sigma =", str(sigma))
229	271		dt = spim.gaussian_filter(input=dt, sigma=sigma)
230	272
231	273		peaks = find_peaks(dt=dt, r_max=r_max)
232		-	print('Initial number of peaks: ', spim.label(peaks)[1])
	274	+	print("Initial number of peaks: ", spim.label(peaks)[1])
233	275		peaks = trim_saddle_points(peaks=peaks, dt=dt, max_iters=500)
234		-	print('Peaks after trimming saddle points: ', spim.label(peaks)[1])
	276	+	print("Peaks after trimming saddle points: ", spim.label(peaks)[1])
235	277		peaks = trim_nearby_peaks(peaks=peaks, dt=dt)
236	278		peaks, N = spim.label(peaks)
237		-	print('Peaks after trimming nearby peaks: ', N)
	279	+	print("Peaks after trimming nearby peaks: ", N)
238	280		tup.peaks = peaks
239	281		if mask:
240	282		mask_solid = im > 0

333	375		combined_dt = 0
334	376		combined_region = 0
335	377		num = [0]
336		-	for i in phases_num:
337		-	print('_' * 60)
	378	+	for i, j in enumerate(phases_num):
	379	+	print("_" * 60)
338	380		if alias is None:
339		-	print('Processing Phase {}'.format(i))
340		-	else:
341		-	print('Processing Phase {}'.format(al[i]))
342		-	phase_snow = snow_partitioning(im == i,
343		-	dt=None, r_max=r_max, sigma=sigma,
344		-	return_all=return_all, mask=mask,
345		-	randomize=randomize)
346		-	if len(phases_num) == 1 and phases_num == 1:
347		-	combined_dt = phase_snow.dt
348		-	combined_region = phase_snow.regions
	381	+	print("Processing Phase {}".format(j))
349	382		else:
350		-	combined_dt += phase_snow.dt
351		-	phase_snow.regions *= phase_snow.im
352		-	phase_snow.regions += num[i - 1]
353		-	phase_ws = phase_snow.regions * phase_snow.im
354		-	phase_ws[phase_ws == num[i - 1]] = 0
355		-	combined_region += phase_ws
	383	+	print("Processing Phase {}".format(al[j]))
	384	+	phase_snow = snow_partitioning(
	385	+	im == j,
	386	+	dt=None,
	387	+	r_max=r_max,
	388	+	sigma=sigma,
	389	+	return_all=return_all,
	390	+	mask=mask,
	391	+	randomize=randomize,
	392	+	)
	393	+	combined_dt += phase_snow.dt
	394	+	phase_snow.regions *= phase_snow.im
	395	+	phase_snow.regions += num[i]
	396	+	phase_ws = phase_snow.regions * phase_snow.im
	397	+	phase_ws[phase_ws == num[i]] = 0
	398	+	combined_region += phase_ws
356	399		num.append(np.amax(combined_region))
357	400		if return_all:
358		-	tup = namedtuple('results', field_names=['im', 'dt', 'phase_max_label',
359		-	'regions'])
	401	+	tup = namedtuple(
	402	+	"results", field_names=["im", "dt", "phase_max_label", "regions"]
	403	+	)
360	404		tup.im = im
361	405		tup.dt = combined_dt
362	406		tup.phase_max_label = num[1:]

399	441		``peaks = peak_local_max(image=dt, min_distance=r, exclude_border=0,
400	442		indices=False)``
401	443
402		-	This automatically uses a square structuring element which is significantly
403		-	faster than using a circular or spherical element.
	444	+	The skimage function automatically uses a square structuring element
	445	+	which is significantly faster than using a circular or spherical element.
404	446		"""
405	447		im = dt > 0
406	448		if im.ndim != im.squeeze().ndim:
407		-	warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
408		-	' Reduce dimensionality with np.squeeze(im) to avoid' +
409		-	' unexpected behavior.')
	449	+	warnings.warn("Input image conains a singleton axis:"
	450	+	+ str(im.shape)
	451	+	+ " Reduce dimensionality with np.squeeze(im) to avoid"
	452	+	+ " unexpected behavior.")
410	453		if footprint is None:
411	454		if im.ndim == 2:
412	455		footprint = disk

447	499		else:
448	500		strel = cube
449	501		markers, N = spim.label(input=peaks, structure=strel(3))
450		-	inds = spim.measurements.center_of_mass(input=peaks,
451		-	labels=markers,
452		-	index=np.arange(1, N+1))
	502	+	inds = spim.measurements.center_of_mass(
	503	+	input=peaks, labels=markers, index=np.arange(1, N + 1)
	504	+	)
453	505		inds = np.floor(inds).astype(int)
454	506		# Centroid may not be on old pixel, so create a new peaks image
455	507		peaks_new = np.zeros_like(peaks, dtype=bool)

457	509		return peaks_new
458	510
459	511
460		-	def trim_saddle_points(peaks, dt, max_iters=10):
	512	+	def trim_saddle_points(peaks, dt, max_iters=10, verbose=1):
461	513		r"""
462	514		Removes peaks that were mistakenly identified because they lied on a
463	515		saddle or ridge in the distance transform that was not actually a true

499	551		slices = spim.find_objects(labels)
500	552		for i in range(N):
501	553		s = extend_slice(s=slices[i], shape=peaks.shape, pad=10)
502		-	peaks_i = labels[s] == i+1
	554	+	peaks_i = labels[s] == i + 1
503	555		dt_i = dt[s]
504	556		im_i = dt_i > 0
505	557		iters = 0
506	558		peaks_dil = np.copy(peaks_i)
507	559		while iters < max_iters:
508	560		iters += 1
509		-	peaks_dil = spim.binary_dilation(input=peaks_dil,
510		-	structure=cube(3))
511		-	peaks_max = peaks_dilnp.amax(dt_ipeaks_dil)
512		-	peaks_extended = (peaks_max == dt_i)*im_i
	561	+	peaks_dil = spim.binary_dilation(input=peaks_dil, structure=cube(3))
	562	+	peaks_max = peaks_dil * np.amax(dt_i * peaks_dil)
	563	+	peaks_extended = (peaks_max == dt_i) * im_i
513	564		if np.all(peaks_extended == peaks_i):
514	565		break # Found a true peak
515		-	elif np.sum(peaks_extended*peaks_i) == 0:
	566	+	elif np.sum(peaks_extended * peaks_i) == 0:
516	567		peaks_i = False
517	568		break # Found a saddle point
518	569		peaks[s] = peaks_i
519		-	if iters >= max_iters:
520		-	print('Maximum number of iterations reached, consider'
521		-	+ 'running again with a larger value of max_iters')
	570	+	if iters >= max_iters and verbose:
	571	+	print(
	572	+	"Maximum number of iterations reached, consider "
	573	+	+ "running again with a larger value of max_iters"
	574	+	)
522	575		return peaks
523	576
524	577

561	614		from skimage.morphology import cube
562	615		peaks, N = spim.label(peaks, structure=cube(3))
563	616		crds = spim.measurements.center_of_mass(peaks, labels=peaks,
564		-	index=np.arange(1, N+1))
	617	+	index=np.arange(1, N + 1))
565	618		crds = np.vstack(crds).astype(int) # Convert to numpy array of ints
566	619		# Get distance between each peak as a distance map
567	620		tree = sptl.cKDTree(data=crds)

583	636		slices = spim.find_objects(input=peaks)
584	637		for s in drop_peaks:
585	638		peaks[slices[s]] = 0
586		-	return (peaks > 0)
	639	+	return peaks > 0
587	640
588	641
589	642		def find_disconnected_voxels(im, conn=None):

649	708		-------
650	709		image : ND-array
651	710		A version of ``im`` but with all the disconnected pores removed.
	711	+	conn : int
	712	+	For 2D the options are 4 and 8 for square and diagonal neighbors, while
	713	+	for the 3D the options are 6 and 26, similarily for square and diagonal
	714	+	neighbors. The default is the maximum option.
652	715
653	716		See Also
654	717		--------

656	719
657	720		"""
658	721		im = np.copy(im)
659		-	holes = find_disconnected_voxels(im)
	722	+	holes = find_disconnected_voxels(im, conn=conn)
660	723		im[holes] = False
661	724		return im
662	725
663	726
664		-	def trim_floating_solid(im):
	727	+	def trim_floating_solid(im, conn=None):
665	728		r"""
666	729		Removes all solid that that is not attached to the edges of the image.
667	730

669	732		----------
670	733		im : ND-array
671	734		The image of the porous material
	735	+	conn : int
	736	+	For 2D the options are 4 and 8 for square and diagonal neighbors, while
	737	+	for the 3D the options are 6 and 26, similarily for square and diagonal
	738	+	neighbors. The default is the maximum option.
672	739
673	740		Returns
674	741		-------

681	748
682	749		"""
683	750		im = np.copy(im)
684		-	holes = find_disconnected_voxels(~im)
	751	+	holes = find_disconnected_voxels(~im, conn=conn)
685	752		im[holes] = True
686	753		return im
687	754
688	755
689		-	def trim_nonpercolating_paths(im, inlet_axis=0, outlet_axis=0):
	756	+	def trim_nonpercolating_paths(im, inlet_axis=0, outlet_axis=0,
	757	+	inlets=None, outlets=None):
690	758		r"""
691	759		Removes all nonpercolating paths between specified edges
692	760

699	767		im : ND-array
700	768		The image of the porous material with ```True`` values indicating the
701	769		phase of interest
702		-
703	770		inlet_axis : int
704	771		Inlet axis of boundary condition. For three dimensional image the
705	772		number ranges from 0 to 2. For two dimensional image the range is
706		-	between 0 to 1.
707		-
	773	+	between 0 to 1. If ``inlets`` is given then this argument is ignored.
708	774		outlet_axis : int
709	775		Outlet axis of boundary condition. For three dimensional image the
710	776		number ranges from 0 to 2. For two dimensional image the range is
711		-	between 0 to 1.
	777	+	between 0 to 1. If ``outlets`` is given then this argument is ignored.
	778	+	inlets : ND-image (optional)
	779	+	A boolean mask indicating locations of inlets. If this argument is
	780	+	supplied then ``inlet_axis`` is ignored.
	781	+	outlets : ND-image (optional)
	782	+	A boolean mask indicating locations of outlets. If this argument is
	783	+	supplied then ``outlet_axis`` is ignored.
712	784
713	785		Returns
714	786		-------

790	866
791	867		"""
792	868		result = im
793		-	if mode in ['maxima', 'extrema']:
794		-	result = reconstruction(seed=im - h, mask=im, method='dilation')
795		-	elif mode in ['minima', 'extrema']:
796		-	result = reconstruction(seed=im + h, mask=im, method='erosion')
	869	+	if mode in ["maxima", "extrema"]:
	870	+	result = reconstruction(seed=im - h, mask=im, method="dilation")
	871	+	elif mode in ["minima", "extrema"]:
	872	+	result = reconstruction(seed=im + h, mask=im, method="erosion")
797	873		return result
798	874
799	875
800		-	@jit(forceobj=True)
801		-	def flood(im, regions=None, mode='max'):
	876	+	def flood(im, regions=None, mode="max"):
802	877		r"""
803	878		Floods/fills each region in an image with a single value based on the
804		-	specific values in that region. The ``mode`` argument is used to
805		-	determine how the value is calculated.
	879	+	specific values in that region.
	880	+
	881	+	The ``mode`` argument is used to determine how the value is calculated.
806	882
807	883		Parameters
808	884		----------
809	885		im : array_like
810		-	An ND image with isolated regions containing 0's elsewhere.
811		-
	886	+	An image with isolated regions with numerical values in each voxel,
	887	+	and 0's elsewhere.
812	888		regions : array_like
813		-	An array the same shape as ``im`` with each region labeled. If None is
814		-	supplied (default) then ``scipy.ndimage.label`` is used with its
815		-	default arguments.
816		-
	889	+	An array the same shape as ``im`` with each region labeled. If
	890	+	``None`` is supplied (default) then ``scipy.ndimage.label`` is used
	891	+	with its default arguments.
817	892		mode : string
818	893		Specifies how to determine which value should be used to flood each
819	894		region. Options are:
820	895
821		-	'max' - Floods each region with the local maximum in that region
	896	+	'maximum' - Floods each region with the local maximum in that region
	897	+
	898	+	'minimum' - Floods each region the local minimum in that region
	899	+
	900	+	'median' - Floods each region the local median in that region
822	901
823		-	'min' - Floods each region the local minimum in that region
	902	+	'mean' - Floods each region the local mean in that region
824	903
825	904		'size' - Floods each region with the size of that region
826	905

841	920		else:
842	921		labels = np.copy(regions)
843	922		N = labels.max()
844		-	I = im.flatten()
845		-	L = labels.flatten()
846		-	if mode.startswith('max'):
847		-	V = np.zeros(shape=N+1, dtype=float)
848		-	for i in range(len(L)):
849		-	if V[L[i]] < I[i]:
850		-	V[L[i]] = I[i]
851		-	elif mode.startswith('min'):
852		-	V = np.ones(shape=N+1, dtype=float)*np.inf
853		-	for i in range(len(L)):
854		-	if V[L[i]] > I[i]:
855		-	V[L[i]] = I[i]
856		-	elif mode.startswith('size'):
857		-	V = np.zeros(shape=N+1, dtype=int)
858		-	for i in range(len(L)):
859		-	V[L[i]] += 1
860		-	im_flooded = np.reshape(V[labels], newshape=im.shape)
861		-	im_flooded = im_flooded*mask
	923	+	mode = "sum" if mode == "size" else mode
	924	+	mode = "maximum" if mode == "max" else mode
	925	+	mode = "minimum" if mode == "min" else mode
	926	+	if mode in ["mean", "median", "maximum", "minimum", "sum"]:
	927	+	f = getattr(spim, mode)
	928	+	vals = f(input=im, labels=labels, index=range(0, N + 1))
	929	+	im_flooded = vals[labels]
	930	+	im_flooded = im_flooded * mask
	931	+	else:
	932	+	raise Exception(mode + " is not a recognized mode")
862	933		return im_flooded
863	934
864	935

886	957		the image. Obviously, voxels with a value of zero have no error.
887	958
888	959		"""
889		-	temp = np.ones(shape=dt.shape)*np.inf
	960	+	temp = np.ones(shape=dt.shape) * np.inf
890	961		for ax in range(dt.ndim):
891	962		dt_lin = distance_transform_lin(np.ones_like(temp, dtype=bool),
892		-	axis=ax, mode='both')
	963	+	axis=ax, mode="both")
893	964		temp = np.minimum(temp, dt_lin)
894	965		result = np.clip(dt - temp, a_min=0, a_max=np.inf)
895	966		return result

913	984		A copy of ``im`` with each voxel value indicating the size of the
914	985		region to which it belongs. This is particularly useful for finding
915	986		chord sizes on the image produced by ``apply_chords``.
	987	+
	988	+	See Also
	989	+	--------
	990	+	flood
916	991		"""
917	992		if im.dtype == bool:
918	993		im = spim.label(im)[0]

962	1037
963	1038		"""
964	1039		if im.ndim != im.squeeze().ndim:
965		-	warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
966		-	' Reduce dimensionality with np.squeeze(im) to avoid' +
967		-	' unexpected behavior.')
	1040	+	warnings.warn(
	1041	+	"Input image conains a singleton axis:"
	1042	+	+ str(im.shape)
	1043	+	+ " Reduce dimensionality with np.squeeze(im) to avoid"
	1044	+	+ " unexpected behavior."
	1045	+	)
968	1046		if spacing < 0:
969		-	raise Exception('Spacing cannot be less than 0')
	1047	+	raise Exception("Spacing cannot be less than 0")
970	1048		if spacing == 0:
971	1049		label = True
972	1050		result = np.zeros(im.shape, dtype=int) # Will receive chords at end
973		-	slxyz = [slice(None, None, spacing*(axis != i) + 1) for i in [0, 1, 2]]
974		-	slices = tuple(slxyz[:im.ndim])
	1051	+	slxyz = [slice(None, None, spacing * (axis != i) + 1) for i in [0, 1, 2]]
	1052	+	slices = tuple(slxyz[: im.ndim])
975	1053		s = [[0, 1, 0], [0, 1, 0], [0, 1, 0]] # Straight-line structuring element
976	1054		if im.ndim == 3: # Make structuring element 3D if necessary
977	1055		s = np.pad(np.atleast_3d(s), pad_width=((0, 0), (0, 0), (1, 1)),
978		-	mode='constant', constant_values=0)
	1056	+	mode="constant", constant_values=0)
979	1057		im = im[slices]
980	1058		s = np.swapaxes(s, 0, axis)
981	1059		chords = spim.label(im, structure=s)[0]

1025	1103
1026	1104		"""
1027	1105		if im.ndim != im.squeeze().ndim:
1028		-	warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
1029		-	' Reduce dimensionality with np.squeeze(im) to avoid' +
1030		-	' unexpected behavior.')
	1106	+	warnings.warn(
	1107	+	"Input image conains a singleton axis:"
	1108	+	+ str(im.shape)
	1109	+	+ " Reduce dimensionality with np.squeeze(im) to avoid"
	1110	+	+ " unexpected behavior."
	1111	+	)
1031	1112		if im.ndim < 3:
1032		-	raise Exception('Must be a 3D image to use this function')
	1113	+	raise Exception("Must be a 3D image to use this function")
1033	1114		if spacing < 0:
1034		-	raise Exception('Spacing cannot be less than 0')
	1115	+	raise Exception("Spacing cannot be less than 0")
1035	1116		ch = np.zeros_like(im, dtype=int)
1036		-	ch[:, ::4+2spacing, ::4+2spacing] = 1 # X-direction
1037		-	ch[::4+2spacing, :, 2::4+2spacing] = 2 # Y-direction
1038		-	ch[2::4+2spacing, 2::4+2spacing, :] = 3 # Z-direction
1039		-	chords = ch*im
	1117	+	ch[:, :: 4 + 2 * spacing, :: 4 + 2 * spacing] = 1 # X-direction
	1118	+	ch[:: 4 + 2 * spacing, :, 2::4 + 2 * spacing] = 2 # Y-direction
	1119	+	ch[2::4 + 2 * spacing, 2::4 + 2 * spacing, :] = 3 # Z-direction
	1120	+	chords = ch * im
1040	1121		if trim_edges:
1041	1122		temp = clear_border(spim.label(chords > 0)[0]) > 0
1042		-	chords = temp*chords
	1123	+	chords = temp * chords
1043	1124		return chords
1044	1125
1045	1126
1046		-	def local_thickness(im, sizes=25, mode='hybrid'):
	1127	+	def local_thickness(im, sizes=25, mode="hybrid", **kwargs):
1047	1128		r"""
1048	1129		For each voxel, this function calculates the radius of the largest sphere
1049	1130		that both engulfs the voxel and fits entirely within the foreground.

1104	1185		function. This is not needed in ``local_thickness`` however.
1105	1186
1106	1187		"""
1107		-	im_new = porosimetry(im=im, sizes=sizes, access_limited=False, mode=mode)
	1188	+	im_new = porosimetry(im=im, sizes=sizes, access_limited=False, mode=mode,
	1189	+	**kwargs)
1108	1190		return im_new
1109	1191
1110	1192
1111		-	def porosimetry(im, sizes=25, inlets=None, access_limited=True,
1112		-	mode='hybrid'):
	1193	+	def porosimetry(im, sizes=25, inlets=None, access_limited=True, mode='hybrid',
	1194	+	fft=True, **kwargs):
1113	1195		r"""
1114		-	Performs a porosimetry simulution on the image
	1196	+	Performs a porosimetry simulution on an image
1115	1197
1116	1198		Parameters
1117	1199		----------
1118	1200		im : ND-array
1119		-	An ND image of the porous material containing True values in the
	1201	+	An ND image of the porous material containing ``True`` values in the
1120	1202		pore space.
1121	1203
1122	1204		sizes : array_like or scalar

1125	1207		the min and max of the distance transform are used.
1126	1208
1127	1209		inlets : ND-array, boolean
1128		-	A boolean mask with True values indicating where the invasion
	1210	+	A boolean mask with ``True`` values indicating where the invasion
1129	1211		enters the image. By default all faces are considered inlets,
1130	1212		akin to a mercury porosimetry experiment. Users can also apply
1131	1213		solid boundaries to their image externally before passing it in,

1134	1216
1135	1217		access_limited : Boolean
1136	1218		This flag indicates if the intrusion should only occur from the
1137		-	surfaces (``access_limited`` is True, which is the default), or
	1219	+	surfaces (``access_limited`` is ``True``, which is the default), or
1138	1220		if the invading phase should be allowed to appear in the core of
1139	1221		the image. The former simulates experimental tools like mercury
1140	1222		intrusion porosimetry, while the latter is useful for comparison

1159	1241		included mostly for comparison purposes. The morphological operations
1160	1242		are done using fft-based method implementations.
1161	1243
	1244	+	fft : boolean (default is ``True``)
	1245	+	Indicates whether to use the ``fftmorphology`` function in
	1246	+	``porespy.filters`` or to use the standard morphology functions in
	1247	+	``scipy.ndimage``. Always use ``fft=True`` unless you have a good
	1248	+	reason not to.
	1249	+
1162	1250		Returns
1163	1251		-------
1164	1252		image : ND-array

1168	1256		capillary pressure by applying a boolean comparison:
1169	1257		``inv_phase = im > r`` where ``r`` is the radius (in voxels) of the
1170	1258		invading sphere. Of course, ``r`` can be converted to capillary
1171		-	pressure using your favorite model.
	1259	+	pressure using a preferred model.
1172	1260
1173	1261		Notes
1174	1262		-----

1179	1267
1180	1268		See Also
1181	1269		--------
1182		-	fftmorphology
1183	1270		local_thickness
1184	1271
1185	1272		"""
1186	1273		if im.ndim != im.squeeze().ndim:
1187		-	warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
1188		-	' Reduce dimensionality with np.squeeze(im) to avoid' +
1189		-	' unexpected behavior.')
	1274	+	warnings.warn("Input image contains a singleton axis:"
	1275	+	+ str(im.shape)
	1276	+	+ " Reduce dimensionality with np.squeeze(im) to avoid"
	1277	+	+ " unexpected behavior.")
1190	1278
1191		-	dt = spim.distance_transform_edt(im > 0)
	1279	+	dt = edt(im > 0)
1192	1280
1193	1281		if inlets is None:
1194		-	inlets = get_border(im.shape, mode='faces')
	1282	+	inlets = get_border(im.shape, mode="faces")
1195	1283
1196	1284		if isinstance(sizes, int):
1197	1285		sizes = np.logspace(start=np.log10(np.amax(dt)), stop=0, num=sizes)

1203	1291		else:
1204	1292		strel = ps_ball
1205	1293
1206		-	if mode == 'mio':
	1294	+	# Parse kwargs for any parallelization arguments
	1295	+	parallel = kwargs.pop('parallel', False)
	1296	+	cores = kwargs.pop('cores', None)
	1297	+	divs = kwargs.pop('divs', 2)
	1298	+
	1299	+	if mode == "mio":
1207	1300		pw = int(np.floor(dt.max()))
1208		-	impad = np.pad(im, mode='symmetric', pad_width=pw)
1209		-	inletspad = np.pad(inlets, mode='symmetric', pad_width=pw)
1210		-	inlets = np.where(inletspad)
1211		-	# sizes = np.unique(np.around(sizes, decimals=0).astype(int))[-1::-1]
	1301	+	impad = np.pad(im, mode="symmetric", pad_width=pw)
	1302	+	inlets = np.pad(inlets, mode="symmetric", pad_width=pw)
	1303	+	# sizes = np.unique(np.around(sizes, decimals=0).astype(int))[-1::-1]
1212	1304		imresults = np.zeros(np.shape(impad))
1213		-	for r in tqdm(sizes):
1214		-	imtemp = fftmorphology(impad, strel(r), mode='erosion')
1215		-	if access_limited:
1216		-	imtemp = trim_disconnected_blobs(imtemp, inlets)
1217		-	imtemp = fftmorphology(imtemp, strel(r), mode='dilation')
1218		-	if np.any(imtemp):
1219		-	imresults[(imresults == 0)*imtemp] = r
	1305	+	with tqdm(sizes) as pbar:
	1306	+	for r in sizes:
	1307	+	pbar.update()
	1308	+	if parallel:
	1309	+	imtemp = chunked_func(func=spim.binary_erosion,
	1310	+	input=impad, structure=strel(r),
	1311	+	overlap=int(2*r) + 1,
	1312	+	cores=cores, divs=divs)
	1313	+	elif fft:
	1314	+	imtemp = fftmorphology(impad, strel(r), mode="erosion")
	1315	+	else:
	1316	+	imtemp = spim.binary_erosion(input=impad,
	1317	+	structure=strel(r))
	1318	+	if access_limited:
	1319	+	imtemp = trim_disconnected_blobs(imtemp, inlets)
	1320	+	if parallel:
	1321	+	imtemp = chunked_func(func=spim.binary_dilation,
	1322	+	input=imtemp, structure=strel(r),
	1323	+	overlap=int(2*r) + 1,
	1324	+	cores=cores, divs=divs)
	1325	+	elif fft:
	1326	+	imtemp = fftmorphology(imtemp, strel(r), mode="dilation")
	1327	+	else:
	1328	+	imtemp = spim.binary_dilation(input=imtemp,
	1329	+	structure=strel(r))
	1330	+	if np.any(imtemp):
	1331	+	imresults[(imresults == 0) * imtemp] = r
1220	1332		imresults = extract_subsection(imresults, shape=im.shape)
1221		-	elif mode == 'dt':
1222		-	inlets = np.where(inlets)
	1333	+	elif mode == "dt":
1223	1334		imresults = np.zeros(np.shape(im))
1224		-	for r in tqdm(sizes):
1225		-	imtemp = dt >= r
1226		-	if access_limited:
1227		-	imtemp = trim_disconnected_blobs(imtemp, inlets)
1228		-	if np.any(imtemp):
1229		-	imtemp = spim.distance_transform_edt(~imtemp) < r
1230		-	imresults[(imresults == 0)*imtemp] = r
1231		-	elif mode == 'hybrid':
1232		-	inlets = np.where(inlets)
	1335	+	with tqdm(sizes) as pbar:
	1336	+	for r in sizes:
	1337	+	pbar.update()
	1338	+	imtemp = dt >= r
	1339	+	if access_limited:
	1340	+	imtemp = trim_disconnected_blobs(imtemp, inlets)
	1341	+	if np.any(imtemp):
	1342	+	imtemp = edt(~imtemp) < r
	1343	+	imresults[(imresults == 0) * imtemp] = r
	1344	+	elif mode == "hybrid":
1233	1345		imresults = np.zeros(np.shape(im))
1234		-	for r in tqdm(sizes):
1235		-	imtemp = dt >= r
1236		-	if access_limited:
1237		-	imtemp = trim_disconnected_blobs(imtemp, inlets)
1238		-	if np.any(imtemp):
1239		-	imtemp = fftconvolve(imtemp, strel(r), mode='same') > 0.0001
1240		-	imresults[(imresults == 0)*imtemp] = r
	1346	+	with tqdm(sizes) as pbar:
	1347	+	for r in sizes:
	1348	+	pbar.update()
	1349	+	imtemp = dt >= r
	1350	+	if access_limited:
	1351	+	imtemp = trim_disconnected_blobs(imtemp, inlets)
	1352	+	if np.any(imtemp):
	1353	+	if parallel:
	1354	+	imtemp = chunked_func(func=spim.binary_dilation,
	1355	+	input=imtemp, structure=strel(r),
	1356	+	overlap=int(2*r) + 1,
	1357	+	cores=cores, divs=divs)
	1358	+	elif fft:
	1359	+	imtemp = fftmorphology(imtemp, strel(r),
	1360	+	mode="dilation")
	1361	+	else:
	1362	+	imtemp = spim.binary_dilation(input=imtemp,
	1363	+	structure=strel(r))
	1364	+	imresults[(imresults == 0) * imtemp] = r
1241	1365		else:
1242		-	raise Exception('Unreckognized mode ' + mode)
	1366	+	raise Exception("Unrecognized mode " + mode)
1243	1367		return imresults
1244	1368
1245	1369
1246		-	def trim_disconnected_blobs(im, inlets):
	1370	+	def trim_disconnected_blobs(im, inlets, strel=None):
1247	1371		r"""
1248	1372		Removes foreground voxels not connected to specified inlets
1249	1373
1250	1374		Parameters
1251	1375		----------
1252	1376		im : ND-array
1253		-	The array to be trimmed
	1377	+	The image containing the blobs to be trimmed
1254	1378		inlets : ND-array or tuple of indices
1255	1379		The locations of the inlets. Can either be a boolean mask the same
1256	1380		shape as ``im``, or a tuple of indices such as that returned by the
1257	1381		``where`` function. Any voxels not connected directly to
1258	1382		the inlets will be trimmed.
	1383	+	strel : ND-array
	1384	+	The neighborhood over which connectivity should be checked. It must
	1385	+	be symmetric and the same dimensionality as the image. It is passed
	1386	+	directly to the ``scipy.ndimage.label`` function as the ``structure``
	1387	+	argument so refer to that docstring for additional info.
1259	1388
1260	1389		Returns
1261	1390		-------