From 7b86728498cbda478e393f1fb3b63d87ec4c10dc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Sat, 10 May 2014 10:12:41 -0400
Subject: [PATCH] Add full test coverage for corner_subpix and fix bug

---
 skimage/feature/corner.py                     | 25 +++++++++++++------
 skimage/feature/tests/test_corner.py          | 21 +++++++++++++++-
 skimage/morphology/greyreconstruct.py         |  4 +--
 .../morphology/tests/test_reconstruction.py   |  8 +++---
 4 files changed, 45 insertions(+), 13 deletions(-)

diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py
index e1dd519b..3b22bbb6 100644
--- a/skimage/feature/corner.py
+++ b/skimage/feature/corner.py
@@ -737,8 +737,12 @@ def corner_subpix(image, corners, window_size=11, alpha=0.99):
         b_edge[:] = byy_y + bxy_x, bxx_x + bxy_y
 
         # estimated positions
-        est_dot = np.linalg.solve(N_dot, b_dot)
-        est_edge = np.linalg.solve(N_edge, b_edge)
+        try:
+            est_dot = np.linalg.solve(N_dot, b_dot)
+            est_edge = np.linalg.solve(N_edge, b_edge)
+        except np.linalg.LinAlgError:
+            corners_subpix[i, :] = np.nan, np.nan
+            continue
 
         # residuals
         ry_dot = y - est_dot[0]
@@ -759,12 +763,19 @@ def corner_subpix(image, corners, window_size=11, alpha=0.99):
                          + winy_winy * rxx_dot)
         var_edge = np.sum(winy_winy * ryy_edge + 2 * winx_winy * rxy_edge \
                           + winx_winx * rxx_edge)
-        # test value (F-distributed)
-        t = var_edge / var_dot
-        # 1 for edge, -1 for dot, 0 for "not classified"
-        corner_class = (t < t_crit_edge) - (t > t_crit_dot)
 
-        if corner_class == - 1:
+        # test value (F-distributed)
+        if var_dot < np.spacing(1) and var_edge < np.spacing(1):
+            t = np.nan
+        elif var_dot == 0:
+            t = np.inf
+        else:
+            t = var_edge / var_dot
+
+        # 1 for edge, -1 for dot, 0 for "not classified"
+        corner_class = int(t < t_crit_edge) - int(t > t_crit_dot)
+
+        if corner_class == -1:
             corners_subpix[i, :] = y0 + est_dot[0], x0 + est_dot[1]
         elif corner_class == 0:
             corners_subpix[i, :] = np.nan, np.nan
diff --git a/skimage/feature/tests/test_corner.py b/skimage/feature/tests/test_corner.py
index 7d2af0df..1c39dd2d 100644
--- a/skimage/feature/tests/test_corner.py
+++ b/skimage/feature/tests/test_corner.py
@@ -188,7 +188,7 @@ def test_rotated_lena():
     assert (np.sort(results[:, 1]) == np.sort(results_rotated[:, 0])).all()
 
 
-def test_subpix():
+def test_subpix_edge():
     img = np.zeros((50, 50))
     img[:25, :25] = 255
     img[25:, 25:] = 255
@@ -197,6 +197,25 @@ def test_subpix():
     assert_array_equal(subpix[0], (24.5, 24.5))
 
 
+def test_subpix_dot():
+    img = np.zeros((50, 50))
+    img[25, 25] = 255
+    corner = peak_local_max(corner_harris(img), num_peaks=1)
+    subpix = corner_subpix(img, corner)
+    assert_array_equal(subpix[0], (25, 25))
+
+
+def test_subpix_no_class():
+    img = np.zeros((50, 50))
+    subpix = corner_subpix(img, np.array([[25, 25]]))
+    assert_array_equal(subpix[0], (np.nan, np.nan))
+
+    img[25, 25] = 1e-10
+    corner = peak_local_max(corner_harris(img), num_peaks=1)
+    subpix = corner_subpix(img, np.array([[25, 25]]))
+    assert_array_equal(subpix[0], (np.nan, np.nan))
+
+
 def test_subpix_border():
     img = np.zeros((50, 50))
     img[1:25,1:25] = 255
diff --git a/skimage/morphology/greyreconstruct.py b/skimage/morphology/greyreconstruct.py
index 3fffd28e..d302c2c0 100644
--- a/skimage/morphology/greyreconstruct.py
+++ b/skimage/morphology/greyreconstruct.py
@@ -131,11 +131,11 @@ def reconstruction(seed, mask, method='dilation', selem=None, offset=None):
     if selem is None:
         selem = np.ones([3] * seed.ndim, dtype=bool)
     else:
-        selem = selem.copy()
+        selem = selem.astype(bool, copy=True)
 
     if offset is None:
         if not all([d % 2 == 1 for d in selem.shape]):
-            ValueError("Footprint dimensions must all be odd")
+            raise ValueError("Footprint dimensions must all be odd")
         offset = np.array([d // 2 for d in selem.shape])
     # Cross out the center of the selem
     selem[[slice(d, d + 1) for d in offset]] = False
diff --git a/skimage/morphology/tests/test_reconstruction.py b/skimage/morphology/tests/test_reconstruction.py
index 5b2763e1..f5678c08 100644
--- a/skimage/morphology/tests/test_reconstruction.py
+++ b/skimage/morphology/tests/test_reconstruction.py
@@ -90,9 +90,11 @@ def test_invalid_seed():
 def test_invalid_selem():
     seed = np.ones((5, 5))
     mask = np.ones((5, 5))
-    selem =
-    assert_raises(ValueError, reconstruction, seed, mask, np.ones((2, 2)))
-    reconstruction(seed, mask, np.ones((3, 3)))
+    assert_raises(ValueError, reconstruction, seed, mask,
+                  selem=np.ones((4, 4)))
+    assert_raises(ValueError, reconstruction, seed, mask,
+                  selem=np.ones((3, 4)))
+    reconstruction(seed, mask, selem=np.ones((3, 3)))
 
 
 if __name__ == '__main__':