Function Reference

Binarization

This module contains functions for thresholding greyscale images into binary images (i.e. True and False) indicating the phases.

porespy.binarization.simple_otsu(im, trim_solid=True)

Uses Otsu’s method to find a threshold, then uses binary opening and closing to remove noise.

Parameters:

im : ND-image

The greyscale image of the porous medium to be binarized.

trim_solid : Boolean

If True (default) then all solid voxels not connected to an image boundary are trimmed.

Returns:

An ND-image the same size as im but with True and False values

indicating the binarized or segmented phases.

Examples

>>> im = ps.generators.blobs([300, 300], porosity=0.5)
>>> im = ps.generators.add_noise(im)
>>> im = spim.gaussian_filter(im, sigma=1)
>>> im = simple_otsu(im)

Filters

This module contains a variety of functions for altering image based on the structural characteristics, such as pore sizes.

porespy.filters.apply_chords(im, spacing=0, axis=0, trim_edges=True)

Adds chords to the void space in the specified direction. The chords are separated by 1 voxel plus the provided spacing.

Parameters:

im : ND-array

An image of the porous material with void marked as True.

spacing : int (default = 0)

Chords are automatically separated by 1 voxel and this argument increases the separation.

axis : int (default = 0)

The axis along which the chords are drawn.

trim_edges : bool (default = True)

Whether or not to remove chords that touch the edges of the image. These chords are artifically shortened, so skew the chord length distribution

Returns:

An ND-array of the same size as `im` with True values indicating the

chords.

See also

apply_chords_3D

porespy.filters.apply_chords_3D(im, spacing=0, trim_edges=True)

Adds chords to the void space in all three principle directions. The chords are seprated by 1 voxel plus the provided spacing. Chords in the X, Y and Z directions are labelled 1, 2 and 3 resepctively.

Parameters:

im : ND-array

A 3D image of the porous material with void space marked as True.

spacing : int (default = 0)

Chords are automatically separed by 1 voxel on all sides, and this argument increases the separation.

trim_edges : bool (default = True)

Whether or not to remove chords that touch the edges of the image. These chords are artifically shortened, so skew the chord length distribution

Returns:

An ND-array of the same size as `im` with values of 1 indicating

x-direction chords, 2 indicating y-direction chords, and 3 indicating

z-direction chords.

See also

apply_chords

Notes

The chords are separated by a spacing of at least 1 voxel so that tools that search for connected components, such as scipy.ndimage.label can detect individual chords.

porespy.filters.local_thickness(im)

For each voxel, this functions calculates the radius of the largest sphere that both engulfs the voxel and fits entirely within the foreground. This is not the same as a simple distance transform, which finds the largest sphere that could be centered on each voxel.

Parameters:

im : array_like

A binary image with the phase of interest set to True

Returns:

An image with the pore size values in each voxel

Notes

The term foreground is used since this function can be applied to both pore space or the solid, whichever is set to True.

porespy.filters.fill_blind_pores(im)

Fills all pores that are not connected to the edges of the image.

Parameters:

im : ND-array

The image of the porous material

Returns:

A version of im but with all the disconnected pores removed.

See also

find_disconnected_voxels

porespy.filters.find_disconnected_voxels(im, conn=None)

This identifies all pore (or solid) voxels that are not connected to the edge of the image. This can be used to find blind pores, or remove artifacts such as solid phase voxels that are floating in space.

Parameters:

im : ND-image

A Boolean image, with True values indicating the phase for which disconnected voxels are sought.

conn : int

For 2D the options are 4 and 8 for square and diagonal neighbors, while for the 3D the options are 6 and 26, similarily for square and diagonal neighbors. The default is max

Returns:

An ND-image the same size as im, with True values indicating voxels of

the phase of interest (i.e. True values in the original image) that are

not connected to the outer edges.

Notes

The returned array (e.g. holes) be used to trim blind pores from im using: im[holes] = False

porespy.filters.flood()

Floods/fills each region in an image with a single value based on the specific values in that region. The mode argument is used to determine how the value is calculated. A region is defined as a connected cluster of voxels surrounded by 0’s for False’s.

Parameters:

im : array_like

An ND image with isolated regions containing 0’s elsewhere.

mode : string

Specifies how to determine which value should be used to flood each region. Options are:

*’max’* : Floods each region with the local maximum in that region

*’min’* : Floods each region the local minimum in that region

*’size’* : Floods each region with the size of that region

Returns:

An ND-array the same size as im with new values placed in each

forground voxel based on the mode.

porespy.filters.local_thickness(im)

For each voxel, this functions calculates the radius of the largest sphere that both engulfs the voxel and fits entirely within the foreground. This is not the same as a simple distance transform, which finds the largest sphere that could be centered on each voxel.

Parameters:

im : array_like

A binary image with the phase of interest set to True

Returns:

An image with the pore size values in each voxel

Notes

The term foreground is used since this function can be applied to both pore space or the solid, whichever is set to True.

porespy.filters.porosimetry(im, npts=25, sizes=None, inlets=None, access_limited=True)

Performs a porosimetry simulution on the image

Parameters:

im : ND-array

An ND image of the porous material containing True values in the pore space.

npts : scalar

The number of invasion points to simulate. Points will be generated spanning the range of sizes in the distance transform. The default is 25 points.

sizes : array_like

The sizes to invade. Use this argument instead of npts for more control of the range and spacing of points.

inlets : ND-array, boolean

