`snow2`#

The SNOW algorithm uses a standard marker-based watershed segementation function to identify pores. The key to the SNOW algorithm is in how the markers are identified. Finding peaks in the distance transform which truly correspond to pore centers requires several careful steps. Each of these steps can be done manually, but the snow2 function was created to combine them all in a single function for ease of use.

In addition to the original SNOW algorithm, snow2 also incorporates a few new developements such as multiphase extractions and parallelized processing. The new snow2 function also handles boundaries better, and returns a different set of geometrical properties. For instance it no longer returns ‘pore.diameter’, which was an opinionated choice that was not appropriate for every case. Instead it returns several different diameters, like ‘pore.equivalent_diameter’ and ‘pore.inscribed_diameter’, and the user must decide which of these they wish to declare as the offical ‘pore.diameter’.

The snow2 function has quite a few arguments, which will be explored and explained in this example.

import porespy as ps
import openpnm as op
import matplotlib.pyplot as plt
import numpy as np
from edt import edt
import scipy.ndimage as spim

/opt/hostedtoolcache/Python/3.8.16/x64/lib/python3.8/site-packages/openpnm/algorithms/_invasion_percolation.py:358: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def _find_trapped_pores(inv_seq, indices, indptr, outlets):  # pragma: no cover

im = ~ps.generators.rsa([400, 400], r=15, clearance=5)
ps.imshow(im);

/tmp/ipykernel_3816/1531839659.py:1: DeprecationWarning: Call to deprecated function (or staticmethod) rsa. (This function will be renamed random_spheres in a future version)
  im = ~ps.generators.rsa([400, 400], r=15, clearance=5)

../../../_images/12c71a6221af6f51f3b91393a2524e28cdfb6b63650f22c050bba8c3742016e4.png

`phases`#

One of the main features of snow2 is the ability to extract networks from multiple phases simultaneously, and have them be interconnected. This is quite useful for modeling electrochemial devices for instance, since transport occurs in the solid and void phases simultaneously, and the processes are coupled. For a powerful example, see Khan et al. where the transient discharage of a Li-ion battery was modelled using PNMs. Li-ion batteries contain 3 distinct phases: liquid electrolyte, solid active material, and carbonaceous binder.

Note that snow2 does not need to operate on mulitphase phases, only that it can if needed. It accomplishes this flexibility by ignoring all 0’s in the input images, and treating all voxels with values of 1, 2, 3, etc. as distinct phases. Therefore supplying an image of 0’s and 1’s will result in a standard single phase extraction, while supplying an image of 1’s and 2’s will create a dual phase network.

A few additional notes:

It should also be mentioned that the image can contain 0’s as well as values (ie. 1, 2, 3), so can have multiple phases and have some regions ignored
The voxels values do not have be integers or sequential. Each unique numerical value in the image is treated as phase.

Let’s explore this:

net1 = ps.networks.snow2(phases=im, boundary_width=0)
net2 = ps.networks.snow2(phases=im.astype(int) + 1, boundary_width=0)

The “object” returned by the snow2 function is a dataclass-like object, meaning it has data attached to it as attributes. These attributes are the segmented image called regions, the OpenPNM compatible network called network, and a copy of the input image called phases which will be padded to be the same size as regions if any boundaries were added.

fig, ax = plt.subplots(1, 2, figsize=[12, 6])
ax[0].imshow(net1.regions, origin='lower', interpolation='none')
ax[1].imshow(net2.regions, origin='lower', interpolation='none');

../../../_images/58c94a12c7cc202e5dfc43e83b7a27b9bcca7321e062093f421c1ee14802b775.png

On the left we can see that the spherical particles are all purple, which indicates 0, because they have not been labelled by the function. On the right the particles each have a unique greyscale value because they have been labelled as distinct regions.

The returned objects (net1 and net2) also contain the pore networks, in OpenPNM format. We can overlay these networks with the above images. Let’s start with the single phase network:

pn1 = op.io.network_from_porespy(net1.network)

ax = op.visualization.plot_connections(pn1, c='k')
ax = op.visualization.plot_coordinates(pn1, markersize=50, ax=ax)
fig, ax = plt.gcf(), plt.gca()
ax.imshow(net1.regions.T)
ax.axis(False);

../../../_images/b2d20eadd0a2b7da5b798b85c446bea5dcc453518b7e9b0c4ce00a3a14e3c490.png

The image needed to be rotated so that the x and y directions lined up with the OpenPNM network. The agreement between the image and network is pretty good. Clearly the the red pore centers lie within each labelled regions and the black throat connections follow the topology quite well.

