Source code for porespy.metrics._funcs

import logging
import numpy as np
import scipy.ndimage as spim
import scipy.spatial as sptl
from scipy import fft as sp_ft
from skimage.measure import regionprops
from deprecated import deprecated
from porespy.tools import extend_slice
from porespy.tools import _check_for_singleton_axes
from porespy.tools import Results
from porespy import settings
from porespy.tools import get_tqdm
from numba import njit


__all__ = [
    "boxcount",
    "chord_counts",
    "chord_length_distribution",
    "find_h",
    "lineal_path_distribution",
    "pore_size_distribution",
    "radial_density_distribution",
    "porosity",
    "porosity_profile",
    "representative_elementary_volume",
    "satn_profile",
    "two_point_correlation",
    "phase_fraction",
    "pc_curve",
    "pc_curve_from_ibip",
    "pc_curve_from_mio",
    "pc_map_to_pc_curve",
]


tqdm = get_tqdm()
logger = logging.getLogger(__name__)


[docs] def boxcount(im, bins=10): r""" Calculates fractal dimension of an image using the tiled box counting method [1]_ Parameters ---------- im : ndarray The image of the porous material. bins : int or array_like, optional The number of box sizes to use. The default is 10 sizes logarithmically spaced between 1 and ``min(im.shape)``. If an array is provided, this is used directly. Returns ------- results An object possessing the following attributes: size : ndarray The box sizes used count : ndarray The number of boxes of each size that contain both solid and void slope : ndarray The gradient of ``count``. This has the same number of elements as ``count`` and References ---------- .. [1] See Boxcounting on `Wikipedia <https://en.wikipedia.org/wiki/Box_counting>`_ Examples -------- `Click here <https://porespy.org/examples/metrics/reference/box_counting.html>`_ to view online example. """ im = np.array(im, dtype=bool) if (len(im.shape) != 2 and len(im.shape) != 3): raise Exception('Image must be 2-dimensional or 3-dimensional') if isinstance(bins, int): Ds = np.unique(np.logspace(1, np.log10(min(im.shape)), bins).astype(int)) else: Ds = np.array(bins).astype(int) N = [] for d in tqdm(Ds, **settings.tqdm): result = 0 for i in range(0, im.shape[0], d): for j in range(0, im.shape[1], d): if len(im.shape) == 2: temp = im[i:i+d, j:j+d] result += np.any(temp) result -= np.all(temp) else: for k in range(0, im.shape[2], d): temp = im[i:i+d, j:j+d, k:k+d] result += np.any(temp) result -= np.all(temp) N.append(result) slope = -1*np.gradient(np.log(np.array(N)), np.log(Ds)) data = Results() data.size = Ds data.count = N data.slope = slope return data
[docs] def representative_elementary_volume(im, npoints=1000): r""" Calculates the porosity of an image as a function subdomain size. This function extracts a specified number of subdomains of random size, then finds their porosity. Parameters ---------- im : ndarray The image of the porous material npoints : int The number of randomly located and sized boxes to sample. The default is 1000. Returns ------- result : Results object A custom object with the following data added as named attributes: 'volume' The total volume of each cubic subdomain tested 'porosity' The porosity of each subdomain tested These attributes can be conveniently plotted by passing the Results object to matplotlib's ``plot`` function using the \* notation: ``plt.plot(\*result, '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. References ---------- .. [1] Bachmat and Bear. On the Concept and Size of a Representative Elementary Volume (Rev), Advances in Transport Phenomena in Porous Media (1987) Examples -------- `Click here <https://porespy.org/examples/metrics/reference/representative_elementary_volume.html>`_ to view online example. """ # TODO: this function is a prime target for parallelization since the # ``npoints`` are calculated independenlty. im_temp = np.zeros_like(im) crds = np.array(np.random.rand(npoints, im.ndim) * im.shape, dtype=int) pads = np.array(np.random.rand(npoints) * np.amin(im.shape) / 2 + 10, dtype=int) im_temp[tuple(crds.T)] = True labels, N = spim.label(input=im_temp) slices = spim.find_objects(input=labels) porosity = np.zeros(shape=(N,), dtype=float) volume = np.zeros(shape=(N,), dtype=int) for i in tqdm(np.arange(0, N), **settings.tqdm): s = slices[i] p = pads[i] new_s = extend_slice(s, shape=im.shape, pad=p) temp = im[new_s] Vp = np.sum(temp, dtype=np.int64) Vt = np.size(temp) porosity[i] = Vp / Vt volume[i] = Vt profile = Results() profile.volume = volume profile.porosity = porosity return profile
[docs] def porosity(im): r""" Calculates the porosity of an image assuming 1's are void space and 0's are solid phase. All other values are ignored, so this can also return the relative fraction of a phase of interest in multiphase images. Parameters ---------- im : ndarray Image of the void space with 1's indicating void phase (or ``True``) and 0's indicating the solid phase (or ``False``). All other values are ignored (see Notes). Returns ------- porosity : float Calculated as the sum of all 1's divided by the sum of all 1's and 0's. See Also -------- phase_fraction find_outer_region 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. Examples -------- `Click here <https://porespy.org/examples/metrics/reference/porosity.html>`_ to view online example. """ im = np.array(im, dtype=np.int64) Vp = np.sum(im == 1, dtype=np.int64) Vs = np.sum(im == 0, dtype=np.int64) e = Vp / (Vs + Vp) return e
[docs] def porosity_profile(im, axis=0): r""" Computes the porosity profile along the specified axis Parameters ---------- im : ndarray The volumetric image for which to calculate the porosity profile. All voxels with a value of 1 (or ``True``) are considered as void. axis : int The axis (0, 1, or 2) along which to calculate the profile. For instance, if `axis` is 0, then the porosity in each YZ plane is calculated and returned as 1D array with 1 value for each X position. Returns ------- result : 1D-array A 1D-array of porosity along the specified axis Examples -------- `Click here <https://porespy.org/examples/metrics/reference/porosity_profile.html>`_ to view online example. """ if axis >= im.ndim: raise Exception('axis out of range') im = np.atleast_3d(im) a = set(range(im.ndim)).difference(set([axis])) a1, a2 = a tmp = np.sum(np.sum(im == 1, axis=a2, dtype=np.int64), axis=a1, dtype=np.int64) prof = tmp / (im.shape[a2] * im.shape[a1]) return prof
[docs] def radial_density_distribution(dt, bins=10, log=False, voxel_size=1): r""" Computes radial density function by analyzing the histogram of voxel values in the distance transform. This function is defined by Torquato [1] as: .. math:: \int_0^\infty P(r)dr = 1.0 where *P(r)dr* is the probability of finding a voxel at a lying at a radial distance between *r* and *dr* from the solid interface. This is equivalent to a probability density function (*pdf*) The cumulative distribution is defined as: .. math:: F(r) = \int_r^\infty P(r)dr which gives the fraction of pore-space with a radius larger than *r*. This is equivalent to the cumulative distribution function (*cdf*). Parameters ---------- dt : ndarray A distance transform of the pore space (the ``edt`` package is recommended). Note that it is recommended to apply ``find_dt_artifacts`` to this image first, and set potentially erroneous values to 0 with ``dt[mask] = 0`` where ``mask = porespy.filters.find_dt_artifaces(dt)``. 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. log : boolean If ``True`` the size data is converted to log (base-10) values before processing. This can help to plot wide size distributions or to better visualize the in the small size region. Note that you should not anti-log the radii values in the retunred ``tuple``, since the binning is performed on the logged radii values. 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 ------- result : Results object A custom object with the following data added as named attributes: ============== ======================================================= Attribute Description ============== ======================================================= *R* or *LogR* Radius, equivalent to ``bin_centers`` *pdf* Probability density function *cdf* Cumulative density function *bin_centers* The center point of each bin *bin_edges* Locations of bin divisions, including 1 more value than the number of bins *bin_widths* Useful for passing to the ``width`` argument of ``matplotlib.pyplot.bar`` ============== ======================================================= Notes ----- Torquato refers to this as the *pore-size density function*, and mentions that it is also known as the *pore-size distribution function*. These terms are avoided here since they have specific connotations in porous media analysis. References ---------- [1] Torquato, S. Random Heterogeneous Materials: Mircostructure and Macroscopic Properties. Springer, New York (2002) - See page 48 & 292 Examples -------- `Click here <https://porespy.org/examples/metrics/reference/radial_density.html>`_ to view online example. """ im = np.copy(dt) x = im[im > 0].flatten() if log: x = np.log10(x) h = np.histogram(x, bins=bins, density=True) h = _parse_histogram(h=h, voxel_size=voxel_size) rdf = Results() rdf[f"{log*'Log' + 'R'}"] = h.bin_centers rdf.pdf = h.pdf rdf.cdf = h.cdf rdf.relfreq = h.relfreq rdf.bin_centers = h.bin_centers rdf.bin_edges = h.bin_edges rdf.bin_widths = h.bin_widths return rdf
[docs] def lineal_path_distribution(im, bins=10, voxel_size=1, log=False): r""" Determines the probability that a point lies within a certain distance of the opposite phase *along a specified direction* This relates directly the radial density function defined by Torquato [1], but instead of reporting the probability of lying within a stated distance to the nearest solid in any direciton, it considers only linear distances along orthogonal directions.The benefit of this is that anisotropy can be detected in materials by performing the analysis in multiple orthogonal directions. Parameters ---------- im : ndarray An image with each voxel containing the distance to the nearest solid along a linear path, as produced by ``distance_transform_lin``. bins : int or array_like The number of bins or a list of specific bins to use voxel_size : scalar The side length of a voxel. This is used to scale the chord lengths into real units. Note this is applied *after* the binning, so ``bins``, if supplied, should be in terms of voxels, not length units. log : boolean If ``True`` (default) the size data is converted to log (base-10) values before processing. This can help to plot wide size distributions or to better visualize data in the small size region. Note that you should not anti-log the radii values in the retunred ``results``, since the binning is performed on the logged radii values. Returns ------- result : Results object A custom object with the following data added as named attributes: *L* or *LogL* Length, equivalent to ``bin_centers`` *pdf* Probability density function *cdf* Cumulative density function *relfreq* Relative frequency chords in each bin. The sum of all bin heights is 1.0. For the cumulative relativce, use *cdf* which is already normalized to 1. *bin_centers* The center point of each bin *bin_edges* Locations of bin divisions, including 1 more value than the number of bins *bin_widths* Useful for passing to the ``width`` argument of ``matplotlib.pyplot.bar`` References ---------- [1] Torquato, S. Random Heterogeneous Materials: Mircostructure and Macroscopic Properties. Springer, New York (2002) Examples -------- `Click here <https://porespy.org/examples/metrics/reference/linearl_path_distribution.html>`_ to view online example. """ x = im[im > 0] if log: x = np.log10(x) h = list(np.histogram(x, bins=bins, density=True)) h = _parse_histogram(h=h, voxel_size=voxel_size) cld = Results() cld[f"{log*'Log' + 'L'}"] = h.bin_centers cld.pdf = h.pdf cld.cdf = h.cdf cld.relfreq = h.relfreq cld.bin_centers = h.bin_centers cld.bin_edges = h.bin_edges cld.bin_widths = h.bin_widths return cld
[docs] def chord_length_distribution(im, bins=10, log=False, voxel_size=1, normalization='count'): r""" Determines the distribution of chord lengths in an image containing chords. Parameters ---------- im : ndarray An image with chords drawn in the pore space, as produced by ``apply_chords`` or ``apply_chords_3d``. ``im`` can be either boolean, in which case each chord will be identified using ``scipy.ndimage.label``, or numerical values in case it is assumed that chords have already been identifed and labeled. In both cases, the size of each chord will be computed as the number of voxels belonging to each labelled region. 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`` (default) the size data is converted to log (base-10) values before processing. This can help to plot wide size distributions or to better visualize the in the small size region. Note that you should not anti-log the radii values in the retunred ``tuple``, since the binning is performed on the logged radii values. normalization : string Indicates how to normalize the bin heights. Options are: *'count' or 'number'* (default) This simply counts the number of chords in each bin in the normal sense of a histogram. This is the rigorous definition according to Torquato [1]. *'length'* This multiplies the number of chords in each bin by the chord length (i.e. bin size). The normalization scheme accounts for the fact that long chords are less frequent than shorert chords, thus giving a more balanced 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 ------- result : Results object A custom object with the following data added as named attributes: *L* or *LogL* Chord length, equivalent to ``bin_centers`` *pdf* Probability density function *cdf* Cumulative density function *relfreq* Relative frequency chords in each bin. The sum of all bin heights is 1.0. For the cumulative relativce, use *cdf* which is already normalized to 1. *bin_centers* The center point of each bin *bin_edges* Locations of bin divisions, including 1 more value than the number of bins *bin_widths* Useful for passing to the ``width`` argument of ``matplotlib.pyplot.bar`` References ---------- [1] Torquato, S. Random Heterogeneous Materials: Mircostructure and Macroscopic Properties. Springer, New York (2002) - See page 45 & 292 Examples -------- `Click here <https://porespy.org/examples/metrics/reference/chord_length_distribution.html>`_ to view online example. """ x = chord_counts(im) if bins is None: bins = np.array(range(0, x.max() + 2)) * voxel_size x = x * voxel_size if log: x = np.log10(x) if normalization == 'length': h = list(np.histogram(x, bins=bins, density=False)) # Scale bin heigths by length h[0] = h[0] * (h[1][1:] + h[1][:-1]) / 2 # Normalize h[0] manually h[0] = h[0] / h[0].sum(dtype=np.int64) / (h[1][1:] - h[1][:-1]) elif normalization in ['number', 'count']: h = np.histogram(x, bins=bins, density=True) else: raise Exception('Unsupported normalization:', normalization) h = _parse_histogram(h) cld = Results() cld[f"{log*'Log' + 'L'}"] = h.bin_centers cld.pdf = h.pdf cld.cdf = h.cdf cld.relfreq = h.relfreq cld.bin_centers = h.bin_centers cld.bin_edges = h.bin_edges cld.bin_widths = h.bin_widths return cld
[docs] def pore_size_distribution(im, bins=10, log=True, voxel_size=1): r""" Calculate a pore-size distribution based on the image produced by the ``porosimetry`` or ``local_thickness`` functions. Parameters ---------- im : ndarray The array of containing the sizes of the largest sphere that overlaps each voxel. Obtained from either ``porosimetry`` or ``local_thickness``. bins : scalar or array_like Either an array of bin sizes to use, or the number of bins that should be automatically generated that span the data range. log : boolean If ``True`` (default) the size data is converted to log (base-10) values before processing. This can help to plot wide size distributions or to better visualize the in the small size region. Note that you should not anti-log the radii values in the retunred ``tuple``, since the binning is performed on the logged radii values. 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 ------- result : Results object A custom object with the following data added as named attributes: *R* or *logR* Radius, equivalent to ``bin_centers`` *pdf* Probability density function *cdf* Cumulative density function *satn* Phase saturation in differential form. For the cumulative saturation, just use *cfd* which is already normalized to 1. *bin_centers* The center point of each bin *bin_edges* Locations of bin divisions, including 1 more value than the number of bins *bin_widths* Useful for passing to the ``width`` argument of ``matplotlib.pyplot.bar`` Notes ----- (1) To ensure the returned values represent actual sizes you can manually scale the input image by the voxel size first (``im *= voxel_size``) plt.bar(psd.R, psd.satn, width=psd.bin_widths, edgecolor='k') Examples -------- `Click here <https://porespy.org/examples/metrics/reference/pore_size_distribution.html>`_ to view online example. """ im = im.flatten() vals = im[im > 0] * voxel_size if log: vals = np.log10(vals) h = _parse_histogram(np.histogram(vals, bins=bins, density=True)) cld = Results() cld[f"{log*'Log' + 'R'}"] = h.bin_centers cld.pdf = h.pdf cld.cdf = h.cdf cld.satn = h.relfreq cld.bin_centers = h.bin_centers cld.bin_edges = h.bin_edges cld.bin_widths = h.bin_widths return cld
def two_point_correlation_bf(im, spacing=10): r""" Calculates the two-point correlation function using brute-force (see Notes) Parameters ---------- im : ndarray 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 ------- result : Results object A custom object with the following data added as named attributes: 'distance' The distance between two points. The distance values are binned as: $$ bins = range(start=0, stop=np.amin(im.shape)/2, stride=spacing) $$ 'probability' The probability that two points of the stated separation distance are within the same phase 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. It is recommended to avoid this. Examples -------- `Click here <https://porespy.org/examples/metrics/reference/two_point_correlation_bf.html>`_ to view online example. """ _check_for_singleton_axes(im) if im.ndim == 2: pts = np.meshgrid(range(0, im.shape[0], spacing), range(0, im.shape[1], spacing)) crds = np.vstack([pts[0].flatten(), pts[1].flatten()]).T elif im.ndim == 3: pts = np.meshgrid(range(0, im.shape[0], spacing), range(0, im.shape[1], spacing), range(0, im.shape[2], spacing)) crds = np.vstack([pts[0].flatten(), pts[1].flatten(), pts[2].flatten()]).T dmat = sptl.distance.cdist(XA=crds, XB=crds) hits = im[tuple(pts)].flatten() dmat = dmat[hits, :] h1 = np.histogram(dmat, bins=range(0, int(np.amin(im.shape) / 2), spacing)) dmat = dmat[:, hits] h2 = np.histogram(dmat, bins=h1[1]) tpcf = Results() tpcf.distance = h2[1][:-1] tpcf.probability = h2[0] / h1[0] return tpcf def _radial_profile(autocorr, bins, pf=None, voxel_size=1): r""" Helper functions to calculate the radial profile of the autocorrelation Masks the image in radial segments from the center and averages the values The distance values are normalized and 100 bins are used as default. Parameters ---------- autocorr : ndarray The image of autocorrelation produced by FFT r_max : int or float The maximum radius in pixels to sum the image over bins : ndarray The edges of the bins to use in summing the radii, ** must be in voxels pf : float the phase fraction (porosity) of the image, used for scaling the normalized autocorrelation down to match the two-point correlation definition as given by Torquato 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 ------- result : tpcf """ if len(autocorr.shape) == 2: adj = np.reshape(autocorr.shape, [2, 1, 1]) # use np.round otherwise with odd image sizes, the mask generated can # be zero, resulting in Div/0 error inds = np.indices(autocorr.shape) - np.round(adj / 2) dt = np.sqrt(inds[0]**2 + inds[1]**2) elif len(autocorr.shape) == 3: adj = np.reshape(autocorr.shape, [3, 1, 1, 1]) # use np.round otherwise with odd image sizes, the mask generated can # be zero, resulting in Div/0 error inds = np.indices(autocorr.shape) - np.round(adj / 2) dt = np.sqrt(inds[0]**2 + inds[1]**2 + inds[2]**2) else: raise Exception('Image dimensions must be 2 or 3') if np.max(bins) > np.max(dt): msg = ( 'Bins specified distances exceeding maximum radial distance for' ' image size. Radial distance cannot exceed distance from center' ' of image to corner.' ) raise Exception(msg) bin_size = bins[1:] - bins[:-1] radial_sum = _get_radial_sum(dt, bins, bin_size, autocorr) # Return normalized bin and radially summed autoc norm_autoc_radial = radial_sum / np.max(autocorr) h = [norm_autoc_radial, bins] h = _parse_histogram(h, voxel_size=1) tpcf = Results() tpcf.distance = h.bin_centers * voxel_size tpcf.bin_centers = h.bin_centers * voxel_size tpcf.bin_edges = h.bin_edges * voxel_size tpcf.bin_widths = h.bin_widths * voxel_size tpcf.probability = norm_autoc_radial tpcf.probability_scaled = norm_autoc_radial * pf tpcf.pdf = h.pdf * pf tpcf.relfreq = h.relfreq return tpcf @njit(parallel=False) def _get_radial_sum(dt, bins, bin_size, autocorr): radial_sum = np.zeros_like(bins[:-1]) for i, r in enumerate(bins[:-1]): mask = (dt <= r) * (dt > (r - bin_size[i])) radial_sum[i] = np.sum(np.ravel(autocorr)[np.ravel(mask)], dtype=np.int64) \ / np.sum(mask) return radial_sum
[docs] def two_point_correlation(im, voxel_size=1, bins=100): r""" Calculate the two-point correlation function using Fourier transforms Parameters ---------- im : ndarray The image of the void space on which the 2-point correlation is desired, in which the phase of interest is labelled as True 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. bins : scalar or array_like Either an array of bin sizes to use, or the number of bins that should be automatically generated that span the data range. The maximum value of the bins, if passed as an array, cannot exceed the distance from the center of the image to the corner. Returns ------- result : tpcf The two-point correlation function object, with named attributes: *distance* The distance between two points, equivalent to bin_centers *bin_centers* The center point of each bin. See distance *bin_edges* Locations of bin divisions, including 1 more value than the number of bins *bin_widths* Useful for passing to the ``width`` argument of ``matplotlib.pyplot.bar`` *probability_normalized* The probability that two points of the stated separation distance are within the same phase normalized to 1 at r = 0 *probability* or *pdf* The probability that two points of the stated separation distance are within the same phase scaled to the phase fraction at r = 0 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 this thesis <https://www.ucl.ac.uk/~ucapikr/projects/KamilaSuankulova_BSc_Project.pdf>`_. Examples -------- `Click here <https://porespy.org/examples/metrics/reference/two_point_correlation.html>`_ to view online example. """ # Get the number of CPUs available to parallel process Fourier transforms cpus = settings.ncores # Get the phase fraction of the image pf = porosity(im) if isinstance(bins, int): # Calculate half lengths of the image r_max = (np.ceil(np.min(np.shape(im))) / 2).astype(int) # Get the bin-size - ensures it will be at least 1 bin_size = int(np.ceil(r_max / bins)) # Calculate the bin divisions, equivalent to bin_edges bins = np.arange(0, r_max + bin_size, bin_size) # set the number of parallel processors to use: with sp_ft.set_workers(cpus): # Fourier Transform and shift image F = sp_ft.ifftshift(sp_ft.rfftn(sp_ft.fftshift(im))) # Compute Power Spectrum P = np.absolute(F**2) # Auto-correlation is inverse of Power Spectrum autoc = np.absolute(sp_ft.ifftshift(sp_ft.irfftn(sp_ft.fftshift(P)))) tpcf = _radial_profile(autoc, bins, pf=pf, voxel_size=voxel_size) return tpcf
def _parse_histogram(h, voxel_size=1, density=True): delta_x = h[1] P = h[0] bin_widths = delta_x[1:] - delta_x[:-1] temp = P * (bin_widths) C = np.cumsum(temp[-1::-1])[-1::-1] S = P * (bin_widths) if not density: P /= np.max(P) temp_sum = np.sum(P * bin_widths) C /= temp_sum S /= temp_sum bin_edges = delta_x * voxel_size bin_widths = (bin_widths) * voxel_size bin_centers = ((delta_x[1:] + delta_x[:-1]) / 2) * voxel_size hist = Results() hist.pdf = P hist.cdf = C hist.relfreq = S hist.bin_centers = bin_centers hist.bin_edges = bin_edges hist.bin_widths = bin_widths return hist
[docs] def chord_counts(im): r""" Find the length of each chord in the supplied image Parameters ---------- im : ndarray An image containing chords drawn in the void space. Returns ------- result : 1D-array 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 ``np.histogram`` to get the histogram data directly. Another useful function is ``np.bincount`` which gives the number of chords of each length in a format suitable for ``plt.plot``. Examples -------- `Click here <https://porespy.org/examples/reference/metrics/chord_counts.html>`_ to view online example. """ labels, N = spim.label(im > 0) props = regionprops(labels) chord_lens = np.array([i.filled_area for i in props], dtype=int) return chord_lens
[docs] def phase_fraction(im, normed=True): r""" Calculate the fraction of each phase in an image Parameters ---------- im : ndarray An ndarray containing integer values normed : boolean If ``True`` (default) the returned values are normalized by the total number of voxels in image, otherwise the voxel count of each phase is returned. Returns ------- result : 1D-array A array of length max(im) with each element containing the number of voxels found with the corresponding label. See Also -------- porosity Examples -------- `Click here <https://porespy.org/examples/metrics/reference/phase_fraction.html>`_ to view online example. """ if im.dtype == bool: im = im.astype(int) labels = np.unique(im) results = {} for label in labels: results[label] = np.sum(im == label, dtype=np.int64) * \ (1 / im.size if normed else 1) return results
[docs] @deprecated("This function is deprecated, use pc_curve instead") def pc_curve_from_ibip(*args, **kwargs): r""" This function is deprecated. Use ``pc_curve`` instead. Note that the ``stepped`` argument is no longer supported since this can be done directly in matplotlib with ``plt.step(...)``. """ return pc_curve(*args, **kwargs)
[docs] @deprecated("This function is deprecated, use pc_curve instead") def pc_curve_from_mio(*args, **kwargs): r""" This function is deprecated. Use ``pc_curve`` instead. """ return pc_curve(*args, **kwargs)
[docs] def pc_curve(im, sizes=None, pc=None, seq=None, sigma=0.072, theta=180, voxel_size=1): r""" Produces a Pc-Snwp curve given a map of meniscus radii or capillary pressures at which each voxel was invaded Parameters ---------- im : ndarray The voxel image of the porous media with ``True`` values indicating the void space sizes : ndarray, optional An image containing the sphere radii at which each voxel was invaded during an invasion experiment. pc : ndarray, optional An image containing the capillary pressures at which each voxel was invaded during an invasion experiment. seq : ndarray, optional An image containing invasion sequence values, such as that returned from the ``ibip`` function. sigma : float, optional The surface tension of the fluid-fluid system of interest. This argument is ignored if ``pc`` are specified, otherwise it is used in the Washburn equation to convert ``sizes`` to capillary pc. theta : float The contact angle measured through the invading phase in degrees. This argument is ignored if ``pc`` are specified, otherwise it is used in the Washburn equation to convert ``sizes`` to capillary pressures. voxel_size : float The voxel resolution of the image. This argument is ignored if ``pc`` are specified, otherwise it is used in the Washburn equation to convert ``sizes`` to capillary pressures. Returns ------- pc_curve : Results object A custom object with the following data added as named attributes: ================== =================================================== Attribute Description ================== =================================================== pc The capillary pressure, either as given in ``pc`` or computed from ``sizes`` (see Notes). snwp The fraction of void space filled by non-wetting phase at each pressure in ``pc`` ================== =================================================== Notes ----- If ``sizes`` is provided, then the Washburn equation is used to convert the radii to capillary pressures, using the given ``sigma`` and ``theta`` values, along with the ``voxel_size`` if the values are in voxel radii. For more control over how capillary pressure model, it can be computed by hand, for example: $$ pc = \frac{-2*0.072*np.cos(np.deg2rad(180))}{sizes \cdot voxel_size} $$ then passed in as the ``pc`` argument. Examples -------- `Click here <https://porespy.org/examples/metrics/reference/pc_curve.html>`_ to view online example. """ tqdm = get_tqdm() if seq is not None: seqs = np.unique(seq)[1:] x = [] y = [] with tqdm(seqs, **settings.tqdm) as pbar: for n in seqs: pbar.update() mask = seq == n if (pc is not None) and (sizes is not None): raise Exception("Only one of pc or sizes can be specified") elif pc is not None: pressure = pc[mask][0] elif sizes is not None: r = sizes[mask][0]*voxel_size pressure = -2*sigma*np.cos(np.deg2rad(theta))/r x.append(pressure) snwp = ((seq <= n)*(seq > 0) * (im == 1)).sum(dtype=np.int64)/im.sum(dtype=np.int64) y.append(snwp) pc_curve = Results() pc_curve.pc = x pc_curve.snwp = y elif sizes is not None: if im is None: im = ~(sizes == 0) sz = np.unique(sizes)[:0:-1] sz = np.hstack((sz[0]*2, sz)) x = [] y = [] with tqdm(sz, **settings.tqdm) as pbar: for n in sz: pbar.update() r = n*voxel_size pc = -2*sigma*np.cos(np.deg2rad(theta))/r x.append(pc) snwp = ((sizes >= n)*(im == 1)).sum(dtype=np.int64)/im.sum(dtype=np.int64) y.append(snwp) pc_curve = Results() pc_curve.pc = x pc_curve.snwp = y elif pc is not None: Ps = np.unique(pc[im]) # Utilize the fact that -inf and +inf will be at locations 0 & -1 in Ps if Ps[-1] == np.inf: Ps[-1] = Ps[-2]*2 if Ps[0] == -np.inf: Ps[0] = Ps[1] - np.abs(Ps[1]/2) else: # Add a point at begining to ensure curve starts a 0, if no residual Ps = np.hstack((Ps[0] - np.abs(Ps[0]/2), Ps)) y = [] Vp = im.sum(dtype=np.int64) temp = pc[im] for p in tqdm(Ps, **settings.tqdm): y.append((temp <= p).sum(dtype=np.int64)/Vp) pc_curve = Results() pc_curve.pc = Ps pc_curve.snwp = y return pc_curve
def pc_map_to_pc_curve(pc, im, seq=None): r""" Converts a pc map into a capillary pressure curve Parameters ---------- pc : ndarray A numpy array with each voxel containing the capillary pressure at which it was invaded. `-inf` indicates voxels which are already filled with non-wetting fluid, and `+inf` indicates voxels that are not invaded by non-wetting fluid (e.g., trapped wetting phase). Solids should be noted by `+inf` but this is also enforced inside the function using `im`. im : ndarray A numpy array with `True` values indicating the void space and `False` elsewhere. This is necessary to define the total void volume of the domain for computing the saturation. seq : ndarray, optional A numpy array with each voxel containing the sequence at which it was invaded. This is required when analyzing results from invasion percolation since the pressures in `pc` do not correspond to the sequence in which they were filled. Returns ------- results : dataclass-like A dataclass like object with the following attributes: ================== ========================================================= Attribute Description ================== ========================================================= pc The capillary pressure snwp The fraction of void space filled by non-wetting phase at each pressure in ``pc`` ================== ========================================================= Notes ----- To use this function with the results of `porosimetry` or `ibip` the sizes map must be converted to a capillary pressure map first. `drainage` and `invasion` both return capillary pressure maps which can be passed directly as `pc`. """ pc[~im] = np.inf # Ensure solid voxels are set to inf invasion pressure if seq is None: pcs, counts = np.unique(pc, return_counts=True) else: vals, index, counts = np.unique(seq, return_index=True, return_counts=True) pcs = pc.flatten()[index] snwp = np.cumsum(counts[pcs < np.inf])/im.sum() pcs = pcs[pcs < np.inf] results = Results() results.pc = pcs results.snwp = snwp return results
[docs] def satn_profile(satn, s=None, im=None, axis=0, span=10, mode='tile'): r""" Computes a saturation profile from an image of fluid invasion Parameters ---------- satn : ndarray An image with each voxel indicating the saturation upon its invasion. 0's are treated as solid and -1's are treated as uninvaded void space. s : scalar The global saturation value for which the profile is desired. If `satn` is a pre-thresholded boolean image then this is ignored, `im` is required. im : ndarray A boolean image with `True` values indicating the void phase. This is used to compute the void volume if `satn` is given as a pre-thresholded boolean mask. axis : int The axis along which to profile should be measured span : int The number of layers to include in the moving average saturation calculation. mode : str How the moving average should be applied. Options are: ======== ============================================================== mode description ======== ============================================================== 'tile' The average is computed for discrete non-overlapping tiles of a size given by ``span`` 'slide' The average is computed in a moving window starting at ``span/2`` and sliding by a single voxel. This method provides more data points but is slower. ======== ============================================================== Returns ------- results : dataclass Results is a custom porespy class with the following attributes: ============= ========================================================= Attribute Description ============= ========================================================= position The position along the given axis at which saturation values are computed. The units are in voxels. saturation The local saturation value at each position. ============= ========================================================= Examples -------- `Click here <https://porespy.org/examples/metrics/reference/satn_profile.html>`_ to view online example. """ span = max(1, span) if s is None: if satn.dtype != bool: msg = 'Must specify a target saturation if saturation map is provided' raise Exception(msg) s = 2 # Will find ALL voxels, then > 0 will limit to only True ones satn = satn.astype(int) satn[satn == 0] = -1 satn[~im] = 0 else: msg = 'The maximum saturation in the image is less than the given threshold' if satn.max() < s: raise Exception(msg) satn = np.swapaxes(satn, 0, axis) if mode == 'tile': y = np.zeros(int(satn.shape[0]/span)) z = np.zeros_like(y) for i in range(int(satn.shape[0]/span)): void = satn[i*span:(i+1)*span, ...] != 0 nwp = (satn[i*span:(i+1)*span, ...] <= s) \ *(satn[i*span:(i+1)*span, ...] > 0) y[i] = nwp.sum(dtype=np.int64)/void.sum(dtype=np.int64) z[i] = i*span + (span-1)/2 if mode == 'slide': y = np.zeros(int(satn.shape[0]-span)) z = np.zeros_like(y) for i in range(int(satn.shape[0]-span)): void = satn[i:i+span, ...] != 0 nwp = (satn[i:i+span, ...] <= s)*(satn[i:i+span, ...] > 0) y[i] = nwp.sum(dtype=np.int64)/void.sum(dtype=np.int64) z[i] = i + (span-1)/2 results = Results() results.position = z results.saturation = y return results
[docs] def find_h(saturation, position=None, srange=[0.01, 0.99]): r""" Given a saturation profile, compute the height between given bounds Parameters ---------- saturation : array_like A list of saturation values as function of ``position`` position : array_like, optional A list of positions corresponding to each saturation. If not provided then each value in ``saturation`` is assumed to be separated by 1 voxel. srange : list The minimum and maximum value of saturation to consider as the start and end of the profile Returns ------- h : scalar The height of the two-phase zone See Also -------- satn_profile Notes ----- The ``satn_profile`` function can be used to obtain the ``saturation`` and ``position`` from an image. Examples -------- `Click here <https://porespy.org/examples/metrics/reference/find_h.html>`_ to view online example. """ r = Results() r.valid = True # First ensure saturation generally descends from left to right if np.mean(saturation[:10]) < np.mean(saturation[-10:]): saturation = np.flip(saturation, axis=0) # Ensure requested saturation limits actually exist if (min(srange) < min(saturation)) or (max(srange) > max(saturation)): srange = max(min(srange), min(saturation)), min(max(srange), max(saturation)) r.valid = False logger.warning(f'The requested saturation range was adjusted to {srange}' ' to accomodate data') # Find zmax x = saturation >= max(srange) zmax = np.where(x)[0][-1] y = saturation <= min(srange) zmin = np.where(y)[0][0] # If position array was given, index into it if position is not None: zmax = position[zmax] zmin = position[zmin] # Add remaining data to results object r.zmax = zmax r.zmin = zmin r.smax = max(srange) r.smin = min(srange) r.h = abs(zmax-zmin) return r