A boolean mask with True values indicating where the invasion enters the image. By default all faces are considered inlets, akin to a mercury porosimetry experiment. Users can also apply solid boundaries to their image externally before passing it in, allowing for complex inlets like circular openings, etc.

access_limited : Boolean

This flag indicates if the intrusion should only occur from the surfaces (access_limited is True, which is the default), or if the invading phase should be allowed to appear in the core of the image. The former simulates experimental tools like mercury intrusion porosimetry, while the latter is useful for comparison to gauge the extent of shielding effects in the sample.

Returns:

An ND-image with voxel values indicating the sphere radius at which it

becomes accessible from the inlets. This image can be used to find

invading fluid configurations as a function of applied capillary pressure

by applying a boolean comparison: inv_phase = im > r where r is

the radius (in voxels) of the invading sphere. Of course, r can be

converted to capillary pressure using your favorite model.

porespy.filters.trim_extrema(im, h, mode='maxima')

This trims local extrema in greyscale values by a specified amount, essentially decapitating peaks or flooding valleys, or both.

Parameters:

im : ND-array

The image whose extrema are to be removed

h : scalar

The height to remove from each peak or fill in each valley

mode : string {‘maxima’ | ‘minima’ | ‘extrema’}

Specifies whether to remove maxima or minima or both

Returns:

A copy of the input image with all the peaks and/or valleys removed.

Notes

This function is referred to as imhmax or imhmin in Mablab.

porespy.filters.trim_floating_solid(im)

Removes all solid that that is not attached to the edges of the image.

Parameters:

im : ND-array

The image of the porous material

Returns:

A version of im but with all the disconnected solid removed.

See also

find_disconnected_voxels

Generators

This module contains a variety of functions for generating artificial images of porous materials, generally for testing, validation, debugging, and illustration purposes.

porespy.generators.add_noise(im, u_void=0.8, u_solid=0.2, s_void=0.15, s_solid=0.15)

Add some normally distributed noise values to the image. This is useful for testing binarization routines.

Parameters:

im : ND-array

The image of the porous media, with 1’s or True’s denoting voids

u_solid : float

The mean greyscale value in the solid space (default is 0.2)

s_solid : float

The standard deviation of the greyscale values in the solid phase (the default is 0.15)

u_void: float

The mean greyscale value in the void space (default is 0.8)

s_void: float

The standard deviation of the greyscale values in the void phase (the default is 0.15)

Returns:

A greyscale image with the same shape as im, but normally distributed

noise values in the void and solid phase. A histogram of the default image

should reveal to distinguishable but overlapping peaks. This can be

adjusted using the mean and standard deviation arguments.

porespy.generators.blobs(shape, porosity=0.5, blobiness=1)

Generates an image containing amorphous blobs

Parameters:

shape : list

The size of the image to generate in [Nx, Ny, Nz] where N is the number of voxels

porosity : float

If specified, this will threshold the image to the specified value prior to returning. If no value is given (the default), then the scalar noise field is returned.

blobiness : array_like (default = 1)

Controls the morphology of the blobs. A higher number results in a larger number of small blobs. If a vector is supplied then the blobs are anisotropic.

Returns:

If porosity is given, then a boolean array with True values denoting

the pore space is returned. If not, then normally distributed and

spatially correlated randomly noise is returned.

See also

norm_to_uniform

porespy.generators.bundle_of_tubes(shape, spacing)

Create a 3D image of a bundle of tubes, in the form of a rectangular plate with randomly sized holes through it.

Parameters:

shape : list

The size the image, with the 3rd dimension indicating the plate thickness. If the 3rd dimension is not given then a thickness of 1 voxel is assumed.

spacing : scalar

The center to center distance of the holes. The hole sizes will be randomly distributed between this values down to 3 voxels.

Returns:

A boolean array with True values denoting the pore space

porespy.generators.circle_pack(shape, radius, offset=0, packing='square')

Generates a 2D packing of circles

Parameters:

shape : list

The size of the image to generate in [Nx, Ny] where N is the number of voxels in each direction

radius : scalar

The radius of circles in the packing (in pixels)

offset : scalar

The amount offset (+ or -) to add between pore centers (in pixels).

packing : string

Specifies the type of cubic packing to create. Options are ‘square’ (default) and ‘triangular’.

Returns:

A boolean array with True values denoting the pore space

porespy.generators.cylinders(shape, radius, nfibers, phi_max=0, theta_max=90)

Generates a binary image of overlapping cylinders. This is a good approximation of a fibrous mat.

Parameters:

phi_max : scalar

A value between 0 and 90 that controls the amount that the fibers lie out of the XY plane, with 0 meaning all fibers lie in the XY plane, and 90 meaning that fibers are randomly oriented out of the plane by as much as +/- 90 degrees.

theta_max : scalar

A value between 0 and 90 that controls the amount rotation in the XY plane, with 0 meaning all fibers point in the X-direction, and 90 meaning they are randomly rotated about the Z axis by as much as +/- 90 degrees.

Returns:

A boolean array with True values denoting the pore space

porespy.generators.insert_shape(im, center, element, value=1)

porespy.generators.line_segment(X0, X1)

Calculate the voxel coordinates of a straight line between the two given end points

Parameters:

X0 and X1 : array_like

The [x, y, z] coordinates of the start and end points of the line.

Returns:

A list of lists containing the X, Y, and Z coordinates of all voxels

that should be drawn between the start and end points to create a solid line.

