porosimetry

This is essentially a drainage algorithm but returns a list of invasion sizes which can be analzyed to extract the pore size distribution.

import matplotlib.pyplot as plt
import numpy as np

import porespy as ps

ps.visualization.set_mpl_style()

The arguments and their defaults are:

import inspect

inspect.signature(ps.filters.porosimetry)

<Signature (im: numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[~_ScalarT]], dt: numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[~_ScalarT]] = None, inlets: numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[~_ScalarT]] = None, sizes: int = None, method: Literal['dsi', 'fft', 'dt'] = 'dt', smooth: bool = True)>

im

The function works on 2D and 3D images:

im = ps.generators.blobs(shape=[200, 200])
sizes = ps.filters.porosimetry(im=im)
fig, ax = plt.subplots(1, 1, figsize=[6, 6])
ax.imshow(sizes / im, origin="lower", interpolation="none")
ax.axis(False);

A cumulative histogram of the voxel values is analogous to a porosimetry experiment:

fig, ax = plt.subplots(1, 1, figsize=[6, 6])
ax.hist(sizes[im], bins=25, cumulative=True, edgecolor="k", density=True)
ax.set_xlabel("pore radius in voxels")
ax.set_ylabel("fraction of pore volume smaller than stated size");

sizes

The number of bins to use, or a list of actual bins to use:

sizes = ps.filters.porosimetry(im=im, sizes=5)
fig, ax = plt.subplots(1, 2, figsize=[12, 6])
ax[0].imshow(sizes / im, origin="lower", interpolation="none")
ax[0].axis(False)
ax[1].hist(sizes[im], bins=25, cumulative=True, edgecolor="k", density=True)
ax[1].set_xlabel("pore radius in voxels")
ax[1].set_ylabel("fraction of pore volume smaller than stated size");

inlets

By default the invading fluid enters from all sizes (typical of a porosimetry experiment), but this can be controlled by specifying inlets:

inlets = np.zeros_like(im)
inlets[0, ...] = True
sizes = ps.filters.porosimetry(im=im, inlets=inlets)
fig, ax = plt.subplots(1, 2, figsize=[12, 6])
ax[0].imshow(sizes / im, origin="lower", interpolation="none")
ax[0].axis(False)
ax[1].hist(sizes[im], bins=25, cumulative=True, edgecolor="k", density=True)
ax[1].set_xlabel("pore radius in voxels")
ax[1].set_ylabel("fraction of pore volume smaller than stated size");

method

This controls which method is used. The default is a ‘dt’ which uses a the threshold of a distance transform to perform an erosion, then additional distance transforms to perform the dilation to generate spheres. Other options are 'conv' which uses convolution for the dilation step, or 'dsi' which direct sphere insertion for the dilation. All results should be exactly the same if they are performed at integer values of sizes. This is necessary since the convolution and direct sphere insertion methods use discrete sphere radii to find insertion sites, while the distance transform method can interpret the fractional distance values as different steps.

sizes = np.arange(50, 1, -1)
sizes1 = ps.filters.porosimetry(im=im, method="dt", sizes=sizes)
sizes2 = ps.filters.porosimetry(im=im, method="conv", sizes=sizes)
sizes3 = ps.filters.porosimetry(im=im, method="dsi", sizes=sizes)
fig, ax = plt.subplots(1, 3, figsize=[16, 6])
ax[0].imshow(sizes1 / im, origin="lower", interpolation="none")
ax[0].axis(False)
ax[1].imshow(sizes2 / im, origin="lower", interpolation="none")
ax[1].axis(False)
ax[2].imshow(sizes3 / im, origin="lower", interpolation="none")
ax[2].axis(False)
assert np.all(sizes1 == sizes2)
assert np.all(sizes2 == sizes3)