Now let’s do the same for the multiphase network:

pn2 = op.io.network_from_porespy(net2.network)

ax = op.visualization.plot_connections(pn2, c='k')
ax = op.visualization.plot_coordinates(pn2, markersize=50, ax=ax)
fig, ax = plt.gcf(), plt.gca()
ax.imshow(net2.regions.T)
ax.axis(False);

../../../_images/18b12d21437596131340efb9e3f6b49bcb40029c5dd130f097e355d662b8ce76.png

The above shows the additional network comprised of the solids (although they don’t actually connect in this 2D example) as well as the void-solid connections. We can improve this plot by coloring the different types of pores and throats differently. If we print pn2 we can see that snow2 has created some labels for us, by assuming phase names of “phase1” and “phase2”.

print(pn2)

══════════════════════════════════════════════════════════════════════════════
net : <openpnm.network.Network at 0x7f11f9188f90>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Properties                                                   Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
throat.conns                                                    717 / 717
pore.coords                                                     251 / 251
pore.region_label                                               251 / 251
pore.phase                                                      251 / 251
throat.phases                                                   717 / 717
pore.region_volume                                              251 / 251
pore.equivalent_diameter                                        251 / 251
pore.local_peak                                                 251 / 251
pore.global_peak                                                251 / 251
pore.geometric_centroid                                         251 / 251
throat.global_peak                                              717 / 717
pore.inscribed_diameter                                         251 / 251
pore.extended_diameter                                          251 / 251
throat.inscribed_diameter                                       717 / 717
throat.total_length                                             717 / 717
throat.direct_length                                            717 / 717
throat.perimeter                                                717 / 717
pore.volume                                                     251 / 251
pore.surface_area                                               251 / 251
throat.cross_sectional_area                                     717 / 717
throat.equivalent_diameter                                      717 / 717
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Labels                                                 Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
pore.all                                                              251
throat.all                                                            717
pore.phase1                                                           102
throat.phase1_phase2                                                  467
pore.phase2                                                           149
throat.phase2_phase1                                                  467
throat.phase2_phase2                                                  250
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

ax = op.visualization.plot_connections(pn2, throats=pn2['throat.phase2_phase2'], c='k')
ax = op.visualization.plot_connections(pn2, throats=pn2['throat.phase2_phase1'], c='grey', ax=ax)
ax = op.visualization.plot_coordinates(pn2, pores=pn2['pore.phase2'], c='r', markersize=50, ax=ax)
ax = op.visualization.plot_coordinates(pn2, pores=pn2['pore.phase1'], c='g', markersize=50, ax=ax)
fig, ax = plt.gcf(), plt.gca()
ax.imshow(net2.regions.T)
ax.axis(False);

../../../_images/059bb969c0abd3b9449c075815d46f94da353304f2d679ed9ee26a86f80180ee.png

Now can see the particles are indicated by a green marker, the void pores by a red marker, the solid-void interconnections are in grey, and the void-void throats are black.

`phase_alias`#

In the previous section the default phase names created by snow2 were the uncreative “phase1” and “phase2”. The snow2 function allows for these to be given more desriptive names. This is done by supplying a dictionary with the voxels value as the key and the desired phase name as the value.

net3 = ps.networks.snow2(phases=im.astype(int) + 1, phase_alias={1: 'solid', 2: 'void'}, boundary_width=0)

pn3 = op.io.network_from_porespy(net3.network)
print(pn3)

══════════════════════════════════════════════════════════════════════════════
net : <openpnm.network.Network at 0x7f11fba96860>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Properties                                                   Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
throat.conns                                                    717 / 717
pore.coords                                                     251 / 251
pore.region_label                                               251 / 251
pore.phase                                                      251 / 251
throat.phases                                                   717 / 717
pore.region_volume                                              251 / 251
pore.equivalent_diameter                                        251 / 251
pore.local_peak                                                 251 / 251
pore.global_peak                                                251 / 251
pore.geometric_centroid                                         251 / 251
throat.global_peak                                              717 / 717
pore.inscribed_diameter                                         251 / 251
pore.extended_diameter                                          251 / 251
throat.inscribed_diameter                                       717 / 717
throat.total_length                                             717 / 717
throat.direct_length                                            717 / 717
throat.perimeter                                                717 / 717
pore.volume                                                     251 / 251
pore.surface_area                                               251 / 251
throat.cross_sectional_area                                     717 / 717
throat.equivalent_diameter                                      717 / 717
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
  #  Labels                                                 Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
