Obtaining the porosity profile along each principle axis

This example illustrates how to use the porosity_profile function to study the distribution of porosity in a sample.

Start by importing the usual packages, and setting the matplotlib style to something more reasonble:

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

[01:04:45] ERROR    PARDISO solver not installed, run `pip install pypardiso`. Otherwise,          _workspace.py:56
                    simulations will be slow. Apple M chips not supported.

Here we generate a test image using the blobs function. We will create an image with different lengths in each direction since this highlights an important point later on when plotting the results.

im = ps.generators.blobs(shape=[200, 300, 400])
fig, ax = plt.subplots()
ax.imshow(im[:, :, 50]);
../../../_images/009f4e5c4c33a7ae8ca20d608c5c994a31c617ca14891ae87f305e99460ba836.png

All results in PoreSpy are reported in voxels, so it is up to the user to scale these values to the physical size of each voxel. Here we’ll assume that each voxel corresponds to 5.9 microns.

voxel_size = 5.9  # microns/voxel

The porosity profile function is straight-forward. It only requires the binary image of the material, and a specification of along which axis to take the profile. Here we’ll do all three directions:

x_profile = ps.metrics.porosity_profile(im, 0)
y_profile = ps.metrics.porosity_profile(im, 1)
z_profile = ps.metrics.porosity_profile(im, 2)

Now finally we can plot the profile in each direction using matplotlib:

fig, ax = plt.subplots(figsize=(5, 4))
ax.plot(np.linspace(0, im.shape[0]*voxel_size, im.shape[0]), x_profile, 'b-', label='yz-plane', alpha=0.5)
ax.plot(np.linspace(0, im.shape[1]*voxel_size, im.shape[1]), y_profile, 'r-', label='xz-plane', alpha=0.5)
ax.plot(np.linspace(0, im.shape[2]*voxel_size, im.shape[2]), z_profile, 'g-', label='xy-plane', alpha=0.5)
ax.set_ylim([0, 1])
ax.set_ylabel('Porosity of slice')
ax.set_xlabel('Position of slice along given axis')
ax.legend();
../../../_images/c919cbce0508f1a24a8a92af60e5e5b198b8a4303d443fda4aeab889fbdaa6f3.png

Note how each line ends at a different position on the x-axis. This is caused by the fact that image had different total length in each direction. We can also do the plots with each axis normalized by their total length:

fig, ax = plt.subplots(figsize=(5, 4))
ax.plot(np.linspace(0, 1, im.shape[0]), x_profile, 'b-', label='yz-plane', alpha=0.5)
ax.plot(np.linspace(0, 1, im.shape[1], im.shape[1]), y_profile, 'r-', label='xz-plane', alpha=0.5)
ax.plot(np.linspace(0, 1, im.shape[2], im.shape[2]), z_profile, 'g-', label='xy-plane', alpha=0.5)
ax.set_ylim([0, 1])
ax.set_ylabel('Porosity of slice')
ax.set_xlabel('Fractional distance along given axis')
ax.legend();
../../../_images/aef5b2f443d74ba6fd6627cf1986c653bc170f92c8f226b20289ef0696d48ab0.png