# coding: utf-8 from __future__ import division from math import sqrt, atan2, pi as PI import itertools from warnings import warn import numpy as np from scipy import ndimage as ndi from ._label import label from . import _moments from functools import wraps __all__ = ['regionprops', 'perimeter'] XY_TO_RC_DEPRECATION_MESSAGE = ( 'regionprops and image moments (including moments, normalized moments, ' 'central moments, and inertia tensor) of 2D images will change from xy ' 'coordinates to rc coordinates in version 0.16.\nSee ' 'http://scikit-image.org/docs/0.14.x/release_notes_and_installation.html#deprecations ' 'for details on how to avoid this message.' ) STREL_4 = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]], dtype=np.uint8) STREL_8 = np.ones((3, 3), dtype=np.uint8) STREL_26_3D = np.ones((3, 3, 3), dtype=np.uint8) PROPS = { 'Area': 'area', 'BoundingBox': 'bbox', 'BoundingBoxArea': 'bbox_area', 'CentralMoments': 'moments_central', 'Centroid': 'centroid', 'ConvexArea': 'convex_area', # 'ConvexHull', 'ConvexImage': 'convex_image', 'Coordinates': 'coords', 'Eccentricity': 'eccentricity', 'EquivDiameter': 'equivalent_diameter', 'EulerNumber': 'euler_number', 'Extent': 'extent', # 'Extrema', 'FilledArea': 'filled_area', 'FilledImage': 'filled_image', 'HuMoments': 'moments_hu', 'Image': 'image', 'Label': 'label', 'MajorAxisLength': 'major_axis_length', 'MaxIntensity': 'max_intensity', 'MeanIntensity': 'mean_intensity', 'MinIntensity': 'min_intensity', 'MinorAxisLength': 'minor_axis_length', 'Moments': 'moments', 'NormalizedMoments': 'moments_normalized', 'Orientation': 'orientation', 'Perimeter': 'perimeter', # 'PixelIdxList', # 'PixelList', 'Solidity': 'solidity', # 'SubarrayIdx' 'WeightedCentralMoments': 'weighted_moments_central', 'WeightedCentroid': 'weighted_centroid', 'WeightedHuMoments': 'weighted_moments_hu', 'WeightedMoments': 'weighted_moments', 'WeightedNormalizedMoments': 'weighted_moments_normalized' } PROP_VALS = set(PROPS.values()) def _cached(f): @wraps(f) def wrapper(obj): cache = obj._cache prop = f.__name__ if not ((prop in cache) and obj._cache_active): cache[prop] = f(obj) return cache[prop] return wrapper def only2d(method): @wraps(method) def func2d(self, *args, **kwargs): if self._ndim > 2: raise NotImplementedError('Property %s is not implemented for ' '3D images' % method.__name__) return method(self, *args, **kwargs) return func2d class _RegionProperties(object): """Please refer to `skimage.measure.regionprops` for more information on the available region properties. """ def __init__(self, slice, label, label_image, intensity_image, cache_active, coordinates): if intensity_image is not None: if not intensity_image.shape == label_image.shape: raise ValueError('Label and intensity image must have the' 'same shape.') self.label = label self._slice = slice self._label_image = label_image self._intensity_image = intensity_image self._cache_active = cache_active self._cache = {} self._ndim = label_image.ndim # Note: in PR 2603, we added support for nD moments in regionprops. # Many properties used xy coordinates, instead of rc. This attribute # helps with the deprecation process and should be removed in 0.16. if label_image.ndim > 2 or coordinates == 'rc': self._use_xy_warning = False self._transpose_moments = False elif coordinates == 'xy': self._use_xy_warning = False # don't warn if 'xy' given explicitly self._transpose_moments = True elif coordinates is None: self._use_xy_warning = True self._transpose_moments = True else: raise ValueError('Incorrect value for regionprops coordinates: %s.' ' Possible values are: "rc", "xy", or None') @_cached def area(self): return np.sum(self.image) def bbox(self): """ Returns ------- A tuple of the bounding box's start coordinates for each dimension, followed by the end coordinates for each dimension """ return tuple([self._slice[i].start for i in range(self._ndim)] + [self._slice[i].stop for i in range(self._ndim)]) def bbox_area(self): return self.image.size def centroid(self): return tuple(self.coords.mean(axis=0)) @_cached def convex_area(self): return np.sum(self.convex_image) @_cached def convex_image(self): from ..morphology.convex_hull import convex_hull_image return convex_hull_image(self.image) def coords(self): indices = np.nonzero(self.image) return np.vstack([indices[i] + self._slice[i].start for i in range(self._ndim)]).T @only2d def eccentricity(self): l1, l2 = self.inertia_tensor_eigvals if l1 == 0: return 0 return sqrt(1 - l2 / l1) def equivalent_diameter(self): if self._ndim == 2: return sqrt(4 * self.area / PI) elif self._ndim == 3: return (6 * self.area / PI) ** (1. / 3) def euler_number(self): euler_array = self.filled_image != self.image _, num = label(euler_array, connectivity=self._ndim, return_num=True, background=0) return -num + 1 def extent(self): return self.area / self.image.size def filled_area(self): return np.sum(self.filled_image) @_cached def filled_image(self): structure = np.ones((3,) * self._ndim) return ndi.binary_fill_holes(self.image, structure) @_cached def image(self): return self._label_image[self._slice] == self.label @_cached def inertia_tensor(self): mu = self.moments_central return _moments.inertia_tensor(self.image, mu) @_cached def inertia_tensor_eigvals(self): return _moments.inertia_tensor_eigvals(self.image, T=self.inertia_tensor) @_cached def intensity_image(self): if self._intensity_image is None: raise AttributeError('No intensity image specified.') return self._intensity_image[self._slice] * self.image def _intensity_image_double(self): return self.intensity_image.astype(np.double) def local_centroid(self): M = self.moments if self._transpose_moments: M = M.T return tuple(M[tuple(np.eye(self._ndim, dtype=int))] / M[(0,) * self._ndim]) def max_intensity(self): return np.max(self.intensity_image[self.image]) def mean_intensity(self): return np.mean(self.intensity_image[self.image]) def min_intensity(self): return np.min(self.intensity_image[self.image]) def major_axis_length(self): l1 = self.inertia_tensor_eigvals[0] return 4 * sqrt(l1) def minor_axis_length(self): l2 = self.inertia_tensor_eigvals[-1] return 4 * sqrt(l2) @_cached def moments(self): M = _moments.moments(self.image.astype(np.uint8), 3) if self._use_xy_warning: warn(XY_TO_RC_DEPRECATION_MESSAGE) if self._transpose_moments: M = M.T return M @_cached def moments_central(self): mu = _moments.moments_central(self.image.astype(np.uint8), self.local_centroid, order=3) if self._use_xy_warning: warn(XY_TO_RC_DEPRECATION_MESSAGE) if self._transpose_moments: mu = mu.T return mu @only2d def moments_hu(self): return _moments.moments_hu(self.moments_normalized) @_cached def moments_normalized(self): return _moments.moments_normalized(self.moments_central, 3) @only2d def orientation(self): a, b, b, c = self.inertia_tensor.flat sign = -1 if self._transpose_moments else 1 if a - c == 0: if b < 0: return -PI / 4. else: return PI / 4. else: return sign * 0.5 * atan2(-2 * b, c - a) @only2d def perimeter(self): return perimeter(self.image, 4) def solidity(self): return self.area / self.convex_area def weighted_centroid(self): ctr = self.weighted_local_centroid return tuple(idx + slc.start for idx, slc in zip(ctr, self._slice)) def weighted_local_centroid(self): M = self.weighted_moments return (M[tuple(np.eye(self._ndim, dtype=int))] / M[(0,) * self._ndim]) @_cached def weighted_moments(self): return _moments.moments(self._intensity_image_double(), 3) @_cached def weighted_moments_central(self): ctr = self.weighted_local_centroid return _moments.moments_central(self._intensity_image_double(), center=ctr, order=3) @only2d def weighted_moments_hu(self): return _moments.moments_hu(self.weighted_moments_normalized) @_cached def weighted_moments_normalized(self): return _moments.moments_normalized(self.weighted_moments_central, 3) def __iter__(self): props = PROP_VALS if self._intensity_image is None: unavailable_props = ('intensity_image', 'max_intensity', 'mean_intensity', 'min_intensity', 'weighted_moments', 'weighted_moments_central', 'weighted_centroid', 'weighted_local_centroid', 'weighted_moments_hu', 'weighted_moments_normalized') props = props.difference(unavailable_props) return iter(sorted(props)) def __getitem__(self, key): value = getattr(self, key, None) if value is not None: return value else: # backwards compatability return getattr(self, PROPS[key]) def __eq__(self, other): if not isinstance(other, _RegionProperties): return False for key in PROP_VALS: try: # so that NaNs are equal np.testing.assert_equal(getattr(self, key, None), getattr(other, key, None)) except AssertionError: return False return True def regionprops(label_image, intensity_image=None, cache=True, coordinates=None): """Measure properties of labeled image regions. Parameters ---------- label_image : (N, M) ndarray Labeled input image. Labels with value 0 are ignored. .. versionchanged:: 0.14.2 Previously, ``label_image`` was processed by ``numpy.squeeze`` and so any number of singleton dimensions was allowed. This resulted in inconsistent handling of images with singleton dimensions. To recover the old behaviour, use ``regionprops(np.squeeze(label_image), ...)``. intensity_image : (N, M) ndarray, optional Intensity (i.e., input) image with same size as labeled image. Default is None. cache : bool, optional Determine whether to cache calculated properties. The computation is much faster for cached properties, whereas the memory consumption increases. coordinates : 'rc' or 'xy', optional Coordinate conventions for 2D images. (Only 'rc' coordinates are supported for 3D images.) Returns ------- properties : list of RegionProperties Each item describes one labeled region, and can be accessed using the attributes listed below. Notes ----- The following properties can be accessed as attributes or keys: **area** : int Number of pixels of region. **bbox** : tuple Bounding box ``(min_row, min_col, max_row, max_col)``. Pixels belonging to the bounding box are in the half-open interval ``[min_row; max_row)`` and ``[min_col; max_col)``. **bbox_area** : int Number of pixels of bounding box. **centroid** : array Centroid coordinate tuple ``(row, col)``. **convex_area** : int Number of pixels of convex hull image. **convex_image** : (H, J) ndarray Binary convex hull image which has the same size as bounding box. **coords** : (N, 2) ndarray Coordinate list ``(row, col)`` of the region. **eccentricity** : float Eccentricity of the ellipse that has the same second-moments as the region. The eccentricity is the ratio of the focal distance (distance between focal points) over the major axis length. The value is in the interval [0, 1). When it is 0, the ellipse becomes a circle. **equivalent_diameter** : float The diameter of a circle with the same area as the region. **euler_number** : int Euler characteristic of region. Computed as number of objects (= 1) subtracted by number of holes (8-connectivity). **extent** : float Ratio of pixels in the region to pixels in the total bounding box. Computed as ``area / (rows * cols)`` **filled_area** : int Number of pixels of filled region. **filled_image** : (H, J) ndarray Binary region image with filled holes which has the same size as bounding box. **image** : (H, J) ndarray Sliced binary region image which has the same size as bounding box. **inertia_tensor** : (2, 2) ndarray Inertia tensor of the region for the rotation around its mass. **inertia_tensor_eigvals** : tuple The two eigen values of the inertia tensor in decreasing order. **intensity_image** : ndarray Image inside region bounding box. **label** : int The label in the labeled input image. **local_centroid** : array Centroid coordinate tuple ``(row, col)``, relative to region bounding box. **major_axis_length** : float The length of the major axis of the ellipse that has the same normalized second central moments as the region. **max_intensity** : float Value with the greatest intensity in the region. **mean_intensity** : float Value with the mean intensity in the region. **min_intensity** : float Value with the least intensity in the region. **minor_axis_length** : float The length of the minor axis of the ellipse that has the same normalized second central moments as the region. **moments** : (3, 3) ndarray Spatial moments up to 3rd order:: m_ji = sum{ array(x, y) * x^j * y^i } where the sum is over the `x`, `y` coordinates of the region. **moments_central** : (3, 3) ndarray Central moments (translation invariant) up to 3rd order:: mu_ji = sum{ array(x, y) * (x - x_c)^j * (y - y_c)^i } where the sum is over the `x`, `y` coordinates of the region, and `x_c` and `y_c` are the coordinates of the region's centroid. **moments_hu** : tuple Hu moments (translation, scale and rotation invariant). **moments_normalized** : (3, 3) ndarray Normalized moments (translation and scale invariant) up to 3rd order:: nu_ji = mu_ji / m_00^[(i+j)/2 + 1] where `m_00` is the zeroth spatial moment. **orientation** : float In 'rc' coordinates, angle between the 0th axis (rows) and the major axis of the ellipse that has the same second moments as the region, ranging from `-pi/2` to `pi/2` counter-clockwise. In `xy` coordinates, as above but the angle is now measured from the "x" or horizontal axis. **perimeter** : float Perimeter of object which approximates the contour as a line through the centers of border pixels using a 4-connectivity. **solidity** : float Ratio of pixels in the region to pixels of the convex hull image. **weighted_centroid** : array Centroid coordinate tuple ``(row, col)`` weighted with intensity image. **weighted_local_centroid** : array Centroid coordinate tuple ``(row, col)``, relative to region bounding box, weighted with intensity image. **weighted_moments** : (3, 3) ndarray Spatial moments of intensity image up to 3rd order:: wm_ji = sum{ array(x, y) * x^j * y^i } where the sum is over the `x`, `y` coordinates of the region. **weighted_moments_central** : (3, 3) ndarray Central moments (translation invariant) of intensity image up to 3rd order:: wmu_ji = sum{ array(x, y) * (x - x_c)^j * (y - y_c)^i } where the sum is over the `x`, `y` coordinates of the region, and `x_c` and `y_c` are the coordinates of the region's weighted centroid. **weighted_moments_hu** : tuple Hu moments (translation, scale and rotation invariant) of intensity image. **weighted_moments_normalized** : (3, 3) ndarray Normalized moments (translation and scale invariant) of intensity image up to 3rd order:: wnu_ji = wmu_ji / wm_00^[(i+j)/2 + 1] where ``wm_00`` is the zeroth spatial moment (intensity-weighted area). Each region also supports iteration, so that you can do:: for prop in region: print(prop, region[prop]) See Also -------- label References ---------- .. [1] Wilhelm Burger, Mark Burge. Principles of Digital Image Processing: Core Algorithms. Springer-Verlag, London, 2009. .. [2] B. Jähne. Digital Image Processing. Springer-Verlag, Berlin-Heidelberg, 6. edition, 2005. .. [3] T. H. Reiss. Recognizing Planar Objects Using Invariant Image Features, from Lecture notes in computer science, p. 676. Springer, Berlin, 1993. .. [4] http://en.wikipedia.org/wiki/Image_moment Examples -------- >>> from skimage import data, util >>> from skimage.measure import label >>> img = util.img_as_ubyte(data.coins()) > 110 >>> label_img = label(img, connectivity=img.ndim) >>> props = regionprops(label_img) >>> # centroid of first labeled object >>> props[0].centroid (22.729879860483141, 81.912285234465827) >>> # centroid of first labeled object >>> props[0]['centroid'] (22.729879860483141, 81.912285234465827) """ if label_image.ndim not in (2, 3): raise TypeError('Only 2-D and 3-D images supported.') if not np.issubdtype(label_image.dtype, np.integer): raise TypeError('Label image must be of integer type.') regions = [] objects = ndi.find_objects(label_image) for i, sl in enumerate(objects): if sl is None: continue label = i + 1 props = _RegionProperties(sl, label, label_image, intensity_image, cache, coordinates=coordinates) regions.append(props) return regions def perimeter(image, neighbourhood=4): """Calculate total perimeter of all objects in binary image. Parameters ---------- image : array Binary image. neighbourhood : 4 or 8, optional Neighborhood connectivity for border pixel determination. Returns ------- perimeter : float Total perimeter of all objects in binary image. References ---------- .. [1] K. Benkrid, D. Crookes. Design and FPGA Implementation of a Perimeter Estimator. The Queen's University of Belfast. http://www.cs.qub.ac.uk/~d.crookes/webpubs/papers/perimeter.doc """ if neighbourhood == 4: strel = STREL_4 else: strel = STREL_8 image = image.astype(np.uint8) eroded_image = ndi.binary_erosion(image, strel, border_value=0) border_image = image - eroded_image perimeter_weights = np.zeros(50, dtype=np.double) perimeter_weights[[5, 7, 15, 17, 25, 27]] = 1 perimeter_weights[[21, 33]] = sqrt(2) perimeter_weights[[13, 23]] = (1 + sqrt(2)) / 2 perimeter_image = ndi.convolve(border_image, np.array([[10, 2, 10], [ 2, 1, 2], [10, 2, 10]]), mode='constant', cval=0) # You can also write # return perimeter_weights[perimeter_image].sum() # but that was measured as taking much longer than bincount + np.dot (5x # as much time) perimeter_histogram = np.bincount(perimeter_image.ravel(), minlength=50) total_perimeter = np.dot(perimeter_histogram, perimeter_weights) return total_perimeter def _parse_docs(): import re import textwrap doc = regionprops.__doc__ matches = re.finditer('\*\*(\w+)\*\* \:.*?\n(.*?)(?=\n [\*\S]+)', doc, flags=re.DOTALL) prop_doc = dict((m.group(1), textwrap.dedent(m.group(2))) for m in matches) return prop_doc def _install_properties_docs(): prop_doc = _parse_docs() for p in [member for member in dir(_RegionProperties) if not member.startswith('_')]: try: getattr(_RegionProperties, p).__doc__ = prop_doc[p] except AttributeError: # For Python 2.x getattr(_RegionProperties, p).im_func.__doc__ = prop_doc[p] setattr(_RegionProperties, p, property(getattr(_RegionProperties, p))) _install_properties_docs()