From 7b0703f663bb3cb062656611e511f58db7dd799d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Mon, 14 May 2012 21:30:41 +0200 Subject: [PATCH] regionprops takes optional intensity images --- skimage/measure/_moments.pyx | 2 +- skimage/measure/_regionprops.py | 122 +++++++++++++++++++--- skimage/measure/tests/test_regionprops.py | 90 ++++++++++++++-- 3 files changed, 192 insertions(+), 22 deletions(-) diff --git a/skimage/measure/_moments.pyx b/skimage/measure/_moments.pyx index d24dbea2..f84e14dd 100644 --- a/skimage/measure/_moments.pyx +++ b/skimage/measure/_moments.pyx @@ -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 diff --git a/skimage/measure/_regionprops.py b/skimage/measure/_regionprops.py index 5c840410..937b80f5 100644 --- a/skimage/measure/_regionprops.py +++ b/skimage/measure/_regionprops.py @@ -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 diff --git a/skimage/measure/tests/test_regionprops.py b/skimage/measure/tests/test_regionprops.py index 780714b9..417311de 100644 --- a/skimage/measure/tests/test_regionprops.py +++ b/skimage/measure/tests/test_regionprops.py @@ -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