From 6fd0cecffdfa1bf3b395ec2b70d7b6cca2b82d31 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Wed, 20 Jan 2016 22:11:23 +0100 Subject: [PATCH] Fix Gaussian kernels in hessian_matrix --- skimage/feature/corner.py | 13 +++++--- skimage/feature/tests/test_corner.py | 50 ++++++++++++++-------------- 2 files changed, 33 insertions(+), 30 deletions(-) diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py index a8fa7d78..f65f30e5 100644 --- a/skimage/feature/corner.py +++ b/skimage/feature/corner.py @@ -149,22 +149,25 @@ def hessian_matrix(image, sigma=1, mode='constant', cval=0): image = _prepare_grayscale_input_2D(image) - # window extent to the left and right, which covers > 99% of the normal - # distribution + # Window extent which covers > 99% of the normal distribution. window_ext = max(1, np.ceil(3 * sigma)) ky, kx = np.mgrid[-window_ext:window_ext + 1, -window_ext:window_ext + 1] - # second derivative Gaussian kernels + # Second derivative Gaussian kernels. gaussian_exp = np.exp(-(kx ** 2 + ky ** 2) / (2 * sigma ** 2)) kernel_xx = 1 / (2 * np.pi * sigma ** 4) * (kx ** 2 / sigma ** 2 - 1) kernel_xx *= gaussian_exp - kernel_xx /= kernel_xx.sum() kernel_xy = 1 / (2 * np.pi * sigma ** 6) * (kx * ky) kernel_xy *= gaussian_exp - kernel_xy /= kernel_xx.sum() kernel_yy = kernel_xx.transpose() + # Remove small kernel values. + eps = np.finfo(kernel_xx.dtype).eps + kernel_xx[np.abs(kernel_xx) < eps * np.abs(kernel_xx).max()] = 0 + kernel_xy[np.abs(kernel_xy) < eps * np.abs(kernel_xy).max()] = 0 + kernel_yy[np.abs(kernel_yy) < eps * np.abs(kernel_yy).max()] = 0 + Hxx = ndi.convolve(image, kernel_xx, mode=mode, cval=cval) Hxy = ndi.convolve(image, kernel_xy, mode=mode, cval=cval) Hyy = ndi.convolve(image, kernel_yy, mode=mode, cval=cval) diff --git a/skimage/feature/tests/test_corner.py b/skimage/feature/tests/test_corner.py index d6c3418f..c1461f08 100644 --- a/skimage/feature/tests/test_corner.py +++ b/skimage/feature/tests/test_corner.py @@ -42,21 +42,21 @@ def test_hessian_matrix(): square = np.zeros((5, 5)) square[2, 2] = 1 Hxx, Hxy, Hyy = hessian_matrix(square, sigma=0.1) - assert_array_equal(Hxx, np.array([[0, 0, 0, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 1, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 0, 0, 0]])) - assert_array_equal(Hxy, np.array([[0, 0, 0, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 0, 0, 0]])) - assert_array_equal(Hyy, np.array([[0, 0, 0, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 1, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 0, 0, 0]])) + assert_almost_equal(Hxx, np.array([[0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, -1591.549431, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]])) + assert_almost_equal(Hxy, np.array([[0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]])) + assert_almost_equal(Hyy, np.array([[0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, -1591.549431, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]])) def test_structure_tensor_eigvals(): @@ -81,16 +81,16 @@ def test_hessian_matrix_eigvals(): square[2, 2] = 1 Hxx, Hxy, Hyy = hessian_matrix(square, sigma=0.1) l1, l2 = hessian_matrix_eigvals(Hxx, Hxy, Hyy) - assert_array_equal(l1, np.array([[0, 0, 0, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 1, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 0, 0, 0]])) - assert_array_equal(l2, np.array([[0, 0, 0, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 1, 0, 0], - [0, 0, 0, 0, 0], - [0, 0, 0, 0, 0]])) + assert_almost_equal(l1, np.array([[0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, -1591.549431, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]])) + assert_almost_equal(l2, np.array([[0, 0, 0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, -1591.549431, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 0]])) @test_parallel()