porespy.generators.noise(shape, porosity=None, octaves=3, frequency=32, mode='simplex')

Generate a field of spatially correlated random noise using the Perlin noise algorithm, or the updated Simplex noise algorithm.

Parameters:

shape : array_like

The size of the image to generate in [Nx, Ny, Nz] where N is the number of voxels.

porosity : float

If specified, this will threshold the image to the specified value prior to returning. If no value is given (the default), then the scalar noise field is returned.

octaves : int

Controls the texture of the noise, with higher octaves giving more complex features over larger length scales.

frequency : array_like

Controls the relative sizes of the features, with higher frequencies giving larger features. A scalar value will apply the same frequency in all directions, given an isotropic field; a vector value will apply the specified values along each axis to create anisotropy.

mode : string

Which noise algorithm to use, either ‘simplex’ (default) or ‘perlin’.

Returns:

If porosity is given, then a boolean array with True values denoting

the pore space is returned. If not, then normally distributed and

spatially correlated randomly noise is returned.

See also

norm_to_uniform

Notes

This method depends the a package called ‘noise’ which must be compiled. It is included in the Anaconda distribution, or a platform specific binary can be downloaded.

porespy.generators.norm_to_uniform(im, scale=None)

Take an image with normally distributed greyscale values and converts it to a uniform (i.e. flat) distribution. It’s also possible to specify the lower and upper limits of the uniform distribution.

Parameters:

im : ND-image

The image containing the normally distributed scalar field

scale : [low, high]

A list or array indicating the lower and upper bounds for the new randomly distributed data. The default is None, which uses the max and min of the original image as the the lower and upper bounds.

Returns:

An ND-image the same size as im with uniformly distributed greyscale

values spanning the specified range, if given.

porespy.generators.overlapping_spheres(shape, radius, porosity)

Generate a packing of overlapping mono-disperse spheres

Parameters:

shape : list

The size of the image to generate in [Nx, Ny, Nz] where Ni is the number of voxels in the *i*th direction.

radius : scalar

The radius of spheres in the packing.

porosity : scalar

The porosity of the final image. This number is approximated by the method so the returned result may not have exactly the specified value.

Returns:

A boolean array with True values denoting the pore space

Notes

This method can also be used to generate a dispersion of hollows by treating porosity as solid volume fraction and inverting the returned image.

porespy.generators.polydisperse_spheres(shape, porosity, dist, nbins=5)

Create an image of spheres with a distribution of radii.

Parameters:

shape : list

The size of the image to generate in [Nx, Ny, Nz] where Ni is the number of voxels in each direction. If shape is only 2D, then an image of polydisperse disks is returns

porosity : scalar

The porosity of the image, defined as the number of void voxels divided by the number of voxels in the image. The specified value is only matched approximately, so it’s suggested to check this value after the image is generated.

dist : scipy.stats distribution object

This should be an initialized distribution the large number of options in the scipy.stats submodule. For instance, a normal distribution with a mean of 20 and a standard deviation of 10 can be obtain with dist = scipy.stats.norm(loc=20, scale=10)

nbins : scalar

The number of discrete sphere sizes that will be used to generate the image. This function generates nbins images of monodisperse spheres that span 0.05 and 0.95 of the possible values produced by the provided distribution, then overlays them to get polydispersivity.

Returns:

A boolean array with True values denoting the pore space

porespy.generators.sphere_pack(shape, radius, offset=0, packing='sc')

Generates a cubic packing of spheres

Parameters:

shape : list

The size of the image to generate in [Nx, Ny, Nz] where N is the number of voxels in each direction.

radius : scalar

The radius of spheres in the packing

offset : scalar

The amount offset (+ or -) to add between pore centers.

packing : string

Specifies the type of cubic packing to create. Options are:

‘sc’ : Simple Cubic (default) ‘fcc’ : Face Centered Cubic ‘bcc’ : Body Centered Cubic

Returns:

A boolean array with True values denoting the pore space

porespy.generators.voronoi_edges(shape, radius, ncells, flat_faces=True)

Create an image of the edges in a Voronoi tessellation

Parameters:

shape : array_like

The size of the image to generate in [Nx, Ny, Nz] where Ni is the number of voxels in each direction.

radius : scalar

The radius to which Voronoi edges should be dilated in the final image.

ncells : scalar

The number of Voronoi cells to include in the tesselation.

flat_faces : Boolean

Whether the Voronoi edges should lie on the boundary of the image (True), or if edges outside the image should be removed (False).

Returns:

A boolean array with True values denoting the pore space

Metrics

A suite of tools for determine key metrics about an image. Typically these are applied to an image after applying a filter, but a few functions can be applied directly to the binary image.

porespy.metrics.chord_length_counts(im)

Determines the length of each chord in the supplied image by looking at its size.

Parameters:

im : ND-array

An image containing chords drawn in the void space.

Returns:

A 1D array with one element for each chord, containing its length.

Notes

The returned array can be passed to plt.hist to plot the histogram,

or to sp.histogram to get the histogram data directly. Another useful

function is sp.bincount which gives the number of chords of each

length in a format suitable for plt.plot.

porespy.metrics.chord_length_distribution(im, bins=25, log=False)

Determines the distribution of chord lengths in a image containing chords.

Parameters:

im : ND-image

An image with chords drawn in the pore space, as produced by apply_chords or apply_chords_3d.

bins : scalar or array_like

If a scalar is given it is interpreted as the number of bins to use, and if an array is given they are used as the bins directly.