pore.all                                                              251
throat.all                                                            717
pore.solid                                                            102
throat.solid_void                                                     467
pore.void                                                             149
throat.void_solid                                                     467
throat.void_void                                                      250
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Now we can see that the phases are named using our preferred nomenclature of “solid” and “void”. This is helpful if you have several phases like “electrolyte”, “binder”, “particles” for instance.

`boundary_width`#

Adding pores (or more generally nodes) on the boundaries of images is helpful when conducting transport simulations. Without boundary pores then boundary condtions must be applied to internal pores which are (a) hard to identify clearly and (b) make it impossible to know the exact length of the domain.

Note that boundary nodes are not physical and thus their sizes are ficticious. This can cause problems when computing the transport conductance for these nodes, so care must be taken.

The boundary_width argument was designed to work analogously to the numpy.pad function. You can specify different paddings along each axis and also along each face. For instance, here are some options (for a 2D image):

[0, 3]: 3 voxels only applied to the y-axis.
[0, [0, 3]]: 3 voxels only applied to the end of y-axis.
[0, [3, 0]]: 3 voxels only applied to the beginning of y-axis.

im = ps.generators.blobs([400, 400], porosity=0.6)
a = ps.networks.snow2(im, boundary_width=[[5, 0], 10])
fig, ax = plt.subplots(figsize=[6, 6])
ax.imshow(a.regions, origin='lower', interpolation='none')
ax.axis(False);

../../../_images/092d65c6c906742431d3f06f413c900b64588a055da1950abb53d0b885d25366.png

A few things to note:

The image contains regions on the left, right and bottom faces, which are added by padding the image after the segmentation was applied, then incrementing labels the padded regions so their labels are greater than any inside the original image.
These new regions are then treated like a normal region when the network topology and geometry is extracted.
The width of each boundary region is nearly the same as the corresponding internal regions, but they are only a few voxels thick (in this case because we specified 5 or 10). This means that the value of ‘pore.inscribed’ sphere for this regions will almost all be 5 or 10 voxels.
The solid touching the boundaries creates gaps in the boundary regions.
Two boundary regions neighboring each other are separated by 2 voxels. This prevents them from being considered as connected during the topology extraction stage.
The top and bottom of the image do not have any boundaries because we set the widths to 0 on this axis.

It is instructive to visualize the network over the image to see how the boundary pores work:

pn = op.io.network_from_porespy(a.network)
ax = op.visualization.plot_connections(pn, c='k')
ax = op.visualization.plot_coordinates(pn, markersize=50, ax=ax)
fig, ax = plt.gcf(), plt.gca()
ax.imshow(a.regions.T)
ax.axis(False);

../../../_images/da8f5c604a287d8d970f99e61b7d26701a4f5602be5bb235158086f913222cc6.png

`accuracy`#

This argument can either take the value of 'high' or 'standard', with standard being the default. It is not recommended to use 'high' since this uses some rather slow functions to measure the surface areas and perimeters. The increase in accuracy is basically futile since (a) the image has already been reduced to a coarse voxelization and (b) the pore network abstraction discards information about the actual microstructure in exchange for speed, so spending ages on the extraction step would undermine this purpose.

`voxel_size`#

This is the voxel size of the image, as in the length of one side of a voxel. Typical microtomography images are in the range of 1-5 um per voxel. NanoCT can be as low as 16 nm per voxel. FIB-SEM might be 4 nm per voxel. Note that all properties returned from the snow2 function, like ‘pore.volume’ will be in the same units are this value. It is strongly recommended to use SI, which is assumed in OpenPNM for most of the fluid property calculations.

`sigma` and `r_max`#

These values are passed on to the gaussian_blur and find_peaks functions.

The gaussian_blur is applied to the distance transfrom, which as discussed in the original snow paper is a very effective way to reduce the number of incorrect peaks that are found in the subsequent steps. A value of sigma=0.4 is the default.

The find_peaks function applies a maximum filter to the distance transform, which replaces each voxel with the maximum it finds with the specified neighborhood. It then looks for voxels that are the same between the filtered image and the original distance transform. Only voxels that are local maximas will be the same between the images.

The following code block illustrates the impact of applying the gaussian_blur filter on finding peaks.

