SNOW network extraction

The SNOW algorithm, published in Physical Review E, uses a marker-based watershed segmentation algorithm to partition an image into regions belonging to each pore. The main contribution of the SNOW algorithm is to find a suitable set of initial markers in the image so that the watershed is not over-segmented. SNOW is an acronym for Sub-Network of an Over-segmented Watershed. This code works on both 2D and 3D images. In this example a 2D image will be segmented using the predefined snow function in PoreSpy.

Start by importing the necessary packages:

import numpy as np
import porespy as ps
import openpnm as op
import matplotlib.pyplot as plt
ps.visualization.set_mpl_style()
np.random.seed(10)
[01:07:01] ERROR    PARDISO solver not installed, run `pip install pypardiso`. Otherwise,          _workspace.py:56
                    simulations will be slow. Apple M chips not supported.                                         

Generate an artificial 2D image for illustration purposes:

im = ps.generators.blobs(shape=[400, 400], porosity=0.6, blobiness=2)
fig, ax = plt.subplots(figsize=(4, 4))
ax.imshow(im);

SNOW is composed of a series of filters, but PoreSpy has a single function that applies all the necessary steps:

snow_output = ps.networks.snow2(im, voxel_size=1)

The snow function returns an object that has a network attribute. This is a dictionary that is suitable for loading into OpenPNM. The best way to get this into OpenPNM is to use the PoreSpy IO class. This splits the data into a network and a geometry:

pn = op.io.network_from_porespy(snow_output.network)

As can be seen by printing the network it contains quite a lot of geometric information:

print(pn)
══════════════════════════════════════════════════════════════════════════════
net : <openpnm.network.Network at 0x7f12529df980>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Properties                                                   Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  throat.conns                                                    283 / 283
  3  pore.coords                                                     262 / 262
  4  pore.region_label                                               262 / 262
  5  pore.phase                                                      262 / 262
  6  throat.phases                                                   283 / 283
  7  pore.region_volume                                              262 / 262
  8  pore.equivalent_diameter                                        262 / 262
  9  pore.local_peak                                                 262 / 262
 10  pore.global_peak                                                262 / 262
 11  pore.geometric_centroid                                         262 / 262
 12  throat.global_peak                                              283 / 283
 13  pore.inscribed_diameter                                         262 / 262
 14  pore.extended_diameter                                          262 / 262
 15  throat.inscribed_diameter                                       283 / 283
 16  throat.total_length                                             283 / 283
 17  throat.direct_length                                            283 / 283
 18  throat.perimeter                                                283 / 283
 19  pore.volume                                                     262 / 262
 20  pore.surface_area                                               262 / 262
 21  throat.cross_sectional_area                                     283 / 283
 22  throat.equivalent_diameter                                      283 / 283
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Labels                                                 Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  2  pore.all                                                              262
  3  throat.all                                                            283
  4  pore.boundary                                                          46
  5  pore.xmin                                                              12
  6  pore.xmax                                                              12
  7  pore.ymin                                                              13
  8  pore.ymax                                                               9
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

You can also overlay the network on the image natively in porespy. Note that you need to transpose the image using im.T, since imshow uses matrix representation, e.g. a (10, 20)-shaped array is shown as 10 pixels in the y-axis, and 20 pixels in the x-axis.

fig, ax = plt.subplots(figsize=[5, 5])
ax.imshow(im.T, cmap=plt.cm.bone);
op.visualization.plot_coordinates(ax=fig,
                                  network=pn,
                                  size_by=pn["pore.inscribed_diameter"],
                                  color_by=pn["pore.inscribed_diameter"],
                                  markersize=200)
op.visualization.plot_connections(network=pn, ax=fig)
ax.axis("off");