log : Boolean

If true, the logarithm of the chord lengths will be used, which can make the data more clear.

Returns:

A tuple containing the chord_length_bins, and four separate pieces of

information: cumulative_chord_count and cumulative_chord_length,

as well as the differenial_chord_count and

differential_chord_length.

porespy.metrics.pore_size_density(im, bins=10, voxel_size=1)

Computes the histogram of the distance transform as an estimator of the pore sizes in the image.

Parameters:

im : ND-array

Either a binary image of the pore space with True indicating the pore phase (or phase of interest), or a pre-calculated distance transform which can save time.

bins : int or array_like

This number of bins (if int) or the location of the bins (if array). This argument is passed directly to Scipy’s histogram function so see that docstring for more information. The default is 10 bins, which reduces produces a relatively smooth distribution.

voxel_size : scalar

The size of a voxel side in preferred units. The default is 1, so the user can apply the scaling to the returned results after the fact.

Returns:

A tuple containing two 1D arrays: ``radius `` is the radius of the voxels,

and count is the number of voxels that are within R of the solid.

Notes

This function should not be taken as a pore size distribution in the explict sense, but rather an indicator of the sizes in the image.

References

[1] Torquato, S. Random Heterogeneous Materials: Mircostructure and Macroscopic Properties. Springer, New York (2002) - See page 292

porespy.metrics.pore_size_distribution(im)

Calculate drainage curve based on the image produced by the porosimetry function.

Parameters:

im : ND-array

The array of entry sizes as produced by the porosimetry function.

Returns:

Rp, Snwp: Two arrays containing (a) the radius of the penetrating

sphere (in voxels) and (b) the volume fraction of pore phase voxels

that are accessible from the specfied inlets.

Notes

This function normalizes the invading phase saturation by total pore volume of the dry image, which is assumed to be all voxels with a value equal to 1. To do porosimetry on images with large outer regions, use the `find_outer_region` function then set these regions to 0 in the input image. In future, this function could be adapted to apply this check by default.

porespy.metrics.porosity(im)

Calculates the porosity of an image assuming 1’s are void space and 0’s are solid phase. All other values are ignored.

Parameters:

im : ND-array

Image of the void space with 1’s indicating void space (or True) and 0’s indicating the solid phase (or False).

Returns:

porosity : float

Calculated as the sum of all 1’s divided by the sum of all 1’s and 0’s.

Notes

This function assumes void is represented by 1 and solid by 0, and all other values are ignored. This is useful, for example, for images of cylindrical cores, where all voxels outside the core are labelled with 2.

Alternatively, images can be processed with find_disconnected_voxels to get an image of only blind pores. This can then be added to the orignal image such that blind pores have a value of 2, thus allowing the calculation of accessible porosity, rather than overall porosity.

porespy.metrics.representative_elementary_volume(im, npoints=1000)

Calculates the porosity of the image as a function subdomain size. This function extracts a specified number of subdomains of random size, then finds their porosity.

Parameters:

im : ND-array

The image of the porous material

npoints : int

The number of randomly located and sized boxes to sample. The default is 1000.

Returns:

A tuple containing the ND-arrays: The subdomain volume and its

porosity. Each of these arrays is npoints long. They can be

conveniently plotted by passing the tuple to matplotlib’s plot function

using the * notation: plt.plot(*the_tuple, 'b.'). The resulting plot

is similar to the sketch given by Bachmat and Bear [1]

Notes

This function is frustratingly slow. Profiling indicates that all the time is spent on scipy’s sum function which is needed to sum the number of void voxels (1’s) in each subdomain.

Also, this function is primed for parallelization since the npoints are calculated independenlty.

References

[1] Bachmat and Bear. On the Concept and Size of a Representative Elementary Volume (Rev), Advances in Transport Phenomena in Porous Media (1987)

porespy.metrics.two_point_correlation_bf(im, spacing=10)

Calculates the two-point correlation function using brute-force (see Notes)

Parameters:

im : ND-array

The image of the void space on which the 2-point correlation is desired

spacing : int

The space between points on the regular grid that is used to generate the correlation (see Notes)

Returns:

A tuple containing the x and y data for plotting the two-point correlation

function, using the *args feature of matplotlib’s plot function. The x

array is the distances between points and the y array is corresponding

probabilities that points of a given distance both lie in the void space.

The distance values are binned as follows:

bins = range(start=0, stop=sp.amin(im.shape)/2, stride=spacing)

Notes

The brute-force approach means overlaying a grid of equally spaced points onto the image, calculating the distance between each and every pair of points, then counting the instances where both pairs lie in the void space.

This approach uses a distance matrix so can consume memory very quickly for large 3D images and/or close spacing.

porespy.metrics.two_point_correlation_fft(im, pad=False)

Calculates the two-point correlation function using fourier transforms

Parameters:

im : ND-array

The image of the void space on which the 2-point correlation is desired

pad : bool

The image is padded with Trues or 1’s depending on dtype around border

Returns:

A tuple containing the x and y data for plotting the two-point correlation

function, using the *args feature of matplotlib’s plot function. The x

array is the distances between points and the y array is corresponding

probabilities that points of a given distance both lie in the void space.

Notes

The fourier transform approach utilizes the fact that the autocorrelation function is the inverse FT of the power spectrum density. For background read the Scipy fftpack docs and for a good explanation see: http://www.ucl.ac.uk/~ucapikr/projects/KamilaSuankulova_BSc_Project.pdf

