PMEAL / porespy

@@ -1,22 +1,21 @@
Loading
1 -
from collections import namedtuple
1 +
import sys
2 +
import dask
3 +
import warnings
2 4
import numpy as np
5 +
from edt import edt
3 6
import operator as op
7 +
from tqdm import tqdm
4 8
import scipy.ndimage as spim
5 9
import scipy.spatial as sptl
6 -
import warnings
7 -
import dask
8 -
from edt import edt
9 -
from dask.diagnostics import ProgressBar
10 +
from collections import namedtuple
10 11
from scipy.signal import fftconvolve
11 -
from tqdm import tqdm
12 -
from numba import jit
13 12
from skimage.segmentation import clear_border
14 13
from skimage.morphology import ball, disk, square, cube, diamond, octahedron
15 14
from skimage.morphology import reconstruction, watershed
16 -
from porespy.tools import randomize_colors, fftmorphology, make_contiguous
15 +
from porespy.tools import randomize_colors, fftmorphology
17 16
from porespy.tools import get_border, extend_slice, extract_subsection
18 -
from porespy.tools import ps_disk, ps_ball
19 17
from porespy.tools import _create_alias_map
18 +
from porespy.tools import ps_disk, ps_ball
20 19
21 20
22 21
def trim_small_clusters(im, size=1):
@@ -44,11 +43,11 @@
Loading
44 43
    elif im.ndims == 3:
45 44
        strel = ball(1)
46 45
    else:
47 -
        raise Exception('Only 2D or 3D images are accepted')
46 +
        raise Exception("Only 2D or 3D images are accepted")
48 47
    filtered_array = np.copy(im)
49 48
    labels, N = spim.label(filtered_array, structure=strel)
50 49
    id_sizes = np.array(spim.sum(im, labels, range(N + 1)))
51 -
    area_mask = (id_sizes <= size)
50 +
    area_mask = id_sizes <= size
52 51
    filtered_array[area_mask[labels]] = 0
53 52
    return filtered_array
54 53
@@ -64,12 +63,12 @@
Loading
64 63
    """
65 64
    A = im
66 65
    B = np.swapaxes(A, axis, -1)
67 -
    updown = np.empty((*B.shape[:-1], B.shape[-1]+1), B.dtype)
66 +
    updown = np.empty((*B.shape[:-1], B.shape[-1] + 1), B.dtype)
68 67
    updown[..., 0], updown[..., -1] = -1, -1
69 68
    np.subtract(B[..., 1:], B[..., :-1], out=updown[..., 1:-1])
70 69
    chnidx = np.where(updown)
71 70
    chng = updown[chnidx]
72 -
    pkidx, = np.where((chng[:-1] > 0) & (chng[1:] < 0) | (chnidx[-1][:-1] == 0))
71 +
    (pkidx,) = np.where((chng[:-1] > 0) & (chng[1:] < 0) | (chnidx[-1][:-1] == 0))
73 72
    pkidx = (*map(op.itemgetter(pkidx), chnidx),)
74 73
    out = np.zeros_like(A)
75 74
    aux = out.swapaxes(axis, -1)
@@ -79,7 +78,7 @@
Loading
79 78
    return result
80 79
81 80
82 -
def distance_transform_lin(im, axis=0, mode='both'):
81 +
def distance_transform_lin(im, axis=0, mode="both"):
83 82
    r"""
84 83
    Replaces each void voxel with the linear distance to the nearest solid
85 84
    voxel along the specified axis.
@@ -113,38 +112,44 @@
Loading
113 112
        the nearest background along the specified axis.
114 113
    """
115 114
    if im.ndim != im.squeeze().ndim:
116 -
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
117 -
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
118 -
                      ' unexpected behavior.')
119 -
    if mode in ['backward', 'reverse']:
115 +
        warnings.warn(
116 +
            "Input image conains a singleton axis:"
117 +
            + str(im.shape)
118 +
            + " Reduce dimensionality with np.squeeze(im) to avoid"
119 +
            + " unexpected behavior."
120 +
        )
121 +
    if mode in ["backward", "reverse"]:
120 122
        im = np.flip(im, axis)
121 -
        im = distance_transform_lin(im=im, axis=axis, mode='forward')
123 +
        im = distance_transform_lin(im=im, axis=axis, mode="forward")
122 124
        im = np.flip(im, axis)
123 125
        return im
124 -
    elif mode in ['both']:
125 -
        im_f = distance_transform_lin(im=im, axis=axis, mode='forward')
126 -
        im_b = distance_transform_lin(im=im, axis=axis, mode='backward')
126 +
    elif mode in ["both"]:
127 +
        im_f = distance_transform_lin(im=im, axis=axis, mode="forward")
128 +
        im_b = distance_transform_lin(im=im, axis=axis, mode="backward")
127 129
        return np.minimum(im_f, im_b)
128 130
    else:
129 131
        b = np.cumsum(im > 0, axis=axis)
130 -
        c = np.diff(b*(im == 0), axis=axis)
132 +
        c = np.diff(b * (im == 0), axis=axis)
131 133
        d = np.minimum.accumulate(c, axis=axis)
132 134
        if im.ndim == 1:
133 -
            e = np.pad(d, pad_width=[1, 0], mode='constant', constant_values=0)
135 +
            e = np.pad(d, pad_width=[1, 0], mode="constant", constant_values=0)
134 136
        elif im.ndim == 2:
135 137
            ax = [[[1, 0], [0, 0]], [[0, 0], [1, 0]]]
136 -
            e = np.pad(d, pad_width=ax[axis], mode='constant', constant_values=0)
138 +
            e = np.pad(d, pad_width=ax[axis], mode="constant", constant_values=0)
137 139
        elif im.ndim == 3:
138 -
            ax = [[[1, 0], [0, 0], [0, 0]],
139 -
                  [[0, 0], [1, 0], [0, 0]],
140 -
                  [[0, 0], [0, 0], [1, 0]]]
141 -
            e = np.pad(d, pad_width=ax[axis], mode='constant', constant_values=0)
142 -
        f = im*(b + e)
140 +
            ax = [
141 +
                [[1, 0], [0, 0], [0, 0]],
142 +
                [[0, 0], [1, 0], [0, 0]],
143 +
                [[0, 0], [0, 0], [1, 0]],
144 +
            ]
145 +
            e = np.pad(d, pad_width=ax[axis], mode="constant", constant_values=0)
146 +
        f = im * (b + e)
143 147
        return f
144 148
145 149
146 -
def snow_partitioning(im, dt=None, r_max=4, sigma=0.4, return_all=False,
147 -
                      mask=True, randomize=True):
150 +
def snow_partitioning(
151 +
    im, dt=None, r_max=4, sigma=0.4, return_all=False, mask=True, randomize=True
152 +
):
148 153
    r"""
