wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-06-30 16:13:58 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #694 from ahojnnes/moments

Browse Source 
Refactor moments functions.
This commit is contained in:
Stefan van der Walt
2013-08-17 12:53:06 -07:00
parent 969927dc2a cc91324625
commit b7274a7f17
5 changed files with 302 additions and 101 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/measure/__init__.py
+6 -1
View File
@@ -2,6 +2,7 @@ from .find_contours import find_contours
from ._regionprops import regionprops, perimeter
from ._structural_similarity import structural_similarity
from ._polygon import approximate_polygon, subdivide_polygon
from ._moments import moments, moments_central, moments_normalized, moments_hu
from .fit import LineModel, CircleModel, EllipseModel, ransac
from .block import block_reduce
@@ -16,4 +17,8 @@ __all__ = ['find_contours',
'CircleModel',
'EllipseModel',
'ransac',
'block_reduce']
'block_reduce',
'moments',
'moments_central',
'moments_normalized',
'moments_hu']
skimage/measure/_moments.pyx
+151 -26
View File
@@ -4,53 +4,178 @@
#cython: wraparound=False
import numpy as np
cimport numpy as cnp
def moments(double[:, :] image, Py_ssize_t order=3):
"""Calculate all raw image moments up to a certain order.
The following properties can be calculated from raw image moments:
* Area as ``m[0, 0]``.
* Centroid as {``m[0, 1] / m[0, 0]``, ``m[1, 0] / m[0, 0]``}.
Note that raw moments are whether translation, scale nor rotation
invariant.
Parameters
----------
image : 2D double array
Rasterized shape as image.
order : int, optional
Maximum order of moments. Default is 3.
Returns
-------
m : (``order + 1``, ``order + 1``) array
Raw image moments.
References
----------
.. [1] Wilhelm Burger, Mark Burge. Principles of Digital Image Processing:
Core Algorithms. Springer-Verlag, London, 2009.
.. [2] B. Jähne. Digital Image Processing. Springer-Verlag,
Berlin-Heidelberg, 6. edition, 2005.
.. [3] T. H. Reiss. Recognizing Planar Objects Using Invariant Image
Features, from Lecture notes in computer science, p. 676. Springer,
Berlin, 1993.
.. [4] http://en.wikipedia.org/wiki/Image_moment
"""
return moments_central(image, 0, 0, order)
def central_moments(cnp.ndarray[cnp.double_t, ndim=2] array, double cr,
double cc, int order):
def moments_central(double[:, :] image, double cr, double cc,
Py_ssize_t order=3):
"""Calculate all central image moments up to a certain order.
Note that central moments are translation invariant but not scale and
rotation invariant.
Parameters
----------
image : 2D double array
Rasterized shape as image.
cr : double
Center row coordinate.
cc : double
Center column coordinate.
order : int, optional
Maximum order of moments. Default is 3.
Returns
-------
mu : (``order + 1``, ``order + 1``) array
Central image moments.
References
----------
.. [1] Wilhelm Burger, Mark Burge. Principles of Digital Image Processing:
Core Algorithms. Springer-Verlag, London, 2009.
.. [2] B. Jähne. Digital Image Processing. Springer-Verlag,
Berlin-Heidelberg, 6. edition, 2005.
.. [3] T. H. Reiss. Recognizing Planar Objects Using Invariant Image
Features, from Lecture notes in computer science, p. 676. Springer,
Berlin, 1993.
.. [4] http://en.wikipedia.org/wiki/Image_moment
"""
cdef Py_ssize_t p, q, r, c
cdef cnp.ndarray[cnp.double_t, ndim=2] mu
mu = np.zeros((order + 1, order + 1), 'double')
cdef double[:, ::1] mu = np.zeros((order + 1, order + 1), dtype=np.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
for r in range(image.shape[0]):
for c in range(image.shape[1]):
mu[p, q] += image[r, c] * (r - cr) ** q * (c - cc) ** p
return np.asarray(mu)
def normalized_moments(cnp.ndarray[cnp.double_t, ndim=2] mu, int order):
def moments_normalized(double[:, :] mu, Py_ssize_t order=3):
"""Calculate all normalized central image moments up to a certain order.
Note that normalized central moments are translation and scale invariant
but not rotation invariant.
Parameters
----------
mu : (M, M) array
Central image moments, where M must be > ``order``.
order : int, optional
Maximum order of moments. Default is 3.
Returns
-------
nu : (``order + 1``, ``order + 1``) array
Normalized central image moments.
References
----------
.. [1] Wilhelm Burger, Mark Burge. Principles of Digital Image Processing:
Core Algorithms. Springer-Verlag, London, 2009.
.. [2] B. Jähne. Digital Image Processing. Springer-Verlag,
Berlin-Heidelberg, 6. edition, 2005.
.. [3] T. H. Reiss. Recognizing Planar Objects Using Invariant Image
Features, from Lecture notes in computer science, p. 676. Springer,
Berlin, 1993.
.. [4] http://en.wikipedia.org/wiki/Image_moment
"""
cdef Py_ssize_t p, q
cdef cnp.ndarray[cnp.double_t, ndim=2] nu
nu = np.zeros((order + 1, order + 1), 'double')
cdef double[:, ::1] nu = np.zeros((order + 1, order + 1), dtype=np.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]**(<double>(p + q) / 2 + 1)
nu[p, q] = mu[p, q] / mu[0, 0] ** (<double>(p + q) / 2 + 1)
else:
nu[p,q] = np.nan
return nu
nu[p, q] = np.nan
return np.asarray(nu)
def hu_moments(cnp.ndarray[cnp.double_t, ndim=2] nu):
cdef cnp.ndarray[cnp.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]
def moments_hu(double[:, :] nu):
"""Calculate Hu's set of image moments.
Note that this set of moments is proofed to be translation, scale and
rotation invariant.
Parameters
----------
nu : (M, M) array
Normalized central image moments, where M must be > 4.
Returns
-------
nu : (7, 1) array
Hu's set of image moments.
References
----------
.. [1] M. K. Hu, "Visual Pattern Recognition by Moment Invariants",
IRE Trans. Info. Theory, vol. IT-8, pp. 179-187, 1962
.. [2] Wilhelm Burger, Mark Burge. Principles of Digital Image Processing:
Core Algorithms. Springer-Verlag, London, 2009.
.. [3] B. Jähne. Digital Image Processing. Springer-Verlag,
Berlin-Heidelberg, 6. edition, 2005.
.. [4] T. H. Reiss. Recognizing Planar Objects Using Invariant Image
Features, from Lecture notes in computer science, p. 676. Springer,
Berlin, 1993.
.. [5] http://en.wikipedia.org/wiki/Image_moment
"""
cdef double[::1] hu = np.zeros((7, ), dtype=np.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]
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[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]
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
return np.asarray(hu)
skimage/measure/_regionprops.py
+61 -62
View File
@@ -18,7 +18,7 @@ STREL_8 = np.ones((3, 3), 'int8')
PROPS = {
'Area': 'area',
'BoundingBox': 'bbox',
'CentralMoments': 'central_moments',
'CentralMoments': 'moments_central',
'Centroid': 'centroid',
'ConvexArea': 'convex_area',
# 'ConvexHull',
@@ -31,7 +31,7 @@ PROPS = {
# 'Extrema',
'FilledArea': 'filled_area',
'FilledImage': 'filled_image',
'HuMoments': 'hu_moments',
'HuMoments': 'moments_hu',
'Image': 'image',
'Label': 'label',
'MajorAxisLength': 'major_axis_length',
@@ -40,7 +40,7 @@ PROPS = {
'MinIntensity': 'min_intensity',
'MinorAxisLength': 'minor_axis_length',
'Moments': 'moments',
'NormalizedMoments': 'normalized_moments',
'NormalizedMoments': 'moments_normalized',
'Orientation': 'orientation',
'Perimeter': 'perimeter',
# 'PixelIdxList',
@@ -49,9 +49,9 @@ PROPS = {
# 'SubarrayIdx'
'WeightedCentralMoments': 'weighted_central_moments',
'WeightedCentroid': 'weighted_centroid',
'WeightedHuMoments': 'weighted_hu_moments',
'WeightedHuMoments': 'weighted_moments_hu',
'WeightedMoments': 'weighted_moments',
'WeightedNormalizedMoments': 'weighted_normalized_moments'
'WeightedNormalizedMoments': 'weighted_moments_normalized'
}
@@ -127,11 +127,6 @@ class _RegionProperties(object):
row, col = self.local_centroid
return row + self._slice[0].start, col + self._slice[1].start
@_cached_property
def central_moments(self):
row, col = self.local_centroid
return _moments.central_moments(self._image_double, row, col, 3)
@_cached_property
def convex_area(self):
return np.sum(self.convex_image)
@@ -176,10 +171,6 @@ class _RegionProperties(object):
def filled_image(self):
return ndimage.binary_fill_holes(self.image, STREL_8)
@_cached_property
def hu_moments(self):
return _moments.hu_moments(self.normalized_moments)
@_cached_property
def image(self):
return self._label_image[self._slice] == self.label
@@ -190,7 +181,7 @@ class _RegionProperties(object):
@_cached_property
def inertia_tensor(self):
mu = self.central_moments
mu = self.moments_central
a = mu[2, 0] / mu[0, 0]
b = -mu[1, 1] / mu[0, 0]
c = mu[0, 2] / mu[0, 0]
@@ -214,10 +205,6 @@ class _RegionProperties(object):
def _intensity_image_double(self):
return self.intensity_image.astype(np.double)
@_cached_property
def moments(self):
return _moments.central_moments(self._image_double, 0, 0, 3)
@_cached_property
def local_centroid(self):
m = self.moments
@@ -248,8 +235,21 @@ class _RegionProperties(object):
return 4 * sqrt(l2)
@_cached_property
def normalized_moments(self):
return _moments.normalized_moments(self.central_moments, 3)
def moments(self):
return _moments.moments(self._image_double, 3)
@_cached_property
def moments_central(self):
row, col = self.local_centroid
return _moments.moments_central(self._image_double, row, col, 3)
@_cached_property
def moments_hu(self):
return _moments.moments_hu(self.moments_normalized)
@_cached_property
def moments_normalized(self):
return _moments.moments_normalized(self.moments_central, 3)
@_cached_property
def orientation(self):
@@ -271,12 +271,6 @@ class _RegionProperties(object):
def solidity(self):
return self.moments[0, 0] / np.sum(self.convex_image)
@_cached_property
def weighted_central_moments(self):
row, col = self.weighted_local_centroid
return _moments.central_moments(self._intensity_image_double,
row, col, 3)
@_cached_property
def weighted_centroid(self):
row, col = self.weighted_local_centroid
@@ -289,17 +283,23 @@ class _RegionProperties(object):
col = m[1, 0] / m[0, 0]
return row, col
@_cached_property
def weighted_hu_moments(self):
return _moments.hu_moments(self.weighted_normalized_moments)
@_cached_property
def weighted_moments(self):
return _moments.central_moments(self._intensity_image_double, 0, 0, 3)
return _moments.moments_central(self._intensity_image_double, 0, 0, 3)
@_cached_property
def weighted_normalized_moments(self):
return _moments.normalized_moments(self.weighted_central_moments, 3)
def weighted_central_moments(self):
row, col = self.weighted_local_centroid
return _moments.moments_central(self._intensity_image_double,
row, col, 3)
@_cached_property
def weighted_moments_hu(self):
return _moments.moments_hu(self.weighted_moments_normalized)
@_cached_property
def weighted_moments_normalized(self):
return _moments.moments_normalized(self.weighted_central_moments, 3)
def __getitem__(self, key):
value = getattr(self, key, None)
@@ -311,7 +311,6 @@ class _RegionProperties(object):
return getattr(self, PROPS[key])
def regionprops(label_image, properties=None,
intensity_image=None, cache=True):
"""Measure properties of labelled image regions.
@@ -346,22 +345,15 @@ def regionprops(label_image, properties=None,
**area** : int
Number of pixels of region.
**bbox** : tuple
Bounding box `(min_row, min_col, max_row, max_col)`
**central_moments** : (3, 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.
Bounding box ``(min_row, min_col, max_row, max_col)``
**centroid** : array
Centroid coordinate tuple `(row, col)`.
Centroid coordinate tuple ``(row, col)``.
**convex_area** : int
Number of pixels of convex hull image.
**convex_image** : (H, J) ndarray
Binary convex hull image which has the same size as bounding box.
**coords** : (N, 2) ndarray
Coordinate list `(row, col)` of the region.
Coordinate list ``(row, col)`` of the region.
**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
@@ -373,22 +365,20 @@ def regionprops(label_image, properties=None,
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)`
Computed as ``area / (rows * cols)``
**filled_area** : int
Number of pixels of filled region.
**filled_image** : (H, J) ndarray
Binary region image with filled holes which has the same size as
bounding box.
**hu_moments** : tuple
Hu moments (translation, scale and rotation invariant).
**image** : (H, J) ndarray
Sliced binary region image which has the same size as bounding box.
**inertia_tensor** : (2, 2) ndarray
Inertia tensor of the region for the rotation around its masss.
Inertia tensor of the region for the rotation around its mass.
**inertia_tensor_eigvals** : tuple
The two eigen values of the inertia tensor in decreasing order.
**label** : int
The label in the labelled input image.
The label in the labeled input image.
**major_axis_length** : float
The length of the major axis of the ellipse that has the same
normalized second central moments as the region.
@@ -407,7 +397,16 @@ def regionprops(label_image, properties=None,
m_ji = sum{ array(x, y) * x^j * y^i }
where the sum is over the `x`, `y` coordinates of the region.
**normalized_moments** : (3, 3) ndarray
**moments_central** : (3, 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.
**moments_hu** : tuple
Hu moments (translation, scale and rotation invariant).
**moments_normalized** : (3, 3) ndarray
Normalized moments (translation and scale invariant) up to 3rd order::
nu_ji = mu_ji / m_00^[(i+j)/2 + 1]
@@ -422,6 +421,15 @@ def regionprops(label_image, properties=None,
through the centers of border pixels using a 4-connectivity.
**solidity** : float
Ratio of pixels in the region to pixels of the convex hull image.
**weighted_centroid** : array
Centroid coordinate tuple ``(row, col)`` weighted with intensity
image.
**weighted_moments** : (3, 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.
**weighted_central_moments** : (3, 3) ndarray
Central moments (translation invariant) of intensity image up to
3rd order::
@@ -430,19 +438,10 @@ def regionprops(label_image, properties=None,
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.
**weighted_centroid** : array
Centroid coordinate tuple `(row, col)` weighted with intensity
image.
**weighted_hu_moments** : tuple
**weighted_moments_hu** : tuple
Hu moments (translation, scale and rotation invariant) of intensity
image.
**weighted_moments** : (3, 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.
**weighted_normalized_moments** : (3, 3) ndarray
**weighted_moments_normalized** : (3, 3) ndarray
Normalized moments (translation and scale invariant) of intensity
image up to 3rd order::
skimage/measure/tests/test_moments.py
+72
View File
@@ -0,0 +1,72 @@
from numpy.testing import assert_equal, assert_almost_equal
import numpy as np
from skimage.measure import (moments, moments_central, moments_normalized,
moments_hu)
def test_moments():
image = np.zeros((20, 20), dtype=np.double)
image[14, 14] = 1
image[15, 15] = 1
image[14, 15] = 0.5
image[15, 14] = 0.5
m = moments(image)
assert_equal(m[0, 0], 3)
assert_almost_equal(m[0, 1] / m[0, 0], 14.5)
assert_almost_equal(m[1, 0] / m[0, 0], 14.5)
def test_moments_central():
image = np.zeros((20, 20), dtype=np.double)
image[14, 14] = 1
image[15, 15] = 1
image[14, 15] = 0.5
image[15, 14] = 0.5
mu = moments_central(image, 14.5, 14.5)
# shift image by dx=2, dy=2
image2 = np.zeros((20, 20), dtype=np.double)
image2[16, 16] = 1
image2[17, 17] = 1
image2[16, 17] = 0.5
image2[17, 16] = 0.5
mu2 = moments_central(image2, 14.5 + 2, 14.5 + 2)
# central moments must be translation invariant
assert_equal(mu, mu2)
def test_moments_normalized():
image = np.zeros((20, 20), dtype=np.double)
image[13:17, 13:17] = 1
mu = moments_central(image, 14.5, 14.5)
nu = moments_normalized(mu)
# shift image by dx=-3, dy=-3 and scale by 0.5
image2 = np.zeros((20, 20), dtype=np.double)
image2[11:13, 11:13] = 1
mu2 = moments_central(image2, 11.5, 11.5)
nu2 = moments_normalized(mu2)
# central moments must be translation and scale invariant
assert_almost_equal(nu, nu2, decimal=1)
def test_moments_hu():
image = np.zeros((20, 20), dtype=np.double)
image[13:15, 13:17] = 1
mu = moments_central(image, 13.5, 14.5)
nu = moments_normalized(mu)
hu = moments_hu(nu)
# shift image by dx=2, dy=3, scale by 0.5 and rotate by 90deg
image2 = np.zeros((20, 20), dtype=np.double)
image2[11, 11:13] = 1
image2 = image2.T
mu2 = moments_central(image2, 11.5, 11)
nu2 = moments_normalized(mu2)
hu2 = moments_hu(nu2)
# central moments must be translation and scale invariant
assert_almost_equal(hu, hu2, decimal=1)
if __name__ == "__main__":
from numpy.testing import run_module_suite
run_module_suite()
skimage/measure/tests/test_regionprops.py
+12 -12
View File
@@ -47,8 +47,8 @@ def test_bbox():
assert_array_almost_equal(bbox, (0, 0, SAMPLE.shape[0], SAMPLE.shape[1]-1))
def test_central_moments():
mu = regionprops(SAMPLE)[0].central_moments
def test_moments_central():
mu = regionprops(SAMPLE)[0].moments_central
# determined with OpenCV
assert_almost_equal(mu[0,2], 436.00000000000045)
# different from OpenCV results, bug in OpenCV
@@ -129,8 +129,8 @@ def test_extent():
assert_almost_equal(extent, 0.4)
def test_hu_moments():
hu = regionprops(SAMPLE)[0].hu_moments
def test_moments_hu():
hu = regionprops(SAMPLE)[0].moments_hu
ref = np.array([
3.27117627e-01,
2.63869194e-02,
@@ -216,8 +216,8 @@ def test_moments():
assert_almost_equal(m[3,0], 95588.0)
def test_normalized_moments():
nu = regionprops(SAMPLE)[0].normalized_moments
def test_moments_normalized():
nu = regionprops(SAMPLE)[0].moments_normalized
# determined with OpenCV
assert_almost_equal(nu[0,2], 0.08410493827160502)
assert_almost_equal(nu[1,1], -0.016846707818929982)
@@ -260,9 +260,9 @@ def test_solidity():
assert_almost_equal(solidity, 0.580645161290323)
def test_weighted_central_moments():
def test_weighted_moments():
wmu = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE
)[0].weighted_central_moments
)[0].weighted_moments_central
ref = np.array(
[[ 7.4000000000e+01, -2.1316282073e-13, 4.7837837838e+02,
-7.5943608473e+02],
@@ -283,9 +283,9 @@ def test_weighted_centroid():
assert_array_almost_equal(centroid, (5.540540540540, 9.445945945945))
def test_weighted_hu_moments():
def test_weighted_moments_hu():
whu = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE
)[0].weighted_hu_moments
)[0].weighted_moments_hu
ref = np.array([
3.1750587329e-01,
2.1417517159e-02,
@@ -314,9 +314,9 @@ def test_weighted_moments():
assert_array_almost_equal(wm, ref)
def test_weighted_normalized_moments():
def test_weighted_moments_normalized():
wnu = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE
)[0].weighted_normalized_moments
)[0].weighted_moments_normalized
ref = np.array(
[[ np.nan, np.nan, 0.0873590903, -0.0161217406],
[ np.nan, -0.0160405109, -0.0031421072, -0.0031376984],