From 84317b6a2f95c16522d2f2c64185d9d22a5b9cda Mon Sep 17 00:00:00 2001
From: Ankit Agrawal <aaaagrawal@iitb.ac.in>
Date: Mon, 26 Aug 2013 22:29:36 +0530
Subject: [PATCH] Speeding up corner_fast_implementation : 2

---
 skimage/feature/corner_cy.pyx        | 49 +++++++++++++++++-----------
 skimage/feature/tests/test_corner.py | 12 ++++---
 2 files changed, 37 insertions(+), 24 deletions(-)

diff --git a/skimage/feature/corner_cy.pyx b/skimage/feature/corner_cy.pyx
index 3ef290a1..2f916d4d 100644
--- a/skimage/feature/corner_cy.pyx
+++ b/skimage/feature/corner_cy.pyx
@@ -167,7 +167,7 @@ def _corner_fast(double[:, ::1] image, char n, double threshold):
     return np.asarray(corner_response)
 
 
-def corner_fast_orientation(image, fast_corners):
+def corner_fast_orientation(image, Py_ssize_t[:, :] fast_corners):
     """Compute the orientation of FAST corners using the first order central
     moment i.e. the center of mass approach. The corner orientation is the
     angle of the vector from the keypoint to the intensity centroid calculated
@@ -199,35 +199,46 @@ def corner_fast_orientation(image, fast_corners):
     if image.ndim != 2:
         raise ValueError("Only 2-D gray-scale images supported.")
 
+    cdef double[:, :] cimage = np.ascontiguousarray(img_as_float(image))
     # Essentially skimage.morphology.octagon(3, 2)
-    cdef char[:, :] circular_mask = np.array([[0, 0, 1, 1, 1, 0, 0],
-                                              [0, 1, 1, 1, 1, 1, 0],
-                                              [1, 1, 1, 1, 1, 1, 1],
-                                              [1, 1, 1, 1, 1, 1, 1],
-                                              [1, 1, 1, 1, 1, 1, 1],
-                                              [0, 1, 1, 1, 1, 1, 0],
-                                              [0, 0, 1, 1, 1, 0, 0]], dtype=np.uint8)
-
-    cdef Py_ssize_t[:, :] cfast_corners = np.ascontiguousarray(fast_corners, dtype=np.intp)
+    cdef char[:, ::1] circular_mask = np.array([[0, 0, 1, 1, 1, 0, 0],
+                                                [0, 1, 1, 1, 1, 1, 0],
+                                                [1, 1, 1, 1, 1, 1, 1],
+                                                [1, 1, 1, 1, 1, 1, 1],
+                                                [1, 1, 1, 1, 1, 1, 1],
+                                                [0, 1, 1, 1, 1, 1, 0],
+                                                [0, 0, 1, 1, 1, 0, 0]], dtype=np.uint8)
 
     cdef Py_ssize_t n_fast_corners = fast_corners.shape[0]
-    cdef Py_ssize_t i, p, q, r, c, x, y
+    cdef Py_ssize_t i, r, c, x_top, y_left
     cdef double[:, ::1] kp_circular_patch, mu
     cdef double[:] kp_orientation = np.zeros(fast_corners.shape[0], dtype=np.double)
+    cdef double m00, m01, m10
 
     for i in range(n_fast_corners):
-        x = cfast_corners[i, 0]
-        y = cfast_corners[i, 1]
+        x_top = fast_corners[i, 0] - 3
+        y_left = fast_corners[i, 1] - 3
 
-        kp_circular_patch = np.ascontiguousarray(image[x - 3:x + 4, y - 3:y + 4])
-        mu = np.zeros((2, 2), dtype=np.double)
+        m00 = 0
+        m01 = 0
+        m10 = 0
+        #kp_circular_patch = np.ascontiguousarray(image[x - 3:x + 4, y - 3:y + 4])
+        #mu = np.zeros((2, 2), dtype=np.double)
         for r in range(7):
             for c in range(7):
                 if circular_mask[r, c]:
-                    for p in range(2):
-                        for q in range(2):
-                            mu[p, q] += kp_circular_patch[r, c] * (r - 3) ** q * (c - 3) ** p
+                    m00 += cimage[x_top + r, y_left + c]
 
-        kp_orientation[i] = atan2(mu[1, 0] / mu[0, 0], mu[0, 1] / mu[0, 0])
+        for r in range(7):
+            for c in range(7):
+                if circular_mask[r, c]:
+                    m01 += cimage[x_top + r, y_left + c] * (c - 3)
+
+        for r in range(7):
+            for c in range(7):
+                if circular_mask[r, c]:
+                    m10 += cimage[x_top + r, y_left + c] * (r - 3)
+
+        kp_orientation[i] = atan2(m10 / m00, m01 / m00)
 
     return np.asarray(kp_orientation)
diff --git a/skimage/feature/tests/test_corner.py b/skimage/feature/tests/test_corner.py
index 56985613..1a520ceb 100644
--- a/skimage/feature/tests/test_corner.py
+++ b/skimage/feature/tests/test_corner.py
@@ -1,5 +1,6 @@
 import numpy as np
-from numpy.testing import assert_array_equal, assert_almost_equal
+from numpy.testing import (assert_array_equal, assert_raises,
+                           assert_almost_equal)
 
 from skimage import data
 from skimage import img_as_float
@@ -181,16 +182,17 @@ def test_corner_fast_lena():
 
 def test_corner_fast_orientation_image_unsupported_error():
     img = np.zeros((20, 20, 3))
-    assert_raises(ValueError, corner_fast_orientation, img, [[7, 7]])
+    assert_raises(ValueError, corner_fast_orientation, img,
+                  np.asarray([[7, 7]]))
 
 
 def test_corner_fast_orientation_lena():
     img = rgb2gray(data.lena())
     corners = corner_peaks(corner_fast(img, 11, 0.35))
-    expected = np.array([-2.79279928, -1.68079274,  2.63070795, -1.81665159,
-                         -2.09631254, -1.41580527])
+    expected = np.array([-1.9195897 , -3.03159624, -1.05991162, -2.89573739,
+                         -2.61607644, 2.98660159])
     actual = corner_fast_orientation(img, corners)
-    assert_array_equal(actual, expected)
+    assert_almost_equal(actual, expected)
 
 
 if __name__ == '__main__':