Tools

This module contains a variety of functions for manipulating images in ways that do NOT return a modified version of the original image.

porespy.tools.extend_slice(s, shape, pad=1)

Adjust slice indices to include additional voxles around the slice. The key to this function is that is does bounds checking to ensure the indices don’t extend outside the image.

Parameters:

s : list of slice objects

A list (or tuple) of N slice objects, where N is the number of dimensions in the image.

shape : array_like

The shape of the image into which the slice objects apply. This is used to check the bounds to prevent indexing beyond the image.

pad : int

The number of voxels to expand in each direction.

Returns:

A list slice objects with the start and stop attributes respectively

incremented and decremented by 1, without extending beyond the image

boundaries.

Examples

>>> from scipy.ndimage import label, find_objects
>>> from porespy.tools import extend_slice
>>> im = sp.array([[1, 0, 0], [1, 0, 0], [0, 0, 1]])
>>> labels = label(im)[0]
>>> s = find_objects(labels)

Using the slices returned by find_objects, set the first label to 3

>>> labels[s[0]] = 3
>>> print(labels)
[[3 0 0]
 [3 0 0]
 [0 0 2]]

Next extend the slice, and use it to set the values to 4

>>> s_ext = extend_slice(s[0], shape=im.shape, pad=1)
>>> labels[s_ext] = 4
>>> print(labels)
[[4 4 0]
 [4 4 0]
 [4 4 2]]

As can be seen by the location of the 4s, the slice was extended by 1, and also handled the extension beyond the boundary correctly.

porespy.tools.extract_subsection(im, shape)

Extracts the middle section of a image

Parameters:

im : ND-array

Image from which to extract the subsection

shape : array_like

Can either specify the size of the extracted section or the fractional size of the image to extact.

Returns:

An ND-array of size given by the shape argument, taken from the center

of the image.

Examples

>>> import scipy as sp
>>> from porespy.tools import extract_subsection
>>> im = sp.array([[1, 1, 1, 1], [1, 2, 2, 2], [1, 2, 3, 3], [1, 2, 3, 4]])
>>> print(im)
[[1 1 1 1]
 [1 2 2 2]
 [1 2 3 3]
 [1 2 3 4]]
>>> im = extract_subsection(im=im, shape=[2, 2])
>>> print(im)
[[2 2]
 [2 3]]

porespy.tools.extract_cylinder(im, r=None, axis=0)

Returns a cylindrical section of the image of specified radius. This is useful for making square images look like cylindrical cores such as those obtained from X-ray tomography.

Parameters:

im : ND-array

The image of the porous material

r : scalr

The radius of the cylinder to extract. If none if given then the default is the largest cylinder that can fit inside the x-y plane.

axis : scalar

The axis along with the cylinder will be oriented.

Returns:

An ND-image the same size im with True values indicating the void space

but with the sample trimmed to a cylindrical section in the center of the

image. The region outside the cylindrical section is labeled with True

values since it is open space.

porespy.tools.find_outer_region(im, r=0)

Finds regions of the image that are outside of the solid matrix. This function uses the rolling ball method to define where the outer region ends and the void space begins.

This function is particularly useful for samples that do not fill the entire rectangular image, such as cylindrical cores or samples with non- parallel faces.

Parameters:

im : ND-array

Image of the porous material with 1’s for void and 0’s for solid

r : scalar

The radius of the rolling ball to use. If not specified then a value is calculated as twice maximum of the distance transform. The image size is padded by this amount in all directions, so the image can become quite large and unwieldy if too large a value is given.

Returns:

A boolean mask the same shape as im, containing True in all voxels

identified as outside the sample.

porespy.tools.get_border(shape, thickness=1, mode='edges')

Creates an array of specified size with corners, edges or faces labelled as True. This can be used as mask to manipulate values laying on the perimeter of an image.

Parameters:

shape : array_like

The shape of the array to return. Can be either 2D or 3D.

thickness : scalar (default is 1)

The number of pixels/voxels to place along perimeter.

mode : string

The type of border to create. Options are ‘faces’, ‘edges’ (default) and ‘corners’. In 2D ‘faces’ and ‘edges’ give the same result.

Returns:

An ND-array of specified shape with True values at the perimeter and False

elsewhere.

Examples

>>> import porespy as ps
>>> import scipy as sp
>>> mask = ps.tools.get_border(shape=[3, 3], mode='corners')
>>> print(mask)
[[ True False  True]
 [False False False]
 [ True False  True]]
>>> mask = ps.tools.get_border(shape=[3, 3], mode='edges')
>>> print(mask)
[[ True  True  True]
 [ True False  True]
 [ True  True  True]]

porespy.tools.get_slice(im, center, size, pad=0)

Given a center location and radius of a feature, returns the slice object into the im that bounds the feature but does not extend beyond the image boundaries.

Parameters:

im : ND-image

The image of the porous media

center : array_like

The coordinates of the center of the feature of interest

size : array_like or scalar

The size of the feature in each direction. If a scalar is supplied, this implies the same size in all directions.

pad : scalar or array_like

The amount to pad onto each side of the slice. The default is 0. A scalar value will increase the slice size equally in all directions, while an array the same shape as im.shape can be passed to pad a specified amount in each direction.

Returns:

A list of slice objects, each indexing into one dimension of the image.

porespy.tools.make_contiguous(im)

