regionprops takes optional intensity images

This commit is contained in:
Johannes Schönberger
2012-05-14 21:30:41 +02:00
parent 58d07c0a05
commit 7b0703f663
3 changed files with 192 additions and 22 deletions
+1 -1
View File
@@ -5,7 +5,7 @@ import numpy as np
cimport numpy as np
def central_moments(np.ndarray[np.uint8_t, ndim=2] array, double cr, double cc,
def central_moments(np.ndarray[np.double_t, ndim=2] array, double cr, double cc,
int order):
cdef int p, q, r, c
cdef np.ndarray[np.double_t, ndim=2] mu
+107 -15
View File
@@ -29,6 +29,9 @@ PROPS = (
'HuMoments',
'Image',
'MajorAxisLength',
'MaxIntensity',
'MeanIntensity',
'MinIntensity',
'MinorAxisLength',
'Moments',
'NormalizedMoments',
@@ -38,19 +41,26 @@ PROPS = (
# 'PixelList',
'Solidity',
# 'SubarrayIdx'
'WeightedCentralMoments',
'WeightedCentroid',
'WeightedHuMoments',
'WeightedMoments',
'WeightedNormalizedMoments'
)
def regionprops(image, properties='all'):
def regionprops(label_image, properties=['Area', 'Centroid'],
intensity_image=None):
"""Measure properties of labelled image regions.
Parameters
----------
image : N x M ndarray
label_image : N x M ndarray
Labelled input image.
properties : {'all', list}
Shape measurements to be determined for each labelled image region.
Default is 'all'. The following properties can be determined:
Default is `['Area', 'Centroid']`. The following properties can be
determined:
* Area : int
Number of pixels of region.
* BoundingBox : tuple
@@ -65,7 +75,7 @@ def regionprops(image, properties='all'):
* ConvexArea : int
Number of pixels of convex hull image.
* ConvexImage : H x J ndarray
Convex hull image which has the same size as bounding box.
Binary convex hull image which has the same size as bounding box.
* Eccentricity : float
Eccentricity of the ellipse that has the same second-moments as the
region. The eccentricity is the ratio of the distance between its
@@ -81,19 +91,25 @@ def regionprops(image, properties='all'):
* FilledArea : int
Number of pixels of filled region.
* FilledImage : H x J ndarray
Region image with filled holes which has the same size as bounding
box.
Binary region image with filled holes which has the same size as
bounding box.
* HuMoments : tuple
Hu moments (translation, scale and rotation invariant).
* Image : H x J ndarray
Sliced region image which has the same size as bounding box.
Sliced binary region image which has the same size as bounding box.
* MajorAxisLength : float
The length of the major axis of the ellipse that has the same
normalized second central moments as the region.
* MaxIntensity: float
Value with the greatest intensity in the region.
* MeanIntensity: float
Value with the mean intensity in the region.
* MinIntensity: float
Value with the least intensity in the region.
* MinorAxisLength : float
The length of the minor axis of the ellipse that has the same
normalized second central moments as the region.
* Moments 3 x 3 ndarray
* Moments : 3 x 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.
@@ -108,6 +124,30 @@ def regionprops(image, properties='all'):
`-pi/2` in counter-clockwise direction.
* Solidity : float
Ratio of pixels in the region to pixels of the convex hull image.
* WeightedCentralMoments : 3 x 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 centroid.
* WeightedCentroid : array
Centroid coordinate tuple `(row, col)` weighted with intensity
image.
* WeightedHuMoments : tuple
Hu moments (translation, scale and rotation invariant) of intensity
image.
* WeightedMoments : 3 x 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.
* WeightedNormalizedMoments : 3 x 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).
intensity_image : N x M ndarray, optional
Intensity image with same size as labelled image. Default is None.
Returns
-------
@@ -134,7 +174,7 @@ def regionprops(image, properties='all'):
>>> props = regionprops(label_img)
>>> props[0]['Centroid'] # centroid of first labelled object
"""
if not np.issubdtype(image.dtype, 'int'):
if not np.issubdtype(label_image.dtype, 'int'):
raise TypeError('labelled image must be of integer dtype')
# determine all properties if nothing specified
@@ -143,7 +183,7 @@ def regionprops(image, properties='all'):
props = []
objects = ndimage.find_objects(image)
objects = ndimage.find_objects(label_image)
for i, sl in enumerate(objects):
label = i + 1
@@ -153,9 +193,7 @@ def regionprops(image, properties='all'):
obj_props['Label'] = label
# binary image of i-th label, converting to uint8 because Cython
# does not have support for bool dtype
array = (image[sl] == label).astype('uint8')
array = (label_image[sl] == label).astype('double')
# upper left corner of object bbox
r0 = sl[0].start
@@ -203,7 +241,10 @@ def regionprops(image, properties='all'):
obj_props['ConvexImage'] = _convex_image
if 'Eccentricity' in properties:
obj_props['Eccentricity'] = sqrt(1 - l2 / l1)
if l1 == 0:
obj_props['Eccentricity'] = 0
else:
obj_props['Eccentricity'] = sqrt(1 - l2 / l1)
if 'EquivDiameter' in properties:
obj_props['EquivDiameter'] = sqrt(4 * m[0,0] / PI)
@@ -251,11 +292,62 @@ def regionprops(image, properties='all'):
obj_props['NormalizedMoments'] = _nu
if 'Orientation' in properties:
obj_props['Orientation'] = - 0.5 * atan(2 * b / (a - c))
if a - c == 0:
obj_props['Orientation'] = PI / 2
else:
obj_props['Orientation'] = - 0.5 * atan(2 * b / (a - c))
if 'Solidity' in properties:
if _convex_image is None:
_convex_image = convex_hull_image(array)
obj_props['Solidity'] = m[0,0] / np.sum(_convex_image)
if intensity_image is not None:
weighted_array = array * intensity_image[sl]
wm = _moments.central_moments(weighted_array, 0, 0, 3)
# weighted centroid
wcr = wm[0,1] / wm[0,0]
wcc = wm[1,0] / wm[0,0]
wmu = _moments.central_moments(weighted_array, wcr, wcc, 3)
# cached results which are used by several properties
_wnu = None
_vals = None
if 'MaxIntensity' in properties:
if _vals is None:
_vals = weighted_array[array.astype('bool')]
obj_props['MaxIntensity'] = np.max(_vals)
if 'MeanIntensity' in properties:
if _vals is None:
_vals = weighted_array[array.astype('bool')]
obj_props['MeanIntensity'] = np.mean(_vals)
if 'MinIntensity' in properties:
if _vals is None:
_vals = weighted_array[array.astype('bool')]
obj_props['MinIntensity'] = np.min(_vals)
if 'WeightedCentralMoments' in properties:
obj_props['WeightedCentralMoments'] = wmu
if 'WeightedCentroid' in properties:
obj_props['WeightedCentroid'] = wcr + r0, wcc + c0
if 'WeightedHuMoments' in properties:
if _wnu is None:
_wnu = _moments.normalized_moments(wmu, 3)
obj_props['WeightedHuMoments'] = _moments.hu_moments(_wnu)
if 'WeightedMoments' in properties:
obj_props['WeightedMoments'] = wm
if 'WeightedNormalizedMoments' in properties:
if _wnu is None:
_wnu = _moments.normalized_moments(wmu, 3)
obj_props['WeightedNormalizedMoments'] = _wnu
return props
+84 -6
View File
@@ -17,6 +17,8 @@ SAMPLE = np.array(
[0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1],
[0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1]]
)
INTENSITY_SAMPLE = SAMPLE.copy()
INTENSITY_SAMPLE[1,9:11] = 2
def test_area():
@@ -120,18 +122,33 @@ def test_filled_area():
area = regionprops(SAMPLE_mod, ['FilledArea'])[0]['FilledArea']
assert area == np.sum(SAMPLE)
def test_minor_axis_length():
length = regionprops(SAMPLE, ['MinorAxisLength'])[0]['MinorAxisLength']
# MATLAB has different interpretation of ellipse than found in literature,
# here implemented as found in literature
assert_almost_equal(length, 9.739302807263)
def test_major_axis_length():
length = regionprops(SAMPLE, ['MajorAxisLength'])[0]['MajorAxisLength']
# MATLAB has different interpretation of ellipse than found in literature,
# here implemented as found in literature
assert_almost_equal(length, 16.7924234999)
def test_max_intensity():
intensity = regionprops(SAMPLE, ['MaxIntensity'], INTENSITY_SAMPLE
)[0]['MaxIntensity']
assert_almost_equal(intensity, 2)
def test_mean_intensity():
intensity = regionprops(SAMPLE, ['MeanIntensity'], INTENSITY_SAMPLE
)[0]['MeanIntensity']
assert_almost_equal(intensity, 1.02777777777777)
def test_min_intensity():
intensity = regionprops(SAMPLE, ['MinIntensity'], INTENSITY_SAMPLE
)[0]['MinIntensity']
assert_almost_equal(intensity, 1)
def test_minor_axis_length():
length = regionprops(SAMPLE, ['MinorAxisLength'])[0]['MinorAxisLength']
# MATLAB has different interpretation of ellipse than found in literature,
# here implemented as found in literature
assert_almost_equal(length, 9.739302807263)
def test_moments():
m = regionprops(SAMPLE, ['Moments'])[0]['Moments']
#: determined with OpenCV
@@ -166,6 +183,67 @@ def test_solidity():
# determined with MATLAB
assert_almost_equal(solidity, 0.580645161290323)
def test_weighted_central_moments():
wmu = regionprops(SAMPLE, ['WeightedCentralMoments'], INTENSITY_SAMPLE
)[0]['WeightedCentralMoments']
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]]
)
np.set_printoptions(precision=10)
print wmu
assert_array_almost_equal(wmu, ref)
def test_weighted_centroid():
centroid = regionprops(SAMPLE, ['WeightedCentroid'], INTENSITY_SAMPLE
)[0]['WeightedCentroid']
assert_array_almost_equal(centroid, (5.540540540540, 9.445945945945))
def test_weighted_hu_moments():
whu = regionprops(SAMPLE, ['WeightedHuMoments'], INTENSITY_SAMPLE
)[0]['WeightedHuMoments']
ref = np.array([
3.1750587329e-01,
2.1417517159e-02,
2.3609322038e-02,
1.2565683360e-03,
8.3014209421e-07,
-3.5073773473e-05,
6.7936409056e-06
])
assert_array_almost_equal(whu, ref)
def test_weighted_moments():
wm = regionprops(SAMPLE, ['WeightedMoments'], INTENSITY_SAMPLE
)[0]['WeightedMoments']
ref = np.array(
[[ 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)
def test_weighted_normalized_moments():
wnu = regionprops(SAMPLE, ['WeightedNormalizedMoments'], INTENSITY_SAMPLE
)[0]['WeightedNormalizedMoments']
ref = np.array(
[[ np.nan, np.nan, 0.0873590903, -0.0161217406],
[ np.nan, -0.0160405109, -0.0031421072, -0.0031376984],
[ 0.230146783, 0.0457932622, 0.0165315478, 0.0043903193],
[-0.0162529732, -0.0104598869, -0.0028544152, -0.0011057191]]
)
assert_array_almost_equal(wnu, ref)
if __name__ == "__main__":
from numpy.testing import run_module_suite