From 32a8442d613665f3c0cffdac0edc00fe2ee3c5db Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Mon, 10 Sep 2012 21:01:03 +0200
Subject: [PATCH] Refactor harris corner detection, so it returns harris
 response image

---
 skimage/feature/interest.py | 103 +++++++++++-------------------------
 1 file changed, 32 insertions(+), 71 deletions(-)

diff --git a/skimage/feature/interest.py b/skimage/feature/interest.py
index ae30a29e..87163247 100644
--- a/skimage/feature/interest.py
+++ b/skimage/feature/interest.py
@@ -1,85 +1,31 @@
-"""
-Harris corner detector
-
-Inspired from Solem's implementation
-http://www.janeriksolem.net/2009/01/harris-corner-detector-in-python.html
-"""
+import numpy as np
 from scipy import ndimage
-
+from skimage.color import rgb2grey
+from skimage.util import img_as_float
 from . import peak
 
 
-def _compute_harris_response(image, eps=1e-6, gaussian_deviation=1):
-    """Compute the Harris corner detector response function
-    for each pixel in the image
+def harris(image, eps=1e-6, gaussian_deviation=1):
+    """Compute Harris response image.
 
     Parameters
     ----------
     image : ndarray of floats
         Input image.
-
     eps : float, optional
         Normalisation factor.
-
-    gaussian_deviation : integer, optional
-        Standard deviation used for the Gaussian kernel.
-
-    Returns
-    --------
-    image : (M, N) ndarray
-        Harris image response
-    """
-    if len(image.shape) == 3:
-        image = image.mean(axis=2)
-
-    # derivatives
-    image = ndimage.gaussian_filter(image, gaussian_deviation)
-    imx = ndimage.sobel(image, axis=0, mode='constant')
-    imy = ndimage.sobel(image, axis=1, mode='constant')
-
-    Wxx = ndimage.gaussian_filter(imx * imx, 1.5, mode='constant')
-    Wxy = ndimage.gaussian_filter(imx * imy, 1.5, mode='constant')
-    Wyy = ndimage.gaussian_filter(imy * imy, 1.5, mode='constant')
-
-    # determinant and trace
-    Wdet = Wxx * Wyy - Wxy**2
-    Wtr = Wxx + Wyy
-    # Alternate formula for Harris response.
-    # Alison Noble, "Descriptions of Image Surfaces", PhD thesis (1989)
-    harris = Wdet / (Wtr + eps)
-
-    return harris
-
-
-def harris(image, min_distance=10, threshold=0.1, eps=1e-6,
-           gaussian_deviation=1):
-    """Return corners from a Harris response image
-
-    Parameters
-    ----------
-    image : ndarray of floats
-        Input image.
-
-    min_distance : int, optional
-        Minimum number of pixels separating interest points and image boundary.
-
-    threshold : float, optional
-        Relative threshold impacting the number of interest points.
-
-    eps : float, optional
-        Normalisation factor.
-
     gaussian_deviation : integer, optional
         Standard deviation used for the Gaussian kernel.
 
     Returns
     -------
     coordinates : (N, 2) array
-        (row, column) coordinates of interest points.
+        `(row, column)` coordinates of interest points.
 
     Examples
     -------
-    >>> square = np.zeros([10,10])
+    >>> from skimage.feature import harris, peak_local_max
+    >>> square = np.zeros([10, 10])
     >>> square[2:8,2:8] = 1
     >>> square
     array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
@@ -92,18 +38,33 @@ def harris(image, min_distance=10, threshold=0.1, eps=1e-6,
            [ 0.,  0.,  1.,  1.,  1.,  1.,  1.,  1.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])
-    >>> harris(square, min_distance=1)
-
-    Corners of the square
-
+    >>> peak_local_max(harris(square), min_distance=1)
     array([[3, 3],
            [3, 6],
            [6, 3],
            [6, 6]])
+
     """
 
-    harrisim = _compute_harris_response(image, eps=eps,
-                                        gaussian_deviation=gaussian_deviation)
-    coordinates = peak.peak_local_max(harrisim, min_distance=min_distance,
-                                      threshold_rel=threshold)
-    return coordinates
+    if image.ndim == 3:
+        image = rgb2grey(image)
+
+    # derivatives
+    image = ndimage.gaussian_filter(image, gaussian_deviation, mode='constant',
+                                    cval=0)
+    imx = ndimage.sobel(image, axis=0, mode='constant', cval=0)
+    imy = ndimage.sobel(image, axis=1, mode='constant', cval=0)
+
+    Wxx = ndimage.gaussian_filter(imx * imx, 1.5, mode='constant', cval=0)
+    Wxy = ndimage.gaussian_filter(imx * imy, 1.5, mode='constant', cval=0)
+    Wyy = ndimage.gaussian_filter(imy * imy, 1.5, mode='constant', cval=0)
+
+    # determinant and trace
+    Wdet = Wxx * Wyy - Wxy**2
+    Wtr = Wxx + Wyy
+
+    # Alternate formula for Harris response.
+    # Alison Noble, "Descriptions of Image Surfaces", PhD thesis (1989)
+    harris = Wdet / (Wtr + eps)
+
+    return harris