np.random.seed(13)
im = ps.generators.overlapping_spheres([100, 100], r=7, porosity=0.7)
dt = edt(im)
mx1 = spim.maximum_filter(dt, footprint=ps.tools.ps_disk(r=4, smooth=False))
pk1 = dt == mx1
dt_blur = spim.gaussian_filter(dt, sigma=0.4)
mx2 = spim.maximum_filter(dt_blur, footprint=ps.tools.ps_disk(r=4, smooth=False))
pk2 = dt_blur == mx2
fig, ax = plt.subplots(3, 2, figsize=[8, 12])
ax[0][0].imshow(dt, origin='lower', interpolation='none')
ax[1][0].imshow(mx1, origin='lower', interpolation='none')
ax[2][0].imshow(pk1/im, origin='lower', interpolation='none')
ax[0][1].imshow(dt_blur, origin='lower', interpolation='none')
ax[1][1].imshow(mx2, origin='lower', interpolation='none')
ax[2][1].imshow(pk2/im, origin='lower', interpolation='none');

../../../_images/b14af3b9a4b2408718bf9ca5364e7f16396fffbf9ed10a8cf92607ba97757aa7.png

`peaks`#

Because finding correct peaks is such an important part marker-based watershed segmentationa, and therefore of snow2, it is possible to supply custom peaks found using any desired process. Obviously, if this is provided then sigma and r_max are ignored.

Let’s use the two sets of peaks found above:

net1 = ps.networks.snow2(phases=im, boundary_width=0, peaks=pk1)
net2 = ps.networks.snow2(phases=im, boundary_width=0, peaks=pk2)

fig, ax = plt.subplots(1, 2, figsize=[8, 4])
ax[0].imshow(net1.regions, origin='lower', interpolation='none')
ax[1].imshow(net2.regions, origin='lower', interpolation='none')
ax[0].axis(False)
ax[1].axis(False);

../../../_images/e9bdf853976a284aaa3fb8c8fd485ab03c6d5ba61adc4d888fae1ac0af1a8611.png

Clearly the image on the left is ‘over segmented’ while the image on the right seems to have identified the pores nearly perfectly, even without subsequent processing of the peaks that normally occur in snow2 like trim_saddle_points and trim_nearby_peaks.

`parallelization`#

The watershed function was identified as the bottleneck of the entire process, especially as images grew in size. Parallelization of the watershed function had not be throughly investigated in the literature, so Khan et al developed a basic protocol for dividing the image into chunks and processing each in parallel. The main challenge with this approach was ensuring sufficient overlap was allowed to prevent any artifacts, while not being overly conservative since large overlaps undermine the speedup obtained by parallelization. The parallelization if accomplished by dask which makes it trivial to dispatch duplicate jobs to different cores.

The snow2 function can optionally parallelize the watershed step. It calls porespy.filters.snow_partitioning_parallel behind the scenes, but this function requires several arguments. Instead of complicating the snow2 function signature, we instead chose to accept a dict of arguments and values which are then passed on. These are:

divs: The number of chunks to divide the image into. A value of 2 means the image is divided in hald along each axis, for a total of 8 chunks in a 3D image. It’s possible to specify a different number of divisions in each direction, like divs=[2, 2, 1]. This could be useful if the image to be processed is very narrow in one direction, such as a piece of fibrous electrode.
overlap: This controls the amount of overlap between each chunk in voxels. Specifying None (the default) means that the algorithm finds the overlap automatically. Khan el al found that using the maximum of the distance tranform was suitable.
cores: This controls the number of cores to use when processing the chunks. If not provided then all cores will be used. This may be problematic on computers without a large amount of RAM since each process will consume some RAM. In these cases in may be beneficial to set cores to a small number and simply wait longer for the result.

If parallelization=None then no parallelization is applied, while if parallelization={} then the default arguments of the snow_partitioning_parallel function are used.

In the following code block, the parallelization option is demonstrated with an overlap that is too small, resulting in a noticable artifacts in the center sectionof the watershed segmentation, where the divsions were located. For comparison, the same image is processed with the default overlap.

np.random.seed(0)
im = ps.generators.overlapping_spheres([400, 400], r=20, porosity=0.8)
net1 = ps.networks.snow2(phases=im, parallelization={'cores': 1, 'divs': [2, 2], 'overlap': 1})
net2 = ps.networks.snow2(phases=im, parallelization={'cores': 1, 'divs': [2, 2], 'overlap': None})

fig, ax = plt.subplots(1, 2, figsize=[8, 4])
temp = ps.tools.randomize_colors(net1.regions)
ax[0].imshow(temp, origin='lower', interpolation='none')
ax[0].axis(False)
temp = ps.tools.randomize_colors(net2.regions)
ax[1].imshow(temp, origin='lower', interpolation='none')
ax[1].axis(False);

../../../_images/43c5302f77486a1556f21ced72558ff8a8b2eb454ba4d05cbe10b12ccbc270a2.png

snow2#

phases#

phase_alias#

boundary_width#

accuracy#

voxel_size#

sigma and r_max#

peaks#

parallelization#