Take an image with arbitrary greyscale values and adjust them to ensure all values fall in a contiguous range starting at 0.

Parameters:

im : array_like

An ND array containing greyscale values

Returns:

An ND-array the same size as im but with all values in contiguous

orders.

porespy.tools.randomize_colors(im, keep_vals=[0])

Takes a greyscale image and randomly shuffles the greyscale values, so that all voxels labeled X will be labelled Y, and all voxels labeled Y will be labeled Z, where X, Y, Z and so on are randomly selected from the values in the input image.

This function is useful for improving the visibility of images with neighboring regions that are only incrementally different from each other, such as that returned by scipy.ndimage.label.

Parameters:

im : array_like

An ND image of greyscale values.

keep_vals : array_like

Indicate which voxel values should NOT be altered. The default is [0] which is useful for leaving the background of the image untouched.

Returns:

An image the same size and type as im but with the greyscale values

reassigned. The unique values in both the input and output images will

be identical.

Notes

If the greyscale values in the input image are not contiguous then the neither will they be in the output.

Examples

>>> import porespy as ps
>>> import scipy as sp
>>> sp.random.seed(0)
>>> im = sp.random.randint(low=0, high=5, size=[4, 4])
>>> print(im)
[[4 0 3 3]
 [3 1 3 2]
 [4 0 0 4]
 [2 1 0 1]]
>>> im_rand = ps.tools.randomize_colors(im)
>>> print(im_rand)
[[2 0 4 4]
 [4 1 4 3]
 [2 0 0 2]
 [3 1 0 1]]

As can be seen, the 2’s have become 3, 3’s have become 4, and 4’s have become 2. 1’s remained 1 by random accident. 0’s remain zeros by default, but this can be controlled using the keep_vals argument.

porespy.tools.subdivide(im, divs=2)

Returns slices into an image describing the specified number of sub-arrays. This function is useful for performing operations on smaller images for memory or speed. Note that most typical operations this will NOT work, since the image borders would cause artifacts (e.g. distance_transform)

Parameters:

im : ND-array

The image of the porous media

divs : scalar or array_like

The number of sub-divisions to create in each axis of the image. If a scalar is given it is assume this value applies in all dimensions.

Returns:

An ND-array containing slice objects for indexing into im that extract

the sub-divided arrays.

Notes

This method uses the array_split package which offers the same functionality as the split method of Numpy’s ND-array, but supports the splitting multidimensional arrays in all dimensions.

Examples

>>> import porespy as ps
>>> import matplotlib.pyplot as plt
>>> im = ps.generators.blobs(shape=[200, 200])
>>> s = ps.tools.subdivide(im, divs=[2, 2])

s contains an array with the shape given by divs. To access the first and last quadrants of im use: >>> print(im[s[0, 0]].shape) (100, 100) >>> print(im[s[1, 1]].shape) (100, 100)

It can be easier to index the array with the slices by applying flatten first: >>> s_flat = s.flatten() >>> for i in s_flat: ... print(im[i].shape) (100, 100) (100, 100) (100, 100) (100, 100)

Network Extraction

This module contains functions for extracting pore networks from images.

porespy.network_extraction.partition_pore_space(im, peaks)

Applies a watershed segmentation to partition the distance transform into discrete pores.

Parameters:

im : ND-array

Either the Boolean array of the pore space or a distance tranform. Passing in a pre-existing distance transform saves time, since one is calculated by the function if a Boolean array is passed.

peaks : ND-array

An array the same size as dt indicating where the pore centers are. If boolean, the peaks are labeled; if numeric then it’s assumed that peaks have already been given desired labels.

Returns:

A ND-array the same size as dt with regions belonging to each peak

labelled. The region number is randomized so that neighboring regions

are contrasting colors for easier visualization.

See also

snow

Notes

Find the appropriate peaks array is the tricky part. In principle, each local maxima in the distance transform is a pore center, but due to loss of fidelity in the voxel representation this leads to many extraneous peaks. Several functions are available to help select peaks which can then be passed to this function.

porespy.network_extraction.align_image_with_openpnm(im)

Rotates an image to agree with the coordinates used in OpenPNM. It is unclear why they are not in agreement to start with. This is necessary for overlaying the image and the network in Paraview.

Parameters:

im : ND-array

The image to be rotated. Can the Boolean image of the pore space or any other image of interest.

Returns:

Returns the image rotated accordingly.

porespy.network_extraction.snow(im, r_max=4, sigma=0.4)

This function partitions the void space into pore regions using a marker-based watershed algorithm. The key to this function is that true local maximum of the distance transform are found by trimming various types of extraneous peaks.

The SNOW network extraction algorithm (Sub-Network of an Over-segmented Watershed) was designed to handle to perculiarities of high porosity materials, but it applies well to other materials as well.

Parameters:

im : array_like

Can be either (a) a boolean image of the domain, with True indicating the pore space and False elsewhere, or (b) a distance transform of the domain calculated externally by the user. Option (b) is faster if a distance transform is already available.

r_max : scalar

The radius of there spherical structuring element to use in the Maximum filter stage that is used to find peaks. The default is 4

sigma : scalar

The standard deviation of the Gaussian filter used in step 1. The default is 0.4. If 0 is given then the filter is not applied, which is useful if a distance transform is supplied as the im argument that has already been processed.

Returns:

A tuple containing several arrays. This is a named tuple meaning that

each array can be accessed using the following attribute names:

