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` 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 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() diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py index beedd379..422db569 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/_moments.pyx b/skimage/measure/_moments.pyx new file mode 100644 index 00000000..f84e14dd --- /dev/null +++ b/skimage/measure/_moments.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.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 + 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.py b/skimage/measure/_regionprops.py new file mode 100644 index 00000000..937b80f5 --- /dev/null +++ b/skimage/measure/_regionprops.py @@ -0,0 +1,353 @@ +# coding: utf-8 +from math import sqrt, atan, pi as PI +import numpy as np +from scipy import ndimage + +from skimage.morphology import convex_hull_image +from . import _moments + + +__all__ = ['regionprops'] + + +STREL_8 = np.ones((3, 3), 'int8') +PROPS = ( + 'Area', + 'BoundingBox', + 'CentralMoments', + 'Centroid', + 'ConvexArea', +# 'ConvexHull', + 'ConvexImage', + 'Eccentricity', + 'EquivDiameter', + 'EulerNumber', + 'Extent', +# 'Extrema', + 'FilledArea', + 'FilledImage', + 'HuMoments', + 'Image', + 'MajorAxisLength', + 'MaxIntensity', + 'MeanIntensity', + 'MinIntensity', + 'MinorAxisLength', + 'Moments', + 'NormalizedMoments', + 'Orientation', +# 'Perimeter', +# 'PixelIdxList', +# 'PixelList', + 'Solidity', +# 'SubarrayIdx' + 'WeightedCentralMoments', + 'WeightedCentroid', + 'WeightedHuMoments', + 'WeightedMoments', + 'WeightedNormalizedMoments' +) + + +def regionprops(label_image, properties=['Area', 'Centroid'], + intensity_image=None): + """Measure properties of labelled image regions. + + Parameters + ---------- + label_image : N x M ndarray + Labelled input image. + properties : {'all', list} + Shape measurements to be determined for each labelled image region. + Default is `['Area', 'Centroid']`. The following properties can be + determined: + * Area : int + Number of pixels of region. + * BoundingBox : tuple + Bounding box `(min_row, min_col, max_row, max_col)` + * CentralMoments : 3 x 3 ndarray + Central moments (translation invariant) up to 3rd order. + 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 + Number of pixels of convex hull image. + * ConvexImage : H x J ndarray + 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 + 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 + 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 : H x J ndarray + 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 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 + 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. + * NormalizedMoments : 3 x 3 ndarray + Normalized moments (translation and scale invariant) up to 3rd + order. + 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 + `-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 + ------- + 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 + ---------- + 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, + from Lecture notes in computer science, p. 676. 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) + >>> props[0]['Centroid'] # centroid of first labelled object + """ + if not np.issubdtype(label_image.dtype, 'int'): + raise TypeError('labelled image must be of integer dtype') + + # determine all properties if nothing specified + if properties == 'all': + properties = PROPS + + props = [] + + objects = ndimage.find_objects(label_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 + + array = (label_image[sl] == label).astype('double') + + # upper left corner of object bbox + r0 = sl[0].start + c0 = sl[1].start + + m = _moments.central_moments(array, 0, 0, 3) + # centroid + cr = m[0,1] / m[0,0] + cc = m[1,0] / m[0,0] + mu = _moments.central_moments(array, cr, cc, 3) + + #: 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 + _convex_image = None + _nu = 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: + 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) + + 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'] = - num + + if 'Extent' in properties: + obj_props['Extent'] = m[0,0] / (array.shape[0] * array.shape[1]) + + if 'HuMoments' in properties: + if _nu is None: + _nu = _moments.normalized_moments(mu, 3) + obj_props['HuMoments'] = _moments.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 'MajorAxisLength' in properties: + obj_props['MajorAxisLength'] = 4 * sqrt(l1) + + if 'MinorAxisLength' in properties: + obj_props['MinorAxisLength'] = 4 * sqrt(l2) + + if 'Moments' in properties: + obj_props['Moments'] = m + + if 'NormalizedMoments' in properties: + if _nu is None: + _nu = _moments.normalized_moments(mu, 3) + obj_props['NormalizedMoments'] = _nu + + if 'Orientation' in properties: + 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/setup.py b/skimage/measure/setup.py index f03b9f3b..4b02a1d1 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(['_moments.pyx'], working_path=base_path) config.add_extension('_find_contours', sources=['_find_contours.c'], include_dirs=[get_numpy_include_dirs()]) + config.add_extension('_moments', sources=['_moments.c'], + include_dirs=[get_numpy_include_dirs()]) return config diff --git a/skimage/measure/tests/test_regionprops.py b/skimage/measure/tests/test_regionprops.py new file mode 100644 index 00000000..417311de --- /dev/null +++ b/skimage/measure/tests/test_regionprops.py @@ -0,0 +1,250 @@ +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]] +) +INTENSITY_SAMPLE = SAMPLE.copy() +INTENSITY_SAMPLE[1,9:11] = 2 + + +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) + # different from OpenCV results, bug in OpenCV + assert_almost_equal(mu[0,3], -737.333333333333) + assert_almost_equal(mu[1,1], -87.33333333333303) + assert_almost_equal(mu[1,2], -127.5555555555593) + assert_almost_equal(mu[2,0], 1259.7777777777774) + assert_almost_equal(mu[2,1], 2000.296296296291) + assert_almost_equal(mu[3,0], -760.0246913580195) + +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'] + assert_almost_equal(eps, 0.814629313427) + +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'] + 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) + +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_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 + 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) + +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 + run_module_suite()