149 154
    Partitions the void space into pore regions using a marker-based watershed
150 155
    algorithm, with specially filtered peaks as markers.
@@ -208,15 +213,15 @@
Loading
208 213
    using marker-based watershed segmenation".  Physical Review E. (2017)
209 214
210 215
    """
211 -
    tup = namedtuple('results', field_names=['im', 'dt', 'peaks', 'regions'])
212 -
    print('_'*60)
216 +
    tup = namedtuple("results", field_names=["im", "dt", "peaks", "regions"])
217 +
    print("-" * 60)
213 218
    print("Beginning SNOW Algorithm")
214 219
    im_shape = np.array(im.shape)
215 220
    if im.dtype is not bool:
216 -
        print('Converting supplied image (im) to boolean')
221 +
        print("Converting supplied image (im) to boolean")
217 222
        im = im > 0
218 223
    if dt is None:
219 -
        print('Peforming Distance Transform')
224 +
        print("Peforming Distance Transform")
220 225
        if np.any(im_shape == 1):
221 226
            ax = np.where(im_shape == 1)[0][0]
222 227
            dt = edt(im.squeeze())
@@ -228,16 +233,16 @@
Loading
228 233
    tup.dt = dt
229 234
230 235
    if sigma > 0:
231 -
        print('Applying Gaussian blur with sigma =', str(sigma))
236 +
        print("Applying Gaussian blur with sigma =", str(sigma))
232 237
        dt = spim.gaussian_filter(input=dt, sigma=sigma)
233 238
234 239
    peaks = find_peaks(dt=dt, r_max=r_max)
235 -
    print('Initial number of peaks: ', spim.label(peaks)[1])
240 +
    print("Initial number of peaks: ", spim.label(peaks)[1])
236 241
    peaks = trim_saddle_points(peaks=peaks, dt=dt, max_iters=500)
237 -
    print('Peaks after trimming saddle points: ', spim.label(peaks)[1])
242 +
    print("Peaks after trimming saddle points: ", spim.label(peaks)[1])
238 243
    peaks = trim_nearby_peaks(peaks=peaks, dt=dt)
239 244
    peaks, N = spim.label(peaks)
240 -
    print('Peaks after trimming nearby peaks: ', N)
245 +
    print("Peaks after trimming nearby peaks: ", N)
241 246
    tup.peaks = peaks
242 247
    if mask:
243 248
        mask_solid = im > 0
@@ -253,8 +258,9 @@
Loading
253 258
        return regions
254 259
255 260
256 -
def snow_partitioning_n(im, r_max=4, sigma=0.4, return_all=True,
257 -
                        mask=True, randomize=False, alias=None):
261 +
def snow_partitioning_n(
262 +
    im, r_max=4, sigma=0.4, return_all=True, mask=True, randomize=False, alias=None
263 +
):
258 264
    r"""
259 265
    This function partitions an imaging oontain an arbitrary number of phases
260 266
    into regions using a marker-based watershed segmentation. Its an extension
@@ -337,15 +343,20 @@
Loading
337 343
    combined_region = 0
338 344
    num = [0]
339 345
    for i in phases_num:
340 -
        print('_' * 60)
346 +
        print("_" * 60)
341 347
        if alias is None:
342 -
            print('Processing Phase {}'.format(i))
348 +
            print("Processing Phase {}".format(i))
343 349
        else:
344 -
            print('Processing Phase {}'.format(al[i]))
345 -
        phase_snow = snow_partitioning(im == i,
346 -
                                       dt=None, r_max=r_max, sigma=sigma,
347 -
                                       return_all=return_all, mask=mask,
348 -
                                       randomize=randomize)
350 +
            print("Processing Phase {}".format(al[i]))
351 +
        phase_snow = snow_partitioning(
352 +
            im == i,
353 +
            dt=None,
354 +
            r_max=r_max,
355 +
            sigma=sigma,
356 +
            return_all=return_all,
357 +
            mask=mask,
358 +
            randomize=randomize,
359 +
        )
349 360
        if len(phases_num) == 1 and phases_num == 1:
350 361
            combined_dt = phase_snow.dt
351 362
            combined_region = phase_snow.regions
@@ -358,8 +369,9 @@
Loading
358 369
            combined_region += phase_ws
359 370
        num.append(np.amax(combined_region))
360 371
    if return_all:
361 -
        tup = namedtuple('results', field_names=['im', 'dt', 'phase_max_label',
362 -
                                                 'regions'])
372 +
        tup = namedtuple(
373 +
            "results", field_names=["im", "dt", "phase_max_label", "regions"]
374 +
        )
363 375
        tup.im = im
364 376
        tup.dt = combined_dt
365 377
        tup.phase_max_label = num[1:]