• im: The binary image of the void space
• dt: The distance transform of the image
• peaks: The peaks of the distance transform after applying the

steps of the SNOW algorithm * regions: The void space partitioned into pores using a marker based watershed with the peaks found by the SNOW algorithm

References

[1] Gostick, J. “A versatile and efficient network extraction algorithm using marker-based watershed segmenation”. Physical Review E. (2017)

porespy.network_extraction.trim_saddle_points(peaks, dt, max_iters=10)

Removes peaks that were mistakenly identified because they lied on a saddle or ridge in the distance transform that was not actually a true local peak.

Parameters:

peaks : ND-array

A boolean image containing True values to mark peaks in the distance transform (dt)

dt : ND-array

The distance transform of the pore space for which the true peaks are sought.

max_iters : int

The maximum number of iterations to run while eroding the saddle points. The default is 10, which is usually not reached; however, a warning is issued if the loop ends prior to removing all saddle points.

Returns:

An image with fewer peaks than was received.

porespy.network_extraction.trim_nearby_peaks(peaks, dt)

Finds pairs of peaks that are nearer to each other than to the solid phase, and removes the peak that is closer to the solid.

Parameters:

peaks : ND-array

A boolean image containing True values to mark peaks in the distance transform (dt)

dt : ND-array

The distance transform of the pore space for which the true peaks are sought.

Returns:

An array the same size as peaks containing a subset of the peaks in

the original image.

Notes

Each pair of peaks is considered simultaneously, so for a triplet of peaks each pair is considered. This ensures that only the single peak that is furthest from the solid is kept. No iteration is required.

porespy.network_extraction.reduce_peaks_to_points(peaks)

Any peaks that are broad or elongated are replaced with a single voxel that is located at the center of mass of the original voxels.

Parameters:

peaks : ND-image

An image containing True values indicating peaks in the distance transform

Returns:

An array with the same number of isolated peaks as the original image, but

fewer total voxels.

Notes

The center of mass of a group of voxels is used as the new single voxel, so if the group has an odd shape (like a horse shoe), the new voxel may not lie on top of the original set.

porespy.network_extraction.find_peaks(dt, r=4, footprint=None)

Returns all local maxima in the distance transform

Parameters:

dt : ND-array

The distance transform of the pore space. This may be calculated and filtered using any means desired.

r : scalar

The size of the structuring element used in the maximum filter. This controls the localness of any maxima. The default is 3 voxels.

footprint : ND-array

Specifies the shape of the structuring element used to define the neighborhood when looking for peaks. If none is specified then a spherical shape is used (or circular in 2D).

Returns:

An ND-array of booleans with True values at the location of any local

maxima.

Notes

It is also possible ot the peak_local_max function from the skimage.feature module as follows:

peaks = peak_local_max(image=dt, min_distance=r, exclude_border=0, indices=False)

This automatically uses a square structuring element which is significantly faster than using a circular or spherical element.

porespy.network_extraction.extract_pore_network(im, dt=None, voxel_size=1)

Analyzes an image that has been partitioned into pore regions and extracts the pore and throat geometry as well as network connectivity.

Parameters:

im : ND-array

An image of the pore space partitioned into individual pore regions. Note that this image must have zeros indicating the solid phase.

dt : ND-array

The distance transform of the pore space. If not given it will be calculated, but it can save time to provide one if available.

voxel_size : scalar

The resolution of the image, expressed as the length of one side of a voxel, so the volume of a voxel would be voxel_size-cubed. The default is 1, which is useful when overlaying the PNM on the original image since the scale of the image is alway 1 unit lenth per voxel.

Returns:

A dictionary containing all the pore and throat size data, as well as the

network topological information. The dictionary names use the OpenPNM

convention (i.e. ‘pore.coords’, ‘throat.conns’) so it may be converted

directly to an OpenPNM network object using the update command.

Simulations

This module contains classes that perform simulations directly on the images to obtain structural information

class porespy.simulations.RandomWalk(im, walkers=1000, max_steps=5000, stride=10, start_frac=0.2)

This class generates objects for performing random walk simulations on.

Parameters:

im: ndarray of bool

2D or 3D image

walkers: int (default = 1000)

The number of walks to initially perform

max_steps: int (default = 5000)

The number of steps to attempt per walk

stride: int (default = 10)

The number of steps taken between each coordinate in path_data

st_frac: float (default = 0.2)

This fraction of the image (in the centre) is searched for a valid starting point for each walk

Examples

Creating a RandomWalk object:

>>> import porespy as ps
>>> im = ps.generators.blobs(300)
>>> rw = ps.simulations.RandomWalk(im, walkers=500, max_steps=3000)

Attributes

max_steps
ndim
path_data
sd_data
shape
start_frac
steps
sterr_data
stride
total_walkers

Methods

find_start_point(start_frac)	Finds a random valid start point in a porous image, searching in the
get_msd([direction])
perform_walks(walkers)	Adds the specified number of walkers to the path_data attribute
save_path(walker, f_name)	Saves the selected walkers path through pore space and free space
show_msd([step_range, direction, showstderr])	Graphs mean squared displacement in pore space and free space vs.
show_path(walker)	Shows the selected walkers path through pore space and free space
walk	This function performs a single random walk through porous image.

find_start_point(start_frac)

Finds a random valid start point in a porous image, searching in the given fraction of the image

Parameters:

start_frac: float

A value between 0 and 1. Determines what fraction of the image is randomly searched for a starting point

Returns:

start_point: tuple

