From 8c1acfb7f51ae76c7dfca5adb0d9dc1722965dba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Tue, 1 May 2012 20:31:26 +0200 Subject: [PATCH 01/25] added regionprops function --- skimage/measure/__init__.py | 1 + skimage/measure/regionprops.pyx | 309 ++++++++++++++++++++++++++++++++ skimage/measure/setup.py | 3 + 3 files changed, 313 insertions(+) create mode 100644 skimage/measure/regionprops.pyx diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py index beedd379..be4a1b43 100755 --- a/skimage/measure/__init__.py +++ b/skimage/measure/__init__.py @@ -1 +1,2 @@ from .find_contours import find_contours +from .regionprops import regionprops diff --git a/skimage/measure/regionprops.pyx b/skimage/measure/regionprops.pyx new file mode 100644 index 00000000..b591b76f --- /dev/null +++ b/skimage/measure/regionprops.pyx @@ -0,0 +1,309 @@ +#cython: boundscheck=False +#cython: wraparound=False +#cython: cdivision=True +from scipy import ndimage +import numpy as np +cimport numpy as np +cimport cython +from libc.math cimport sqrt, atan2, fabs, fmin, fmax + +from skimage.morphology import convex_hull_image + + +__all__ = ['regionprops'] + + +STREL_8 = np.ones((3, 3), 'int8') +cdef float PI = 3.14159265 +cdef tuple PROPS = ( + 'Area', + 'BoundingBox', + 'CentralMoments', + 'Centroid', + 'ConvexArea', +# 'ConvexHull', + 'ConvexImage', + 'Eccentricity', + 'EquivDiameter', + 'EulerNumber', + 'Extent', +# 'Extrema', + 'FilledArea', + 'FilledImage', + 'HuMoments', + 'Image', + 'MajorAxisLength', + 'MinorAxisLength', + 'Moments', + 'NormalizedMoments', + 'Orientation', +# 'Perimeter', +# 'PixelIdxList', +# 'PixelList', + 'Solidity', +# 'SubarrayIdx' +) + + +def _moments(np.ndarray[np.uint8_t, ndim=2] array, int order): + cdef int p, q, r, c + cdef np.ndarray[np.double_t, ndim=2] m + m = np.zeros((order+1, order+1), 'double') + for p in range(order+1): + for q in range(order+1): + for r in range(array.shape[0]): + for c in range(array.shape[1]): + m[p,q] += array[r,c] * r**q * c**p + return m + +def _central_moments(np.ndarray[np.uint8_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 + mu = np.zeros((order+1, order+1), 'double') + for p in range(order+1): + for q in range(order+1): + for r in range(array.shape[0]): + for c in range(array.shape[1]): + mu[p,q] += array[r,c] * (r-cr)**q * (c-cc)**p + return mu + +def _normalized_moments(np.ndarray[np.double_t, ndim=2] mu, int order): + cdef int p, q + cdef np.ndarray[np.double_t, ndim=2] nu + nu = np.zeros((order+1, order+1), 'double') + for p in range(order+1): + for q in range(order+1): + if p + q >= 2: + nu[p,q] = mu[p,q] / mu[0,0]**((p+q) / 2 + 1) + else: + nu[p,q] = np.nan + return nu + +def _hu_moments(np.ndarray[np.double_t, ndim=2] nu): + cdef np.ndarray[np.double_t, ndim=1] hu = np.zeros((7,), 'double') + cdef double t0 = nu[3,0] + nu[1,2] + cdef double t1 = nu[2,1] + nu[0,3] + cdef double q0 = t0 * t0 + cdef double q1 = t1 * t1 + cdef double n4 = 4 * nu[1,1] + cdef double s = nu[2,0] + nu[0,2] + cdef double d = nu[2,0] - nu[0,2] + hu[0] = s + hu[1] = d * d + n4 * nu[1,1] + hu[3] = q0 + q1 + hu[5] = d * (q0 - q1) + n4 * t0 * t1 + t0 *= q0 - 3 * q1 + t1 *= 3 * q0 - q1 + q0 = nu[3,0]- 3 * nu[1,2] + q1 = 3 * nu[2,1] - nu[0,3] + hu[2] = q0 * q0 + q1 * q1 + hu[4] = q0 * t0 + q1 * t1 + hu[6] = q1 * t0 - q0 * t1 + return hu + +def regionprops(image, properties='all'): + """Measure properties of labeled image regions. + + Parameters + ---------- + image : NxM ndarray + Labelled input image. + properties : {'all', list, tuple} + Shape measurements to be determined for each labeled image region. + Default is 'all'. The following properties can be determined: + * Area : int + Number of pixels of region. + * BoundingBox : tuple + Bounding box `(minr, minc, maxr, maxc)` + * CentralMoments : 3x3 ndarray + Central moments (translation invariant) Mu_pq up to 3rd order. + * Centroid : array + Centroid coordinate tuple `(r, c)`. + * ConvexArea : int + Number of pixels of convex hull image. + * ConvexImage : HxJ ndarray + Convex hull image which has the same size as bounding box. + * Eccentricity : float + Linear eccentricity of the ellipse that has the same second-moments + as the region (0 <= eccentricity <= 1). + * EquivDiameter : float + The diameter of a circle with the same area as the region. + * EulerNumber : int + Euler number 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)` + * FilledArea : int + Number of pixels of filled region. + * FilledImage : HxJ ndarray + Region image with filled holes which has the same size as bounding + box. + * HuMoments : tuple + Hu moments (translation, scale and rotation invariant). + * Image : HxJ ndarray + Sliced 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. + * MinorAxisLength : float + The length of the minor axis of the ellipse that has the same + normalized second central moments as the region. + * Moments 3x3 ndarray + Spatial moments Mu_pq up to 3rd order. + * NormalizedMoments : 3x3 ndarray + Normalized moments (translation and scale invariant) Nu_pq up to 3rd + order. + * Orientation : float + Angle between the X-axis and the major axis of the ellipse that has + the same second-moments as the region. Ranging from `-pi/2` to + `-pi/2` in counter-clockwise direction. + * Solidity : float + Ratio of pixels in the region to pixels of the convex hull image. + + Returns + ------- + properties : list of dicts + List containing a property dict for each region. The property dicts + contain all the specified properties plus a 'Label' field. + + References + ---------- + B. Jähne. Digitale Bildverarbeitung. Springer-Verlag, + Berlin-Heidelberg, 6. edition, 2005. + T. H. Reiss. Recognizing Planar Objects Using Invariant Image Features, + Bd. 676 von Lecture notes in computer science. Springer, Berlin, 1993. + http://en.wikipedia.org/wiki/Image_moment + + Examples + -------- + >>> from skimage.data import coins + >>> from skimage.morphology import label + >>> img = coins() > 110 + >>> label_img = label(img) + >>> props = regionprops(label_img) + """ + cdef int i, r0, c0, label + cdef np.ndarray[np.double_t, ndim=2] m, mu, nu + cdef double cr, cc, a, b, c + + # determine all properties if nothing specified + if properties == 'all': + properties = PROPS + + props = [] + + objects = ndimage.find_objects(image) + for i, sl in enumerate(objects): + label = i + 1 + + # create property dict for current label + obj_props = {} + props.append(obj_props) + + 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') + + # upper left corner of object bbox + r0 = sl[0].start + c0 = sl[1].start + + m = _moments(array, 3) + # centroid + cr = m[0,1] / m[0,0] + cc = m[1,0] / m[0,0] + mu = _central_moments(array, cr, cc, 3) + nu = _normalized_moments(mu, 3) + + # elements of second order central moment covariance matrix + a = mu[2,0] / mu[0,0] + b = mu[1,1] / mu[0,0] + c = mu[0,2] / mu[0,0] + # eigenvalues of covariance matrix + l1 = fabs(0.5*(a+c-sqrt((a-c)**2+4*b**2))) + l2 = fabs(0.5*(a+c+sqrt((a-c)**2+4*b**2))) + + # cached results which are used by several properties + _filled_image = None + _convex_image = None + + if 'Area' in properties: + obj_props['Area'] = m[0,0] + + if 'BoundingBox' in properties: + obj_props['BoundingBox'] = (r0, c0, sl[0].stop, sl[1].stop) + + if 'Centroid' in properties: + obj_props['Centroid'] = cr+r0, cc+c0 + + if 'CentralMoments' in properties: + obj_props['CentralMoments'] = mu + + if 'ConvexArea' in properties: + if _convex_image is None: + _convex_image = convex_hull_image(array) + obj_props['ConvexArea'] = np.sum(_convex_image) + + if 'ConvexImage' in properties: + if _convex_image is None: + _convex_image = convex_hull_image(array) + obj_props['ConvexImage'] = _convex_image + + if 'Eccentricity' in properties: + # linear eccentricity of ellipse + obj_props['Eccentricity'] = sqrt(1-(fmin(l1, l2)/fmax(l1, l2))**2) + + if 'EquivDiameter' in properties: + obj_props['EquivDiameter'] = sqrt(4 * m[0,0] / PI) + + if 'EulerNumber' in properties: + if _filled_image is None: + _filled_image = ndimage.binary_fill_holes(array, STREL_8) + euler_array = _filled_image != array + _, num = ndimage.label(euler_array, STREL_8) + obj_props['EulerNumber'] = 1 - num + + if 'Extent' in properties: + obj_props['Extent'] = m[0,0] / (array.shape[0] * array.shape[1]) + + if 'HuMoments' in properties: + obj_props['HuMoments'] = _hu_moments(nu) + + if 'Image' in properties: + obj_props['Image'] = array + + if 'FilledArea' in properties: + if _filled_image is None: + _filled_image = ndimage.binary_fill_holes(array, STREL_8) + obj_props['FilledArea'] = np.sum(_filled_image) + + if 'FilledImage' in properties: + if _filled_image is None: + _filled_image = ndimage.binary_fill_holes(array, STREL_8) + obj_props['FilledImage'] = _filled_image + + if 'MinorAxisLength' in properties: + obj_props['MinorAxisLength'] = fmin(l1, l2) + + if 'MajorAxisLength' in properties: + obj_props['MajorAxisLength'] = fmax(l1, l2) + + if 'Moments' in properties: + obj_props['Moments'] = m + + if 'NormalizedMoments' in properties: + obj_props['NormalizedMoments'] = nu + + if 'Orientation' in properties: + obj_props['Orientation'] = - 0.5 * atan2(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) + + return props diff --git a/skimage/measure/setup.py b/skimage/measure/setup.py index f03b9f3b..c16e7a13 100644 --- a/skimage/measure/setup.py +++ b/skimage/measure/setup.py @@ -12,9 +12,12 @@ def configuration(parent_package='', top_path=None): config.add_data_dir('tests') cython(['_find_contours.pyx'], working_path=base_path) + cython(['regionprops.pyx'], working_path=base_path) config.add_extension('_find_contours', sources=['_find_contours.c'], include_dirs=[get_numpy_include_dirs()]) + config.add_extension('regionprops', sources=['regionprops.c'], + include_dirs=[get_numpy_include_dirs()]) return config From 6905603b9810ccf140b1553ef567737639026a9c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Wed, 2 May 2012 22:26:42 +0200 Subject: [PATCH 02/25] adapted to PEP8 whitespace conventions --- skimage/measure/regionprops.pyx | 35 +++++++++++++++++---------------- 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/skimage/measure/regionprops.pyx b/skimage/measure/regionprops.pyx index b591b76f..c0e23d05 100644 --- a/skimage/measure/regionprops.pyx +++ b/skimage/measure/regionprops.pyx @@ -48,34 +48,34 @@ cdef tuple PROPS = ( def _moments(np.ndarray[np.uint8_t, ndim=2] array, int order): cdef int p, q, r, c cdef np.ndarray[np.double_t, ndim=2] m - m = np.zeros((order+1, order+1), 'double') - for p in range(order+1): - for q in range(order+1): + m = np.zeros((order + 1, order + 1), 'double') + for p in range(order + 1): + for q in range(order + 1): for r in range(array.shape[0]): for c in range(array.shape[1]): - m[p,q] += array[r,c] * r**q * c**p + m[p,q] += array[r,c] * r ** q * c ** p return m def _central_moments(np.ndarray[np.uint8_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 - mu = np.zeros((order+1, order+1), 'double') - for p in range(order+1): - for q in range(order+1): + mu = np.zeros((order + 1, order + 1), 'double') + for p in range(order + 1): + for q in range(order + 1): for r in range(array.shape[0]): for c in range(array.shape[1]): - mu[p,q] += array[r,c] * (r-cr)**q * (c-cc)**p + mu[p,q] += array[r,c] * (r - cr) ** q * (c - cc) ** p return mu def _normalized_moments(np.ndarray[np.double_t, ndim=2] mu, int order): cdef int p, q cdef np.ndarray[np.double_t, ndim=2] nu - nu = np.zeros((order+1, order+1), 'double') - for p in range(order+1): - for q in range(order+1): + nu = np.zeros((order + 1, order + 1), 'double') + for p in range(order + 1): + for q in range(order + 1): if p + q >= 2: - nu[p,q] = mu[p,q] / mu[0,0]**((p+q) / 2 + 1) + nu[p,q] = mu[p,q] / mu[0,0]**((p + q) / 2 + 1) else: nu[p,q] = np.nan return nu @@ -224,8 +224,8 @@ def regionprops(image, properties='all'): b = mu[1,1] / mu[0,0] c = mu[0,2] / mu[0,0] # eigenvalues of covariance matrix - l1 = fabs(0.5*(a+c-sqrt((a-c)**2+4*b**2))) - l2 = fabs(0.5*(a+c+sqrt((a-c)**2+4*b**2))) + l1 = fabs(0.5 * (a + c - sqrt((a - c) ** 2 + 4 * b ** 2))) + l2 = fabs(0.5 * (a + c + sqrt((a - c) ** 2 + 4 * b ** 2))) # cached results which are used by several properties _filled_image = None @@ -238,7 +238,7 @@ def regionprops(image, properties='all'): obj_props['BoundingBox'] = (r0, c0, sl[0].stop, sl[1].stop) if 'Centroid' in properties: - obj_props['Centroid'] = cr+r0, cc+c0 + obj_props['Centroid'] = cr + r0, cc + c0 if 'CentralMoments' in properties: obj_props['CentralMoments'] = mu @@ -255,7 +255,8 @@ def regionprops(image, properties='all'): if 'Eccentricity' in properties: # linear eccentricity of ellipse - obj_props['Eccentricity'] = sqrt(1-(fmin(l1, l2)/fmax(l1, l2))**2) + obj_props['Eccentricity'] = \ + sqrt(1 - (fmin(l1, l2) / fmax(l1, l2)) ** 2) if 'EquivDiameter' in properties: obj_props['EquivDiameter'] = sqrt(4 * m[0,0] / PI) @@ -299,7 +300,7 @@ def regionprops(image, properties='all'): obj_props['NormalizedMoments'] = nu if 'Orientation' in properties: - obj_props['Orientation'] = - 0.5 * atan2(2*b, a-c) + obj_props['Orientation'] = - 0.5 * atan2(2 * b, a - c) if 'Solidity' in properties: if _convex_image is None: From cb673ba5164525a50024a3e3ac1199478a659934 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Wed, 2 May 2012 22:44:28 +0200 Subject: [PATCH 03/25] extended and improved description of properties in doc string --- skimage/measure/regionprops.pyx | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/skimage/measure/regionprops.pyx b/skimage/measure/regionprops.pyx index c0e23d05..189cb877 100644 --- a/skimage/measure/regionprops.pyx +++ b/skimage/measure/regionprops.pyx @@ -115,11 +115,14 @@ def regionprops(image, properties='all'): * Area : int Number of pixels of region. * BoundingBox : tuple - Bounding box `(minr, minc, maxr, maxc)` + Bounding box `(min_row, min_col, max_row, max_col)` * CentralMoments : 3x3 ndarray - Central moments (translation invariant) Mu_pq up to 3rd order. + Central moments (translation invariant) up to 3rd order. + .. math:: + \texttt{mu} _{ji} = \sum _{x,y} \left (\texttt{array} (x,y) \\ + \cdot (x - \bar{x} )^j \cdot (y - \bar{y} )^i \right) * Centroid : array - Centroid coordinate tuple `(r, c)`. + Centroid coordinate tuple `(row, col)`. * ConvexArea : int Number of pixels of convex hull image. * ConvexImage : HxJ ndarray @@ -152,9 +155,15 @@ def regionprops(image, properties='all'): normalized second central moments as the region. * Moments 3x3 ndarray Spatial moments Mu_pq up to 3rd order. + .. math:: + \texttt{m} _{ji}= \sum _{x,y} \left (\texttt{array} (x,y) \\ + \cdot x^j \cdot y^i \right) * NormalizedMoments : 3x3 ndarray Normalized moments (translation and scale invariant) Nu_pq up to 3rd order. + .. math:: + \texttt{nu} _{ji} = \\ + \frac{\texttt{mu}_{ji}}{\texttt{m}_{00}^{(i+j)/2+1}} * Orientation : float Angle between the X-axis and the major axis of the ellipse that has the same second-moments as the region. Ranging from `-pi/2` to @@ -254,7 +263,6 @@ def regionprops(image, properties='all'): obj_props['ConvexImage'] = _convex_image if 'Eccentricity' in properties: - # linear eccentricity of ellipse obj_props['Eccentricity'] = \ sqrt(1 - (fmin(l1, l2) / fmax(l1, l2)) ** 2) From e845dc77462d142b8e8d4918bb4ceb112092dec1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Wed, 2 May 2012 22:49:14 +0200 Subject: [PATCH 04/25] moments and central moments no longer use separate functions --- skimage/measure/regionprops.pyx | 13 +------------ 1 file changed, 1 insertion(+), 12 deletions(-) diff --git a/skimage/measure/regionprops.pyx b/skimage/measure/regionprops.pyx index 189cb877..0da3b7c3 100644 --- a/skimage/measure/regionprops.pyx +++ b/skimage/measure/regionprops.pyx @@ -45,17 +45,6 @@ cdef tuple PROPS = ( ) -def _moments(np.ndarray[np.uint8_t, ndim=2] array, int order): - cdef int p, q, r, c - cdef np.ndarray[np.double_t, ndim=2] m - m = np.zeros((order + 1, order + 1), 'double') - for p in range(order + 1): - for q in range(order + 1): - for r in range(array.shape[0]): - for c in range(array.shape[1]): - m[p,q] += array[r,c] * r ** q * c ** p - return m - def _central_moments(np.ndarray[np.uint8_t, ndim=2] array, double cr, double cc, int order): cdef int p, q, r, c @@ -221,7 +210,7 @@ def regionprops(image, properties='all'): r0 = sl[0].start c0 = sl[1].start - m = _moments(array, 3) + m = _central_moments(array, 0, 0, 3) # centroid cr = m[0,1] / m[0,0] cc = m[1,0] / m[0,0] From 658dbec381e3a96b418354446dbd346aae271842 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Wed, 2 May 2012 22:51:06 +0200 Subject: [PATCH 05/25] removed old description of moments in doc string --- skimage/measure/regionprops.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/measure/regionprops.pyx b/skimage/measure/regionprops.pyx index 0da3b7c3..41317971 100644 --- a/skimage/measure/regionprops.pyx +++ b/skimage/measure/regionprops.pyx @@ -143,12 +143,12 @@ def regionprops(image, properties='all'): The length of the minor axis of the ellipse that has the same normalized second central moments as the region. * Moments 3x3 ndarray - Spatial moments Mu_pq up to 3rd order. + Spatial moments up to 3rd order. .. math:: \texttt{m} _{ji}= \sum _{x,y} \left (\texttt{array} (x,y) \\ \cdot x^j \cdot y^i \right) * NormalizedMoments : 3x3 ndarray - Normalized moments (translation and scale invariant) Nu_pq up to 3rd + Normalized moments (translation and scale invariant) up to 3rd order. .. math:: \texttt{nu} _{ji} = \\ From 6b0c0e194850456fc4491b8080987b485762925a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Sat, 5 May 2012 10:08:03 +0200 Subject: [PATCH 06/25] changed reference notes in doc string --- skimage/measure/regionprops.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/measure/regionprops.pyx b/skimage/measure/regionprops.pyx index 41317971..a58ba90b 100644 --- a/skimage/measure/regionprops.pyx +++ b/skimage/measure/regionprops.pyx @@ -168,10 +168,10 @@ def regionprops(image, properties='all'): References ---------- - B. Jähne. Digitale Bildverarbeitung. Springer-Verlag, + B. Jähne. Digital Image Processing. Springer-Verlag, Berlin-Heidelberg, 6. edition, 2005. T. H. Reiss. Recognizing Planar Objects Using Invariant Image Features, - Bd. 676 von Lecture notes in computer science. Springer, Berlin, 1993. + LNICS, p. 676. Springer, Berlin, 1993. http://en.wikipedia.org/wiki/Image_moment Examples From 0027b91e840e5bd3d79b469194959141c8bf015f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Sat, 5 May 2012 13:39:21 +0200 Subject: [PATCH 07/25] added test: labelled image must be of integer type --- skimage/measure/regionprops.pyx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/skimage/measure/regionprops.pyx b/skimage/measure/regionprops.pyx index a58ba90b..25cb1e48 100644 --- a/skimage/measure/regionprops.pyx +++ b/skimage/measure/regionprops.pyx @@ -186,6 +186,9 @@ def regionprops(image, properties='all'): cdef np.ndarray[np.double_t, ndim=2] m, mu, nu cdef double cr, cc, a, b, c + if not np.issubdtype(image.dtype, 'int'): + raise TypeError('labelled image must be of integer dtype') + # determine all properties if nothing specified if properties == 'all': properties = PROPS From e99312432980209be34ad1d503d09abe03bc7270 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Sat, 5 May 2012 15:05:36 +0200 Subject: [PATCH 08/25] most functionality of regionprops now pure python except for moments --- skimage/measure/_regionprops.pyx | 52 +++++++++++ .../{regionprops.pyx => regionprops.py} | 92 +++++-------------- skimage/measure/setup.py | 4 +- 3 files changed, 75 insertions(+), 73 deletions(-) create mode 100644 skimage/measure/_regionprops.pyx rename skimage/measure/{regionprops.pyx => regionprops.py} (75%) diff --git a/skimage/measure/_regionprops.pyx b/skimage/measure/_regionprops.pyx new file mode 100644 index 00000000..d24dbea2 --- /dev/null +++ b/skimage/measure/_regionprops.pyx @@ -0,0 +1,52 @@ +#cython: boundscheck=False +#cython: wraparound=False +#cython: cdivision=True +import numpy as np +cimport numpy as np + + +def central_moments(np.ndarray[np.uint8_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 + mu = np.zeros((order + 1, order + 1), 'double') + for p in range(order + 1): + for q in range(order + 1): + for r in range(array.shape[0]): + for c in range(array.shape[1]): + mu[p,q] += array[r,c] * (r - cr) ** q * (c - cc) ** p + return mu + +def normalized_moments(np.ndarray[np.double_t, ndim=2] mu, int order): + cdef int p, q + cdef np.ndarray[np.double_t, ndim=2] nu + nu = np.zeros((order + 1, order + 1), 'double') + for p in range(order + 1): + for q in range(order + 1): + if p + q >= 2: + nu[p,q] = mu[p,q] / mu[0,0]**((p + q) / 2 + 1) + else: + nu[p,q] = np.nan + return nu + +def hu_moments(np.ndarray[np.double_t, ndim=2] nu): + cdef np.ndarray[np.double_t, ndim=1] hu = np.zeros((7,), 'double') + cdef double t0 = nu[3,0] + nu[1,2] + cdef double t1 = nu[2,1] + nu[0,3] + cdef double q0 = t0 * t0 + cdef double q1 = t1 * t1 + cdef double n4 = 4 * nu[1,1] + cdef double s = nu[2,0] + nu[0,2] + cdef double d = nu[2,0] - nu[0,2] + hu[0] = s + hu[1] = d * d + n4 * nu[1,1] + hu[3] = q0 + q1 + hu[5] = d * (q0 - q1) + n4 * t0 * t1 + t0 *= q0 - 3 * q1 + t1 *= 3 * q0 - q1 + q0 = nu[3,0]- 3 * nu[1,2] + q1 = 3 * nu[2,1] - nu[0,3] + hu[2] = q0 * q0 + q1 * q1 + hu[4] = q0 * t0 + q1 * t1 + hu[6] = q1 * t0 - q0 * t1 + return hu diff --git a/skimage/measure/regionprops.pyx b/skimage/measure/regionprops.py similarity index 75% rename from skimage/measure/regionprops.pyx rename to skimage/measure/regionprops.py index 25cb1e48..5b2286d4 100644 --- a/skimage/measure/regionprops.pyx +++ b/skimage/measure/regionprops.py @@ -1,21 +1,17 @@ -#cython: boundscheck=False -#cython: wraparound=False -#cython: cdivision=True -from scipy import ndimage +# coding: utf-8 +import math import numpy as np -cimport numpy as np -cimport cython -from libc.math cimport sqrt, atan2, fabs, fmin, fmax +from scipy import ndimage from skimage.morphology import convex_hull_image +from . import _regionprops __all__ = ['regionprops'] STREL_8 = np.ones((3, 3), 'int8') -cdef float PI = 3.14159265 -cdef tuple PROPS = ( +PROPS = ( 'Area', 'BoundingBox', 'CentralMoments', @@ -45,52 +41,6 @@ cdef tuple PROPS = ( ) -def _central_moments(np.ndarray[np.uint8_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 - mu = np.zeros((order + 1, order + 1), 'double') - for p in range(order + 1): - for q in range(order + 1): - for r in range(array.shape[0]): - for c in range(array.shape[1]): - mu[p,q] += array[r,c] * (r - cr) ** q * (c - cc) ** p - return mu - -def _normalized_moments(np.ndarray[np.double_t, ndim=2] mu, int order): - cdef int p, q - cdef np.ndarray[np.double_t, ndim=2] nu - nu = np.zeros((order + 1, order + 1), 'double') - for p in range(order + 1): - for q in range(order + 1): - if p + q >= 2: - nu[p,q] = mu[p,q] / mu[0,0]**((p + q) / 2 + 1) - else: - nu[p,q] = np.nan - return nu - -def _hu_moments(np.ndarray[np.double_t, ndim=2] nu): - cdef np.ndarray[np.double_t, ndim=1] hu = np.zeros((7,), 'double') - cdef double t0 = nu[3,0] + nu[1,2] - cdef double t1 = nu[2,1] + nu[0,3] - cdef double q0 = t0 * t0 - cdef double q1 = t1 * t1 - cdef double n4 = 4 * nu[1,1] - cdef double s = nu[2,0] + nu[0,2] - cdef double d = nu[2,0] - nu[0,2] - hu[0] = s - hu[1] = d * d + n4 * nu[1,1] - hu[3] = q0 + q1 - hu[5] = d * (q0 - q1) + n4 * t0 * t1 - t0 *= q0 - 3 * q1 - t1 *= 3 * q0 - q1 - q0 = nu[3,0]- 3 * nu[1,2] - q1 = 3 * nu[2,1] - nu[0,3] - hu[2] = q0 * q0 + q1 * q1 - hu[4] = q0 * t0 + q1 * t1 - hu[6] = q1 * t0 - q0 * t1 - return hu - def regionprops(image, properties='all'): """Measure properties of labeled image regions. @@ -182,10 +132,6 @@ def regionprops(image, properties='all'): >>> label_img = label(img) >>> props = regionprops(label_img) """ - cdef int i, r0, c0, label - cdef np.ndarray[np.double_t, ndim=2] m, mu, nu - cdef double cr, cc, a, b, c - if not np.issubdtype(image.dtype, 'int'): raise TypeError('labelled image must be of integer dtype') @@ -213,24 +159,24 @@ def regionprops(image, properties='all'): r0 = sl[0].start c0 = sl[1].start - m = _central_moments(array, 0, 0, 3) + m = _regionprops._central_moments(array, 0, 0, 3) # centroid cr = m[0,1] / m[0,0] cc = m[1,0] / m[0,0] - mu = _central_moments(array, cr, cc, 3) - nu = _normalized_moments(mu, 3) + mu = _regionprops.central_moments(array, cr, cc, 3) # elements of second order central moment covariance matrix a = mu[2,0] / mu[0,0] b = mu[1,1] / mu[0,0] c = mu[0,2] / mu[0,0] # eigenvalues of covariance matrix - l1 = fabs(0.5 * (a + c - sqrt((a - c) ** 2 + 4 * b ** 2))) - l2 = fabs(0.5 * (a + c + sqrt((a - c) ** 2 + 4 * b ** 2))) + l1 = abs(0.5 * (a + c - math.sqrt((a - c) ** 2 + 4 * b ** 2))) + l2 = abs(0.5 * (a + c + math.sqrt((a - c) ** 2 + 4 * b ** 2))) # cached results which are used by several properties _filled_image = None _convex_image = None + _nu = None if 'Area' in properties: obj_props['Area'] = m[0,0] @@ -256,10 +202,10 @@ def regionprops(image, properties='all'): if 'Eccentricity' in properties: obj_props['Eccentricity'] = \ - sqrt(1 - (fmin(l1, l2) / fmax(l1, l2)) ** 2) + math.sqrt(1 - (min(l1, l2) / max(l1, l2)) ** 2) if 'EquivDiameter' in properties: - obj_props['EquivDiameter'] = sqrt(4 * m[0,0] / PI) + obj_props['EquivDiameter'] = math.sqrt(4 * m[0,0] / math.pi) if 'EulerNumber' in properties: if _filled_image is None: @@ -272,7 +218,9 @@ def regionprops(image, properties='all'): obj_props['Extent'] = m[0,0] / (array.shape[0] * array.shape[1]) if 'HuMoments' in properties: - obj_props['HuMoments'] = _hu_moments(nu) + if _nu is None: + _nu = _regionprops.normalized_moments(mu, 3) + obj_props['HuMoments'] = _regionprops.hu_moments(_nu) if 'Image' in properties: obj_props['Image'] = array @@ -288,19 +236,21 @@ def regionprops(image, properties='all'): obj_props['FilledImage'] = _filled_image if 'MinorAxisLength' in properties: - obj_props['MinorAxisLength'] = fmin(l1, l2) + obj_props['MinorAxisLength'] = min(l1, l2) if 'MajorAxisLength' in properties: - obj_props['MajorAxisLength'] = fmax(l1, l2) + obj_props['MajorAxisLength'] = max(l1, l2) if 'Moments' in properties: obj_props['Moments'] = m if 'NormalizedMoments' in properties: - obj_props['NormalizedMoments'] = nu + if _nu is None: + _nu = _regionprops.normalized_moments(mu, 3) + obj_props['NormalizedMoments'] = _nu if 'Orientation' in properties: - obj_props['Orientation'] = - 0.5 * atan2(2 * b, a - c) + obj_props['Orientation'] = - 0.5 * math.atan2(2 * b, a - c) if 'Solidity' in properties: if _convex_image is None: diff --git a/skimage/measure/setup.py b/skimage/measure/setup.py index c16e7a13..267487c1 100644 --- a/skimage/measure/setup.py +++ b/skimage/measure/setup.py @@ -12,11 +12,11 @@ def configuration(parent_package='', top_path=None): config.add_data_dir('tests') cython(['_find_contours.pyx'], working_path=base_path) - cython(['regionprops.pyx'], working_path=base_path) + cython(['_regionprops.pyx'], working_path=base_path) config.add_extension('_find_contours', sources=['_find_contours.c'], include_dirs=[get_numpy_include_dirs()]) - config.add_extension('regionprops', sources=['regionprops.c'], + config.add_extension('_regionprops', sources=['_regionprops.c'], include_dirs=[get_numpy_include_dirs()]) return config From 1b2bed75016d05146df490af9f2b20a582c60bf8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Fri, 11 May 2012 21:30:18 +0200 Subject: [PATCH 09/25] renamed regionprops source files for better separation between function name and source file --- skimage/measure/__init__.py | 2 +- skimage/measure/{_regionprops.pyx => _moments.pyx} | 0 skimage/measure/{regionprops.py => _regionprops.py} | 12 ++++++------ skimage/measure/setup.py | 4 ++-- 4 files changed, 9 insertions(+), 9 deletions(-) rename skimage/measure/{_regionprops.pyx => _moments.pyx} (100%) rename skimage/measure/{regionprops.py => _regionprops.py} (96%) diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py index be4a1b43..422db569 100755 --- a/skimage/measure/__init__.py +++ b/skimage/measure/__init__.py @@ -1,2 +1,2 @@ from .find_contours import find_contours -from .regionprops import regionprops +from ._regionprops import regionprops diff --git a/skimage/measure/_regionprops.pyx b/skimage/measure/_moments.pyx similarity index 100% rename from skimage/measure/_regionprops.pyx rename to skimage/measure/_moments.pyx diff --git a/skimage/measure/regionprops.py b/skimage/measure/_regionprops.py similarity index 96% rename from skimage/measure/regionprops.py rename to skimage/measure/_regionprops.py index 5b2286d4..d5c1c211 100644 --- a/skimage/measure/regionprops.py +++ b/skimage/measure/_regionprops.py @@ -4,7 +4,7 @@ import numpy as np from scipy import ndimage from skimage.morphology import convex_hull_image -from . import _regionprops +from . import _moments __all__ = ['regionprops'] @@ -159,11 +159,11 @@ def regionprops(image, properties='all'): r0 = sl[0].start c0 = sl[1].start - m = _regionprops._central_moments(array, 0, 0, 3) + m = _moments.central_moments(array, 0, 0, 3) # centroid cr = m[0,1] / m[0,0] cc = m[1,0] / m[0,0] - mu = _regionprops.central_moments(array, cr, cc, 3) + mu = _moments.central_moments(array, cr, cc, 3) # elements of second order central moment covariance matrix a = mu[2,0] / mu[0,0] @@ -219,8 +219,8 @@ def regionprops(image, properties='all'): if 'HuMoments' in properties: if _nu is None: - _nu = _regionprops.normalized_moments(mu, 3) - obj_props['HuMoments'] = _regionprops.hu_moments(_nu) + _nu = _moments.normalized_moments(mu, 3) + obj_props['HuMoments'] = _moments.hu_moments(_nu) if 'Image' in properties: obj_props['Image'] = array @@ -246,7 +246,7 @@ def regionprops(image, properties='all'): if 'NormalizedMoments' in properties: if _nu is None: - _nu = _regionprops.normalized_moments(mu, 3) + _nu = _moments.normalized_moments(mu, 3) obj_props['NormalizedMoments'] = _nu if 'Orientation' in properties: diff --git a/skimage/measure/setup.py b/skimage/measure/setup.py index 267487c1..4b02a1d1 100644 --- a/skimage/measure/setup.py +++ b/skimage/measure/setup.py @@ -12,11 +12,11 @@ def configuration(parent_package='', top_path=None): config.add_data_dir('tests') cython(['_find_contours.pyx'], working_path=base_path) - cython(['_regionprops.pyx'], working_path=base_path) + cython(['_moments.pyx'], working_path=base_path) config.add_extension('_find_contours', sources=['_find_contours.c'], include_dirs=[get_numpy_include_dirs()]) - config.add_extension('_regionprops', sources=['_regionprops.c'], + config.add_extension('_moments', sources=['_moments.c'], include_dirs=[get_numpy_include_dirs()]) return config From 4eb3028a64b0d5ca7124de7f9453d277eb796918 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Fri, 11 May 2012 21:33:59 +0200 Subject: [PATCH 10/25] fixed typo and improved example in doc string --- skimage/measure/_regionprops.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/skimage/measure/_regionprops.py b/skimage/measure/_regionprops.py index d5c1c211..9cfa1ac9 100644 --- a/skimage/measure/_regionprops.py +++ b/skimage/measure/_regionprops.py @@ -42,14 +42,14 @@ PROPS = ( def regionprops(image, properties='all'): - """Measure properties of labeled image regions. + """Measure properties of labelled image regions. Parameters ---------- image : NxM ndarray Labelled input image. properties : {'all', list, tuple} - Shape measurements to be determined for each labeled image region. + Shape measurements to be determined for each labelled image region. Default is 'all'. The following properties can be determined: * Area : int Number of pixels of region. @@ -131,6 +131,7 @@ def regionprops(image, properties='all'): >>> img = coins() > 110 >>> label_img = label(img) >>> props = regionprops(label_img) + >>> props[0]['Centroid'] # centroid of first labelled object """ if not np.issubdtype(image.dtype, 'int'): raise TypeError('labelled image must be of integer dtype') From 19614cc8d745e41ad79efd6460168328079544b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Fri, 11 May 2012 23:00:08 +0200 Subject: [PATCH 11/25] euler number now behaves like MATLAB implementation --- skimage/measure/_regionprops.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/measure/_regionprops.py b/skimage/measure/_regionprops.py index 9cfa1ac9..f3af9e47 100644 --- a/skimage/measure/_regionprops.py +++ b/skimage/measure/_regionprops.py @@ -213,7 +213,7 @@ def regionprops(image, properties='all'): _filled_image = ndimage.binary_fill_holes(array, STREL_8) euler_array = _filled_image != array _, num = ndimage.label(euler_array, STREL_8) - obj_props['EulerNumber'] = 1 - num + obj_props['EulerNumber'] = - num if 'Extent' in properties: obj_props['Extent'] = m[0,0] / (array.shape[0] * array.shape[1]) From 04e444d90f686d45cf57226280ef785fa7e5e1e4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Fri, 11 May 2012 23:00:33 +0200 Subject: [PATCH 12/25] added test cases for regionprops --- skimage/measure/tests/test_regionprops.py | 177 ++++++++++++++++++++++ 1 file changed, 177 insertions(+) create mode 100644 skimage/measure/tests/test_regionprops.py diff --git a/skimage/measure/tests/test_regionprops.py b/skimage/measure/tests/test_regionprops.py new file mode 100644 index 00000000..df526ee2 --- /dev/null +++ b/skimage/measure/tests/test_regionprops.py @@ -0,0 +1,177 @@ +from numpy.testing import assert_array_equal, assert_almost_equal, \ + assert_array_almost_equal +import numpy as np + +from skimage.measure import regionprops + + +SAMPLE = np.array( + [[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0], + [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]] +) + + +def test_area(): + area = regionprops(SAMPLE, ['Area'])[0]['Area'] + assert area == np.sum(SAMPLE) + +def test_bbox(): + bbox = regionprops(SAMPLE, ['BoundingBox'])[0]['BoundingBox'] + assert_array_almost_equal(bbox, (0, 0, SAMPLE.shape[0], SAMPLE.shape[1])) + + SAMPLE_mod = SAMPLE.copy() + SAMPLE_mod[:,-1] = 0 + bbox = regionprops(SAMPLE_mod, ['BoundingBox'])[0]['BoundingBox'] + assert_array_almost_equal(bbox, (0, 0, SAMPLE.shape[0], SAMPLE.shape[1]-1)) + +def test_central_moments(): + mu = regionprops(SAMPLE, ['CentralMoments'])[0]['CentralMoments'] + #: determined with OpenCV + assert_almost_equal(mu[0,2], 436.00000000000045) + + # bug in OpenCV and Wikipedia? + assert_almost_equal(mu[0,3], -0.016762262088312843) + + 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(): + centroid = regionprops(SAMPLE, ['Centroid'])[0]['Centroid'] + # determined with MATLAB + assert_array_almost_equal(centroid, (5.66666666666666, 9.444444444444444)) + +def test_convex_area(): + area = regionprops(SAMPLE, ['ConvexArea'])[0]['ConvexArea'] + # determined with MATLAB + assert area == 124 + +def test_convex_image(): + img = regionprops(SAMPLE, ['ConvexImage'])[0]['ConvexImage'] + # determined with MATLAB + ref = np.array( + [[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0], + [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]] + ) + assert_array_equal(img, ref) + +def test_eccentricity(): + eps = regionprops(SAMPLE, ['Eccentricity'])[0]['Eccentricity'] + # MATLAB has different interpretation of ellipse than found in literature, + # here implemented as found in literature + assert_almost_equal(eps, 0.941726665966) + +def test_equiv_diameter(): + diameter = regionprops(SAMPLE, ['EquivDiameter'])[0]['EquivDiameter'] + # determined with MATLAB + assert_almost_equal(diameter, 9.57461472963) + +def test_euler_number(): + en = regionprops(SAMPLE, ['EulerNumber'])[0]['EulerNumber'] + assert en == 0 + + SAMPLE_mod = SAMPLE.copy() + SAMPLE_mod[7,-3] = 0 + en = regionprops(SAMPLE_mod, ['EulerNumber'])[0]['EulerNumber'] + assert en == -1 + +def test_extent(): + extent = regionprops(SAMPLE, ['Extent'])[0]['Extent'] + assert_almost_equal(extent, 0.4) + +def test_hu_moments(): + hu = regionprops(SAMPLE, ['HuMoments'])[0]['HuMoments'] + # determined with OpenCV + ref = np.array( + [[ 3.27117627e-01], + [ 2.63869194e-02], + [ 1.86845507e-02], + [ 2.47503247e-03], + [-6.25438993e-06], + [-2.02071583e-04], + [ 1.56259006e-05]] + ) + # bug in OpenCV caused in Central Moments calculation? + assert_array_almost_equal(hu, ref) + +def test_image(): + img = regionprops(SAMPLE, ['Image'])[0]['Image'] + assert_array_equal(img, SAMPLE) + +def test_filled_area(): + area = regionprops(SAMPLE, ['FilledArea'])[0]['FilledArea'] + assert area == np.sum(SAMPLE) + + SAMPLE_mod = SAMPLE.copy() + SAMPLE_mod[7,-3] = 0 + 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, 5.92837619822) + +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, 17.6240929376) + +def test_moments(): + m = regionprops(SAMPLE, ['Moments'])[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) + +def test_normalized_moments(): + nu = regionprops(SAMPLE, ['NormalizedMoments'])[0]['NormalizedMoments'] + #: 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) + +def test_orientation(): + orientation = regionprops(SAMPLE, ['Orientation'])[0]['Orientation'] + # determined with MATLAB + assert_almost_equal(orientation, 0.10446844651921) + +def test_solidity(): + solidity = regionprops(SAMPLE, ['Solidity'])[0]['Solidity'] + # determined with MATLAB + assert_almost_equal(solidity, 0.580645161290323) + + +if __name__ == "__main__": + from numpy.testing import run_module_suite + run_module_suite() From 87ddbbeabe719077b629ee454953bcfb35f83644 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Sat, 12 May 2012 10:05:14 +0200 Subject: [PATCH 13/25] more constistent way of determining ellipse parameters --- skimage/measure/_regionprops.py | 59 ++++++++++++++++------- skimage/measure/tests/test_regionprops.py | 8 ++- 2 files changed, 44 insertions(+), 23 deletions(-) diff --git a/skimage/measure/_regionprops.py b/skimage/measure/_regionprops.py index f3af9e47..765f68de 100644 --- a/skimage/measure/_regionprops.py +++ b/skimage/measure/_regionprops.py @@ -1,5 +1,5 @@ # coding: utf-8 -import math +from math import sqrt, atan, pi as PI import numpy as np from scipy import ndimage @@ -67,8 +67,15 @@ def regionprops(image, properties='all'): * ConvexImage : HxJ ndarray Convex hull image which has the same size as bounding box. * Eccentricity : float - Linear eccentricity of the ellipse that has the same second-moments - as the region (0 <= eccentricity <= 1). + Eccentricity of the ellipse that has the same second-moments as the + region (1 <= Eccentricity < Inf), where Eccentricity = 1 corresponds + to a circular disk and elongated regions have Eccentricity > 1. + .. math:: + a_1 = \mu _{2,0} + \mu _{0,2} + \sqrt{(\mu _{2,0} - \\ + \mu _{0,2})^2 + 4 {\mu _{1,1}}^2} + a_2 = \mu _{2,0} + \mu _{0,2} - \sqrt{(\mu _{2,0} - \\ + \mu _{0,2})^2 + 4 {\mu _{1,1}}^2} + Ecc = \frac{a_1}{a_2} * EquivDiameter : float The diameter of a circle with the same area as the region. * EulerNumber : int @@ -89,9 +96,17 @@ def regionprops(image, properties='all'): * MajorAxisLength : float The length of the major axis of the ellipse that has the same normalized second central moments as the region. + .. math:: + a_1 = \mu _{2,0} + \mu _{0,2} + \sqrt{(\mu _{2,0} - \\ + \mu _{0,2})^2 + 4 {\mu _{1,1}}^2} + A = \sqrt{\frac{2 a_1}{\mu _{0,0}}} * MinorAxisLength : float The length of the minor axis of the ellipse that has the same normalized second central moments as the region. + .. math:: + a_2 = \mu _{2,0} + \mu _{0,2} - \sqrt{(\mu _{2,0} - \\ + \mu _{0,2})^2 + 4 {\mu _{1,1}}^2} + A = \sqrt{\frac{2 a_2}{\mu _{0,0}}} * Moments 3x3 ndarray Spatial moments up to 3rd order. .. math:: @@ -118,6 +133,8 @@ def regionprops(image, properties='all'): References ---------- + Wilhelm Burger, Mark Burge. Principles of Digital Image Processing: Core + Algorithms. Springer-Verlag, London, 2009. B. Jähne. Digital Image Processing. Springer-Verlag, Berlin-Heidelberg, 6. edition, 2005. T. H. Reiss. Recognizing Planar Objects Using Invariant Image Features, @@ -166,18 +183,17 @@ def regionprops(image, properties='all'): cc = m[1,0] / m[0,0] mu = _moments.central_moments(array, cr, cc, 3) - # elements of second order central moment covariance matrix - a = mu[2,0] / mu[0,0] - b = mu[1,1] / mu[0,0] - c = mu[0,2] / mu[0,0] - # eigenvalues of covariance matrix - l1 = abs(0.5 * (a + c - math.sqrt((a - c) ** 2 + 4 * b ** 2))) - l2 = abs(0.5 * (a + c + math.sqrt((a - c) ** 2 + 4 * b ** 2))) + # elements of the inertia tensor + a = mu[2,0] + b = mu[1,1] + c = mu[0,2] # cached results which are used by several properties _filled_image = None _convex_image = None _nu = None + _a1 = None + _a2 = None if 'Area' in properties: obj_props['Area'] = m[0,0] @@ -202,11 +218,14 @@ def regionprops(image, properties='all'): obj_props['ConvexImage'] = _convex_image if 'Eccentricity' in properties: - obj_props['Eccentricity'] = \ - math.sqrt(1 - (min(l1, l2) / max(l1, l2)) ** 2) + if _a2 is None: + _a2 = a + c - sqrt((a - c)**2 + 4 * b ** 2) + if _a1 is None: + _a1 = a + c + sqrt((a - c)**2 + 4 * b ** 2) + obj_props['Eccentricity'] = _a1 / _a2 if 'EquivDiameter' in properties: - obj_props['EquivDiameter'] = math.sqrt(4 * m[0,0] / math.pi) + obj_props['EquivDiameter'] = sqrt(4 * m[0,0] / PI) if 'EulerNumber' in properties: if _filled_image is None: @@ -236,11 +255,15 @@ def regionprops(image, properties='all'): _filled_image = ndimage.binary_fill_holes(array, STREL_8) obj_props['FilledImage'] = _filled_image - if 'MinorAxisLength' in properties: - obj_props['MinorAxisLength'] = min(l1, l2) - if 'MajorAxisLength' in properties: - obj_props['MajorAxisLength'] = max(l1, l2) + if _a1 is None: + _a1 = a + c + sqrt((a - c)**2 + 4 * b ** 2) + obj_props['MajorAxisLength'] = sqrt(2 * _a1 / m[0,0]) + + if 'MinorAxisLength' in properties: + if _a2 is None: + _a2 = a1 = a + c - sqrt((a - c)**2 + 4 * b ** 2) + obj_props['MinorAxisLength'] = sqrt(2 * _a2 / m[0,0]) if 'Moments' in properties: obj_props['Moments'] = m @@ -251,7 +274,7 @@ def regionprops(image, properties='all'): obj_props['NormalizedMoments'] = _nu if 'Orientation' in properties: - obj_props['Orientation'] = - 0.5 * math.atan2(2 * b, a - c) + obj_props['Orientation'] = - 0.5 * atan(2 * b / (a - c)) if 'Solidity' in properties: if _convex_image is None: diff --git a/skimage/measure/tests/test_regionprops.py b/skimage/measure/tests/test_regionprops.py index df526ee2..313f9c70 100644 --- a/skimage/measure/tests/test_regionprops.py +++ b/skimage/measure/tests/test_regionprops.py @@ -75,9 +75,7 @@ def test_convex_image(): def test_eccentricity(): eps = regionprops(SAMPLE, ['Eccentricity'])[0]['Eccentricity'] - # MATLAB has different interpretation of ellipse than found in literature, - # here implemented as found in literature - assert_almost_equal(eps, 0.941726665966) + assert_almost_equal(eps, 2.9728364645382) def test_equiv_diameter(): diameter = regionprops(SAMPLE, ['EquivDiameter'])[0]['EquivDiameter'] @@ -129,13 +127,13 @@ 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, 5.92837619822) + assert_almost_equal(length, 4.869651403631) 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, 17.6240929376) + assert_almost_equal(length, 8.396211749968) def test_moments(): m = regionprops(SAMPLE, ['Moments'])[0]['Moments'] From 8ec52869c7b7c632c353abe0bc08329e3ad4a907 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Sat, 12 May 2012 11:16:53 +0200 Subject: [PATCH 14/25] fix wrong indentation --- skimage/measure/tests/test_regionprops.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/skimage/measure/tests/test_regionprops.py b/skimage/measure/tests/test_regionprops.py index 313f9c70..ef70ce2d 100644 --- a/skimage/measure/tests/test_regionprops.py +++ b/skimage/measure/tests/test_regionprops.py @@ -160,9 +160,9 @@ def test_normalized_moments(): assert_almost_equal(nu[3,0], -0.017278118992041805) def test_orientation(): - orientation = regionprops(SAMPLE, ['Orientation'])[0]['Orientation'] - # determined with MATLAB - assert_almost_equal(orientation, 0.10446844651921) + orientation = regionprops(SAMPLE, ['Orientation'])[0]['Orientation'] + # determined with MATLAB + assert_almost_equal(orientation, 0.10446844651921) def test_solidity(): solidity = regionprops(SAMPLE, ['Solidity'])[0]['Solidity'] From b04138c39cbd291ef0e44f971f9064d73893d848 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Sun, 13 May 2012 18:23:20 +0200 Subject: [PATCH 15/25] parameters of ellipse match Matlab results --- skimage/measure/_regionprops.py | 43 +++++++---------------- skimage/measure/tests/test_regionprops.py | 6 ++-- 2 files changed, 15 insertions(+), 34 deletions(-) diff --git a/skimage/measure/_regionprops.py b/skimage/measure/_regionprops.py index 765f68de..d3d68745 100644 --- a/skimage/measure/_regionprops.py +++ b/skimage/measure/_regionprops.py @@ -68,14 +68,8 @@ def regionprops(image, properties='all'): 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 (1 <= Eccentricity < Inf), where Eccentricity = 1 corresponds - to a circular disk and elongated regions have Eccentricity > 1. - .. math:: - a_1 = \mu _{2,0} + \mu _{0,2} + \sqrt{(\mu _{2,0} - \\ - \mu _{0,2})^2 + 4 {\mu _{1,1}}^2} - a_2 = \mu _{2,0} + \mu _{0,2} - \sqrt{(\mu _{2,0} - \\ - \mu _{0,2})^2 + 4 {\mu _{1,1}}^2} - Ecc = \frac{a_1}{a_2} + region. The eccentricity is the ratio of the distance between its + minor and major axis length. The value is between 0 and 1. * EquivDiameter : float The diameter of a circle with the same area as the region. * EulerNumber : int @@ -96,17 +90,9 @@ def regionprops(image, properties='all'): * MajorAxisLength : float The length of the major axis of the ellipse that has the same normalized second central moments as the region. - .. math:: - a_1 = \mu _{2,0} + \mu _{0,2} + \sqrt{(\mu _{2,0} - \\ - \mu _{0,2})^2 + 4 {\mu _{1,1}}^2} - A = \sqrt{\frac{2 a_1}{\mu _{0,0}}} * MinorAxisLength : float The length of the minor axis of the ellipse that has the same normalized second central moments as the region. - .. math:: - a_2 = \mu _{2,0} + \mu _{0,2} - \sqrt{(\mu _{2,0} - \\ - \mu _{0,2})^2 + 4 {\mu _{1,1}}^2} - A = \sqrt{\frac{2 a_2}{\mu _{0,0}}} * Moments 3x3 ndarray Spatial moments up to 3rd order. .. math:: @@ -183,10 +169,13 @@ def regionprops(image, properties='all'): cc = m[1,0] / m[0,0] mu = _moments.central_moments(array, cr, cc, 3) - # elements of the inertia tensor - a = mu[2,0] - b = mu[1,1] - c = mu[0,2] + #: elements of the inertia tensor [a b; b c] + a = mu[2,0] / mu[0,0] + b = mu[1,1] / mu[0,0] + c = mu[0,2] / mu[0,0] + #: eigen values of inertia tensor + l1 = (a + c) / 2 + sqrt(4 * b ** 2 + (a - c) ** 2) / 2 + l2 = (a + c) / 2 - sqrt(4 * b ** 2 + (a - c) ** 2) / 2 # cached results which are used by several properties _filled_image = None @@ -218,11 +207,7 @@ def regionprops(image, properties='all'): obj_props['ConvexImage'] = _convex_image if 'Eccentricity' in properties: - if _a2 is None: - _a2 = a + c - sqrt((a - c)**2 + 4 * b ** 2) - if _a1 is None: - _a1 = a + c + sqrt((a - c)**2 + 4 * b ** 2) - obj_props['Eccentricity'] = _a1 / _a2 + obj_props['Eccentricity'] = sqrt(1 - l2 / l1) if 'EquivDiameter' in properties: obj_props['EquivDiameter'] = sqrt(4 * m[0,0] / PI) @@ -256,14 +241,10 @@ def regionprops(image, properties='all'): obj_props['FilledImage'] = _filled_image if 'MajorAxisLength' in properties: - if _a1 is None: - _a1 = a + c + sqrt((a - c)**2 + 4 * b ** 2) - obj_props['MajorAxisLength'] = sqrt(2 * _a1 / m[0,0]) + obj_props['MajorAxisLength'] = 4 * sqrt(l1) if 'MinorAxisLength' in properties: - if _a2 is None: - _a2 = a1 = a + c - sqrt((a - c)**2 + 4 * b ** 2) - obj_props['MinorAxisLength'] = sqrt(2 * _a2 / m[0,0]) + obj_props['MinorAxisLength'] = 4 * sqrt(l2) if 'Moments' in properties: obj_props['Moments'] = m diff --git a/skimage/measure/tests/test_regionprops.py b/skimage/measure/tests/test_regionprops.py index ef70ce2d..cb9f3035 100644 --- a/skimage/measure/tests/test_regionprops.py +++ b/skimage/measure/tests/test_regionprops.py @@ -75,7 +75,7 @@ def test_convex_image(): def test_eccentricity(): eps = regionprops(SAMPLE, ['Eccentricity'])[0]['Eccentricity'] - assert_almost_equal(eps, 2.9728364645382) + assert_almost_equal(eps, 0.814629313427) def test_equiv_diameter(): diameter = regionprops(SAMPLE, ['EquivDiameter'])[0]['EquivDiameter'] @@ -127,13 +127,13 @@ 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, 4.869651403631) + 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, 8.396211749968) + assert_almost_equal(length, 16.7924234999) def test_moments(): m = regionprops(SAMPLE, ['Moments'])[0]['Moments'] From 06759d2bdfc5d5861578ade18db8b114359cab47 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Sun, 13 May 2012 18:27:25 +0200 Subject: [PATCH 16/25] make all test cases of regionprops work --- skimage/measure/tests/test_regionprops.py | 25 ++++++++++------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/skimage/measure/tests/test_regionprops.py b/skimage/measure/tests/test_regionprops.py index cb9f3035..780714b9 100644 --- a/skimage/measure/tests/test_regionprops.py +++ b/skimage/measure/tests/test_regionprops.py @@ -36,10 +36,8 @@ def test_central_moments(): mu = regionprops(SAMPLE, ['CentralMoments'])[0]['CentralMoments'] #: determined with OpenCV assert_almost_equal(mu[0,2], 436.00000000000045) - - # bug in OpenCV and Wikipedia? - assert_almost_equal(mu[0,3], -0.016762262088312843) - + # 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) @@ -97,16 +95,15 @@ def test_extent(): def test_hu_moments(): hu = regionprops(SAMPLE, ['HuMoments'])[0]['HuMoments'] - # determined with OpenCV - ref = np.array( - [[ 3.27117627e-01], - [ 2.63869194e-02], - [ 1.86845507e-02], - [ 2.47503247e-03], - [-6.25438993e-06], - [-2.02071583e-04], - [ 1.56259006e-05]] - ) + ref = np.array([ + 3.27117627e-01, + 2.63869194e-02, + 2.35390060e-02, + 1.23151193e-03, + 1.38882330e-06, + -2.72586158e-05, + 6.48350653e-06 + ]) # bug in OpenCV caused in Central Moments calculation? assert_array_almost_equal(hu, ref) From c414d8df074746fe44f8ea1d81fee5d0700dfcbf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Sun, 13 May 2012 19:46:41 +0200 Subject: [PATCH 17/25] update contribution for regionprops --- CONTRIBUTORS.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 9907ba6f..baaf064b 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -106,3 +106,4 @@ - Johannes Schönberger Polygon, circle and ellipse drawing functions Adaptive thresholding + Implementation of Matlab's `regionprops` From 8d1f2dc38faf9f7fdb9624aed76f27496d52dc8f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Mon, 14 May 2012 09:59:02 +0200 Subject: [PATCH 18/25] remove unused variables --- skimage/measure/_regionprops.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/skimage/measure/_regionprops.py b/skimage/measure/_regionprops.py index d3d68745..9afe3174 100644 --- a/skimage/measure/_regionprops.py +++ b/skimage/measure/_regionprops.py @@ -181,8 +181,6 @@ def regionprops(image, properties='all'): _filled_image = None _convex_image = None _nu = None - _a1 = None - _a2 = None if 'Area' in properties: obj_props['Area'] = m[0,0] From e30fcfc493adaa21c360905070afb6867c8b3579 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Mon, 14 May 2012 10:01:15 +0200 Subject: [PATCH 19/25] make matrix shape more readable in doc string --- skimage/measure/_regionprops.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/skimage/measure/_regionprops.py b/skimage/measure/_regionprops.py index 9afe3174..dc9d2ff1 100644 --- a/skimage/measure/_regionprops.py +++ b/skimage/measure/_regionprops.py @@ -46,7 +46,7 @@ def regionprops(image, properties='all'): Parameters ---------- - image : NxM ndarray + image : N x M ndarray Labelled input image. properties : {'all', list, tuple} Shape measurements to be determined for each labelled image region. @@ -64,7 +64,7 @@ def regionprops(image, properties='all'): Centroid coordinate tuple `(row, col)`. * ConvexArea : int Number of pixels of convex hull image. - * ConvexImage : HxJ ndarray + * ConvexImage : H x J ndarray 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 @@ -80,12 +80,12 @@ def regionprops(image, properties='all'): Computed as `Area / (rows*cols)` * FilledArea : int Number of pixels of filled region. - * FilledImage : HxJ ndarray + * FilledImage : H x J ndarray Region image with filled holes which has the same size as bounding box. * HuMoments : tuple Hu moments (translation, scale and rotation invariant). - * Image : HxJ ndarray + * Image : H x J ndarray Sliced 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 From 627ac3cbb9a37fc05754b62a9281a7b6f1be8fc2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Mon, 14 May 2012 10:02:59 +0200 Subject: [PATCH 20/25] change reference note in doc string --- skimage/measure/_regionprops.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/measure/_regionprops.py b/skimage/measure/_regionprops.py index dc9d2ff1..b6199524 100644 --- a/skimage/measure/_regionprops.py +++ b/skimage/measure/_regionprops.py @@ -124,7 +124,7 @@ def regionprops(image, properties='all'): B. Jähne. Digital Image Processing. Springer-Verlag, Berlin-Heidelberg, 6. edition, 2005. T. H. Reiss. Recognizing Planar Objects Using Invariant Image Features, - LNICS, p. 676. Springer, Berlin, 1993. + from Lecture notes in computer science, p. 676. Springer, Berlin, 1993. http://en.wikipedia.org/wiki/Image_moment Examples From 0c225219162f224ede69719b487abafed0ea90c1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Mon, 14 May 2012 10:05:19 +0200 Subject: [PATCH 21/25] reduce parameter choice of properties --- skimage/measure/_regionprops.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/measure/_regionprops.py b/skimage/measure/_regionprops.py index b6199524..78779192 100644 --- a/skimage/measure/_regionprops.py +++ b/skimage/measure/_regionprops.py @@ -48,7 +48,7 @@ def regionprops(image, properties='all'): ---------- image : N x M ndarray Labelled input image. - properties : {'all', list, tuple} + properties : {'all', list} Shape measurements to be determined for each labelled image region. Default is 'all'. The following properties can be determined: * Area : int From 58d07c0a05bb654d16538d09120ba50b1c2268e5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Mon, 14 May 2012 16:22:14 +0200 Subject: [PATCH 22/25] replace inline latex equations in doc string with plain text --- skimage/measure/_regionprops.py | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/skimage/measure/_regionprops.py b/skimage/measure/_regionprops.py index 78779192..5c840410 100644 --- a/skimage/measure/_regionprops.py +++ b/skimage/measure/_regionprops.py @@ -55,11 +55,11 @@ def regionprops(image, properties='all'): Number of pixels of region. * BoundingBox : tuple Bounding box `(min_row, min_col, max_row, max_col)` - * CentralMoments : 3x3 ndarray + * CentralMoments : 3 x 3 ndarray Central moments (translation invariant) up to 3rd order. - .. math:: - \texttt{mu} _{ji} = \sum _{x,y} \left (\texttt{array} (x,y) \\ - \cdot (x - \bar{x} )^j \cdot (y - \bar{y} )^i \right) + 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. * Centroid : array Centroid coordinate tuple `(row, col)`. * ConvexArea : int @@ -93,17 +93,15 @@ def regionprops(image, properties='all'): * MinorAxisLength : float The length of the minor axis of the ellipse that has the same normalized second central moments as the region. - * Moments 3x3 ndarray + * Moments 3 x 3 ndarray Spatial moments up to 3rd order. - .. math:: - \texttt{m} _{ji}= \sum _{x,y} \left (\texttt{array} (x,y) \\ - \cdot x^j \cdot y^i \right) - * NormalizedMoments : 3x3 ndarray + m_ji = sum{ array(x, y) * x^j * y^i } + where the sum is over the `x`, `y` coordinates of the region. + * NormalizedMoments : 3 x 3 ndarray Normalized moments (translation and scale invariant) up to 3rd order. - .. math:: - \texttt{nu} _{ji} = \\ - \frac{\texttt{mu}_{ji}}{\texttt{m}_{00}^{(i+j)/2+1}} + nu_ji = mu_ji / m_00^[(i+j)/2 + 1] + where `m_00` is the zeroth spatial moment. * Orientation : float Angle between the X-axis and the major axis of the ellipse that has the same second-moments as the region. Ranging from `-pi/2` to 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 23/25] 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 From 12a878dd590af2228039043d06b4c663e9e34418 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Wed, 16 May 2012 23:53:55 +0200 Subject: [PATCH 24/25] add example script for regionprops function to documentation --- doc/examples/plot_regionprops.py | 64 ++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) create mode 100644 doc/examples/plot_regionprops.py diff --git a/doc/examples/plot_regionprops.py b/doc/examples/plot_regionprops.py new file mode 100644 index 00000000..90b40a89 --- /dev/null +++ b/doc/examples/plot_regionprops.py @@ -0,0 +1,64 @@ +""" +========================= +Measure region properties +========================= + +This example shows how to measure properties of labelled image regions. +""" + +import math +import matplotlib.pyplot as plt +import numpy as np + +from skimage.draw import ellipse +from skimage.morphology import label +from skimage.measure import regionprops +from scipy.ndimage import geometric_transform + + +ANGLE = 0.2 + +def rotate(xy): + x, y = xy + out_x = math.cos(ANGLE) * x - math.sin(ANGLE) * y + out_y = math.sin(ANGLE) * x + math.cos(ANGLE) * y + return (out_x, out_y) + +image = np.zeros((600, 600), 'int') + +rr, cc = ellipse(300, 350, 100, 220) +image[rr,cc] = 1 + +image = geometric_transform(image, rotate) + +label_img = label(image) +props = regionprops(label_img, [ + 'BoundingBox', + 'Centroid', + 'Orientation', + 'MajorAxisLength', + 'MinorAxisLength' +]) + +plt.imshow(image) + +for prop in props: + x0 = prop['Centroid'][1] + y0 = prop['Centroid'][0] + x1 = x0 + math.cos(prop['Orientation']) * 0.5 * prop['MajorAxisLength'] + y1 = y0 - math.sin(prop['Orientation']) * 0.5 * prop['MajorAxisLength'] + x2 = x0 - math.sin(prop['Orientation']) * 0.5 * prop['MinorAxisLength'] + y2 = y0 - math.cos(prop['Orientation']) * 0.5 * prop['MinorAxisLength'] + + plt.plot((x0, x1), (y0, y1), '-r', linewidth=2.5) + plt.plot((x0, x2), (y0, y2), '-r', linewidth=2.5) + plt.plot(x0, y0, '.g', markersize=15) + + minr, minc, maxr, maxc = prop['BoundingBox'] + bx = (minc, maxc, maxc, minc, minc) + by = (minr, minr, maxr, maxr, minr) + plt.plot(bx, by, '-b', linewidth=2.5) + +plt.gray() +plt.axis((0, 600, 600, 0)) +plt.show() From 8e09617516028aa482cddd55ccc0277034e024e2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Scho=CC=88nberger?= Date: Tue, 22 May 2012 23:30:37 +0200 Subject: [PATCH 25/25] add moments Cython extension to bento.info --- bento.info | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/bento.info b/bento.info index 07809644..fe9cf164 100644 --- a/bento.info +++ b/bento.info @@ -4,10 +4,10 @@ Summary: Image processing routines for SciPy Url: http://scikits-image.org DownloadUrl: http://github.com/scikits-image/scikits-image Description: Image Processing SciKit - + Image processing algorithms for SciPy, including IO, morphology, filtering, warping, color manipulation, object detection, etc. - + Please refer to the online documentation at http://scikits-image.org/ Maintainer: Stefan van der Walt @@ -52,6 +52,9 @@ Library: Extension: skimage.measure._find_contours Sources: skimage/measure/_find_contours.pyx + Extension: skimage.measure._moments + Sources: + skimage/measure/_moments.pyx Extension: skimage.graph._mcp Sources: skimage/graph/_mcp.pyx