@@ -407,9 +419,12 @@
Loading
407 419
    """
408 420
    im = dt > 0
409 421
    if im.ndim != im.squeeze().ndim:
410 -
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
411 -
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
412 -
                      ' unexpected behavior.')
422 +
        warnings.warn(
423 +
            "Input image conains a singleton axis:"
424 +
            + str(im.shape)
425 +
            + " Reduce dimensionality with np.squeeze(im) to avoid"
426 +
            + " unexpected behavior."
427 +
        )
413 428
    if footprint is None:
414 429
        if im.ndim == 2:
415 430
            footprint = disk
@@ -417,8 +432,8 @@
Loading
417 432
            footprint = ball
418 433
        else:
419 434
            raise Exception("only 2-d and 3-d images are supported")
420 -
    mx = spim.maximum_filter(dt + 2*(~im), footprint=footprint(r_max))
421 -
    peaks = (dt == mx)*im
435 +
    mx = spim.maximum_filter(dt + 2 * (~im), footprint=footprint(r_max))
436 +
    peaks = (dt == mx) * im
422 437
    return peaks
423 438
424 439
@@ -450,9 +465,9 @@
Loading
450 465
    else:
451 466
        strel = cube
452 467
    markers, N = spim.label(input=peaks, structure=strel(3))
453 -
    inds = spim.measurements.center_of_mass(input=peaks,
454 -
                                            labels=markers,
455 -
                                            index=np.arange(1, N+1))
468 +
    inds = spim.measurements.center_of_mass(
469 +
        input=peaks, labels=markers, index=np.arange(1, N + 1)
470 +
    )
456 471
    inds = np.floor(inds).astype(int)
457 472
    # Centroid may not be on old pixel, so create a new peaks image
458 473
    peaks_new = np.zeros_like(peaks, dtype=bool)
@@ -502,26 +517,27 @@
Loading
502 517
    slices = spim.find_objects(labels)
503 518
    for i in range(N):
504 519
        s = extend_slice(s=slices[i], shape=peaks.shape, pad=10)
505 -
        peaks_i = labels[s] == i+1
520 +
        peaks_i = labels[s] == i + 1
506 521
        dt_i = dt[s]
507 522
        im_i = dt_i > 0
508 523
        iters = 0
509 524
        peaks_dil = np.copy(peaks_i)
510 525
        while iters < max_iters:
511 526
            iters += 1
512 -
            peaks_dil = spim.binary_dilation(input=peaks_dil,
513 -
                                             structure=cube(3))
514 -
            peaks_max = peaks_dil*np.amax(dt_i*peaks_dil)
515 -
            peaks_extended = (peaks_max == dt_i)*im_i
527 +
            peaks_dil = spim.binary_dilation(input=peaks_dil, structure=cube(3))
528 +
            peaks_max = peaks_dil * np.amax(dt_i * peaks_dil)
529 +
            peaks_extended = (peaks_max == dt_i) * im_i
516 530
            if np.all(peaks_extended == peaks_i):
517 531
                break  # Found a true peak
518 -
            elif np.sum(peaks_extended*peaks_i) == 0:
532 +
            elif np.sum(peaks_extended * peaks_i) == 0:
519 533
                peaks_i = False
520 534
                break  # Found a saddle point
521 535
        peaks[s] = peaks_i
522 536
        if iters >= max_iters:
523 -
            print('Maximum number of iterations reached, consider'
524 -
                  + 'running again with a larger value of max_iters')
537 +
            print(
538 +
                "Maximum number of iterations reached, consider "
539 +
                + "running again with a larger value of max_iters"
540 +
            )
525 541
    return peaks
526 542
527 543
@@ -563,8 +579,9 @@
Loading
563 579
    else:
564 580
        from skimage.morphology import cube
565 581
    peaks, N = spim.label(peaks, structure=cube(3))
566 -
    crds = spim.measurements.center_of_mass(peaks, labels=peaks,
567 -
                                            index=np.arange(1, N+1))
582 +
    crds = spim.measurements.center_of_mass(
583 +
        peaks, labels=peaks, index=np.arange(1, N + 1)
584 +
    )
568 585
    crds = np.vstack(crds).astype(int)  # Convert to numpy array of ints
569 586
    # Get distance between each peak as a distance map
570 587
    tree = sptl.cKDTree(data=crds)
@@ -586,7 +603,7 @@
Loading
586 603
    slices = spim.find_objects(input=peaks)
587 604
    for s in drop_peaks:
588 605
        peaks[slices[s]] = 0
589 -
    return (peaks > 0)
606 +
    return peaks > 0
590 607
591 608
592 609
def find_disconnected_voxels(im, conn=None):
@@ -620,23 +637,26 @@
Loading
620 637
621 638
    """
622 639
    if im.ndim != im.squeeze().ndim:
623 -
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
624 -
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
625 -
                      ' unexpected behavior.')
640 +
        warnings.warn(
641 +
            "Input image conains a singleton axis:"
642 +
            + str(im.shape)
643 +
            + " Reduce dimensionality with np.squeeze(im) to avoid"
644 +
            + " unexpected behavior."
645 +
        )
626 646
    if im.ndim == 2:
627 647
        if conn == 4:
628 648
            strel = disk(1)
629 649
        elif conn in [None, 8]:
630 650
            strel = square(3)
631 651
        else:
632 -
            raise Exception('Received conn is not valid')
652 +
            raise Exception("Received conn is not valid")
633 653
    elif im.ndim == 3:
634 654
        if conn == 6:
635 655
            strel = ball(1)
636 656
        elif conn in [None, 26]:
637 657
            strel = cube(3)
638 658
        else:
639 -
            raise Exception('Received conn is not valid')
659 +
            raise Exception("Received conn is not valid")
640 660
    labels, N = spim.label(input=im, structure=strel)
641 661
    holes = clear_border(labels=labels) > 0
642 662
    return holes
@@ -700,8 +720,9 @@
Loading
700 720
    return im
701 721
702 722
703 -
def trim_nonpercolating_paths(im, inlet_axis=0, outlet_axis=0,
704 -
                              inlets=None, outlets=None):
723 +
def trim_nonpercolating_paths(
724 +
    im, inlet_axis=0, outlet_axis=0, inlets=None, outlets=None
725 +
):
705 726
    r"""
706 727
    Removes all nonpercolating paths between specified edges
707 728
@@ -742,9 +763,12 @@
Loading
742 763
743 764
    """
