Merge pull request scikit-image#3303 from jni/fix-inertia-tensor

Fix incorrect inertia tensor

skimage/measure/_moments.py

@@ -410,7 +410,7 @@ def inertia_tensor(image, mu=None):
            Scientific Applications. (Chapter 8: Tensor Methods) Springer, 1993.
     """
     if mu is None:
-        mu = moments_central(image)
+        mu = moments_central(image, order=2)  # don't need higher-order moments
     mu0 = mu[(0,) * image.ndim]
     result = np.zeros((image.ndim, image.ndim))
 
@@ -419,7 +419,12 @@ def inertia_tensor(image, mu=None):
     corners2 = tuple(2 * np.eye(image.ndim, dtype=int))
     d = np.diag(result)
     d.flags.writeable = True
-    d[:] = mu[corners2] / mu0
+    # See https://ocw.mit.edu/courses/aeronautics-and-astronautics/
+    #             16-07-dynamics-fall-2009/lecture-notes/MIT16_07F09_Lec26.pdf
+    # Iii is the sum of second-order moments of every axis *except* i, not the
+    # second order moment of axis i.
+    # See also https://github.com/scikit-image/scikit-image/issues/3229
+    d[:] = (np.sum(mu[corners2]) - mu[corners2]) / mu0
 
     for dims in itertools.combinations(range(image.ndim), 2):
         mu_index = np.zeros(image.ndim, dtype=int)

skimage/measure/_regionprops.py

@@ -273,13 +273,14 @@ def moments_normalized(self):
     @only2d
     def orientation(self):
         a, b, b, c = self.inertia_tensor.flat
+        sign = -1 if self._transpose_moments else 1
         if a - c == 0:
             if b < 0:
                 return -PI / 4.
             else:
                 return PI / 4.
         else:
-            return -0.5 * atan2(-2 * b, (a - c))
+            return sign * 0.5 * atan2(-2 * b, c - a)
 
     @only2d
     def perimeter(self):
@@ -475,9 +476,12 @@ def regionprops(label_image, intensity_image=None, cache=True,
 
         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.
+        In 'rc' coordinates, angle between the 0th axis (rows) and the major
+        axis of the ellipse that has the same second moments as the region,
+        ranging from `-pi/2` to `pi/2` counter-clockwise.
+
+        In `xy` coordinates, as above but the angle is now measured from the
+        "x" or horizontal axis.
     **perimeter** : float
         Perimeter of object which approximates the contour as a line
         through the centers of border pixels using a 4-connectivity.

skimage/measure/tests/test_moments.py

@@ -1,8 +1,10 @@
 import numpy as np
+from scipy import ndimage as ndi
 from skimage import draw
 from skimage.measure import (moments, moments_central, moments_coords,
                              moments_coords_central, moments_normalized,
-                             moments_hu, centroid)
+                             moments_hu, centroid, inertia_tensor,
+                             inertia_tensor_eigvals)
 
 from skimage._shared import testing
 from skimage._shared.testing import (assert_equal, assert_almost_equal,
@@ -150,3 +152,37 @@ def test_centroid():
     image[15, 14:16] = 1/3
     image_centroid = centroid(image)
     assert_allclose(image_centroid, (14.25, 14.5))
+
+
+def test_inertia_tensor_2d():
+    image = np.zeros((40, 40))
+    image[15:25, 5:35] = 1  # big horizontal rectangle (aligned with axis 1)
+    T = inertia_tensor(image)
+    assert T[0, 0] > T[1, 1]
+    np.testing.assert_allclose(T[0, 1], 0)
+    v0, v1 = inertia_tensor_eigvals(image, T=T)
+    np.testing.assert_allclose(np.sqrt(v0/v1), 3, rtol=0.01, atol=0.05)
+
+
+def test_inertia_tensor_3d():
+    image = draw.ellipsoid(10, 5, 3)
+    T0 = inertia_tensor(image)
+    eig0, V0 = np.linalg.eig(T0)
+    # principal axis of ellipse = eigenvector of smallest eigenvalue
+    v0 = V0[:, np.argmin(eig0)]
+
+    assert np.allclose(v0, [1, 0, 0]) or np.allclose(-v0, [1, 0, 0])
+
+    imrot = ndi.rotate(image.astype(float), 30, axes=(0, 1), order=1)
+    Tr = inertia_tensor(imrot)
+    eigr, Vr = np.linalg.eig(Tr)
+    vr = Vr[:, np.argmin(eigr)]
+
+    # Check that axis has rotated by expected amount
+    pi, cos, sin = np.pi, np.cos, np.sin
+    R = np.array([[ cos(pi/6), -sin(pi/6), 0],
+                  [ sin(pi/6),  cos(pi/6), 0],
+                  [         0,          0, 1]])
+    expected_vr = R @ v0
+    assert (np.allclose(vr, expected_vr, atol=1e-3, rtol=0.01) or
+            np.allclose(-vr, expected_vr, atol=1e-3, rtol=0.01))

skimage/measure/tests/test_regionprops.py

@@ -305,15 +305,15 @@ def test_moments_normalized():
 
 
 def test_orientation():
-    orientation = regionprops(SAMPLE.T)[0].orientation
+    orientation = regionprops(SAMPLE, coordinates='xy')[0].orientation
     # determined with MATLAB
     assert_almost_equal(orientation, 0.10446844651921)
     # test correct quadrant determination
-    orientation2 = regionprops(SAMPLE)[0].orientation
-    assert_almost_equal(orientation2, math.pi / 2 - orientation)
+    orientation2 = regionprops(SAMPLE, coordinates='rc')[0].orientation
+    assert_almost_equal(orientation2, -math.pi / 2 + orientation)
     # test diagonal regions
     diag = np.eye(10, dtype=int)
-    orientation_diag = regionprops(diag)[0].orientation
+    orientation_diag = regionprops(diag, coordinates='xy')[0].orientation
     assert_almost_equal(orientation_diag, -math.pi / 4)
     orientation_diag = regionprops(np.flipud(diag))[0].orientation
     assert_almost_equal(orientation_diag, math.pi / 4)