Merge pull request #1045 from PMEAL/play_w_blobs

Added `periodic` argument to `blobs` which makes a periodic image that also has uniform density near edges #api

src/porespy/generators/_imgen.py

@@ -608,7 +608,8 @@ def polydisperse_spheres(
     dist,
     nbins: int = 5,
     r_min: int = 5,
-    seed=None):
+    seed=None,
+):
     r"""
     Create an image of randomly placed, overlapping spheres with a
     distribution of radii.
@@ -1027,7 +1028,8 @@ def blobs(
     porosity: float = 0.5,
     blobiness: int = 1,
     divs: int = 1,
-    seed=None,
+    seed: int = None,
+    periodic: bool = True,
 ):
     """
     Generates an image containing amorphous blobs
@@ -1042,7 +1044,7 @@ def blobs(
         prior to returning.  If ``None`` is specified, then the scalar
         noise field is converted to a uniform distribution and returned
         without thresholding.
-    blobiness : int or list of ints(default = 1)
+    blobiness : int or list of ints (default = 1)
         Controls the morphology of the blobs.  A higher number results in
         a larger number of small blobs.  If a list is supplied then the
         blobs are anisotropic.
@@ -1052,10 +1054,14 @@ def blobs(
         equivalent to ``[2, 2, 2]`` for a 3D image.  The number of cores
         used is specified in ``porespy.settings.ncores`` and defaults to
         all cores.
-    seed : int, optional, default = `None`
+    seed : int, default = `None`
         Initializes numpy's random number generator to the specified state. If not
         provided, the current global value is used. This means calls to
         ``np.random.state(seed)`` prior to calling this function will be respected.
+    periodic : bool, default = `True`
+        If `True` the blobs will be periodic, meaning that the image can be tiled
+        and the phases will be continuous. `False` will provide the "legacy" version
+        of an image, which has high-porosity artifacts at the image boundaries.
 
     Returns
     -------
@@ -1093,6 +1099,7 @@ def blobs(
     if isinstance(blobiness, int):
         blobiness = [blobiness]*len(shape)
     blobiness = np.array(blobiness)
+    mode = 'wrap' if periodic else 'reflect'
     parallel = False
     if isinstance(divs, int):
         divs = [divs]*len(shape)
@@ -1107,7 +1114,7 @@ def blobs(
                           input=im, sigma=sigma,
                           divs=divs, overlap=overlap)
     else:
-        im = spim.gaussian_filter(im, sigma=sigma)
+        im = spim.gaussian_filter(im, sigma=sigma, mode=mode)
     im = all_to_uniform(im, scale=[0, 1])
     if porosity:
         im = im < porosity

src/porespy/tools/_funcs.py

@@ -2,6 +2,7 @@
 
 import numpy as np
 import scipy.ndimage as spim
+from scipy.stats import rankdata
 from numba import boolean, njit
 from skimage.morphology import ball, disk
 from skimage.segmentation import relabel_sequential
@@ -1097,6 +1098,9 @@ def all_to_uniform(im, scale=None):
     """
     if scale is None:
         scale = [im.min(), im.max()]
+    # Alternative, might be faster
+    # im2 = rankdata(im).reshape(im.shape)
+    # im = (im2 - im2.min())/(im2.max() - im2.min())*(scale[1] - scale[0]) + scale[0]
     aargsort_im = np.argsort(np.argsort(im.flatten()))  # Twice for the inverse permutation
     linspace_im = np.linspace(scale[0], scale[1], len(aargsort_im), endpoint=True)
     uniform_flatten_im = linspace_im[aargsort_im]

test/integration/test_drainage.py

@@ -13,6 +13,7 @@ def test_drainage(plot=False):
         porosity=0.708328,
         blobiness=1.5,
         seed=6,
+        periodic=False,
     )
     inlets = np.zeros_like(im)
     inlets[0, :] = True

test/integration/test_ibip_example_script.py

@@ -9,7 +9,12 @@ def test_ibip():
     np.random.seed(0)
 
     # Generate or load a test image
-    im = ps.generators.blobs(shape=[200, 200], porosity=0.605475, blobiness=2)
+    im = ps.generators.blobs(
+        shape=[200, 200],
+        porosity=0.605475,
+        blobiness=2,
+        periodic=False,
+    )
 
     bd = np.zeros_like(im)
     bd[:, 0] = True

test/unit/test_dns.py

@@ -27,12 +27,13 @@ def test_tortuosity_different_solvers(self):
         np.testing.assert_allclose(t.tortuosity, 1.35995, rtol=1e-4)
 
     def test_exception_if_no_pores_remain_after_trimming_floating_pores(self):
-        im = ps.generators.blobs(shape=[200, 200], porosity=0.05)
+        im = ps.generators.blobs(shape=[200, 200], porosity=0.05, periodic=False,)
         with pytest.raises(Exception):
             _ = ps.simulations.tortuosity_fd(im=im, axis=1)
 
     def test_flux(self):
-        im = ps.generators.blobs(shape=[15, 20, 25], porosity=0.85, blobiness=1.5)
+        im = ps.generators.blobs(
+            shape=[15, 20, 25], porosity=0.85, blobiness=1.5, periodic=False,)
         for axis in range(3):
             out = ps.simulations.tortuosity_fd(im, axis=axis)
             c = out["im_conc"]
@@ -43,7 +44,8 @@ def test_flux(self):
             np.testing.assert_allclose(rate, rate[0], rtol=1e-5)
 
     def test_tau_from_cmap(self):
-        im = ps.generators.blobs(shape=[15, 20, 25], porosity=0.85, blobiness=1.5)
+        im = ps.generators.blobs(
+            shape=[15, 20, 25], porosity=0.85, blobiness=1.5, periodic=False,)
         for axis in range(3):
             out = ps.simulations.tortuosity_fd(im, axis=axis)
             c = out["im_conc"]

test/unit/test_filters.py

@@ -19,7 +19,8 @@ def setup_class(self):
         self.im = ps.generators.blobs(shape=[100, 100, 100],
                                       blobiness=2,
                                       seed=0,
-                                      porosity=0.499829)
+                                      porosity=0.499829,
+                                      periodic=False,)
         # Ensure that im was generated as expected
         assert self.im.sum()/self.im.size == 0.499829
         self.im_dt = edt(self.im)
@@ -129,7 +130,8 @@ def test_find_disconnected_voxels_3d_conn6(self):
 
     def test_trim_nonpercolating_paths_2d_axis0(self):
         np.random.seed(0)
-        im = ps.generators.blobs(shape=[200, 200], porosity=0.55875, blobiness=2)
+        im = ps.generators.blobs(
+            shape=[200, 200], porosity=0.55875, blobiness=2, periodic=False,)
         assert im.sum()/im.size == 0.55875
         inlets = np.zeros_like(im)
         inlets[0, :] = 1
@@ -143,7 +145,8 @@ def test_trim_nonpercolating_paths_2d_axis0(self):
 
     def test_trim_nonpercolating_paths_2d_axis1(self):
         np.random.seed(0)
-        im = ps.generators.blobs(shape=[200, 200], porosity=0.55875, blobiness=2)
+        im = ps.generators.blobs(
+            shape=[200, 200], porosity=0.55875, blobiness=2, periodic=False,)
         assert im.sum()/im.size == 0.55875
         inlets = np.zeros_like(im)
         inlets[:, 0] = 1
@@ -157,7 +160,8 @@ def test_trim_nonpercolating_paths_2d_axis1(self):
 
     def test_trim_nonpercolating_paths_no_paths(self):
         np.random.seed(0)
-        im = ps.generators.blobs([200, 200], porosity=0.2535, blobiness=2)
+        im = ps.generators.blobs(
+            shape=[200, 200], porosity=0.2535, blobiness=2, periodic=False,)
         assert im.sum()/im.size == 0.2535
         inlets = np.zeros_like(im)
         inlets[:, 0] = 1
@@ -171,7 +175,8 @@ def test_trim_nonpercolating_paths_no_paths(self):
 
     def test_trim_nonpercolating_paths_3d_axis2(self):
         np.random.seed(0)
-        im = ps.generators.blobs([100, 100, 100], porosity=0.550191, blobiness=2)
+        im = ps.generators.blobs(
+            shape=[100, 100, 100], porosity=0.550191, blobiness=2, periodic=False,)
         assert im.sum()/im.size == 0.550191
         inlets = np.zeros_like(im)
         inlets[..., 0] = 1
@@ -185,7 +190,8 @@ def test_trim_nonpercolating_paths_3d_axis2(self):
 
     def test_trim_nonpercolating_paths_3d_axis1(self):
         np.random.seed(0)
-        im = ps.generators.blobs([100, 100, 100], porosity=0.550191, blobiness=2)
+        im = ps.generators.blobs(
+            shape=[100, 100, 100], porosity=0.550191, blobiness=2, periodic=False,)
         assert im.sum()/im.size == 0.550191
         inlets = np.zeros_like(im)
         inlets[:, 0, :] = 1
@@ -199,7 +205,8 @@ def test_trim_nonpercolating_paths_3d_axis1(self):
 
     def test_trim_nonpercolating_paths_3d_axis0(self):
         np.random.seed(0)
-        im = ps.generators.blobs([100, 100, 100], porosity=0.550191, blobiness=2)
+        im = ps.generators.blobs(
+            shape=[100, 100, 100], porosity=0.550191, blobiness=2, periodic=False,)
         assert im.sum()/im.size == 0.550191
         inlets = np.zeros_like(im)
         inlets[0, ...] = 1
@@ -214,7 +221,8 @@ def test_trim_nonpercolating_paths_3d_axis0(self):
     def test_trim_disconnected_blobs(self):
         s = disk(1)
         np.random.seed(0)
-        im = ps.generators.blobs([200, 200], porosity=0.55875, blobiness=2)
+        im = ps.generators.blobs(
+            shape=[200, 200], porosity=0.55875, blobiness=2, periodic=False,)
         assert im.sum()/im.size == 0.55875
         inlets = np.zeros_like(im)
         inlets[0, ...] = 1
@@ -239,7 +247,8 @@ def test_fill_blind_pores_w_surface(self):
         assert im3.sum() == 0
 
     def test_fill_blind_pores_surface_blobs_2D(self):
-        im = ps.generators.blobs([100, 100], porosity=0.6021, seed=0)
+        im = ps.generators.blobs(
+            shape=[100, 100], porosity=0.6021, seed=0, periodic=False,)
         assert im.sum()/im.size == 0.6021
         im2 = ps.filters.fill_blind_pores(im)
         assert im.sum() == 6021
@@ -248,7 +257,8 @@ def test_fill_blind_pores_surface_blobs_2D(self):
         assert im3.sum() < im2.sum()
 
     def test_fill_blind_pores_surface_blobs_3D(self):
-        im = ps.generators.blobs([100, 100, 100], porosity=0.497569, seed=0)
+        im = ps.generators.blobs(
+            shape=[100, 100, 100], porosity=0.497569, seed=0, periodic=False,)
         assert im.sum()/im.size == 0.497569
         im2 = ps.filters.fill_blind_pores(im, surface=True)
         labels, N = spim.label(im2, ps.tools.ps_rect(3, ndim=3))
@@ -411,7 +421,8 @@ def test_find_dt_artifacts(self):
         assert np.all(dt[inds] - ar[inds] == 1)
 
     def test_snow_partitioning_n_2D(self):
-        im = ps.generators.blobs([500, 500], porosity=0.494604, blobiness=1, seed=0)
+        im = ps.generators.blobs(
+            shape=[500, 500], porosity=0.494604, blobiness=1, seed=0, periodic=False,)
         assert im.sum()/im.size == 0.494604
         snow = ps.filters.snow_partitioning_n(im + 1, r_max=4, sigma=0.4)
         assert np.amax(snow.regions) == 136
@@ -423,7 +434,8 @@ def test_snow_partitioning_n_3D(self):
         im = ps.generators.blobs(shape=[100, 100, 100],
                                  porosity=0.495157,
                                  blobiness=0.75,
-                                 seed=0)
+                                 seed=0,
+                                 periodic=False,)
         assert im.sum()/im.size == 0.495157
         snow = ps.filters.snow_partitioning_n(im + 1, r_max=4, sigma=0.4)
         assert np.amax(snow.regions) == 620
@@ -478,7 +490,8 @@ def test_chunked_func_3d_w_strel(self):
 
     def test_chunked_func_w_ill_defined_filter(self):
         import scipy.signal as spsg
-        im = ps.generators.blobs(shape=[100, 100, 100], porosity=0.497569, seed=0)
+        im = ps.generators.blobs(
+            shape=[100, 100, 100], porosity=0.497569, seed=0, periodic=False,)
         assert im.sum()/im.size == 0.497569
         with pytest.raises(IndexError):
             ps.filters.chunked_func(func=spsg.convolve,
@@ -497,7 +510,7 @@ def test_prune_branches(self):
 
     def test_prune_branches_n2(self):
         from skimage.morphology import skeletonize
-        im = ps.generators.random_spheres([100, 100, 100], r=4, seed=0)
+        im = ps.generators.random_spheres(shape=[100, 100, 100], r=4, seed=0)
         skel1 = skeletonize(im)
         skel2 = ps.filters.prune_branches(skel1, iterations=1)
         skel3 = ps.filters.prune_branches(skel1, iterations=2)
@@ -508,7 +521,8 @@ def test_prune_branches_n2(self):
 
     def test_apply_padded(self):
         from skimage.morphology import skeletonize
-        im = ps.generators.blobs(shape=[100, 100])
+        im = ps.generators.blobs(
+            shape=[100, 100], periodic=False,)
         skel1 = skeletonize(im)
         skel2 = ps.filters.apply_padded(
             im=im,
@@ -522,7 +536,8 @@ def test_trim_small_clusters(self):
         im = ps.generators.blobs(shape=[100, 100],
                                  blobiness=2,
                                  porosity=0.4028,
-                                 seed=0)
+                                 seed=0,
+                                 periodic=False,)
         assert im.sum()/im.size == 0.4028
         im5 = ps.filters.trim_small_clusters(im=im, size=5)
         im10 = ps.filters.trim_small_clusters(im=im, size=10)
@@ -551,7 +566,8 @@ def test_nl_means_layered(self):
         im = ps.generators.blobs(shape=[50, 50, 50],
                                  porosity=0.492664,
                                  blobiness=0.5,
-                                 seed=0)
+                                 seed=0,
+                                 periodic=False,)
         assert im.sum()/im.size == 0.492664
         np.random.seed(0)
         im2 = random_noise(im)
@@ -564,7 +580,8 @@ def test_trim_nearby_peaks(self):
         im = ps.generators.blobs(shape=[400, 400],
                                  blobiness=[2, 1],
                                  porosity=0.5916375,
-                                 seed=0)
+                                 seed=0,
+                                 periodic=False,)
         assert im.sum()/im.size == 0.5916375
         im_dt = edt(im)
         dt = spim.gaussian_filter(input=im_dt, sigma=0.4)
@@ -582,7 +599,8 @@ def test_trim_nearby_peaks_threshold(self):
         im = ps.generators.blobs(shape=[400, 400],
                                  blobiness=[2, 1],
                                  porosity=0.5916375,
-                                 seed=0)
+                                 seed=0,
+                                 periodic=False,)
         assert im.sum()/im.size == 0.5916375
         im_dt = edt(im)
         dt = spim.gaussian_filter(input=im_dt, sigma=0.4)
@@ -594,20 +612,23 @@ def test_trim_nearby_peaks_threshold(self):
         assert num_peaks_after_far_trim <= num_peaks_after_close_trim
 
     def test_regions_size(self):
-        im = ps.generators.blobs([50, 50], porosity=0.1164, seed=0)
+        im = ps.generators.blobs(
+            shape=[50, 50], porosity=0.1164, seed=0, periodic=False,)
         assert im.sum()/im.size == 0.1164
         s = ps.filters.region_size(im)
         hits = [1, 2, 3, 4, 5, 6, 8, 9, 18, 23, 24, 26, 28, 31]
         assert np.all(hits == np.unique(s)[1:])
         np.random.seed(0)
-        im = ps.generators.blobs([20, 20, 20], porosity=0.121, seed=0)
+        im = ps.generators.blobs(
+            shape=[20, 20, 20], porosity=0.121, seed=0, periodic=False,)
         assert im.sum()/im.size == 0.121
         s = ps.filters.region_size(im)
         hits = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 15, 16, 17, 19, 31, 32, 37]
         assert np.all(hits == np.unique(s)[1:])
 
     def test_find_trapped_regions_return_mask_side_outlet(self):
-        im = ps.generators.blobs(shape=[100, 100], porosity=0.6, seed=7)
+        im = ps.generators.blobs(
+            shape=[100, 100], porosity=0.6, seed=7, periodic=False,)
         inlets = np.zeros_like(im)
         inlets[0, :] = True
         outlets = np.zeros_like(im)
@@ -631,7 +652,8 @@ def test_find_trapped_regions_return_mask_side_outlet(self):
         assert np.all(trp1 == trp2)
 
     def test_find_trapped_regions_return_mask_top_outlet(self):
-        im = ps.generators.blobs(shape=[100, 100], porosity=0.6, seed=7)
+        im = ps.generators.blobs(
+            shape=[100, 100], porosity=0.6, seed=7, periodic=False,)
         inlets = np.zeros_like(im)
         inlets[0, :] = True
         outlets = np.zeros_like(im)
@@ -655,7 +677,8 @@ def test_find_trapped_regions_return_mask_top_outlet(self):
         assert np.all(trp1 == trp2)
 
     def test_find_trapped_regions_top_outlet(self):
-        im = ps.generators.blobs(shape=[100, 100], porosity=0.6, seed=7)
+        im = ps.generators.blobs(
+            shape=[100, 100], porosity=0.6, seed=7, periodic=False,)
         inlets = np.zeros_like(im)
         inlets[0, :] = True
         outlets = np.zeros_like(im)

test/unit/test_filters_size_seq_satn.py

@@ -17,11 +17,13 @@ def setup_class(self):
         self.bd = bd
         self.im2D = ps.generators.blobs(shape=[51, 51],
                                         seed=0,
-                                        porosity=0.48212226066897346)
+                                        porosity=0.48212226066897346,
+                                        periodic=False,)
         assert self.im2D.sum()/self.im2D.size == 0.48212226066897346
         self.im3D = ps.generators.blobs(shape=[51, 51, 51],
                                         seed=0,
-                                        porosity=0.49954391599007925)
+                                        porosity=0.49954391599007925,
+                                        periodic=False,)
         assert self.im3D.sum()/self.im3D.size == 0.49954391599007925
 
     def test_satn_to_seq(self):
@@ -192,7 +194,8 @@ def test_size_to_satn_uninvaded(self):
         assert satn[0, 2] == 0.95
 
     def test_compare_size_and_seq_to_satn(self):
-        im = ps.generators.blobs(shape=[250, 250], seed=0, porosity=0.496064)
+        im = ps.generators.blobs(
+            shape=[250, 250], seed=0, porosity=0.496064, periodic=False,)
         assert im.sum()/im.size == 0.496064
         dt = edt(im)
         sizes = np.arange(int(dt.max())+1, 0, -1)