744 765
    if im.ndim != im.squeeze().ndim:
745 -
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
746 -
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
747 -
                      ' unexpected behavior.')
766 +
        warnings.warn(
767 +
            "Input image conains a singleton axis:"
768 +
            + str(im.shape)
769 +
            + " Reduce dimensionality with np.squeeze(im) to avoid"
770 +
            + " unexpected behavior."
771 +
        )
748 772
    im = trim_floating_solid(~im)
749 773
    labels = spim.label(~im)[0]
750 774
    if inlets is None:
@@ -775,14 +799,14 @@
Loading
775 799
                outlets[-1, :] = True
776 800
            elif outlet_axis == 1:
777 801
                outlets[:, -1] = True
778 -
    IN = np.unique(labels*inlets)
779 -
    OUT = np.unique(labels*outlets)
802 +
    IN = np.unique(labels * inlets)
803 +
    OUT = np.unique(labels * outlets)
780 804
    new_im = np.isin(labels, list(set(IN) ^ set(OUT)), invert=True)
781 805
    im[new_im == 0] = True
782 806
    return ~im
783 807
784 808
785 -
def trim_extrema(im, h, mode='maxima'):
809 +
def trim_extrema(im, h, mode="maxima"):
786 810
    r"""
787 811
    Trims local extrema in greyscale values by a specified amount.
788 812
@@ -810,14 +834,14 @@
Loading
810 834
811 835
    """
812 836
    result = im
813 -
    if mode in ['maxima', 'extrema']:
814 -
        result = reconstruction(seed=im - h, mask=im, method='dilation')
815 -
    elif mode in ['minima', 'extrema']:
816 -
        result = reconstruction(seed=im + h, mask=im, method='erosion')
837 +
    if mode in ["maxima", "extrema"]:
838 +
        result = reconstruction(seed=im - h, mask=im, method="dilation")
839 +
    elif mode in ["minima", "extrema"]:
840 +
        result = reconstruction(seed=im + h, mask=im, method="erosion")
817 841
    return result
818 842
819 843
820 -
def flood(im, regions=None, mode='max'):
844 +
def flood(im, regions=None, mode="max"):
821 845
    r"""
822 846
    Floods/fills each region in an image with a single value based on the
823 847
    specific values in that region.
@@ -864,16 +888,16 @@
Loading
864 888
    else:
865 889
        labels = np.copy(regions)
866 890
        N = labels.max()
867 -
    mode = 'sum' if mode == 'size' else mode
868 -
    mode = 'maximum' if mode == 'max' else mode
869 -
    mode = 'minimum' if mode == 'min' else mode
870 -
    if mode in ['mean', 'median', 'maximum', 'minimum', 'sum']:
891 +
    mode = "sum" if mode == "size" else mode
892 +
    mode = "maximum" if mode == "max" else mode
893 +
    mode = "minimum" if mode == "min" else mode
894 +
    if mode in ["mean", "median", "maximum", "minimum", "sum"]:
871 895
        f = getattr(spim, mode)
872 -
        vals = f(input=im, labels=labels, index=range(0, N+1))
896 +
        vals = f(input=im, labels=labels, index=range(0, N + 1))
873 897
        im_flooded = vals[labels]
874 -
        im_flooded = im_flooded*mask
898 +
        im_flooded = im_flooded * mask
875 899
    else:
876 -
        raise Exception(mode + ' is not a recognized mode')
900 +
        raise Exception(mode + " is not a recognized mode")
877 901
    return im_flooded
878 902
879 903
@@ -901,10 +925,11 @@
Loading
901 925
        the image.  Obviously, voxels with a value of zero have no error.
902 926
903 927
    """
904 -
    temp = np.ones(shape=dt.shape)*np.inf
928 +
    temp = np.ones(shape=dt.shape) * np.inf
905 929
    for ax in range(dt.ndim):
906 -
        dt_lin = distance_transform_lin(np.ones_like(temp, dtype=bool),
907 -
                                        axis=ax, mode='both')
930 +
        dt_lin = distance_transform_lin(
931 +
            np.ones_like(temp, dtype=bool), axis=ax, mode="both"
932 +
        )
908 933
        temp = np.minimum(temp, dt_lin)
909 934
    result = np.clip(dt - temp, a_min=0, a_max=np.inf)
910 935
    return result
@@ -977,20 +1002,27 @@
Loading
977 1002
978 1003
    """
979 1004
    if im.ndim != im.squeeze().ndim:
980 -
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
981 -
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
982 -
                      ' unexpected behavior.')
1005 +
        warnings.warn(
1006 +
            "Input image conains a singleton axis:"
1007 +
            + str(im.shape)
1008 +
            + " Reduce dimensionality with np.squeeze(im) to avoid"
1009 +
            + " unexpected behavior."
1010 +
        )
983 1011
    if spacing < 0:
984 -
        raise Exception('Spacing cannot be less than 0')
1012 +
        raise Exception("Spacing cannot be less than 0")
985 1013
    if spacing == 0:
986 1014
        label = True
987 1015
    result = np.zeros(im.shape, dtype=int)  # Will receive chords at end
988 -
    slxyz = [slice(None, None, spacing*(axis != i) + 1) for i in [0, 1, 2]]
989 -
    slices = tuple(slxyz[:im.ndim])
1016 +
    slxyz = [slice(None, None, spacing * (axis != i) + 1) for i in [0, 1, 2]]
1017 +
    slices = tuple(slxyz[: im.ndim])
990 1018
    s = [[0, 1, 0], [0, 1, 0], [0, 1, 0]]  # Straight-line structuring element
991 1019
    if im.ndim == 3:  # Make structuring element 3D if necessary
992 -
        s = np.pad(np.atleast_3d(s), pad_width=((0, 0), (0, 0), (1, 1)),
993 -
                   mode='constant', constant_values=0)
1020 +
        s = np.pad(
1021 +
            np.atleast_3d(s),
1022 +
            pad_width=((0, 0), (0, 0), (1, 1)),
1023 +
            mode="constant",
1024 +
            constant_values=0,
1025 +
        )
994 1026
    im = im[slices]
995 1027
    s = np.swapaxes(s, 0, axis)
996 1028
    chords = spim.label(im, structure=s)[0]
@@ -1040,25 +1072,28 @@
Loading
1040 1072
1041 1073
    """
