Repository: pmeal/porespy
Issue: metrics region_volumes gives the same or negative values for different regions

The problem is when choosing 'voxel' method, the returned volumes of regions is the same for all pore regions (this can be fixed by changing:
https://github.com/PMEAL/porespy/blob/76e0c3132e69d07eeb115d872302ad12203a287a/porespy/metrics/_meshtools.py#L49
to:
```python
vols[i]=region.sum()
```
Another problem is when choosing 'marching cube' method, the returned volumes of regions are all negative. I will be checking the root cause and updating.

example of the problem:
```python
import matplotlib.pyplot as plt
import numpy as np
import porespy as ps
np.random.seed(10)
im = ps.generators.blobs(shape=[50,50,50])
snow = ps.filters.snow_partitioning(im)
regions = snow.regions
vols = ps.metrics.region_volumes(regions=regions)
print(vols[:5])
#[ -17.74722219 -103.5931484   -18.03024671 -231.43682392 -689.57165046]
vols = ps.metrics.region_volumes(regions=regions, mode='voxel')
print(vols[:5])
#[7474405. 7474405. 7474405. 7474405. 7474405.]
```


Changed files:
test/integration/test_drainage.py, test/integration/test_variable_Bo_study.py, test/unit/test_io.py, porespy/metrics/_funcs.py, test/integration/test_drainage_from_top.py, test/integration/test_inverse_Bo_study.py, test/integration/test_snow_example_script.py, porespy/metrics/_meshtools.py, porespy/visualization/_funcs.py, test/integration/test_ibip_example_script.py, test/unit/test_metrics.py, test/unit/test_visualization.py

diff --git a/porespy/metrics/_funcs.py b/porespy/metrics/_funcs.py
index ee03abf3..5f7c69ad 100644
--- a/porespy/metrics/_funcs.py
+++ b/porespy/metrics/_funcs.py
@@ -787,6 +787,7 @@ def phase_fraction(im, normed=True):
         results[label] = np.sum(im == label) * (1 / im.size if normed else 1)
     return results
 
+
 @deprecated("This function is deprecated, use pc_curve instead")
 def pc_curve_from_ibip(*args, **kwargs):
     r"""
diff --git a/porespy/metrics/_meshtools.py b/porespy/metrics/_meshtools.py
index 52c68cab..b5a4a4a5 100644
--- a/porespy/metrics/_meshtools.py
+++ b/porespy/metrics/_meshtools.py
@@ -46,7 +46,7 @@ def region_volumes(regions, mode='marching_cubes'):
         if mode == 'marching_cubes':
             vols[i] = mesh_volume(region)
         elif mode.startswith('voxel'):
-            vols[i] = regions.sum()
+            vols[i] = region.sum()
     return vols
 
 
@@ -70,7 +70,7 @@ def mesh_volume(region):
     mc = mesh_region(region > 0)
     m = Trimesh(vertices=mc.verts, faces=mc.faces, vertex_normals=mc.norm)
     if m.is_watertight:
-        vol = m.volume
+        vol = np.abs(m.volume)
     else:
         vol = np.nan
     return vol
diff --git a/porespy/visualization/_funcs.py b/porespy/visualization/_funcs.py
index 6d740375..274bac11 100644
--- a/porespy/visualization/_funcs.py
+++ b/porespy/visualization/_funcs.py
@@ -176,7 +176,7 @@ def prep_for_imshow(im, mask=None, axis=0, slice=None):
     ----------
     im : ndimage
         The image to show. If ``im`` includes ``+inf`` or ``-inf`` values,
-        they are converted to 1 above or below the maximum and minimum finite
+        they are converted to 1 above or below the minimum and maximum finite
         values in ``im``, respectively.
     mask : ndimage, optional
         An image of the porous material with ``True`` indicating voxels of
@@ -205,8 +205,10 @@ def prep_for_imshow(im, mask=None, axis=0, slice=None):
         'X'             The adjusted image with ``+inf`` replaced by
                         ``vmax + 1``, and all solid voxels replacd by
                         ``np.nan`` to show as white in ``imshow``
-        'vmax'          The maximum of ``values`` not including ``+inf`` or values in
-        'vmin'          The minimum of ``values`` not including ``-inf``
+        'vmax'          The maximum of ``values`` not including ``+inf`` or
+                        values in ``False`` voxels in ``mask``.
+        'vmin'          The minimum of ``values`` not including ``-inf`` or
+                        values in ``False`` voxels in ``mask``.
         'interpolation' Set to 'none' to avoid artifacts in ``imshow``
         'origin'        Set to 'lower' to put (0, 0) on the bottom-left corner
         =============== =======================================================
diff --git a/test/integration/test_drainage.py b/test/integration/test_drainage.py
index 2be26190..d80ed86c 100644
--- a/test/integration/test_drainage.py
+++ b/test/integration/test_drainage.py
@@ -1,16 +1,17 @@
 import numpy as np
 import porespy as ps
 import matplotlib.pyplot as plt
-np.random.seed(6)
 
 
 def test_drainage():
+    np.random.seed(6)
     im = ps.generators.blobs(shape=[500, 500], porosity=0.7, blobiness=1.5)
     inlets = np.zeros_like(im)
     inlets[0, :] = True
     outlets = np.zeros_like(im)
     outlets[-1, :] = True
-    im = ps.filters.trim_nonpercolating_paths(im=im, inlets=inlets, outlets=outlets)
+    im = ps.filters.trim_nonpercolating_paths(im=im, inlets=inlets,
+                                              outlets=outlets)
     pc = None
     lt = ps.filters.local_thickness(im)
     residual = lt > 25
@@ -22,7 +23,8 @@ def test_drainage():
     g = 0
     bg = 'grey'
 
-    drn1 = ps.simulations.drainage(im=im, voxel_size=voxel_size, inlets=inlets, g=g)
+    drn1 = ps.simulations.drainage(im=im, voxel_size=voxel_size,
+                                   inlets=inlets, g=g)
     drn2 = ps.simulations.drainage(im=im,
                                    voxel_size=voxel_size,
                                    inlets=inlets,
@@ -70,38 +72,38 @@ def test_drainage():
         ax = plt.axes()
         ax.set_facecolor(bg)
         plt.step(np.log10(drn1.pc), drn1.snwp, 'b-o', where='post',
-                label="No trapping, no residual")
+                 label="No trapping, no residual")
         plt.step(np.log10(drn2.pc), drn2.snwp, 'r--o', where='post',
-                label="With trapping, no residual")
+                 label="With trapping, no residual")
         plt.step(np.log10(drn3.pc), drn3.snwp, 'g--o', where='post',
-                label="No trapping, with residual")
+                 label="No trapping, with residual")
         plt.step(np.log10(drn4.pc), drn4.snwp, 'm--o', where='post',
-                label="With trapping, with residual")
+                 label="With trapping, with residual")
         plt.legend()
 
     # %% Now repeat with some gravity
     g = 9.81
 
     drn1 = ps.simulations.drainage(im=im,
-                                voxel_size=voxel_size,
-                                inlets=inlets,
-                                g=g)
+                                   voxel_size=voxel_size,
+                                   inlets=inlets,
+                                   g=g)
     drn2 = ps.simulations.drainage(im=im,
-                                voxel_size=voxel_size,
-                                inlets=inlets,
-                                outlets=outlets,
-                                g=g)
+                                   voxel_size=voxel_size,
+                                   inlets=inlets,
+                                   outlets=outlets,
+                                   g=g)
     drn3 = ps.simulations.drainage(im=im,
-                                voxel_size=voxel_size,
-                                inlets=inlets,
-                                residual=residual,
-                                g=g)
+                                   voxel_size=voxel_size,
+                                   inlets=inlets,
+                                   residual=residual,
+                                   g=g)
     drn4 = ps.simulations.drainage(im=im,
-                                voxel_size=voxel_size,
-                                inlets=inlets,
-                                outlets=outlets,
-                                residual=residual,
-                                g=g)
+                                   voxel_size=voxel_size,
+                                   inlets=inlets,
+                                   outlets=outlets,
+                                   residual=residual,
+                                   g=g)
 
     # Ensure initial saturations correspond to amount of residual present
     assert drn1.snwp[0] == 0
diff --git a/test/integration/test_drainage_from_top.py b/test/integration/test_drainage_from_top.py
index 89960ad3..fd597d83 100644
--- a/test/integration/test_drainage_from_top.py
+++ b/test/integration/test_drainage_from_top.py
@@ -1,16 +1,17 @@
 import numpy as np
 import porespy as ps
 import matplotlib.pyplot as plt
-np.random.seed(6)
 
 
 def test_drainage_from_top():
+    np.random.seed(6)
     im = ps.generators.blobs(shape=[300, 300], porosity=0.7, blobiness=1.5)
     inlets = np.zeros_like(im)
     inlets[-1, :] = True
     outlets = np.zeros_like(im)
     outlets[0, :] = True
-    im = ps.filters.trim_nonpercolating_paths(im=im, inlets=inlets, outlets=outlets)
+    im = ps.filters.trim_nonpercolating_paths(im=im, inlets=inlets,
+                                              outlets=outlets)
     pc = None
     lt = ps.filters.local_thickness(im)
     residual = lt > 25
@@ -90,11 +91,11 @@ def test_drainage_from_top():
         ax = plt.axes()
         ax.set_facecolor(bg)
         plt.step(np.log10(drn1.pc), drn1.snwp, 'b-o', where='post',
-                label="No trapping, no residual")
+                 label="No trapping, no residual")
         plt.step(np.log10(drn2.pc), drn2.snwp, 'r--o', where='post',
-                label="With trapping, no residual")
+                 label="With trapping, no residual")
         plt.step(np.log10(drn3.pc), drn3.snwp, 'g--o', where='post',
-                label="No trapping, with residual")
+                 label="No trapping, with residual")
         plt.step(np.log10(drn4.pc), drn4.snwp, 'm--o', where='post',
-                label="With trapping, with residual")
+                 label="With trapping, with residual")
         plt.legend()
diff --git a/test/integration/test_ibip_example_script.py b/test/integration/test_ibip_example_script.py
index 76ccf951..d7fb3864 100644
--- a/test/integration/test_ibip_example_script.py
+++ b/test/integration/test_ibip_example_script.py
@@ -25,5 +25,5 @@ def test_ibip():
         cmap = copy(plt.cm.viridis)
         cmap.set_under(color='black')
         plt.imshow(inv_seq, cmap=cmap, vmin=1e-3,
-                interpolation='none', origin='lower')
+                   interpolation='none', origin='lower')
         mov = ps.visualization.satn_to_movie(im, inv_satn)
diff --git a/test/integration/test_inverse_Bo_study.py b/test/integration/test_inverse_Bo_study.py
index 33e33475..40154a5f 100644
--- a/test/integration/test_inverse_Bo_study.py
+++ b/test/integration/test_inverse_Bo_study.py
@@ -26,22 +26,26 @@ def test_inverse_Bo_study():
         Bo = 1/inv_Bo[i]
         delta_rho = Bo*sigma/(g*a**2)  # delta_rho is found given the Bo
         sim1[i] = ps.simulations.drainage(im=im,
-                                        voxel_size=vx,
-                                        inlets=inlets,
-                                        delta_rho=delta_rho,
-                                        sigma=sigma,
-                                        g=g,
-                                        bins=25)
+                                          voxel_size=vx,
+                                          inlets=inlets,
+                                          delta_rho=delta_rho,
+                                          sigma=sigma,
+                                          g=g,
+                                          bins=25)
 
     # %%  Process data to make 1/Bo vs H plot
     data = []
     smin, smax = 0.1, 0.90
     for h in range(len(inv_Bo)):
         for s in np.arange(0.2, 1.0, 0.1):
-            prof = ps.metrics.satn_profile(satn=sim1[h].im_satn, s=s, span=1, mode='slide')
+            prof = ps.metrics.satn_profile(satn=sim1[h].im_satn, s=s, span=1,
+                                           mode='slide')
             if 0:
-                plt.plot(prof.position, prof.saturation)
-            H = ps.metrics.find_h(prof.saturation, prof.position, srange=[smin, smax])
+                plt.plot(prof.position,
+                         prof.saturation)
+            H = ps.metrics.find_h(prof.saturation,
+                                  prof.position,
+                                  srange=[smin, smax])
             data.append((inv_Bo[h], s, H.zmin, H.zmax, H.h))
     df = pd.DataFrame(data, columns=['1/Bo', 's', 'zmin', 'zmax', 'h'])
 
diff --git a/test/integration/test_snow_example_script.py b/test/integration/test_snow_example_script.py
index defcdf87..5f08dec4 100644
--- a/test/integration/test_snow_example_script.py
+++ b/test/integration/test_snow_example_script.py
@@ -2,11 +2,11 @@ import numpy as np
 import porespy as ps
 import matplotlib.pyplot as plt
 import openpnm as op
-np.random.seed(0)
 
 
 def test_snow_example_script():
     plot = False
+    np.random.seed(0)
 
     im1 = ps.generators.blobs(shape=[600, 400], porosity=None, blobiness=1) < 0.4
     im2 = ps.generators.blobs(shape=[600, 400], porosity=None, blobiness=1) < 0.7
@@ -14,22 +14,20 @@ def test_snow_example_script():
     # phases = phases > 0
 
     snow_n = ps.networks.snow2(phases,
-                            phase_alias={1: 'solid', 2: 'void'},
-                            boundary_width=5,
-                            accuracy='high',
-                            parallelization=None)
+                               phase_alias={1: 'solid', 2: 'void'},
+                               boundary_width=5,
+                               parallelization=None)
 
-    assert snow_n.regions.max() == 211
+    assert snow_n.regions.max() == 210
     # Remove all but 1 pixel-width of boundary regions
     temp = ps.tools.extract_subsection(im=snow_n.regions,
-                                    shape=np.array(snow_n.regions.shape)-8)
-    assert temp.max() == 211
+                                       shape=np.array(snow_n.regions.shape)-8)
+    assert temp.max() == 210
     # Remove complete boundary region
     temp = ps.tools.extract_subsection(im=snow_n.regions,
-                                    shape=np.array(snow_n.regions.shape)-10)
-    assert temp.max() == 164
+                                       shape=np.array(snow_n.regions.shape)-10)
+    assert temp.max() == 163
 
-    # %% Plot the final extraction overlayed with snow segmentation
     if plot:
         fig, ax = plt.subplots(1, 1)
         ax.imshow(ps.tools.randomize_colors(snow_n.regions.T))
diff --git a/test/integration/test_variable_Bo_study.py b/test/integration/test_variable_Bo_study.py
index 03f855f5..22e2e227 100644
--- a/test/integration/test_variable_Bo_study.py
+++ b/test/integration/test_variable_Bo_study.py
@@ -74,7 +74,8 @@ def test_variable_Bo_study():
         vmax = np.amax(temp)
         temp[temp == 0] = vmin - 1
         temp[im == 0] = vmax + 1
-        ax.imshow(temp, vmax=vmax, vmin=vmin, cmap=cmap, origin='lower', interpolation='none')
+        ax.imshow(temp, vmax=vmax, vmin=vmin, cmap=cmap, origin='lower',
+                  interpolation='none')
         ax.axis('off')
         plt.colorbar(ax.imshow(temp, vmax=vmax, vmin=vmin, cmap=cmap, origin='lower'))
 
diff --git a/test/unit/test_io.py b/test/unit/test_io.py
index bdb65425..7da40404 100644
--- a/test/unit/test_io.py
+++ b/test/unit/test_io.py
@@ -47,36 +47,6 @@ class ExportTest():
         assert a < b
         os.remove('dictvtk.vti')
 
-    def test_openpnm_to_im(self):
-        net = op.network.Cubic(shape=[5, 5, 5])
-        spacing = op.topotools.get_spacing(net)
-        shape = op.topotools.get_shape(net)
-        # FIXME: get rid of try/except once openpnm v3 is out
-        try:
-            geom = op.geometry.SpheresAndCylinders(network=net,
-                                                   pores=net.Ps,
-                                                   throats=net.Ts)
-        except AttributeError:
-            geom = op.geometry.StickAndBall(network=net,
-                                            pores=net.Ps,
-                                            throats=net.Ts)
-
-        geom.add_model(propname="pore.volume",
-                       model=op.models.geometry.pore_volume.cube)
-        geom.add_model(propname="throat.volume",
-                       model=op.models.geometry.throat_volume.cylinder)
-        geom.regenerate_models()
-
-        im = ps.io.openpnm_to_im(network=net, pore_shape="cube",
-                                 throat_shape="cylinder", rtol=0.01)
-        porosity_actual = im.astype(bool).sum() / np.prod(im.shape)
-
-        volume_void = net["pore.volume"].sum() + net["throat.volume"].sum()
-        volume_total = np.prod(spacing * shape)
-        porosity_desired = volume_void / volume_total
-
-        assert_allclose(actual=porosity_actual, desired=porosity_desired, rtol=0.1)
-
     def test_to_stl(self):
         im = ps.generators.blobs(shape=[50, 50, 50])
         ps.io.to_stl(im, filename="im2stl")
diff --git a/test/unit/test_metrics.py b/test/unit/test_metrics.py
index 10598b54..60301cf0 100644
--- a/test/unit/test_metrics.py
+++ b/test/unit/test_metrics.py
@@ -6,6 +6,7 @@ from skimage import io
 from edt import edt
 from pathlib import Path
 import scipy.ndimage as spim
+from skimage.morphology import ball
 from numpy.testing import assert_allclose
 ps.settings.tqdm['disable'] = True
 
@@ -151,6 +152,23 @@ class MetricsTest():
         # assert np.all(ia.conns[0] == [2, 19])
         # assert np.around(ia.area[0], decimals=2) == 3.59
 
+    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')
+        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.])
+        assert_allclose(np.mean(vols_march), 1907.8062788852674)
+        assert_allclose(np.mean(vols_vox), 1952.125)
+
+    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')
+        assert_allclose(vol_march, 4102.28678846)
+        assert_allclose(vol_vox, 4169.)
+
     def test_phase_fraction(self):
         im = np.reshape(np.random.randint(0, 10, 1000), [10, 10, 10])
         labels = np.unique(im, return_counts=True)[1]
diff --git a/test/unit/test_visualization.py b/test/unit/test_visualization.py
index b31d1588..b1549436 100644
--- a/test/unit/test_visualization.py
+++ b/test/unit/test_visualization.py
@@ -85,7 +85,7 @@ class VisualizationTest():
     def test_prep_for_imshow_3D(self):
         a = ps.visualization.prep_for_imshow(self.lt, self.im)
         assert a['X'].shape == (51, 51)
-        assert a['vmin'] == 0
+        assert a['vmin'] == 1.0
         b = ps.visualization.prep_for_imshow(self.lt, self.im, axis=None)
         assert b['X'].shape == (51, 51, 51)