Repository: PMEAL/porespy

Issue:
Use keyword `method` instead of `mode`
I have noticed that some numpy/scipy functions have (started) using this.  I think I like it and it's slightly more descriptive.  Is it worth making the change?

Discussion:
@jgostick: `'mode'` would be used when a different result is returned, while `'method'` would be a different way to compute the same result.  Example would be `method='tessellation'` vs `method='distance_matrix'`.  @jgostick: We should do a "find all" for the `mode'` in the function definitions, then create a list of all functions where this change might make sense here.  Most function will *probably* remain with `mode`.@jgostick: For example, `find_trapped_regions` currently uses cluster labeling, which works fine/fast for pressure based drainage but for invasion percolation it is too slow, so we have created `find_trapped_regions2`. This should be accessible via `find_trapped_regions` by switching the `method`.  

Solution patch:
diff --git a/src/porespy/filters/_funcs.py b/src/porespy/filters/_funcs.py
index eed041908..e6efc9964 100644
--- a/src/porespy/filters/_funcs.py
+++ b/src/porespy/filters/_funcs.py
@@ -303,11 +303,13 @@ def find_disconnected_voxels(im, conn: int = None, surface: bool = False):
     if not surface:
         holes = clear_border(labels=labels) > 0
     else:
-        keep = set(np.unique(labels))
+        keep = set()
         for ax in range(labels.ndim):
             labels = np.swapaxes(labels, 0, ax)
-            keep.intersection_update(set(np.unique(labels[0, ...])))
-            keep.intersection_update(set(np.unique(labels[-1, ...])))
+            s1 = set(np.unique(labels[0, ...]))
+            s2 = set(np.unique(labels[-1, ...]))
+            tmp = s1.intersection(s2)
+            keep.update(tmp)
             labels = np.swapaxes(labels, 0, ax)
         holes = np.isin(labels, list(keep), invert=True)
     return holes
diff --git a/src/porespy/filters/_size_seq_satn.py b/src/porespy/filters/_size_seq_satn.py
index 956cabdef..dad427a76 100644
--- a/src/porespy/filters/_size_seq_satn.py
+++ b/src/porespy/filters/_size_seq_satn.py
@@ -494,7 +494,7 @@ def f(r, sigma, theta, voxel_size):
     tmp = np.copy(inv.im_seq)
     tmp[tmp == 1] = -1  # Set some to uninvaded/trapped
     tmp[tmp > 0] = tmp.max() - tmp[tmp > 0] + 1  # Invert sequence numbers
-    satn = ps.filters.seq_to_satn(seq=tmp, im=im, mode='imb')
+    satn = ps.filters.seq_to_satn(seq=tmp, im=im, mode='imbibition')
     assert (1 - (satn > 0).sum()/im.sum()) == satn[satn > 0].min()
     # Now convert it back!
     seq = satn_to_seq(satn, im=im, mode='imbibition')
