Local thickness

This example explains how to use the local_thickness filter to get information about the pore size distribution from an image. The local thickness is probably the closest you can get to an actual pore size distribution. Unlike porosimetry experiments or simulations it is unaffected by artifacts such as edge effects. The implementation in PoreSpy is slightly different than the common approach done in ImageJ, as will be explained below.

These notebooks were generated with a specific version of Python, PoreSpy, Numpy, etc. For reference to future viewers, the version information is given below. This notebook may or may not work with later versions, but we can assert it works with the version listed below:

Import the needed packages from the Scipy stack:

import numpy as np
import porespy as ps
import scipy.ndimage as spim
import matplotlib.pyplot as plt
ps.visualization.set_mpl_style()

Generate Test Image

Start by generating an image. We’ll use the RSA generator for fun:

im = np.zeros([300, 300])
im = ps.generators.RSA(im, r=20, volume_fraction=0.2)
im = ps.generators.RSA(im, r=15, volume_fraction=0.4)
im = ps.generators.RSA(im, r=10, volume_fraction=0.6)
im = im == 0
fig, ax = plt.subplots()
ax.imshow(im);

Apply Local Thickness Filter

The local thickness filter is called by simply passing in the image. Like all filters in PoreSpy it is applied to the foreground, indicated by 1’s or True:

thk = ps.filters.local_thickness(im, mode='dt')
fig, ax = plt.subplots()
ax.imshow(thk, cmap=plt.cm.jet);
import tqdm.notebook as tqdm
from time import sleep
for i in tqdm.tqdm(range(10)):
    sleep(0.1)

from tqdm.notebook import tqdm, tqdm_notebook
for i in tqdm_notebook(range(10)):
    sleep(0.1)

Extracting PSD as a Histogram

Obtaining pore size distribution information from this image requires obtaining a histogram of voxel values. A function in the metrics module does this for us:

psd = ps.metrics.pore_size_distribution(im=thk)

The result returned into psd is a “named-tuple”, which is a list of arrays that can be accessed by location (i.e. psd[0]), but has the added benefit of accessing arrays by name so you know what you’re getting. You can print a list of available arrays as follows:

print(psd._fields)

('LogR', 'pdf', 'cdf', 'satn', 'bin_centers', 'bin_edges', 'bin_widths')

Let’s plot a pore-size distribution histogram:

fig, ax = plt.subplots(figsize=(5, 4))
ax.set_xlabel('log(Pore Radius) [voxels]')
ax.set_ylabel('Normalized Volume Fraction')
ax.bar(x=psd.LogR, height=psd.pdf, width=psd.bin_widths, edgecolor='k');

PoreSpy Implementation

The local_thickness filter in PoreSpy is implemented differently that the normal approach such as the ImageJ plugin, though the end result is comparible though not identical.

In our approach, we use a form of image dilation and erosion. We start with a large spherical structuring element, and note all the places where this fits in the pore space. This gives a result like that below for a structuring element of radius R=10:

R = 10
strel = ps.tools.ps_disk(R)
im_temp = spim.binary_opening(im, structure=strel)
fig, ax = plt.subplots()
ax.imshow(im_temp*2.0 + ~im);

The key is to make a master array containing the numerical value of the largest sphere that covers each voxel. We’ll initialize a new array with the current locations where R=10 fits:

im_result = im_temp*R

Now this is repeated for a range of decreasing structuring element sizes. For illustration, do R = 8:

R = 8
strel = ps.tools.ps_disk(R)
im_temp = spim.binary_opening(im, structure=strel)

This new image must be added to the im_result array, but only in places that were not filled at any larger radius. This is done using boolean logic as follows:

im_result[(im_result == 0)*im_temp] = R

There are now 2 values in the im_results array indicating the locations where the structuring element of size 10 fits, and where size 8 fit on the subsequent step:

fig, ax = plt.subplots()
ax.imshow(im_result + ~im);

The procedure is then repeated for smaller structuring elements down to R = 1. It’s possible to specify which sizes are used, but by default all integers between \(R_{max}\) and 1. This yields the image showed above:

fig, ax = plt.subplots()
ax.imshow(thk, cmap=plt.cm.jet);