1042 1074
    if im.ndim != im.squeeze().ndim:
1043 -
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
1044 -
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
1045 -
                      ' unexpected behavior.')
1075 +
        warnings.warn(
1076 +
            "Input image conains a singleton axis:"
1077 +
            + str(im.shape)
1078 +
            + " Reduce dimensionality with np.squeeze(im) to avoid"
1079 +
            + " unexpected behavior."
1080 +
        )
1046 1081
    if im.ndim < 3:
1047 -
        raise Exception('Must be a 3D image to use this function')
1082 +
        raise Exception("Must be a 3D image to use this function")
1048 1083
    if spacing < 0:
1049 -
        raise Exception('Spacing cannot be less than 0')
1084 +
        raise Exception("Spacing cannot be less than 0")
1050 1085
    ch = np.zeros_like(im, dtype=int)
1051 -
    ch[:, ::4+2*spacing, ::4+2*spacing] = 1  # X-direction
1052 -
    ch[::4+2*spacing, :, 2::4+2*spacing] = 2  # Y-direction
1053 -
    ch[2::4+2*spacing, 2::4+2*spacing, :] = 3  # Z-direction
1054 -
    chords = ch*im
1086 +
    ch[:, :: 4 + 2 * spacing, :: 4 + 2 * spacing] = 1       # X-direction
1087 +
    ch[:: 4 + 2 * spacing, :, 2::4 + 2 * spacing] = 2     # Y-direction
1088 +
    ch[2::4 + 2 * spacing, 2::4 + 2 * spacing, :] = 3   # Z-direction
1089 +
    chords = ch * im
1055 1090
    if trim_edges:
1056 1091
        temp = clear_border(spim.label(chords > 0)[0]) > 0
1057 -
        chords = temp*chords
1092 +
        chords = temp * chords
1058 1093
    return chords
1059 1094
1060 1095
1061 -
def local_thickness(im, sizes=25, mode='hybrid'):
1096 +
def local_thickness(im, sizes=25, mode="hybrid"):
1062 1097
    r"""
1063 1098
    For each voxel, this function calculates the radius of the largest sphere
1064 1099
    that both engulfs the voxel and fits entirely within the foreground.
@@ -1123,8 +1158,7 @@
Loading
1123 1158
    return im_new
1124 1159
1125 1160
1126 -
def porosimetry(im, sizes=25, inlets=None, access_limited=True,
1127 -
                mode='hybrid'):
1161 +
def porosimetry(im, sizes=25, inlets=None, access_limited=True, mode="hybrid"):
1128 1162
    r"""
1129 1163
    Performs a porosimetry simulution on the image
1130 1164
@@ -1199,14 +1233,17 @@
Loading
1199 1233
1200 1234
    """
1201 1235
    if im.ndim != im.squeeze().ndim:
1202 -
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
1203 -
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
1204 -
                      ' unexpected behavior.')
1236 +
        warnings.warn(
1237 +
            "Input image conains a singleton axis:"
1238 +
            + str(im.shape)
1239 +
            + " Reduce dimensionality with np.squeeze(im) to avoid"
1240 +
            + " unexpected behavior."
1241 +
        )
1205 1242
1206 1243
    dt = edt(im > 0)
1207 1244
1208 1245
    if inlets is None:
1209 -
        inlets = get_border(im.shape, mode='faces')
1246 +
        inlets = get_border(im.shape, mode="faces")
1210 1247
1211 1248
    if isinstance(sizes, int):
1212 1249
        sizes = np.logspace(start=np.log10(np.amax(dt)), stop=0, num=sizes)
@@ -1218,40 +1255,40 @@
Loading
1218 1255
    else:
1219 1256
        strel = ps_ball
1220 1257
1221 -
    if mode == 'mio':
1258 +
    if mode == "mio":
1222 1259
        pw = int(np.floor(dt.max()))
1223 -
        impad = np.pad(im, mode='symmetric', pad_width=pw)
1224 -
        inlets = np.pad(inlets, mode='symmetric', pad_width=pw)
1225 -
#        sizes = np.unique(np.around(sizes, decimals=0).astype(int))[-1::-1]
1260 +
        impad = np.pad(im, mode="symmetric", pad_width=pw)
1261 +
        inlets = np.pad(inlets, mode="symmetric", pad_width=pw)
1262 +
        # sizes = np.unique(np.around(sizes, decimals=0).astype(int))[-1::-1]
1226 1263
        imresults = np.zeros(np.shape(impad))
1227 -
        for r in tqdm(sizes):
1228 -
            imtemp = fftmorphology(impad, strel(r), mode='erosion')
1264 +
        for r in tqdm(sizes, file=sys.stdout):
1265 +
            imtemp = fftmorphology(impad, strel(r), mode="erosion")
1229 1266
            if access_limited:
1230 1267
                imtemp = trim_disconnected_blobs(imtemp, inlets)
1231 -
            imtemp = fftmorphology(imtemp, strel(r), mode='dilation')
1268 +
            imtemp = fftmorphology(imtemp, strel(r), mode="dilation")
1232 1269
            if np.any(imtemp):
1233 -
                imresults[(imresults == 0)*imtemp] = r
1270 +
                imresults[(imresults == 0) * imtemp] = r
1234 1271
        imresults = extract_subsection(imresults, shape=im.shape)
1235 -
    elif mode == 'dt':
1272 +
    elif mode == "dt":
1236 1273
        imresults = np.zeros(np.shape(im))
1237 -
        for r in tqdm(sizes):
1274 +
        for r in tqdm(sizes, file=sys.stdout):
1238 1275
            imtemp = dt >= r
1239 1276
            if access_limited:
1240 1277
                imtemp = trim_disconnected_blobs(imtemp, inlets)
1241 1278
            if np.any(imtemp):
1242 1279
                imtemp = edt(~imtemp) < r
1243 -
                imresults[(imresults == 0)*imtemp] = r
1244 -
    elif mode == 'hybrid':
1280 +
                imresults[(imresults == 0) * imtemp] = r
1281 +
    elif mode == "hybrid":
1245 1282
        imresults = np.zeros(np.shape(im))
1246 -
        for r in tqdm(sizes):
1283 +
        for r in tqdm(sizes, file=sys.stdout):
1247 1284
            imtemp = dt >= r
1248 1285
            if access_limited:
1249 1286
                imtemp = trim_disconnected_blobs(imtemp, inlets)
1250 1287
            if np.any(imtemp):
1251 -
                imtemp = fftconvolve(imtemp, strel(r), mode='same') > 0.0001
1252 -
                imresults[(imresults == 0)*imtemp] = r
1288 +
                imtemp = fftconvolve(imtemp, strel(r), mode="same") > 0.0001
1289 +
                imresults[(imresults == 0) * imtemp] = r
1253 1290
    else:
1254 -
        raise Exception('Unreckognized mode ' + mode)
1291 +
        raise Exception("Unreckognized mode " + mode)
1255 1292
    return imresults
1256 1293
1257 1294
@@ -1287,8 +1324,9 @@
Loading
1287 1324
    elif (inlets.shape == im.shape) and (inlets.max() == 1):
1288 1325
        inlets = inlets.astype(bool)
1289 1326
    else:
1290 -
        raise Exception('inlets not valid, refer to docstring for info')
1327 +
        raise Exception("inlets not valid, refer to docstring for info")
1291 1328
    from skimage.morphology import square
1329 +
1292 1330
    if im.ndim == 3:
1293 1331
        square = cube
1294 1332
    else:
@@ -1300,15 +1338,15 @@
Loading
1300 1338
        im2 = np.reshape(np.in1d(labels, keep), newshape=im.shape)
1301 1339
    else:
1302 1340
        im2 = np.zeros_like(im)
1303 -
    im2 = im2*im
1341 +
    im2 = im2 * im
1304 1342
    return im2
1305 1343
1306 1344
1307 1345
def _get_axial_shifts(ndim=2, include_diagonals=False):
1308 -
    r'''
1346 +
    r"""