diff --git a/src/porespy/metrics/_funcs.py b/src/porespy/metrics/_funcs.py
index 992ede86e..67b4f669a 100644
--- a/src/porespy/metrics/_funcs.py
+++ b/src/porespy/metrics/_funcs.py
@@ -12,6 +12,7 @@
 from porespy.filters import (
     local_thickness,
     pc_to_seq,
+    fill_blind_pores,
 )
 from porespy.tools import (
     Results,
@@ -201,14 +202,11 @@ def representative_elementary_volume(im, npoints=1000):
     return profile
 
 
-def porosity(im, mask=None):
+def porosity(im, mask=None, fill_hidden=False, fill_surface=False):
     r"""
     Calculates the porosity of an image assuming 1's are void space and 0's
     are solid phase.
 
-    All other values are ignored, so this can also return the relative
-    fraction of a phase of interest in multiphase images.
-
     Parameters
     ----------
     im : ndarray
@@ -221,6 +219,12 @@ def porosity(im, mask=None):
         entire array, like cylindrical cores.  Note that setting values in `im`
         to 2 will also exclude them from consideration so provides the same effect
         as `mask`, but providing a `mask` is usually much easier.
+    fill_hideen : bool (default = `False`)
+        A flag to indicate if hidden pores (not connected to any image boundary)
+        should be filled or not before computing the porosity.
+    fill_surface : bool (default = `False`)
+        A flag to indicate if blind pores connected only to one surface should be
+        filled or not befor computing the porosity.
 
     Returns
     -------
@@ -238,11 +242,6 @@ def porosity(im, mask=None):
     other values are ignored.  This is useful, for example, for images of
     cylindrical cores, where all voxels outside the core are labelled with 2.
 
-    Alternatively, images can be processed with ``find_disconnected_voxels``
-    to get an image of only blind pores.  This can then be added to the orignal
-    image such that blind pores have a value of 2, thus allowing the
-    calculation of accessible porosity, rather than overall porosity.
-
     Examples
     --------
     `Click here
@@ -250,8 +249,16 @@ def porosity(im, mask=None):
     to view online example.
 
     """
+    im = im.copy()
     if mask is not None:
         im = np.array(im, dtype=np.int64)*mask
+    if fill_hidden or fill_surface:
+        hidden_pores = im * ~fill_blind_pores(im, surface=False)
+        surface_pores = im * ~fill_blind_pores(im, surface=True) * ~hidden_pores
+        if fill_hidden:
+            im[hidden_pores] = 0
+        if fill_surface:
+            im[surface_pores] = 0
     Vp = np.sum(im == 1, dtype=np.int64)
     Vs = np.sum(im == 0, dtype=np.int64)
     e = Vp / (Vs + Vp)
@@ -1277,8 +1284,19 @@ def find_h(saturation, position=None, srange=[0.01, 0.99]):
 
     Returns
     -------
-    h : scalar
-        The height of the two-phase zone
+    A dataclass-like object with the following attributes:
+
+        =========== ================================================================
+        Attribute   Description
+        =========== ================================================================
+        `zmax`      The position where the saturation first exceeds `smax`
+        `zmin`      The position where the saturation first exceeds `smin`
+        `smax`      The value defining the start of the saturation profile
+        `smin`      The value defining the end of the saturation profile
+        `h`         The total distance in voxels between `zmax` and `zmin`
+        `valid`     A flag indicating whether the requested saturation difference
+                    (between `smin` and `smax`) was found.
+        =========== ================================================================
 
     See Also
     --------
@@ -1286,8 +1304,9 @@ def find_h(saturation, position=None, srange=[0.01, 0.99]):
 
     Notes
     -----
-    The ``satn_profile`` function can be used to obtain the ``saturation``
-    and ``position`` from an image.
+    The `satn_profile` function can be used to obtain the ``saturation``
+    and `position` from an image, such as a displacement map produced by
+    `drainage` or `imbibition`.
 
     Examples
     --------
diff --git a/src/porespy/metrics/_meshtools.py b/src/porespy/metrics/_meshtools.py
index 22bc5336d..3874c832c 100644
--- a/src/porespy/metrics/_meshtools.py
+++ b/src/porespy/metrics/_meshtools.py
@@ -27,7 +27,7 @@
 logger = logging.getLogger(__name__)
 
 
-def region_volumes(regions, mode='marching_cubes', voxel_size=(1, 1, 1)):
+def region_volumes(regions, method='marching_cubes', voxel_size=(1, 1, 1)):
     r"""
     Compute volume of each labelled region in an image
 
@@ -35,17 +35,18 @@ def region_volumes(regions, mode='marching_cubes', voxel_size=(1, 1, 1)):
     ----------
     regions : ndarray
         An image with labelled regions
-    mode : string
+    method : string
         Controls the method used. Options are:
 
-        'marching_cubes' (default)
-            Finds a mesh for each region using the marching cubes algorithm
-            from ``scikit-image``, then finds the volume of the mesh using the
-            ``trimesh`` package.
-
-        'voxel'
-            Calculates the region volume as the sum of voxels within each
-            region.
+        ================ ===========================================================
+        Option           Description
+        ================ ===========================================================
+        `marching_cubes` (default) Finds a mesh for each region using the marching
+                         cubes algorithm from ``scikit-image``, then finds the
+                         volume of the mesh using the ``trimesh`` package.
+        `voxel`          Calculates the region volume as the sum of voxels within
+                         each region.
+        ================ ===========================================================
 
     Returns
     -------
@@ -65,9 +66,9 @@ def region_volumes(regions, mode='marching_cubes', voxel_size=(1, 1, 1)):
     msg = "Computing region volumes".ljust(60)
     for i, s in enumerate(tqdm(slices, desc=msg, **settings.tqdm)):
         region = regions[s] == (i + 1)
-        if mode == 'marching_cubes':
+        if method == 'marching_cubes':
             vols[i] = mesh_volume(region, voxel_size=voxel_size)
-        elif mode.startswith('voxel'):
+        elif method.startswith('voxel'):
             vols[i] = region.sum(dtype=np.int64)
     return vols
 
diff --git a/src/porespy/networks/_getnet_orig.py b/src/porespy/networks/_getnet_orig.py
index 22c43888d..792759c60 100644
--- a/src/porespy/networks/_getnet_orig.py
+++ b/src/porespy/networks/_getnet_orig.py
@@ -257,7 +257,7 @@ def regions_to_network(
         logger.warning(msg)
         accuracy = 'standard'
     if (accuracy == 'high'):
-        net['pore.volume'] = region_volumes(regions=im, mode='marching_cubes')
+        net['pore.volume'] = region_volumes(regions=im, method='marching_cubes')
         areas = region_surface_areas(regions=im, voxel_size=voxel_size)
         net['pore.surface_area'] = areas
         interface_area = region_interface_areas(regions=im, areas=areas,
diff --git a/src/porespy/tools/_funcs.py b/src/porespy/tools/_funcs.py
index 52a01f051..76d3303a0 100644
--- a/src/porespy/tools/_funcs.py
+++ b/src/porespy/tools/_funcs.py
@@ -51,12 +51,37 @@
     'ps_rect',
     'ps_round',
     'subdivide',
+    'tilde',
     'unpad',
     'jit_extend_slice',
     'pad',
 ]
 
 
+def tilde(im):
+    r"""
+    Inverts an image.
+
+    This is to replace the `~` operator which is being deprecated. It also ensures
+    `im` is actually a `boolean` array since inverting `1` and `0` gives different
+    results than `True` and `False`.
+
+    Parameters
+    ----------
+    im : ndarray
+        The image whose boolean values are to be inverted. If image is not a boolean
+        it will be converted to one.
+
+    Returns
+    -------
+    inv : ndarray
+        The result of calling `np.logical_not` on `im`.
+    """
+    im = im.astype(bool)
+    inv = np.logical_not(im)
+    return inv
+
+
 def im_to_slabs(im, axis=0, span=50, step=None, mode='tile'):
     r"""
     Generates a list of slice objects which can be used to obtain slabs of an image


Test patch:
diff --git a/test/unit/test_metrics.py b/test/unit/test_metrics.py
index e11f27fd6..b0d206798 100644
--- a/test/unit/test_metrics.py
+++ b/test/unit/test_metrics.py
@@ -43,8 +43,15 @@ def setup_class(self):
         self.regions = np.array(io.imread(path))
 
     def test_porosity(self):
-        phi = ps.metrics.porosity(im=self.im2D)
-        assert np.allclose(phi, 0.66856)
+        im = ps.generators.blobs([300, 300], porosity=0.6, seed=1)
+        phi = ps.metrics.porosity(im)
+        assert np.allclose(phi, 0.6)
+        phi = ps.metrics.porosity(im, fill_surface=True, fill_hidden=False)
+        assert np.allclose(phi, 0.3963111111111111)
+        phi = ps.metrics.porosity(im, fill_surface=False, fill_hidden=True)
+        assert np.allclose(phi, 0.5996222222222222)
+        phi = ps.metrics.porosity(im, fill_surface=True, fill_hidden=True)
+        assert np.allclose(phi, 0.39593333333333336)
 
     def test_tpcf_fft_2d(self):
         tpcf_fft_1 = ps.metrics.two_point_correlation(self.im2D)
@@ -175,7 +182,7 @@ def test_region_interface_areas(self):
     def test_region_volumes(self):
         regions = self.regions[:50, :50, :50]
         vols_march = ps.metrics.region_volumes(regions=regions)
-        vols_vox = ps.metrics.region_volumes(regions=regions, mode='voxel')
+        vols_vox = ps.metrics.region_volumes(regions=regions, method='voxel')
         assert_allclose(vols_march[:5], [1498.85320453, 2597.90798652,
                                          2158.34548652, 1281.17978573, 1172.39853573])
         assert_allclose(vols_vox[:5], [1540., 2648., 2206., 1320., 1210.])
@@ -185,7 +192,7 @@ def test_region_volumes(self):
     def test_region_volumes_for_sphere(self):
         region = ball(10)
         vol_march = ps.metrics.region_volumes(regions=region)
-        vol_vox = ps.metrics.region_volumes(region, mode='voxel')
+        vol_vox = ps.metrics.region_volumes(region, method='voxel')
         assert_allclose(vol_march, 4102.28678846)
         assert_allclose(vol_vox, 4169.)