diff --git a/skimage/measure/_regionprops.py b/skimage/measure/_regionprops.py index dc01fc24..297a6390 100644 --- a/skimage/measure/_regionprops.py +++ b/skimage/measure/_regionprops.py @@ -1,4 +1,5 @@ # coding: utf-8 +from __future__ import division from math import sqrt, atan2, pi as PI import numpy as np from scipy import ndimage as ndi @@ -17,6 +18,7 @@ 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', @@ -73,6 +75,16 @@ def _cached(f): 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. @@ -80,11 +92,11 @@ class _RegionProperties(object): def __init__(self, slice, label, label_image, intensity_image, cache_active): - + 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 @@ -93,31 +105,38 @@ class _RegionProperties(object): self._cache_active = cache_active self._cache = {} + self._ndim = label_image.ndim + @_cached def area(self): - return self.moments[0, 0] + return np.sum(self.image) def bbox(self): - return (self._slice[0].start, self._slice[1].start, - self._slice[0].stop, self._slice[1].stop) + return tuple([self._slice[i].start for i in range(self._ndim)] + + [self._slice[i].stop for i in range(self._ndim)]) def centroid(self): - row, col = self.local_centroid - return row + self._slice[0].start, col + self._slice[1].start + centroid_coords = self.local_centroid + return tuple(centroid_coords + + np.array([self._slice[i].start + for i in range(self._ndim)])) + @only2d def convex_area(self): return np.sum(self.convex_image) @_cached + @only2d def convex_image(self): from ..morphology.convex_hull import convex_hull_image return convex_hull_image(self.image) def coords(self): - rr, cc = np.nonzero(self.image) - return np.vstack((rr + self._slice[0].start, - cc + self._slice[1].start)).T + 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: @@ -125,8 +144,12 @@ class _RegionProperties(object): return sqrt(1 - l2 / l1) def equivalent_diameter(self): - return sqrt(4 * self.moments[0, 0] / PI) + if self._ndim == 2: + return sqrt(4 * self.area / PI) + elif self._ndim == 3: + return (6 * self.area / PI) ** (1. / 3) + @only2d def euler_number(self): euler_array = self.filled_image != self.image _, num = label(euler_array, neighbors=8, return_num=True, @@ -134,21 +157,22 @@ class _RegionProperties(object): return -num + 1 def extent(self): - rows, cols = self.image.shape - return self.moments[0, 0] / (rows * cols) + return self.area / self.image.size def filled_area(self): return np.sum(self.filled_image) @_cached def filled_image(self): - return ndi.binary_fill_holes(self.image, STREL_8) + structure = STREL_8 if self._ndim == 2 else STREL_26_3D + return ndi.binary_fill_holes(self.image, structure) @_cached def image(self): return self._label_image[self._slice] == self.label @_cached + @only2d def inertia_tensor(self): mu = self.moments_central a = mu[2, 0] / mu[0, 0] @@ -157,6 +181,7 @@ class _RegionProperties(object): return np.array([[a, b], [b, c]]) @_cached + @only2d def inertia_tensor_eigvals(self): a, b, b, c = self.inertia_tensor.flat # eigen values of inertia tensor @@ -173,6 +198,7 @@ class _RegionProperties(object): def _intensity_image_double(self): return self.intensity_image.astype(np.double) + @only2d def local_centroid(self): m = self.moments row = m[0, 1] / m[0, 0] @@ -188,31 +214,38 @@ class _RegionProperties(object): def min_intensity(self): return np.min(self.intensity_image[self.image]) + @only2d def major_axis_length(self): l1, _ = self.inertia_tensor_eigvals return 4 * sqrt(l1) + @only2d def minor_axis_length(self): _, l2 = self.inertia_tensor_eigvals return 4 * sqrt(l2) @_cached + @only2d def moments(self): return _moments.moments(self.image.astype(np.uint8), 3) @_cached + @only2d def moments_central(self): row, col = self.local_centroid return _moments.moments_central(self.image.astype(np.uint8), row, col, 3) + @only2d def moments_hu(self): return _moments.moments_hu(self.moments_normalized) @_cached + @only2d def moments_normalized(self): return _moments.moments_normalized(self.moments_central, 3) + @only2d def orientation(self): a, b, b, c = self.inertia_tensor.flat b = -b @@ -224,16 +257,20 @@ class _RegionProperties(object): else: return - 0.5 * atan2(2 * b, (a - c)) + @only2d def perimeter(self): return perimeter(self.image, 4) + @only2d def solidity(self): return self.moments[0, 0] / np.sum(self.convex_image) + @only2d def weighted_centroid(self): row, col = self.weighted_local_centroid return row + self._slice[0].start, col + self._slice[1].start + @only2d def weighted_local_centroid(self): m = self.weighted_moments row = m[0, 1] / m[0, 0] @@ -241,20 +278,23 @@ class _RegionProperties(object): return row, col @_cached + @only2d def weighted_moments(self): - return _moments.moments_central(self._intensity_image_double(), - 0, 0, 3) + return _moments.moments_central(self._intensity_image_double(), 0, 0, 3) @_cached + @only2d def weighted_moments_central(self): row, col = self.weighted_local_centroid return _moments.moments_central(self._intensity_image_double(), row, col, 3) + @only2d def weighted_moments_hu(self): return _moments.moments_hu(self.weighted_moments_normalized) @_cached + @only2d def weighted_moments_normalized(self): return _moments.moments_normalized(self.weighted_moments_central, 3) @@ -476,8 +516,8 @@ def regionprops(label_image, intensity_image=None, cache=True): label_image = np.squeeze(label_image) - if label_image.ndim != 2: - raise TypeError('Only 2-D images supported.') + 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 integral type.') diff --git a/skimage/measure/tests/test_regionprops.py b/skimage/measure/tests/test_regionprops.py index 01330cdf..f724577d 100644 --- a/skimage/measure/tests/test_regionprops.py +++ b/skimage/measure/tests/test_regionprops.py @@ -22,6 +22,10 @@ SAMPLE = np.array( INTENSITY_SAMPLE = SAMPLE.copy() INTENSITY_SAMPLE[1, 9:11] = 2 +SAMPLE_3D = np.zeros((6, 6, 6), dtype=np.uint8) +SAMPLE_3D[1:3, 1:3, 1:3] = 1 +SAMPLE_3D[3, 2, 2] = 1 +INTENSITY_SAMPLE_3D = SAMPLE_3D.copy() def test_all_props(): region = regionprops(SAMPLE, INTENSITY_SAMPLE)[0] @@ -29,6 +33,14 @@ def test_all_props(): assert_almost_equal(region[prop], getattr(region, PROPS[prop])) +def test_all_props_3d(): + region = regionprops(SAMPLE_3D, INTENSITY_SAMPLE_3D)[0] + for prop in PROPS: + try: + assert_almost_equal(region[prop], getattr(region, PROPS[prop])) + except NotImplementedError: + pass + def test_dtype(): regionprops(np.zeros((10, 10), dtype=np.int)) regionprops(np.zeros((10, 10), dtype=np.uint)) @@ -42,12 +54,15 @@ def test_ndim(): regionprops(np.zeros((10, 10), dtype=np.int)) regionprops(np.zeros((10, 10, 1), dtype=np.int)) regionprops(np.zeros((10, 10, 1, 1), dtype=np.int)) - assert_raises(TypeError, regionprops, np.zeros((10, 10, 2), dtype=np.int)) + regionprops(np.zeros((10, 10, 10), dtype=np.int)) + assert_raises(TypeError, regionprops, np.zeros((10, 10, 10, 2), dtype=np.int)) def test_area(): area = regionprops(SAMPLE)[0].area assert area == np.sum(SAMPLE) + area = regionprops(SAMPLE_3D)[0].area + assert area == np.sum(SAMPLE_3D) def test_bbox(): @@ -59,18 +74,21 @@ def test_bbox(): bbox = regionprops(SAMPLE_mod)[0].bbox assert_array_almost_equal(bbox, (0, 0, SAMPLE.shape[0], SAMPLE.shape[1]-1)) + bbox = regionprops(SAMPLE_3D)[0].bbox + assert_array_almost_equal(bbox, (1, 1, 1, 4, 3, 3)) + def test_moments_central(): mu = regionprops(SAMPLE)[0].moments_central # determined with OpenCV - assert_almost_equal(mu[0, 2], 436.00000000000045) + assert_almost_equal(mu[0,2], 436.00000000000045) # different from OpenCV results, bug in OpenCV - assert_almost_equal(mu[0, 3], -737.333333333333) - assert_almost_equal(mu[1, 1], -87.33333333333303) - assert_almost_equal(mu[1, 2], -127.5555555555593) - assert_almost_equal(mu[2, 0], 1259.7777777777774) - assert_almost_equal(mu[2, 1], 2000.296296296291) - assert_almost_equal(mu[3, 0], -760.0246913580195) + assert_almost_equal(mu[0,3], -737.333333333333) + assert_almost_equal(mu[1,1], -87.33333333333303) + assert_almost_equal(mu[1,2], -127.5555555555593) + assert_almost_equal(mu[2,0], 1259.7777777777774) + assert_almost_equal(mu[2,1], 2000.296296296291) + assert_almost_equal(mu[3,0], -760.0246913580195) def test_centroid(): @@ -110,6 +128,11 @@ def test_coordinates(): prop_coords = regionprops(sample)[0].coords assert_array_equal(prop_coords, coords) + sample = np.zeros((6, 6, 6), dtype=np.int8) + coords = np.array([[1, 1, 1], [1, 2, 1], [1, 3, 1]]) + sample[coords[:, 0], coords[:, 1], coords[:, 2]] = 1 + prop_coords = regionprops(sample)[0].coords + assert_array_equal(prop_coords, coords) def test_eccentricity(): eps = regionprops(SAMPLE)[0].eccentricity @@ -161,11 +184,17 @@ def test_image(): img = regionprops(SAMPLE)[0].image assert_array_equal(img, SAMPLE) + img = regionprops(SAMPLE_3D)[0].image + assert_array_equal(img, SAMPLE_3D[1:4, 1:3, 1:3]) + def test_label(): label = regionprops(SAMPLE)[0].label assert_array_equal(label, 1) + label = regionprops(SAMPLE_3D)[0].label + assert_array_equal(label, 1) + def test_filled_area(): area = regionprops(SAMPLE)[0].filled_area @@ -217,27 +246,27 @@ def test_minor_axis_length(): def test_moments(): m = regionprops(SAMPLE)[0].moments # determined with OpenCV - assert_almost_equal(m[0, 0], 72.0) - assert_almost_equal(m[0, 1], 408.0) - assert_almost_equal(m[0, 2], 2748.0) - assert_almost_equal(m[0, 3], 19776.0) - assert_almost_equal(m[1, 0], 680.0) - assert_almost_equal(m[1, 1], 3766.0) - assert_almost_equal(m[1, 2], 24836.0) - assert_almost_equal(m[2, 0], 7682.0) - assert_almost_equal(m[2, 1], 43882.0) - assert_almost_equal(m[3, 0], 95588.0) + assert_almost_equal(m[0,0], 72.0) + assert_almost_equal(m[0,1], 408.0) + assert_almost_equal(m[0,2], 2748.0) + assert_almost_equal(m[0,3], 19776.0) + assert_almost_equal(m[1,0], 680.0) + assert_almost_equal(m[1,1], 3766.0) + assert_almost_equal(m[1,2], 24836.0) + assert_almost_equal(m[2,0], 7682.0) + assert_almost_equal(m[2,1], 43882.0) + assert_almost_equal(m[3,0], 95588.0) def test_moments_normalized(): nu = regionprops(SAMPLE)[0].moments_normalized # determined with OpenCV - assert_almost_equal(nu[0, 2], 0.08410493827160502) - assert_almost_equal(nu[1, 1], -0.016846707818929982) - assert_almost_equal(nu[1, 2], -0.002899800614433943) - assert_almost_equal(nu[2, 0], 0.24301268861454037) - assert_almost_equal(nu[2, 1], 0.045473992910668816) - assert_almost_equal(nu[3, 0], -0.017278118992041805) + assert_almost_equal(nu[0,2], 0.08410493827160502) + assert_almost_equal(nu[1,1], -0.016846707818929982) + assert_almost_equal(nu[1,2], -0.002899800614433943) + assert_almost_equal(nu[2,0], 0.24301268861454037) + assert_almost_equal(nu[2,1], 0.045473992910668816) + assert_almost_equal(nu[3,0], -0.017278118992041805) def test_orientation(): @@ -277,14 +306,14 @@ def test_weighted_moments_central(): wmu = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE )[0].weighted_moments_central ref = np.array( - [[ 7.4000000000e+01, -2.1316282073e-13, - 4.7837837838e+02, -7.5943608473e+02], - [ 3.7303493627e-14, -8.7837837838e+01, - -1.4801314828e+02, -1.2714707125e+03], - [ 1.2602837838e+03, 2.1571526662e+03, - 6.6989799420e+03, 1.5304076361e+04], - [-7.6561796932e+02, -4.2385971907e+03, - -9.9501164076e+03, -3.3156729271e+04]] + [[ 7.4000000000e+01, -2.1316282073e-13, 4.7837837838e+02, + -7.5943608473e+02], + [ 3.7303493627e-14, -8.7837837838e+01, -1.4801314828e+02, + -1.2714707125e+03], + [ 1.2602837838e+03, 2.1571526662e+03, 6.6989799420e+03, + 1.5304076361e+04], + [ -7.6561796932e+02, -4.2385971907e+03, -9.9501164076e+03, + -3.3156729271e+04]] ) np.set_printoptions(precision=10) assert_array_almost_equal(wmu, ref) @@ -315,10 +344,14 @@ def test_weighted_moments(): wm = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE )[0].weighted_moments ref = np.array( - [[7.4000000e+01, 4.1000000e+02, 2.7500000e+03, 1.9778000e+04], - [6.9900000e+02, 3.7850000e+03, 2.4855000e+04, 1.7500100e+05], - [7.8630000e+03, 4.4063000e+04, 2.9347700e+05, 2.0810510e+06], - [9.7317000e+04, 5.7256700e+05, 3.9007170e+06, 2.8078871e+07]] + [[ 7.4000000000e+01, 4.1000000000e+02, 2.7500000000e+03, + 1.9778000000e+04], + [ 6.9900000000e+02, 3.7850000000e+03, 2.4855000000e+04, + 1.7500100000e+05], + [ 7.8630000000e+03, 4.4063000000e+04, 2.9347700000e+05, + 2.0810510000e+06], + [ 9.7317000000e+04, 5.7256700000e+05, 3.9007170000e+06, + 2.8078871000e+07]] ) assert_array_almost_equal(wm, ref)