wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-02 17:40:31 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1505 from emmanuelle/regionprops3d

Browse Source 
Add preliminary and incomplete 3D support to regionprops.
This commit is contained in:
Juan Nunez-Iglesias
2015-12-29 11:59:48 +11:00
parent 99871d77b4 6c9325aadf
commit 2e812e0b12
2 changed files with 127 additions and 54 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/measure/_regionprops.py
+58 -18
View File
@@ -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.')
skimage/measure/tests/test_regionprops.py
+69 -36
View File
@@ -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)