Source code for porespy.generators._imgen

import numpy as np
import inspect as insp
from edt import edt
import porespy as ps
from numba import njit
import scipy.spatial as sptl
import scipy.ndimage as spim
from porespy.tools import norm_to_uniform, ps_ball, ps_disk, get_border
from porespy import settings
from typing import List
from loguru import logger
tqdm = ps.tools.get_tqdm()


[docs]def insert_shape(im, element, center=None, corner=None, value=1, mode="overwrite"): r""" Inserts sub-image into a larger image at the specified location. If the inserted image extends beyond the boundaries of the image it will be cropped accordingly. Parameters ---------- im : ndarray The image into which the sub-image will be inserted element : ndarray The sub-image to insert center : tuple Coordinates indicating the position in the main image where the inserted imaged will be centered. If ``center`` is given then ``corner`` cannot be specified. Note that ``center`` can only be used if all dimensions of ``element`` are odd, otherwise the meaning of center is not defined. corner : tuple Coordinates indicating the position in the main image where the lower corner (i.e. [0, 0, 0]) of the inserted image should be anchored. If ``center`` is given then ``corner`` cannot be specified. value : scalar A scalar value to apply to the sub-image. The default is 1. mode : string If 'overwrite' (default) the inserted image replaces the values in the main image. If 'overlay' the inserted image is added to the main image. In both cases the inserted image is multiplied by ``value`` first. Returns ------- im : ndarray A copy of ``im`` with the supplied element inserted. Examples -------- `Click here <https://porespy.org/examples/generators/howtos/insert_shape.html>`_ to view online example. """ im = im.copy() if im.ndim != element.ndim: msg = f"Image shape {im.shape} and element shape {element.shape} do not match" raise Exception(msg) s_im = [] s_el = [] if (center is not None) and (corner is None): for dim in range(im.ndim): r, d = np.divmod(element.shape[dim], 2) if d == 0: raise Exception( "Cannot specify center point when element " + "has one or more even dimension" ) lower_im = np.amax((center[dim] - r, 0)) upper_im = np.amin((center[dim] + r + 1, im.shape[dim])) s_im.append(slice(lower_im, upper_im)) lower_el = np.amax((lower_im - center[dim] + r, 0)) upper_el = np.amin((upper_im - center[dim] + r, element.shape[dim])) s_el.append(slice(lower_el, upper_el)) elif (corner is not None) and (center is None): for dim in range(im.ndim): L = int(element.shape[dim]) lower_im = np.amax((corner[dim], 0)) upper_im = np.amin((corner[dim] + L, im.shape[dim])) s_im.append(slice(lower_im, upper_im)) lower_el = np.amax((lower_im - corner[dim], 0)) upper_el = np.amin((upper_im - corner[dim], element.shape[dim])) s_el.append(slice(min(lower_el, upper_el), upper_el)) else: raise Exception("Cannot specify both corner and center") if mode == "overlay": im[tuple(s_im)] = im[tuple(s_im)] + element[tuple(s_el)] * value elif mode == "overwrite": im[tuple(s_im)] = element[tuple(s_el)] * value else: raise Exception("Invalid mode " + mode) return im
[docs]def RSA(im_or_shape: np.array, r: int, volume_fraction: int = 1, clearance: int = 0, n_max: int = 100000, mode: str = "contained", return_spheres: bool = False, smooth: bool = True): r""" Generates a sphere or disk packing using Random Sequential Addition Parameters ---------- im_or_shape : ndarray or list To provide flexibility, this argument accepts either an image into which the spheres are inserted, or a shape which is used to create an empty image. In both cases the spheres are added as ``True`` values to the background. Since ``True`` is considered the pore space, then the added spheres represent holes. r : int The radius of the disk or sphere to insert. volume_fraction : scalar (default is 1.0) The fraction of the image that should be filled with spheres. The spheres are added as ``True``'s, so each sphere addition increases the ``volume_fraction`` until the specified limit is reached. Note that if ``n_max`` is reached first, then ``volume_fraction`` will not be acheived. Also, ``volume_fraction`` is not counted correctly if the ``mode`` is ``'extended'``. clearance : int (optional, default = 0) The amount of space to put between each sphere. Negative values are acceptable to create overlaps, so long as ``abs(clearance) < r``. n_max : int (default is 100,000) The maximum number of spheres to add. Using a low value may halt the addition process prior to reaching the specified ``volume_fraction``. If ``None`` is given, then no limit is used. mode : string (default is 'contained') Controls how the edges of the image are handled. Options are: 'contained' Spheres are all completely within the image 'extended' Spheres are allowed to extend beyond the edge of the image. In this mode the volume fraction will be less that requested since some spheres extend beyond the image, but their entire volume is counted as added for computational efficiency. return_spheres : bool If ``True`` then an image containing only the spheres is returned rather than the input image with the spheres added, which is the default behavior. smooth : bool Indicates whether balls should have smooth faces (``True``) or should include the bumps on the extremities (``False``). Returns ------- image : ndarray An image with spheres of specified radius *added* to the background. See Also -------- pseudo_gravity_packing pseudo_electrostatic_packing Notes ----- This algorithm ensures that spheres do not overlap but does not guarantee they are tightly packed. This function adds spheres to the background of the received ``im``, which allows iteratively adding spheres of different radii to the unfilled space by repeatedly passing in the result of previous calls to RSA. References ---------- [1] Random Heterogeneous Materials, S. Torquato (2001) Examples -------- `Click here <https://porespy.org/examples/generators/howtos/RSA.html>`_ to view online example. """ logger.debug(f"RSA: Adding spheres of size {r}") if len(im_or_shape) <= 3: im = np.zeros(shape=im_or_shape, dtype=bool) else: im = im_or_shape im = im.astype(bool) if return_spheres: im_temp = np.copy(im) if n_max is None: n_max = np.inf vf_final = volume_fraction vf_start = im.sum() / im.size logger.debug(f"Initial volume fraction: {vf_start}") if im.ndim == 2: template_lg = ps_disk((r + clearance) * 2) template_sm = ps_disk(r, smooth=smooth) else: template_lg = ps_ball((r + clearance) * 2) template_sm = ps_ball(r, smooth=smooth) vf_template = template_sm.sum() / im.size # Pad image by the radius of large template to enable insertion near edges im = np.pad(im, pad_width=2 * (r + clearance), mode="edge") # Depending on mode, adjust mask to remove options around edge if mode == "contained": border = get_border(im.shape, thickness=2 * (r + clearance), mode="faces") elif mode == "extended": border = get_border(im.shape, thickness=(r + clearance) + 1, mode="faces") else: raise Exception("Unrecognized mode: ", mode) # Remove border pixels im[border] = True # Dilate existing objects by strel to remove pixels near them # from consideration for sphere placement logger.trace("Dilating foreground features by sphere radius") dt = edt(im == 0) options_im = dt >= r # Begin inserting the spheres vf = vf_start free_sites = np.flatnonzero(options_im) i = 0 while (vf <= vf_final) and (i < n_max) and (len(free_sites) > 0): c, count = _make_choice(options_im, free_sites=free_sites) # The 100 below is arbitrary and may change performance if count > 100: # Regenerate list of free_sites logger.debug(f"Regenerating `free_sites` after {i} iterations") free_sites = np.flatnonzero(options_im) if all(np.array(c) == -1): break s_sm = tuple([slice(x - r, x + r + 1, None) for x in c]) s_lg = tuple([slice(x - 2 * (r + clearance), x + 2 * (r + clearance) + 1, None) for x in c]) im[s_sm] += template_sm # Add ball to image options_im[s_lg][template_lg] = False # Update extended region vf += vf_template i += 1 logger.trace(f"Number of spheres inserted: {i}") # Get slice into returned image to retain original size s = tuple([slice(2 * (r + clearance), d - 2 * (r + clearance), None) for d in im.shape]) im = im[s] vf = im.sum() / im.size logger.debug("Final volume fraction:", vf) if return_spheres: im = im * (~im_temp) return im
@njit def _make_choice(options_im, free_sites): r""" This function is called by _begin_inserting to find valid insertion points. Parameters ---------- options_im : ndarray An array with ``True`` at all valid locations and ``False`` at all locations where a sphere already exists PLUS a region of radius R around each sphere since these points are also invalid insertion points. free_sites : array_like A 1D array containing valid insertion indices. This list is used to select insertion points from a limited which occasionally gets smaller. Returns ------- coords : list The XY or XYZ coordinates of the next insertion point count : int The number of attempts that were needed to find valid point. If this value gets too high, a short list of ``free_sites`` should be generated in the calling function. """ choice = False count = 0 upper_limit = len(free_sites) maxiter = upper_limit * 20 if options_im.ndim == 2: coords = [-1, -1] Nx, Ny = options_im.shape while not choice: if count >= maxiter: coords = [-1, -1] break ind = np.random.randint(0, upper_limit) # This numpy function is not supported by numba yet # c1, c2 = np.unravel_index(free_sites[ind], options_im.shape) # So using manual unraveling coords[1] = free_sites[ind] % Ny coords[0] = (free_sites[ind] // Ny) % Nx choice = options_im[coords[0], coords[1]] count += 1 if options_im.ndim == 3: coords = [-1, -1, -1] Nx, Ny, Nz = options_im.shape while not choice: if count >= maxiter: coords = [-1, -1, -1] break ind = np.random.randint(0, upper_limit) # This numpy function is not supported by numba yet # c1, c2, c3 = np.unravel_index(free_sites[ind], options_im.shape) # So using manual unraveling coords[2] = free_sites[ind] % Nz coords[1] = (free_sites[ind] // Nz) % Ny coords[0] = (free_sites[ind] // (Nz * Ny)) % Nx choice = options_im[coords[0], coords[1], coords[2]] count += 1 return coords, count
[docs]def bundle_of_tubes(shape: List[int], spacing: int): r""" 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 : int The center to center distance of the holes. The hole sizes will be randomly distributed between this values down to 3 voxels. Returns ------- image : ndarray A boolean array with ``True`` values denoting the pore space Examples -------- `Click here <https://porespy.org/examples/generators/howtos/bundle_of_tubes.html>`_ to view online example. """ shape = np.array(shape) if np.size(shape) == 1: shape = np.full((3,), int(shape)) if np.size(shape) == 2: shape = np.hstack((shape, [1])) temp = np.zeros(shape=shape[:2]) Xi = np.ceil( np.linspace(spacing / 2, shape[0] - (spacing / 2) - 1, int(shape[0] / spacing)) ) Xi = np.array(Xi, dtype=int) Yi = np.ceil( np.linspace(spacing / 2, shape[1] - (spacing / 2) - 1, int(shape[1] / spacing)) ) Yi = np.array(Yi, dtype=int) temp[tuple(np.meshgrid(Xi, Yi))] = 1 inds = np.where(temp) for i in range(len(inds[0])): r = np.random.randint(1, (spacing / 2)) try: s1 = slice(inds[0][i] - r, inds[0][i] + r + 1) s2 = slice(inds[1][i] - r, inds[1][i] + r + 1) temp[s1, s2] = ps_disk(r) except ValueError: odd_shape = np.shape(temp[s1, s2]) temp[s1, s2] = ps_disk(r)[: odd_shape[0], : odd_shape[1]] im = np.broadcast_to(array=np.atleast_3d(temp), shape=shape) return im
[docs]def polydisperse_spheres( shape: List[int], porosity: float, dist, nbins: int = 5, r_min: int = 5 ): r""" Create an image of randomly placed, overlapping 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 : float 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 chosen from 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 obtained with ``dist = scipy.stats.norm(loc=20, scale=10)`` nbins : int 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 ------- image : ndarray A boolean array with ``True`` values denoting the pore space Examples -------- `Click here <https://porespy.org/examples/generators/howtos/polydisperse_spheres.html>`_ to view online example. """ shape = np.array(shape) if np.size(shape) == 1: shape = np.full((3,), int(shape)) Rs = dist.interval(np.linspace(0.05, 0.95, nbins)) Rs = np.vstack(Rs).T Rs = (Rs[:-1] + Rs[1:]) / 2 Rs = np.clip(Rs.flatten(), a_min=r_min, a_max=None) phi_desired = 1 - (1 - porosity) / (len(Rs)) im = np.ones(shape, dtype=bool) for r in Rs: phi_im = im.sum() / np.prod(shape) phi_corrected = 1 - (1 - phi_desired) / phi_im temp = overlapping_spheres(shape=shape, r=r, porosity=phi_corrected) im = im * temp return im
[docs]def voronoi_edges(shape: List[int], r: int, ncells: int, flat_faces: bool = True): r""" Create an image from the edges of 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 : int The radius to which Voronoi edges should be dilated in the final image. ncells : int The number of Voronoi cells to include in the tesselation. flat_faces : bool Whether the Voronoi edges should lie on the boundary of the image (``True``), or if edges outside the image should be removed (``False``). Returns ------- image : ndarray A boolean array with ``True`` values denoting the pore space Examples -------- `Click here <https://porespy.org/examples/generators/howtos/voronoi_edges.html>`_ to view online example. """ logger.trace(f"Generating {ncells} cells") shape = np.array(shape) if np.size(shape) == 1: shape = np.full((3,), int(shape)) im = np.zeros(shape, dtype=bool) base_pts = np.random.rand(ncells, 3) * shape if flat_faces: # Reflect base points Nx, Ny, Nz = shape orig_pts = base_pts base_pts = np.vstack((base_pts, [-1, 1, 1] * orig_pts + [2.0 * Nx, 0, 0])) base_pts = np.vstack((base_pts, [1, -1, 1] * orig_pts + [0, 2.0 * Ny, 0])) base_pts = np.vstack((base_pts, [1, 1, -1] * orig_pts + [0, 0, 2.0 * Nz])) base_pts = np.vstack((base_pts, [-1, 1, 1] * orig_pts)) base_pts = np.vstack((base_pts, [1, -1, 1] * orig_pts)) base_pts = np.vstack((base_pts, [1, 1, -1] * orig_pts)) vor = sptl.Voronoi(points=base_pts) vor.vertices = np.around(vor.vertices) vor.vertices *= (np.array(im.shape) - 1) / np.array(im.shape) vor.edges = _get_Voronoi_edges(vor) for row in vor.edges: pts = vor.vertices[row].astype(int) if np.all(pts >= 0) and np.all(pts < im.shape): line_pts = line_segment(pts[0], pts[1]) im[tuple(line_pts)] = True im = edt(~im) > r return im
def _get_Voronoi_edges(vor): r""" Given a Voronoi object as produced by the scipy.spatial.Voronoi class, this function calculates the start and end points of eeach edge in the Voronoi diagram, in terms of the vertex indices used by the received Voronoi object. Parameters ---------- vor : scipy.spatial.Voronoi object Returns ------- A 2-by-N array of vertex indices, indicating the start and end points of each vertex in the Voronoi diagram. These vertex indices can be used to index straight into the ``vor.vertices`` array to get spatial positions. """ edges = [[], []] for facet in vor.ridge_vertices: # Create a closed cycle of vertices that define the facet edges[0].extend(facet[:-1] + [facet[-1]]) edges[1].extend(facet[1:] + [facet[0]]) edges = np.vstack(edges).T # Convert to scipy-friendly format mask = np.any(edges == -1, axis=1) # Identify edges at infinity edges = edges[~mask] # Remove edges at infinity edges = np.sort(edges, axis=1) # Move all points to upper triangle # Remove duplicate pairs edges = edges[:, 0] + 1j * edges[:, 1] # Convert to imaginary edges = np.unique(edges) # Remove duplicates edges = np.vstack((np.real(edges), np.imag(edges))).T # Back to real edges = np.array(edges, dtype=int) return edges
[docs]def lattice_spheres(shape: List[int], r: int, spacing: int = None, offset: int = None, smooth: bool = True, lattice: str = "sc"): r""" Generate a cubic packing of spheres in a specified lattice arrangement. Parameters ---------- shape : list The size of the image to generate in [Nx, Ny, Nz] where N is the number of voxels in each direction. For a 2D image, use [Nx, Ny]. radius : int The radius of spheres (circles) in the packing. spacing : int or List[int] The spacing between unit cells. If the spacing is too small then spheres may overlap. If an ``int`` is given it will be applied in all directions, while a list of ``int`` will be interpreted to apply along each axis. offset : int or List[int] The amount offset to add between sphere centers and the edges of the image. A single ``int`` will be applied in all directions, while a list of ``int`` will be interpreted to apply along each axis. smooth : bool, default=True If ``True`` (default) the outer extremities of the sphere will not have the little bumps on each face. lattice : str Specifies the type of lattice to create. Options are: 'sc' Simple Cubic (default) 'fcc' Face Centered Cubic 'bcc' Body Centered Cubic Returns ------- image : ndarray A boolean array with ``True`` values denoting the pore space. Notes ----- For 2D images, 'sc' gives a square lattice and both 'fcc' and 'bcc' give a triangular lattice. Examples -------- `Click here <https://porespy.org/examples/generators/howtos/lattice_spheres.html>`_ to view online example. """ logger.debug(f"Generating {lattice} lattice") shape = np.array(shape) im = np.zeros(shape, dtype=bool) # Parse lattice type lattice = lattice.lower() if im.ndim == 2: if lattice in ['sc', 'square', 'cubic', 'simple cubic']: lattice = 'sq' elif lattice in ['tri', 'triangular']: lattice = 'tri' else: raise Exception(f'Unrecognized mode: {lattice}') else: if lattice in ['sc', 'cubic', 'simple cubic']: lattice = 'sc' elif lattice in ['bcc', 'body centered cubic']: lattice = 'bcc' elif lattice in ['fcc', 'face centered cubic']: lattice = 'fcc' else: raise Exception(f'Unrecognized mode: {lattice}') # Parse offset and spacing args if spacing is None: spacing = 2*r if isinstance(spacing, int): spacing = [spacing]*im.ndim if offset is None: offset = r if isinstance(offset, int): offset = [offset]*im.ndim if lattice == 'sq': im[offset[0]::spacing[0], offset[1]::spacing[1]] = True elif lattice == 'tri': im[offset[0]::spacing[0], offset[1]::spacing[1]] = True im[offset[0]+int(spacing[0]/2)::spacing[0], offset[1]+int(spacing[1]/2)::spacing[1]] = True elif lattice == 'sc': im[offset[0]::spacing[0], offset[1]::spacing[1], offset[2]::spacing[2]] = True elif lattice == 'bcc': im[offset[0]::spacing[0], offset[1]::spacing[1], offset[2]::spacing[2]] = True im[offset[0]+int(spacing[0]/2)::spacing[0], offset[1]+int(spacing[1]/2)::spacing[1], offset[2]+int(spacing[2]/2)::spacing[2]] = True elif lattice == 'fcc': im[offset[0]::spacing[0], offset[1]::spacing[1], offset[2]::spacing[2]] = True # xy-plane im[offset[0]+int(spacing[0]/2)::spacing[0], offset[1]+int(spacing[1]/2)::spacing[1], offset[2]::spacing[2]] = True # xz-plane im[offset[0]+int(spacing[0]/2)::spacing[0], offset[1]::spacing[1], offset[2]+int(spacing[2]/2)::spacing[2]] = True # yz-plane im[offset[0]::spacing[0], offset[1]+int(spacing[1]/2)::spacing[1], offset[2]+int(spacing[2]/2)::spacing[2]] = True if smooth: im = ~(edt(~im) < r) else: im = ~(edt(~im) <= r) return im
[docs]def overlapping_spheres(shape: List[int], r: int, porosity: float, maxiter: int = 10, tol: float = 0.01): r""" 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. r : scalar The radius of spheres in the packing. porosity : scalar The porosity of the final image, accurate to the given tolerance. maxiter : int Maximum number of iterations for the iterative algorithm that improves the porosity of the final image to match the given value. tol : float Tolerance for porosity of the final image compared to the given value. Returns ------- image : ndarray 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. Examples -------- `Click here <https://porespy.org/examples/generators/howtos/overlapping_spheres.html>`_ to view online example. """ shape = np.array(shape) if np.size(shape) == 1: shape = np.full((3, ), int(shape)) ndim = (shape != 1).sum() s_vol = ps_disk(r).sum() if ndim == 2 else ps_ball(r).sum() bulk_vol = np.prod(shape) N = int(np.ceil((1 - porosity) * bulk_vol / s_vol)) im = np.random.random(size=shape) # Helper functions for calculating porosity: phi = g(f(N)) def f(N): return edt(im > N / bulk_vol) < r def g(im): r"""Returns fraction of 0s, given a binary image""" return 1 - im.sum() / np.prod(shape) # # Newton's method for getting image porosity match the given # w = 1.0 # Damping factor # dN = 5 if ndim == 2 else 25 # Perturbation # for i in range(maxiter): # err = g(f(N)) - porosity # d_err = (g(f(N+dN)) - g(f(N))) / dN # if d_err == 0: # break # if abs(err) <= tol: # break # N2 = N - int(err/d_err) # xnew = xold - f/df # N = w * N2 + (1-w) * N # Bisection search: N is always undershoot (bc. of overlaps) N_low, N_high = N, 4 * N for i in range(maxiter): N = np.mean([N_high, N_low], dtype=int) err = g(f(N)) - porosity if err > 0: N_low = N else: N_high = N if abs(err) <= tol: break return ~f(N)
[docs]def blobs(shape: List[int], porosity: float = 0.5, blobiness: int = 1, divs: int = 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 ``None`` is specified, then the scalar noise field is converted to a uniform distribution and returned without thresholding. blobiness : int or list of ints(default = 1) Controls the morphology of the blobs. A higher number results in a larger number of small blobs. If a list is supplied then the blobs are anisotropic. divs : int or array_like The number of times to divide the image for parallel processing. If ``1`` then parallel processing does not occur. ``2`` is equivalent to ``[2, 2, 2]`` for a 3D image. The number of cores used is specified in ``porespy.settings.ncores`` and defaults to all cores. Returns ------- image : ndarray A boolean array with ``True`` values denoting the pore space See Also -------- norm_to_uniform Notes ----- This function generates random noise, the applies a gaussian blur to the noise with a sigma controlled by the blobiness argument as: $$ np.mean(shape) / (40 * blobiness) $$ The value of 40 was chosen so that a ``blobiness`` of 1 gave a reasonable result. Examples -------- `Click here <https://porespy.org/examples/generators/howtos/blobs.html>`_ to view online example. """ if isinstance(shape, int): shape = [shape]*3 if len(shape) == 1: shape = [shape[0]]*3 shape = np.array(shape) if isinstance(blobiness, int): blobiness = [blobiness]*len(shape) blobiness = np.array(blobiness) parallel = False if isinstance(divs, int): divs = [divs]*len(shape) if max(divs) > 1: parallel = True logger.info(f'Performing {insp.currentframe().f_code.co_name} in parallel') sigma = np.mean(shape) / (40 * blobiness) im = np.random.random(shape) if parallel: overlap = max([int(s*4) for s in sigma]) im = ps.filters.chunked_func(func=spim.gaussian_filter, input=im, sigma=sigma, divs=divs, overlap=overlap) else: im = spim.gaussian_filter(im, sigma=sigma) im = norm_to_uniform(im, scale=[0, 1]) if porosity: im = im < porosity return im
def _cylinders(shape: List[int], r: int, ncylinders: int, phi_max: float = 0, theta_max: float = 90, length: float = None, verbose: bool = True): r""" Generates a binary image of overlapping cylinders. This is a good approximation of a fibrous mat. Parameters ---------- shape : list The size of the image to generate in [Nx, Ny, Nz] where N is the number of voxels. 2D images are not permitted. r : int The radius of the cylinders in voxels ncylinders : int The number of cylinders to add to the domain. Adjust this value to control the final porosity, which is not easily specified since cylinders overlap and intersect different fractions of the domain. phi_max : int A value between 0 and 90 that controls the amount that the cylinders lie *out of* the XY plane, with 0 meaning all cylinders lie in the XY plane, and 90 meaning that cylinders are randomly oriented out of the plane by as much as +/- 90 degrees. theta_max : int A value between 0 and 90 that controls the amount of rotation *in the* XY plane, with 0 meaning all cylinders point in the X-direction, and 90 meaning they are randomly rotated about the Z axis by as much as +/- 90 degrees. length : int The length of the cylinders to add. If ``None`` (default) then the cylinders will extend beyond the domain in both directions so no ends will exist. If a scalar value is given it will be interpreted as the Euclidean distance between the two ends of the cylinder. Note that one or both of the ends *may* still lie outside the domain, depending on the randomly chosen center point of the cylinder. Returns ------- image : ndarray A boolean array with ``True`` values denoting the pore space Examples -------- `Click here <https://porespy.org/examples/generators/cylinders.html>`_ to view online example. """ shape = np.array(shape) if np.size(shape) == 1: shape = np.full((3, ), int(shape)) elif np.size(shape) == 2: raise Exception("2D cylinders don't make sense") # Find hypotenuse of domain from [0,0,0] to [Nx,Ny,Nz] H = np.sqrt(np.sum(np.square(shape))).astype(int) if length is None: # Assume cylinders span domain if length not given length = 2 * H R = min(int(length / 2), 2 * H) # Trim given length to 2H if too long # Adjust max angles to be between 0 and 90 if (phi_max > 90) or (phi_max < 0): raise Exception('phi_max must be betwen 0 and 90') if (theta_max > 90) or (theta_max < 0): raise Exception('theta_max must be betwen 0 and 90') # Create empty image for inserting into im = np.zeros(shape, dtype=bool) n = 0 L = min(H, R) # Disable tqdm if called from another tqdm to prevent double pbars tqdm_settings = settings.tqdm.copy() if not settings.tqdm["disable"]: tqdm_settings = {**settings.tqdm, **{'disable': not verbose}} with tqdm(ncylinders, **tqdm_settings) as pbar: while n < ncylinders: # Choose a random starting point in domain x = np.random.rand(3) * (shape + 2 * L) # Chose a random phi and theta within given ranges phi = (np.pi / 2 - np.pi * np.random.rand()) * phi_max / 90 theta = (np.pi / 2 - np.pi * np.random.rand()) * theta_max / 90 X0 = R * np.array([np.cos(phi) * np.cos(theta), np.cos(phi) * np.sin(theta), np.sin(phi)]) [X0, X1] = [x + X0, x - X0] crds = line_segment(X0, X1) lower = ~np.any(np.vstack(crds).T < [L, L, L], axis=1) upper = ~np.any(np.vstack(crds).T >= shape + L, axis=1) valid = upper * lower if np.any(valid): im[crds[0][valid] - L, crds[1][valid] - L, crds[2][valid] - L] = 1 n += 1 pbar.update() im = np.array(im, dtype=bool) dt = edt(~im) < r return ~dt
[docs]def cylinders(shape: List[int], r: int, ncylinders: int = None, porosity: float = None, phi_max: float = 0, theta_max: float = 90, length: float = None, maxiter: int = 3): r""" Generates a binary image of overlapping cylinders given porosity OR number of cylinders. This is a good approximation of a fibrous mat. Parameters ---------- shape : list The size of the image to generate in [Nx, Ny, Nz] where N is the number of voxels. 2D images are not permitted. r : scalar The radius of the cylinders in voxels ncylinders : int The number of cylinders to add to the domain. Adjust this value to control the final porosity, which is not easily specified since cylinders overlap and intersect different fractions of the domain. porosity : float The targeted value for the porosity of the generated mat. The function uses an algorithm for predicted the number of required number of cylinder, and refines this over a certain number of fractional insertions (according to the 'iterations' input). phi_max : int A value between 0 and 90 that controls the amount that the cylinders lie *out of* the XY plane, with 0 meaning all cylinders lie in the XY plane, and 90 meaning that cylinders are randomly oriented out of the plane by as much as +/- 90 degrees. theta_max : int A value between 0 and 90 that controls the amount of rotation *in the* XY plane, with 0 meaning all cylinders point in the X-direction, and 90 meaning they are randomly rotated about the Z axis by as much as +/- 90 degrees. length : int The length of the cylinders to add. If ``None`` (default) then the cylinders will extend beyond the domain in both directions so no ends will exist. If a scalar value is given it will be interpreted as the Euclidean distance between the two ends of the cylinder. Note that one or both of the ends *may* still lie outside the domain, depending on the randomly chosen center point of the cylinder. maxiter : int The number of fractional fiber insertions used to target the requested porosity. By default a value of 3 is used (and this is typically effective in getting very close to the targeted porosity), but a greater number can be input to improve the achieved porosity. Returns ------- image : ndarray A boolean array with ``True`` values denoting the pore space Notes ----- The cylinders_porosity function works by estimating the number of cylinders needed to be inserted into the domain by estimating cylinder length, and exploiting the fact that, when inserting any potentially overlapping objects randomly into a volume v_total (which has units of pixels and is equal to dimx x dimy x dimz, for example), such that the total volume of objects added to the volume is v_added (and includes any volume that was inserted but overlapped with already occupied space), the resulting porosity will be equal to exp(-v_added/v_total). After intially estimating the cylinder number and inserting a small fraction of the estimated number, the true cylinder volume is calculated, the estimate refined, and a larger fraction of cylinders inserted. This is repeated a number of times according to the ``maxiter`` argument, yielding an image with a porosity close to the goal. Examples -------- `Click here <https://porespy.org/examples/generators/howtos/cylinders.html>`_ to view online example. """ if ncylinders is not None: im = _cylinders( shape=shape, r=r, ncylinders=ncylinders, phi_max=phi_max, theta_max=theta_max, length=length, ) return im if porosity is None: raise Exception("'ncylinders' and 'porosity' can't be both None") # if maxiter < 3: # raise Exception("Iterations must be greater than or equal to 3") vol_total = float(np.prod(shape)) def get_num_pixels(porosity): r""" Helper method to calculate number of pixels given a porosity """ return -np.log(porosity) * vol_total # Crudely estimate fiber length as cube root of product of dims length_estimate = vol_total ** (1 / 3) if length is None else length # Rough fiber volume estimate vol_fiber = length_estimate * np.pi * r * r n_pixels_to_add = get_num_pixels(porosity) # Rough estimate of n_fibers n_fibers_added = 0 # Calculate fraction of fibers to be added in each iteration. subdif = 0.8 / np.sum(np.arange(1, maxiter) ** 2) fractions = [0.2] for i in range(1, maxiter): fractions.append(fractions[i - 1] + (maxiter - i) ** 2 * subdif) im = np.ones(shape, dtype=bool) for frac in tqdm(fractions, **settings.tqdm): n_fibers_total = n_pixels_to_add / vol_fiber n_fibers = int(np.ceil(frac * n_fibers_total) - n_fibers_added) if n_fibers > 0: im = im & _cylinders(shape, r, n_fibers, phi_max, theta_max, length, verbose=False) n_fibers_added += n_fibers # Update parameters for next iteration porosity = ps.metrics.porosity(im) vol_added = get_num_pixels(porosity) vol_fiber = vol_added / n_fibers_added logger.debug(f"{n_fibers_added} fibers added to reach target porosity.") return im
[docs]def line_segment(X0, X1): r""" Calculate the voxel coordinates of a straight line between the two given end points Parameters ---------- X0 and X1 : array_like The [x, y] or [x, y, z] coordinates of the start and end points of the line. Returns ------- coords : list of lists 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. Examples -------- `Click here <https://porespy.org/examples/generators/howtos/line_segment.html>`_ to view online example. """ X0 = np.around(X0).astype(int) X1 = np.around(X1).astype(int) if len(X0) == 3: L = np.amax(np.absolute([[X1[0] - X0[0]], [X1[1] - X0[1]], [X1[2] - X0[2]]])) + 1 x = np.rint(np.linspace(X0[0], X1[0], L)).astype(int) y = np.rint(np.linspace(X0[1], X1[1], L)).astype(int) z = np.rint(np.linspace(X0[2], X1[2], L)).astype(int) return [x, y, z] else: L = np.amax(np.absolute([[X1[0] - X0[0]], [X1[1] - X0[1]]])) + 1 x = np.rint(np.linspace(X0[0], X1[0], L)).astype(int) y = np.rint(np.linspace(X0[1], X1[1], L)).astype(int) return [x, y]