1309 1347
    Helper function to generate the axial shifts that will be performed on
1310 1348
    the image to identify bordering pixels/voxels
1311 -
    '''
1349 +
    """
1312 1350
    if ndim == 2:
1313 1351
        if include_diagonals:
1314 1352
            neighbors = square(3)
@@ -1333,40 +1371,37 @@
Loading
1333 1371
1334 1372
1335 1373
def _make_stack(im, include_diagonals=False):
1336 -
    r'''
1374 +
    r"""
1337 1375
    Creates a stack of images with one extra dimension to the input image
1338 1376
    with length equal to the number of borders to search + 1.
1339 1377
    Image is rolled along the axial shifts so that the border pixel is
1340 1378
    overlapping the original pixel. First image in stack is the original.
1341 1379
    Stacking makes direct vectorized array comparisons possible.
1342 -
    '''
1380 +
    """
1343 1381
    ndim = len(np.shape(im))
1344 1382
    axial_shift = _get_axial_shifts(ndim, include_diagonals)
1345 1383
    if ndim == 2:
1346 -
        stack = np.zeros([np.shape(im)[0],
1347 -
                          np.shape(im)[1],
1348 -
                          len(axial_shift)+1])
1384 +
        stack = np.zeros([np.shape(im)[0], np.shape(im)[1], len(axial_shift) + 1])
1349 1385
        stack[:, :, 0] = im
1350 1386
        for i in range(len(axial_shift)):
1351 1387
            ax0, ax1 = axial_shift[i]
1352 1388
            temp = np.roll(np.roll(im, ax0, 0), ax1, 1)
1353 -
            stack[:, :, i+1] = temp
1389 +
            stack[:, :, i + 1] = temp
1354 1390
        return stack
1355 1391
    elif ndim == 3:
1356 -
        stack = np.zeros([np.shape(im)[0],
1357 -
                          np.shape(im)[1],
1358 -
                          np.shape(im)[2],
1359 -
                          len(axial_shift)+1])
1392 +
        stack = np.zeros(
1393 +
            [np.shape(im)[0], np.shape(im)[1], np.shape(im)[2], len(axial_shift) + 1]
1394 +
        )
1360 1395
        stack[:, :, :, 0] = im
1361 1396
        for i in range(len(axial_shift)):
1362 1397
            ax0, ax1, ax2 = axial_shift[i]
1363 1398
            temp = np.roll(np.roll(np.roll(im, ax0, 0), ax1, 1), ax2, 2)
1364 -
            stack[:, :, :, i+1] = temp
1399 +
            stack[:, :, :, i + 1] = temp
1365 1400
        return stack
1366 1401
1367 1402
1368 1403
def nphase_border(im, include_diagonals=False):
1369 -
    r'''
1404 +
    r"""
1370 1405
    Identifies the voxels in regions that border *N* other regions.
1371 1406
1372 1407
    Useful for finding triple-phase boundaries.
@@ -1387,17 +1422,20 @@
Loading
1387 1422
    image : ND-array
1388 1423
        A copy of ``im`` with voxel values equal to the number of uniquely
1389 1424
        different bordering values
1390 -
    '''