A tuple containing the index of a valid start point. If the image is 2D, start_point will have 0 as the first coord

get_msd(direction=None)

perform_walks(walkers)

Adds the specified number of walkers to the path_data attribute

Parameters:

walkers: int

The number of walks to perform. This many rows will be added to each of the path_data arrays (path_data.pore_space and path_data.free_space)

Notes

This function will continue to use the max_steps and stride value stored on the object

save_path(walker, f_name)

Saves the selected walkers path through pore space and free space

Parameters:

walker: int

The walker whose path will be saved

f_name: string

A path to a filename. Note that this function will add the suffix pore_space and free_space to differentiate between paths

show_msd(step_range=None, direction=None, showstderr=True)

Graphs mean squared displacement in pore space and free space vs. number of steps taken

Parameters:

step_range: array_like (optional)

A sequence with the min and max range of steps to plot. If no argument is given, the range will be zero to max_steps

direction: int (optional)

The axis along which to show mean squared displacement in. If no argument is given, total msd is shown

showstderr: bool (optional)

If set to True, the graph will show error bars

Returns:

out: namedtuple

show_path(walker)

Shows the selected walkers path through pore space and free space using a matplotlib figure

Parameters:

walker: int

The walker whose path will be shown

walk

This function performs a single random walk through porous image. It returns an array containing the walker path in the image, and the walker path in free space.

Parameters:

start_point: array_like

A sequence containing the indices of a valid starting point. If the image is 2D, the first coordinate needs to be zero

max_steps: int

The number of steps to attempt per walk

stride: int

Number of steps taken between each returned coordinate

Returns:

paths: namedtuple

A tuple containing 2 arrays showing the paths taken by the walker: path.pore_space and path.free_space

Notes

This function only performs a single walk through the image. Use RandomWalk.-perform_walks() to add more

class porespy.simulations.Porosimetry(im, voxel_size=1)

Simulates a porosimetry experiment on a binary image. This function is equivalent to the morphological image opening and/or the full morphology approaches.

Parameters:

im : ND-array

The Boolean image of the porous material, with True values indicating the pore space

Attributes

result

Methods

get_drainage_data()	Calculate drainage curve based on the image produced by the porosimetry function.
plot_drainage_curve([fig])
plot_size_histogram([fig])
run([npts, sizes, inlets, access_limited])	Performs a porosimetry simulution on the image

get_drainage_data()

Calculate drainage curve based on the image produced by the porosimetry function.

Returns:

Rp, Snwp: Two arrays containing (a) the radius of the penetrating

sphere (in voxels) and (b) the volume fraction of pore phase voxels

that are accessible from the specfied inlets.

Notes

This function normalizes the invading phase saturation by total pore volume of the dry image, which is assumed to be all voxels with a value equal to 1. To do porosimetry on images with large outer regions, use the `find_outer_region` function then set these regions to 0 in the input image. In future, this function could be adapted to apply this check by default.

plot_drainage_curve(fig=None, **kwargs)

plot_size_histogram(fig=None, **kwargs)

run(npts=25, sizes=None, inlets=None, access_limited=True)

Performs a porosimetry simulution on the image

Parameters:

im : ND-array

An ND image of the porous material containing True values in the pore space.

npts : scalar

The number of invasion points to simulate. Points will be generated spanning the range of sizes in the distance transform. The default is 25 points.

sizes : array_like

The sizes to invade. Use this argument instead of npts for more control of the range and spacing of points.

inlets : ND-array, boolean

A boolean mask with True values indicating where the invasion enters the image. By default all faces are considered inlets, akin to a mercury porosimetry experiment. Users can also apply solid boundaries to their image externally before passing it in, allowing for complex inlets like circular openings, etc.

access_limited : Boolean

This flag indicates if the intrusion should only occur from the surfaces (access_limited is True, which is the default), or if the invading phase should be allowed to appear in the core of the image. The former simulates experimental tools like mercury intrusion porosimetry, while the latter is useful for comparison to gauge the extent of shielding effects in the sample.

Returns:

An ND-image with voxel values indicating the sphere radius at which it

becomes accessible from the inlets. This image can be used to find

invading fluid configurations as a function of applied capillary pressure

by applying a boolean comparison: inv_phase = im > r where r is

the radius (in voxels) of the invading sphere. Of course, r can be

converted to capillary pressure using your favorite model.

Visualization

This module contains functions for quickly visualizing 3D images in 2D views.

porespy.visualization.sem(im, direction='X')

Simulates an SEM photograph looking into the porous material in the specified direction. Features are colored according to their depth into the image, so darker features are further away.

Parameters:

im : array_like

ND-image of the porous material with the solid phase marked as 1 or True

direction : string

Specify the axis along which the camera will point. Options are ‘X’, ‘Y’, and ‘Z’.

Returns:

A 2D greyscale image suitable for use in matplotlib’s `imshow`

function.

porespy.visualization.xray(im, direction='X')

Simulates an X-ray radiograph looking through the porouls material in the specfied direction. The resulting image is colored according to the amount of attenuation an X-ray would experience, so regions with more solid will appear darker.

Parameters:

im : array_like

ND-image of the porous material with the solid phase marked as 1 or True

direction : string

Specify the axis along which the camera will point. Options are ‘X’, ‘Y’, and ‘Z’.

Returns:

A 2D greyscale image suitable for use in matplotlib’s `imshow`

function.

porespy.visualization.set_mpl_style()