import logging
import numpy as np
import scipy.ndimage as spim
from pandas import DataFrame
from skimage.measure import mesh_surface_area, regionprops
from skimage.measure._regionprops import RegionProperties
from skimage.morphology import ball, skeletonize_3d
from porespy.tools import bbox_to_slices, extract_subsection
try:
from skimage.measure import marching_cubes
except ImportError:
from skimage.measure import marching_cubes_lewiner as marching_cubes
try:
from pyedt import edt
except ModuleNotFoundError:
from edt import edt
__all__ = [
"regionprops_3D",
"props_to_DataFrame",
"prop_to_image",
]
logger = logging.getLogger(__name__)
[docs]
def props_to_DataFrame(regionprops):
r"""
Create a ``pandas`` DataFrame containing all the scalar metrics for each
region, such as volume, sphericity, and so on, calculated by
``regionprops_3D``.
Parameters
----------
regionprops : list
This is a list of properties for each region that is computed
by ``regionprops_3D``. Because ``regionprops_3D`` returns data in
the same ``list`` format as the ``regionprops`` function in **Skimage**
you can pass in either.
Returns
-------
DataFrame : Pandas DataFrame
A Pandas DataFrame with each region corresponding to a row and each
column corresponding to a key metric. All the values for a given
property (e.g. 'sphericity') can be obtained as
``val = df['sphericity']``. Conversely, all the key metrics for a
given region can be found with ``df.iloc[1]``.
See Also
--------
prop_to_image
regionprops_3d
Examples
--------
`Click here
<https://porespy.org/examples/metrics/reference/props_to_DataFrame.html>`_
to view online example.
"""
# Parse the regionprops list and pull out all props with scalar values
metrics = []
reg = regionprops[0]
for item in reg.__dir__():
if not item.startswith('_'):
try:
if np.shape(getattr(reg, item)) == ():
metrics.append(item)
except (TypeError, NotImplementedError, AttributeError):
pass
# Create a dictionary of all metrics that are simple scalar properties
d = {}
for i, k in enumerate(metrics):
logger.info("Processing {k}")
try:
d[k] = np.array([r[k] for r in regionprops])
except ValueError: # pragma: no cover
logger.error(f'Error encountered evaluating {k} so skipping it')
# Create pandas data frame an return
df = DataFrame(d)
return df
[docs]
def prop_to_image(regionprops, shape, prop):
r"""
Create an image with each region colored according the specified ``prop``,
as obtained by ``regionprops_3d``.
Parameters
----------
regionprops : list
This is a list of properties for each region that is computed
by PoreSpy's ``regionprops_3D`` or Skimage's ``regionsprops``.
shape : array_like
The shape of the original image for which ``regionprops`` was obtained.
prop : string
The region property of interest. Can be a scalar item such as 'volume'
in which case the the regions will be colored by their respective
volumes, or can be an image-type property such as 'border' or
'convex_image', which will return an image composed of the sub-images.
Returns
-------
image : ndarray
An ndarray the same size as the original image, with each region
represented by the values specified in ``prop``.
See Also
--------
props_to_DataFrame
regionprops_3d
Examples
--------
`Click here
<https://porespy.org/examples/metrics/reference/prop_to_image.html>`_
to view online example.
"""
im = np.zeros(shape=shape)
for r in regionprops:
if prop == 'convex':
mask = r.convex_image
else:
mask = r.image
temp = mask * r[prop]
s = bbox_to_slices(r.bbox)
im[s] += temp
return im
[docs]
def regionprops_3D(im):
r"""
Calculates various metrics for each labeled region in a 3D image.
This functions offers a few extras for 3D images that are not provided
by the ``regionprops`` function in ``scikit-image``.
Parameters
----------
im : array_like
An image containing at least one labeled region. If a boolean image
is received than the ``True`` voxels are treated as a single region
labeled ``1``. Regions labeled 0 are ignored in all cases.
Returns
-------
props : list
An augmented version of the list returned by skimage's ``regionprops``.
Information, such as ``volume``, can be found for region A using the
following syntax: ``result[A-1].volume``.
The returned list contains all the metrics normally returned by
**skimage.measure.regionprops** plus the following:
'slices'
Slice indices into the image that can be used to extract the
region
'volume'
Volume of the region in number of voxels.
'bbox_volume'
Volume of the bounding box that contains the region.
'border'
The edges of the region, found as the locations where the distance
transform is 1.
'inscribed_sphere'
An image containing the largest sphere can can fit entirely inside
the region.
'surface_mesh_vertices'
Obtained by applying the marching cubes algorithm on the region,
AFTER first blurring the voxel image. This allows marching cubes
more freedom to fit the surface contours.
See also ``surface_mesh_simplices``
'surface_mesh_simplices'
This accompanies ``surface_mesh_vertices`` and together they can
be used to define the region as a mesh.
'surface_area'
Calculated using the mesh obtained as described above, using the
``porespy.metrics.mesh_surface_area`` method.
'sphericity'
Defined as the ratio of the area of a sphere with the same volume
as the region to the actual surface area of the region.
'skeleton'
The medial axis of the region obtained using the ``skeletonize_3D``
method from **skimage**.
'convex_volume'
Same as convex_area, but translated to a more meaningful name.
See Also
--------
snow_partitioning
Notes
-----
Regions can be identified using a watershed algorithm, which can be a bit
tricky to obtain desired results. *PoreSpy* includes the SNOW algorithm,
which may be helpful.
Examples
--------
`Click here
<https://porespy.org/examples/metrics/reference/regionprops_3D.html>`_
to view online example.
"""
results = regionprops(im)
for i, obj in enumerate(results):
a = results[i]
b = RegionPropertiesPS(a.slice,
a.label,
a._label_image,
a._intensity_image,
a._cache_active)
results[i] = b
return results
class RegionPropertiesPS(RegionProperties):
@property
def mask(self):
return self.image
@property
def slices(self):
return self._slice
@property
def volume(self):
return self.area
@property
def bbox_volume(self):
mask = self.mask
return np.prod(mask.shape)
@property
def border(self):
return self.dt == 1
@property
def dt(self):
mask = self.mask
mask_padded = np.pad(mask, pad_width=1, mode='constant')
temp = edt(mask_padded)
return extract_subsection(temp, shape=mask.shape)
@property
def inscribed_sphere(self):
dt = self.dt
r = dt.max()
inv_dt = edt(dt < r)
return inv_dt < r
@property
def sphericity(self):
vol = self.volume
r = (3 / 4 / np.pi * vol)**(1 / 3)
a_equiv = 4 * np.pi * r**2
a_region = self.surface_area
return a_equiv / a_region
@property
def skeleton(self):
return skeletonize_3d(self.mask)
@property
def surface_area(self):
mask = self.mask
tmp = np.pad(np.atleast_3d(mask), pad_width=1, mode='constant')
tmp = spim.convolve(tmp, weights=ball(1)) / 5
verts, faces, norms, vals = marching_cubes(volume=tmp, level=0)
self._surface_mesh_vertices = verts
self._surface_mesh_simplices = faces
area = mesh_surface_area(verts, faces)
return area
@property
def surface_mesh_vertices(self):
if not hasattr(self, '_surface_mesh_vertices'):
_ = self.surface_area
return self._surface_mesh_vertices
@property
def surface_mesh_simplices(self):
if not hasattr(self, '_surface_mesh_simplices'):
_ = self.surface_area
return self._surface_mesh_simplices
@property
def convex_volume(self):
return self.convex_area