1425 +
    """
1391 1426
    if im.ndim != im.squeeze().ndim:
1392 -
        warnings.warn('Input image conains a singleton axis:' + str(im.shape) +
1393 -
                      ' Reduce dimensionality with np.squeeze(im) to avoid' +
1394 -
                      ' unexpected behavior.')
1427 +
        warnings.warn(
1428 +
            "Input image conains a singleton axis:"
1429 +
            + str(im.shape)
1430 +
            + " Reduce dimensionality with np.squeeze(im) to avoid"
1431 +
            + " unexpected behavior."
1432 +
        )
1395 1433
    # Get dimension of image
1396 1434
    ndim = len(np.shape(im))
1397 1435
    if ndim not in [2, 3]:
1398 1436
        raise NotImplementedError("Function only works for 2d and 3d images")
1399 1437
    # Pad image to handle edges
1400 -
    im = np.pad(im, pad_width=1, mode='edge')
1438 +
    im = np.pad(im, pad_width=1, mode="edge")
1401 1439
    # Stack rolled images for each neighbor to be inspected
1402 1440
    stack = _make_stack(im, include_diagonals)
1403 1441
    # Sort the stack along the last axis
@@ -1407,9 +1445,9 @@
Loading
1407 1445
    # Number of changes is number of unique bordering regions
1408 1446
    for k in range(np.shape(stack)[ndim])[1:]:
1409 1447
        if ndim == 2:
1410 -
            mask = stack[:, :, k] != stack[:, :, k-1]
1448 +
            mask = stack[:, :, k] != stack[:, :, k - 1]
1411 1449
        elif ndim == 3:
1412 -
            mask = stack[:, :, :, k] != stack[:, :, :, k-1]
1450 +
            mask = stack[:, :, :, k] != stack[:, :, :, k - 1]
1413 1451
        out += mask
1414 1452
    # Un-pad
1415 1453
    if ndim == 2:
@@ -1447,12 +1485,12 @@
Loading
1447 1485
    im_result = np.zeros_like(skel)
1448 1486
    # If branch points are not supplied, attempt to find them
1449 1487
    if branch_points is None:
1450 -
        branch_points = spim.convolve(skel*1.0, weights=cube(3)) > 3
1451 -
        branch_points = branch_points*skel
1488 +
        branch_points = spim.convolve(skel * 1.0, weights=cube(3)) > 3
1489 +
        branch_points = branch_points * skel
1452 1490
    # Store original branch points before dilating
1453 1491
    pts_orig = branch_points
1454 1492
    # Find arcs of skeleton by deleting branch points
1455 -
    arcs = skel*(~branch_points)
1493 +
    arcs = skel * (~branch_points)
1456 1494
    # Label arcs
1457 1495
    arc_labels = spim.label(arcs, structure=cube(3))[0]
1458 1496
    # Dilate branch points so they overlap with the arcs
@@ -1464,26 +1502,32 @@
Loading
1464 1502
    for s in slices:
1465 1503
        label_num += 1
1466 1504
        # Find branch point labels the overlap current arc
1467 -
        hits = pts_labels[s]*(arc_labels[s] == label_num)
1505 +
        hits = pts_labels[s] * (arc_labels[s] == label_num)
1468 1506
        # If image contains 2 branch points, then it's not a tail.
1469 1507
        if len(np.unique(hits)) == 3:
1470 1508
            im_result[s] += arc_labels[s] == label_num
1471 1509
    # Add missing branch points back to arc image to make complete skeleton
1472 -
    im_result += skel*pts_orig
1510 +
    im_result += skel * pts_orig
1473 1511
    if iterations > 1:
1474 1512
        iterations -= 1
1475 1513
        im_temp = np.copy(im_result)
1476 -
        im_result = prune_branches(skel=im_result,
1477 -
                                   branch_points=None,
1478 -
                                   iterations=iterations)
1514 +
        im_result = prune_branches(
1515 +
            skel=im_result, branch_points=None, iterations=iterations
1516 +
        )
1479 1517
        if np.all(im_temp == im_result):
1480 1518
            iterations = 0
1481 1519
    return im_result
1482 1520
1483 1521
1484 -
def chunked_func(func, overlap=None, divs=2, cores=None,
1485 -
                 im_arg=['input', 'image', 'im'],
1486 -
                 strel_arg=['strel', 'structure', 'footprint'], **kwargs):
1522 +
def chunked_func(
1523 +
    func,
1524 +
    overlap=None,
1525 +
    divs=2,
1526 +
    cores=None,
1527 +
    im_arg=["input", "image", "im"],
1528 +
    strel_arg=["strel", "structure", "footprint"],
1529 +
    **kwargs
1530 +
):
1487 1531
    r"""
1488 1532
    Performs the specfied operation "chunk-wise" in parallel
1489 1533
@@ -1564,12 +1608,15 @@
Loading
1564 1608
    True
1565 1609
1566 1610
    """
1611 +
1567 1612
    @dask.delayed
1568 1613
    def apply_func(func, **kwargs):
1569 1614
        # Apply function on sub-slice of overall image
1570 1615
        return func(**kwargs)
1616 +
1571 1617
    # Import the array_split methods
1572 1618
    from array_split import shape_split, ARRAY_BOUNDS
1619 +
1573 1620
    # Determine the value for im_arg
1574 1621
    if type(im_arg) == str:
1575 1622
        im_arg = [im_arg]
@@ -1581,10 +1628,10 @@
Loading
1581 1628
    # Fetch image from the kwargs dict
1582 1629
    im = kwargs[im_arg]
1583 1630
    # Determine the number of divisions to create
1584 -
    divs = np.ones((im.ndim, ), dtype=int)*np.array(divs)
1631 +
    divs = np.ones((im.ndim,), dtype=int) * np.array(divs)
1585 1632
    # If overlap given then use it, otherwise search for strel in kwargs
1586 1633
    if overlap is not None:
1587 -
        halo = overlap*(divs > 1)
1634 +
        halo = overlap * (divs > 1)
1588 1635
    else:
1589 1636
        if type(strel_arg) == str:
1590 1637
            strel_arg = [strel_arg]
@@ -1593,9 +1640,11 @@
Loading
1593 1640
                strel = kwargs[item]
1594 1641
                break
1595 1642
        halo = np.array(strel.shape) * (divs > 1)
1596 -
    slices = np.ravel(shape_split(im.shape, axis=divs,
1597 -
                                  halo=halo.tolist(),
1598 -
                                  tile_bounds_policy=ARRAY_BOUNDS))
1643 +
    slices = np.ravel(
1644 +
        shape_split(
1645 +
            im.shape, axis=divs, halo=halo.tolist(), tile_bounds_policy=ARRAY_BOUNDS
1646 +
        )
1647 +
    )
1599 1648
    # Apply func to each subsection of the image
1600 1649
    res = []
1601 1650
    # print('Image will be broken into the following chunks:')

@@ -1,3 +1,4 @@
Loading
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 @@
Loading
35 36
    directly to an OpenPNM network object using the ``update`` command.
36 37
37 38
    """
