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()

/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

The arguments and their defaults are:

import inspect
inspect.signature(ps.filters.porosimetry)

<Signature (im, sizes=25, inlets=None, access_limited=True, mode='hybrid', divs=1)>

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);
../../../_images/88a876b73bb210047c78feecce17b21ddaeeedb3fdfad0355af404f7274fb617.png

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');
../../../_images/48bab3441a4ba7fd9905e43ee3c821f51c37806797c6d2baa6752500aafc2ad2.png

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');
../../../_images/669d72cee0e51cd377a5b41425a22f0de9aa4e9129dfe58d44b3c379398baaac.png

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');
../../../_images/abdbb23932828f129e2cfb452bc0f839eb05c7ccd867aa8afabb6bcd016026e6.png

access limited#

If True (default) then the invading fluid can only enter from the edges of the image (or inlets). If set to False then fluid will enter into the largest possible voids first regardless of physical connectivity. Note that access_limited=False is identical to use the local_thickness filter.

sizes1 = ps.filters.porosimetry(im=im, access_limited=False)
sizes2 = ps.filters.porosimetry(im=im, access_limited=True)
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(sizes2[im], bins=25, cumulative=True, edgecolor='k', density=True, label='access limited')
ax[1].hist(sizes1[im], bins=25, cumulative=True, edgecolor='k', density=True, label='not access limited')
ax[1].set_xlabel('pore radius in voxels')
ax[1].set_ylabel('fraction of pore volume smaller than stated size')
ax[1].legend();
../../../_images/7aa34072ce1dcb5d39e6ac3e1f79330a6cc888d98f1b7743cd158ddfa7c81903.png

As can be seen the case with access limitations appears to have smaller pores, which is because the invading fluid cannot reach the large pores that are shielded by smaller ones, so their volume gets attributed to the smaller voids.

mode#

This controls which method is used. The default is a ‘hybrid’ which uses a the threshold of a distance transform to perform an erosion, then fft-based dilation to generate spheres. Other options are dt which uses a distance transform for both steps, or ‘mio’ which uses morphological operations for both steps. The selected method will affect the speed, but this depends on the computer being used. All results should be exactly the same, and this is ensured in a unit test:

sizes1 = ps.filters.porosimetry(im=im, mode='dt')
sizes2 = ps.filters.porosimetry(im=im, mode='hybrid')
sizes3 = ps.filters.porosimetry(im=im, mode='mio')
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);
../../../_images/13ff28df93eab748a1a59e424e8b7709c0763ec3b0243aab5f37858312b1cc8c.png

divs#

Allows the fucntion to run in parallel by subdiving the image into the given number of division, either a scalar or a list for each axis. The result will be identical but should be faster on computers with many cores.

sizes1 = ps.filters.porosimetry(im=im, divs=1)
sizes2 = ps.filters.porosimetry(im=im, divs=[2, 1])
sizes3 = ps.filters.porosimetry(im=im, divs=3)
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);
../../../_images/13ff28df93eab748a1a59e424e8b7709c0763ec3b0243aab5f37858312b1cc8c.png