Source code for porespy.networks._getnet

import logging
import numpy as np
import scipy.ndimage as spim
from skimage.morphology import disk, ball
from edt import edt
from porespy.tools import extend_slice
from porespy import settings
from porespy.tools import get_tqdm, make_contiguous
from porespy.metrics import region_surface_areas, region_interface_areas
from porespy.metrics import region_volumes


__all__ = [
    "regions_to_network",
]


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


[docs] def regions_to_network(regions, phases=None, voxel_size=1, accuracy='standard'): r""" Analyzes an image that has been partitioned into pore regions and extracts the pore and throat geometry as well as network connectivity. Parameters ---------- regions : ndarray An image of the material partitioned into individual regions. Zeros in this image are ignored. phases : ndarray, optional An image indicating to which phase each voxel belongs. The returned network contains a 'pore.phase' array with the corresponding value. If not given a value of 1 is assigned to every pore. voxel_size : scalar (default = 1) 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. accuracy : string Controls how accurately certain properties are calculated. Options are: 'standard' (default) Computes the surface areas and perimeters by simply counting voxels. This is *much* faster but does not properly account for the rough, voxelated nature of the surfaces. 'high' Computes surface areas using the marching cube method, and perimeters using the fast marching method. These are substantially slower but better account for the voxelated nature of the images. Returns ------- net : dict 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'). Notes ----- The meaning of each of the values returned in ``net`` are outlined below: 'pore.region_label' The region labels corresponding to the watershed extraction. The pore indices and regions labels will be offset by 1, so pore 0 will be region 1. 'throat.conns' An *Nt-by-2* array indicating which pores are connected to each other 'pore.region_label' Mapping of regions in the watershed segmentation to pores in the network 'pore.local_peak' The coordinates of the location of the maxima of the distance transform performed on the pore region in isolation 'pore.global_peak' The coordinates of the location of the maxima of the distance transform performed on the full image 'pore.geometric_centroid' The center of mass of the pore region as calculated by ``skimage.measure.center_of_mass`` 'throat.global_peak' The coordinates of the location of the maxima of the distance transform performed on the full image 'pore.region_volume' The volume of the pore region computed by summing the voxels 'pore.volume' The volume of the pore found by as volume of a mesh obtained from the marching cubes algorithm 'pore.surface_area' The surface area of the pore region as calculated by either counting voxels or using the fast marching method to generate a tetramesh (if ``accuracy`` is set to ``'high'``.) 'throat.cross_sectional_area' The cross-sectional area of the throat found by either counting voxels or using the fast marching method to generate a tetramesh (if ``accuracy`` is set to ``'high'``.) 'throat.perimeter' The perimeter of the throat found by counting voxels on the edge of the region defined by the intersection of two regions. 'pore.inscribed_diameter' The diameter of the largest sphere inscribed in the pore region. This is found as the maximum of the distance transform on the region in isolation. 'pore.extended_diameter' The diamter of the largest sphere inscribed in overal image, which can extend outside the pore region. This is found as the local maximum of the distance transform on the full image. 'throat.inscribed_diameter' The diameter of the largest sphere inscribed in the throat. This is found as the local maximum of the distance transform in the area where to regions meet. 'throat.total_length' The length between pore centered via the throat center 'throat.direct_length' The length between two pore centers on a straight line between them that does not pass through the throat centroid. Examples -------- `Click here <https://porespy.org/examples/networks/reference/regions_to_network.html>`_ to view online example. """ logger.info('Extracting pore/throat information') im = make_contiguous(regions) struc_elem = disk if im.ndim == 2 else ball voxel_size = float(voxel_size) if phases is None: phases = (im > 0).astype(int) if im.size != phases.size: raise Exception('regions and phase are different sizes, probably ' + 'because boundary regions were not added to phases') dt = np.zeros_like(phases, dtype="float32") # since edt returns float32 for i in np.unique(phases[phases.nonzero()]): dt += edt(phases == i) # Get 'slices' into im for each pore region slices = spim.find_objects(im) # Initialize arrays Ps = np.arange(1, np.amax(im)+1) Np = np.size(Ps) p_coords_cm = np.zeros((Np, im.ndim), dtype=float) p_coords_dt = np.zeros((Np, im.ndim), dtype=float) p_coords_dt_global = np.zeros((Np, im.ndim), dtype=float) p_volume = np.zeros((Np, ), dtype=float) p_dia_local = np.zeros((Np, ), dtype=float) p_dia_global = np.zeros((Np, ), dtype=float) p_label = np.zeros((Np, ), dtype=int) p_area_surf = np.zeros((Np, ), dtype=int) p_phase = np.zeros((Np, ), dtype=int) # The number of throats is not known at the start, so lists are used # which can be dynamically resized more easily. t_conns = [] t_dia_inscribed = [] t_area = [] t_perimeter = [] t_coords = [] # Start extracting size information for pores and throats msg = "Extracting pore and throat properties" for i in tqdm(Ps, desc=msg, **settings.tqdm): pore = i - 1 if slices[pore] is None: continue s = extend_slice(slices[pore], im.shape) sub_im = im[s] sub_dt = dt[s] pore_im = sub_im == i padded_mask = np.pad(pore_im, pad_width=1, mode='constant') pore_dt = edt(padded_mask) s_offset = np.array([i.start for i in s]) p_label[pore] = i p_coords_cm[pore, :] = spim.center_of_mass(pore_im) + s_offset temp = np.vstack(np.where(pore_dt == pore_dt.max()))[:, 0] p_coords_dt[pore, :] = temp + s_offset p_phase[pore] = (phases[s]*pore_im).max() temp = np.vstack(np.where(sub_dt == sub_dt.max()))[:, 0] p_coords_dt_global[pore, :] = temp + s_offset p_volume[pore] = np.sum(pore_im, dtype=np.int64) p_dia_local[pore] = 2*np.amax(pore_dt) p_dia_global[pore] = 2*np.amax(sub_dt) # The following is overwritten if accuracy is set to 'high' p_area_surf[pore] = np.sum(pore_dt == 1, dtype=np.int64) im_w_throats = spim.binary_dilation(input=pore_im, structure=struc_elem(1)) im_w_throats = im_w_throats*sub_im Pn = np.unique(im_w_throats)[1:] - 1 for j in Pn: if j > pore: t_conns.append([pore, j]) vx = np.where(im_w_throats == (j + 1)) t_dia_inscribed.append(2*np.amax(sub_dt[vx])) # The following is overwritten if accuracy is set to 'high' t_perimeter.append(np.sum(sub_dt[vx] < 2, dtype=np.int64)) # The following is overwritten if accuracy is set to 'high' t_area.append(np.size(vx[0])) p_area_surf[pore] -= np.size(vx[0]) t_inds = tuple([i+j for i, j in zip(vx, s_offset)]) temp = np.where(dt[t_inds] == np.amax(dt[t_inds]))[0][0] t_coords.append(tuple([t_inds[k][temp] for k in range(im.ndim)])) # Clean up values p_coords = p_coords_cm Nt = len(t_dia_inscribed) # Get number of throats if im.ndim == 2: # If 2D, add 0's in 3rd dimension p_coords = np.vstack((p_coords_cm.T, np.zeros((Np, )))).T t_coords = np.vstack((np.array(t_coords).T, np.zeros((Nt, )))).T net = {} ND = im.ndim # Define all the fundamental stuff net['throat.conns'] = np.array(t_conns) net['pore.coords'] = np.array(p_coords)*voxel_size net['pore.all'] = np.ones_like(net['pore.coords'][:, 0], dtype=bool) net['throat.all'] = np.ones_like(net['throat.conns'][:, 0], dtype=bool) net['pore.region_label'] = np.array(p_label) net['pore.phase'] = np.array(p_phase, dtype=int) net['throat.phases'] = net['pore.phase'][net['throat.conns']] V = np.copy(p_volume)*(voxel_size**ND) net['pore.region_volume'] = V # This will be an area if image is 2D f = 3/4 if ND == 3 else 1.0 net['pore.equivalent_diameter'] = 2*(V/np.pi * f)**(1/ND) # Extract the geometric stuff net['pore.local_peak'] = np.copy(p_coords_dt)*voxel_size net['pore.global_peak'] = np.copy(p_coords_dt_global)*voxel_size net['pore.geometric_centroid'] = np.copy(p_coords_cm)*voxel_size net['throat.global_peak'] = np.array(t_coords)*voxel_size net['pore.inscribed_diameter'] = np.copy(p_dia_local)*voxel_size net['pore.extended_diameter'] = np.copy(p_dia_global)*voxel_size net['throat.inscribed_diameter'] = np.array(t_dia_inscribed)*voxel_size P12 = net['throat.conns'] PT1 = np.sqrt(np.sum(((p_coords[P12[:, 0]]-t_coords)*voxel_size)**2, axis=1)) PT2 = np.sqrt(np.sum(((p_coords[P12[:, 1]]-t_coords)*voxel_size)**2, axis=1)) net['throat.total_length'] = PT1 + PT2 PT1 = PT1-p_dia_local[P12[:, 0]]/2*voxel_size PT2 = PT2-p_dia_local[P12[:, 1]]/2*voxel_size dist = (p_coords[P12[:, 0]] - p_coords[P12[:, 1]])*voxel_size net['throat.direct_length'] = np.sqrt(np.sum(dist**2, axis=1, dtype=np.int64)) net['throat.perimeter'] = np.array(t_perimeter)*voxel_size if (accuracy == 'high') and (im.ndim == 2): msg = "accuracy='high' only available in 3D, reverting to 'standard'" logger.warning(msg) accuracy = 'standard' if (accuracy == 'high'): net['pore.volume'] = region_volumes(regions=im, mode='marching_cubes') areas = region_surface_areas(regions=im, voxel_size=voxel_size) net['pore.surface_area'] = areas interface_area = region_interface_areas(regions=im, areas=areas, voxel_size=voxel_size) A = interface_area.area net['throat.cross_sectional_area'] = A net['throat.equivalent_diameter'] = (4*A/np.pi)**(1/2) else: net['pore.volume'] = np.copy(p_volume)*(voxel_size**ND) net['pore.surface_area'] = np.copy(p_area_surf)*(voxel_size**2) A = np.array(t_area)*(voxel_size**2) net['throat.cross_sectional_area'] = A net['throat.equivalent_diameter'] = (4*A/np.pi)**(1/2) return net