38 -
    print('_'*60)
39 +
    print('-'*60)
39 40
    print('Extracting pore and throat information from image')
40 41
    from skimage.morphology import disk, ball
41 42
    struc_elem = disk if im.ndim == 2 else ball
@@ -67,7 +68,7 @@
Loading
67 68
    # dt_shape = np.array(dt.shape)
68 69
69 70
    # Start extracting size information for pores and throats
70 -
    for i in tqdm(Ps):
71 +
    for i in tqdm(Ps, file=sys.stdout):
71 72
        pore = i - 1
72 73
        if slices[pore] is None:
73 74
            continue

@@ -3,9 +3,9 @@
Loading
3 3
4 4
def set_mpl_style():
5 5
6 -
    sfont = 25
7 -
    mfont = 35
8 -
    lfont = 45
6 +
    sfont = 15
7 +
    mfont = 15
8 +
    lfont = 15
9 9
10 10
    line_props = {'linewidth': 4,
11 11
                  'markersize': 10}

@@ -1,10 +1,11 @@
Loading
1 +
import sys
1 2
import numpy as np
2 -
import scipy.ndimage as spim
3 3
from tqdm import tqdm
4 +
import scipy.ndimage as spim
4 5
from porespy.tools import extract_subsection, bbox_to_slices
5 -
from skimage.measure import regionprops
6 6
from skimage.measure import mesh_surface_area, marching_cubes_lewiner
7 7
from skimage.morphology import skeletonize_3d, ball
8 +
from skimage.measure import regionprops
8 9
from pandas import DataFrame
9 10
10 11
@@ -176,11 +177,11 @@
Loading
176 177
    which may be helpful.
177 178
178 179
    """
179 -
    print('_'*60)
180 +
    print('-'*60)
180 181
    print('Calculating regionprops')
181 182
182 183
    results = regionprops(im)
183 -
    for i in tqdm(range(len(results))):
184 +
    for i in tqdm(range(len(results)), file=sys.stdout):
184 185
        mask = results[i].image
185 186
        mask_padded = np.pad(mask, pad_width=1, mode='constant')
186 187
        temp = spim.distance_transform_edt(mask_padded)

@@ -1,8 +1,9 @@
Loading
1 1
import numpy as np
2 +
from edt import edt
3 +
import porespy as ps
2 4
import scipy.spatial as sptl
3 5
import scipy.ndimage as spim
4 6
from numba import njit, jit
5 -
from edt import edt
6 7
from porespy.tools import norm_to_uniform, ps_ball, ps_disk, get_border
7 8
from porespy.tools import insert_sphere, fftmorphology
8 9
from typing import List
@@ -60,8 +61,8 @@
Loading
60 61
        for dim in range(im.ndim):
61 62
            r, d = np.divmod(element.shape[dim], 2)
62 63
            if d == 0:
63 -
                raise Exception('Cannot specify center point when element ' +
64 -
                                'has one or more even dimension')
64 +
                raise Exception('Cannot specify center point when element '
65 +
                                + 'has one or more even dimension')
65 66
            lower_im = np.amax((center[dim] - r, 0))
66 67
            upper_im = np.amin((center[dim] + r + 1, im.shape[dim]))
67 68
            s_im.append(slice(lower_im, upper_im))
@@ -153,8 +154,8 @@
Loading
153 154
    [1] Random Heterogeneous Materials, S. Torquato (2001)
154 155
155 156
    """
156 -
    print(78*'―')
157 -
    print('RSA: Adding spheres of size ' + str(radius))
157 +
    print(80*'-')
158 +
    print(f'RSA: Adding spheres of size {radius}')
158 159
    im = im.astype(bool)
159 160
    if n_max is None:
160 161
        n_max = 10000

@@ -651,7 +651,7 @@
Loading
651 651
        area shared by regions 0 and 5.
652 652
653 653
    """
654 -
    print('_'*60)
654 +
    print('-'*60)
655 655
    print('Finding interfacial areas between each region')
656 656
    from skimage.morphology import disk, ball
657 657
    im = regions.copy()
@@ -735,7 +735,7 @@
Loading
735 735
        that the surface area of region 1 is stored in element 0 of the list.
736 736
737 737
    """
738 -
    print('_'*60)
738 +
    print('-'*60)
739 739
    print('Finding surface area of each region')
740 740
    im = regions.copy()
741 741
    # Get 'slices' into im for each pore region
Files Coverage
porespy 87.9%
docs/conf.py 100.0%
examples.py 100.0%
Project Totals (17 files) 88.1%
codecov-umbrella
Build #73179761 -
unittests
1
codecov:
2
  branch: master
3

4
coverage:
5
  precision: 2
6
  round: down
7
  range: "50...100"
8

9
  status:
10
    project:
11
      default:
12
        target: auto
13
        threshold: null
14
        branches: null
15

16
    patch:
17
      default:
18
        target: auto
19
        branches: null
20

21
comment:
22
  layout: "header, diff, changes, sunburst, uncovered"
23
  branches: null
24
  behavior: default
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading