Merge pull request #77 from PMEAL/SNOW_bug

Bug resolved and added extra functionality.

porespy/network_extraction/__funcs__.py

@@ -68,17 +68,17 @@ def add_boundary_regions(regions=None, faces=['front', 'back', 'left',
             regions = sp.pad(regions, 2, 'edge')
 
             # Remove unselected faces
-            if 'top' not in faces:
-                regions = regions[:, 3:, :]
-            if 'bottom' not in faces:
-                regions = regions[:, :-3, :]
             if 'front' not in faces:
-                regions = regions[3:, :, :]
+                regions = regions[:, 3:, :]  # y
             if 'back' not in faces:
-                regions = regions[:-3, :, :]
+                regions = regions[:, :-3, :]
             if 'left' not in faces:
-                regions = regions[:, :, 3:]
+                regions = regions[3:, :, :]  # x
             if 'right' not in faces:
+                regions = regions[:-3, :, :]
+            if 'bottom' not in faces:
+                regions = regions[:, :, 3:]  # z
+            if 'top' not in faces:
                 regions = regions[:, :, :-3]
 
         elif regions.ndim == 2:
@@ -99,13 +99,13 @@ def add_boundary_regions(regions=None, faces=['front', 'back', 'left',
             regions = sp.pad(regions, 2, 'edge')
 
             # Remove unselected faces
-            if 'top' not in faces:
-                regions = regions[3:, :]
-            if 'bottom' not in faces:
-                regions = regions[:-3, :]
             if 'left' not in faces:
-                regions = regions[:, 3:]
+                regions = regions[3:, :]  # x
             if 'right' not in faces:
+                regions = regions[:-3, :]
+            if 'front' not in faces and 'bottom' not in faces:
+                regions = regions[:, 3:]  # y
+            if 'back' not in faces and 'top' not in faces:
                 regions = regions[:, :-3]
         else:
             print('add_boundary_regions works only on 2D and 3D images')

porespy/network_extraction/__snow__.py

@@ -1,15 +1,94 @@
 from porespy.network_extraction import regions_to_network, add_boundary_regions
 from porespy.filters import snow_partitioning
+from porespy.tools import make_contiguous
+import scipy as sp
 
 
 def snow(im, voxel_size=1, boundary_faces=['top', 'bottom', 'left',
                                            'right', 'front', 'back']):
     r"""
+    Analyzes an image that has been partitioned into void and solid regions
+    and extracts the void and solid phase geometry as well as network
+    connectivity.
+
+    Parameters
+    ----------
+    im : ND-array
+        Binary image in the Boolean form with True’s as void phase and False’s
+        as solid phase.
+
+    voxel_size : scalar
+        The resolution of the image, expressed as the length of one side of a
+        voxel, so the volume of a voxel would be **voxel_size**-cubed.  The
+        default is 1, which is useful when overlaying the PNM on the original
+        image since the scale of the image is alway 1 unit lenth per voxel.
+
+    boundary_faces : list of strings
+        Boundary faces labels are provided to assign hypothetical boundary
+        nodes having zero resistance to transport process. For cubical
+        geometry, the user can choose ‘left’, ‘right’, ‘top’, ‘bottom’,
+        ‘front’ and ‘back’ face labels to assign boundary nodes. If no label is
+        assigned then all six faces will be selected as boundary nodes
+        automatically which can be trimmed later on based on user requirements.
+
+    Returns
+    -------
+    A dictionary containing the void phase size data, as well as the network
+    topological information.  The dictionary names use the OpenPNM
+    convention (i.e. 'pore.coords', 'throat.conns') so it may be converted
+    directly to an OpenPNM network object using the ``update`` command.
     """
+
     regions = snow_partitioning(im=im, return_all=True)
     im = regions.im
     dt = regions.dt
     regions = regions.regions
+    b_num = sp.amax(regions)
     regions = add_boundary_regions(regions=regions, faces=boundary_faces)
-    net = regions_to_network(im=regions*im, dt=dt, voxel_size=voxel_size)
+    f = boundary_faces
+    if f is not None:
+        if im.ndim == 2:
+            faces = [(int('left' in f)*3, int('right' in f)*3),
+                     (int(('front') in f)*3 or int(('bottom') in f)*3,
+                      int(('back') in f)*3 or int(('top') in f)*3)]
+
+        if im.ndim == 3:
+            faces = [(int('left' in f)*3, int('right' in f)*3),
+                     (int('front' in f)*3, int('back' in f)*3),
+                     (int('top' in f)*3, int('bottom' in f)*3)]
+        dt = sp.pad(dt, pad_width=faces, mode='edge')
+        im = sp.pad(im, pad_width=faces, mode='edge')
+    else:
+        dt = dt
+    regions = regions*im
+    regions = make_contiguous(regions)
+    net = regions_to_network(im=regions, dt=dt, voxel_size=voxel_size)
+    boundary_labels = net['pore.label'] > b_num
+    loc1 = net['throat.conns'][:, 0] < b_num
+    loc2 = net['throat.conns'][:, 1] >= b_num
+    pore_labels = net['pore.label'] <= b_num
+    loc3 = net['throat.conns'][:, 0] < b_num
+    loc4 = net['throat.conns'][:, 1] < b_num
+    net['pore.boundary'] = boundary_labels
+    net['throat.boundary'] = loc1 * loc2
+    net['pore.internal'] = pore_labels
+    net['throat.internal'] = loc3 * loc4
+    # -------------------------------------------------------------------------
+    # label boundary pore faces
+    if f is not None:
+        coords = net['pore.coords']
+        condition = coords[net['pore.internal']]
+        dic = {'left': 0, 'right': 0, 'front': 1, 'back': 1,
+               'top': 2, 'bottom': 2}
+        if all(coords[:, 2] == 0):
+            dic['top'] = 1
+            dic['bottom'] = 1
+        print(dic)
+        for i in f:
+            if i in ['left', 'front', 'bottom']:
+                net['pore.{}'.format(i)] = (coords[:, dic[i]] <
+                                            min(condition[:, dic[i]]))
+            elif i in ['right', 'back', 'top']:
+                net['pore.{}'.format(i)] = (coords[:, dic[i]] >
+                                            max(condition[:, dic[i]]))
     return net

porespy/network_extraction/__snow_dual__.py

@@ -69,10 +69,14 @@ def snow_dual(im, voxel_size=1, boundary_faces=['top', 'bottom', 'left',
     # Padding distance transform to extract geometrical properties
     f = boundary_faces
     if f is not None:
-        faces = [(int('top' in f)*3, int('bottom' in f)*3),
-                 (int('left' in f)*3, int('right' in f)*3)]
+        if im.ndim == 2:
+            faces = [(int('left' in f)*3, int('right' in f)*3),
+                     (int(('front') in f)*3 or int(('bottom') in f)*3,
+                      int(('back') in f)*3 or int(('top') in f)*3)]
         if im.ndim == 3:
-            faces = [(int('front' in f)*3, int('back' in f)*3)] + faces
+            faces = [(int('left' in f)*3, int('right' in f)*3),
+                     (int('front' in f)*3, int('back' in f)*3),
+                     (int('top' in f)*3, int('bottom' in f)*3)]
         dt = sp.pad(dt, pad_width=faces, mode='edge')
     else:
         dt = dt

test/unit/test_network_extraction.py

@@ -38,6 +38,14 @@ def test_extract_pore_network_3d(self):
                 found_nans = True
         assert found_nans is False
 
+    def test_snow(self):
+        net = ps.network_extraction.snow(self.im3d)
+        found_nans = False
+        for key in net.keys():
+            if np.any(np.isnan(net[key])):
+                found_nans = True
+        assert found_nans is False
+
 #    def test_snow_dual_2d(self):
 #        net = ps.network_extraction.snow_dual(self.im)
 #        found_nans = False
@@ -59,27 +67,30 @@ def test_add_bounadary_regions_2D(self):
         regions = ps.filters.snow_partitioning(im)
         f = ['left', 'right']
         bd = ps.network_extraction.add_boundary_regions(regions, faces=f)
-        assert bd.shape[1] > regions.shape[1]
-        f = ['bottom', 'top']
-        bd = ps.network_extraction.add_boundary_regions(regions, faces=f)
         assert bd.shape[0] > regions.shape[0]
-        f = ['bottom', 'top', 'left', 'right']
+        f = ['bottom', 'top']
         bd = ps.network_extraction.add_boundary_regions(regions, faces=f)
-        assert bd.shape[0] > regions.shape[0]
         assert bd.shape[1] > regions.shape[1]
         f = ['front', 'back']
         bd = ps.network_extraction.add_boundary_regions(regions, faces=f)
-        assert bd.shape == regions.shape
+        assert bd.shape[1] > regions.shape[1]
+        f = ['bottom', 'top', 'left', 'right', 'front', 'back']
+        bd = ps.network_extraction.add_boundary_regions(regions, faces=f)
+        assert bd.shape[0] > regions.shape[0]
+        assert bd.shape[1] > regions.shape[1]
 
     def test_add_bounadary_regions_3D(self):
         im = self.im3d
         regions = ps.filters.snow_partitioning(im)
         f = ['left', 'right']
         bd = ps.network_extraction.add_boundary_regions(regions, faces=f)
-        assert bd.shape[2] > regions.shape[2]
-        f = ['bottom', 'top']
+        assert bd.shape[0] > regions.shape[0]
+        f = ['front', 'back']
         bd = ps.network_extraction.add_boundary_regions(regions, faces=f)
         assert bd.shape[1] > regions.shape[1]
+        f = ['bottom', 'top']
+        bd = ps.network_extraction.add_boundary_regions(regions, faces=f)
+        assert bd.shape[2] > regions.shape[2]
         f = ['bottom', 'top', 'left', 'right', 'front', 'back']
         bd = ps.network_extraction.add_boundary_regions(regions, faces=f)
         assert bd.shape[0] > regions.shape[0]