Expanded on Hessian eigenvalue test

[mdmueller]

This is a start on addressing #1393. Unfortunately I'm not familiar with computing the Hessian via convolution for an image, so I can't think of a good way to test the hessian_matrix function in more detail. I did try testing the hessian_matrix_det function by matching against the determinant of the matrix returned by hessian_matrix, but the values weren't even close, and varied considerably based on sigma. For example,

In [3]: im = np.random.rand(5, 5)

In [4]: im
Out[4]: 
array([[ 0.29763742,  0.05449981,  0.71780469,  0.51500638,  0.63904432],
       [ 0.67701664,  0.00113305,  0.05631246,  0.004884  ,  0.03438824],
       [ 0.03267068,  0.26552206,  0.31953294,  0.8475709 ,  0.14632559],
       [ 0.18858282,  0.44127783,  0.7114495 ,  0.45750935,  0.89106171],
       [ 0.52487659,  0.74190839,  0.49202619,  0.30676331,  0.61698863]])

In [5]: Hxx, Hxy, Hyy = hessian_matrix(im, sigma=5)

In [6]: det = hessian_matrix_det(im, sigma=5)

In [7]: h0 = [[Hxx[0][0], Hxy[0][0]], [Hxy[0][0], Hyy[0][0]]]

In [8]: h0
Out[8]: 
[[1.8388767163921964, 0.00033013261040972804],
 [0.00033013261040972804, 1.7744890649800114]]

In [9]: det[0][0], h0[0][0]*h0[1][1] - h0[0][1]*h0[1][0]
Out[9]: (0.0031705952509722057, 3.263066516096762)

I'm probably missing something here, so if anyone could clear that up I'd be happy to write tests for hessian_matrix_det and hessian_matrix.

[coveralls]

Coverage increased (+0.01%) to 92.71% when pulling 879ce83 on mdmueller:test-hessian into e91dcba on scikit-image:master.

[ahojnnes]

Thanks also here for investigating. However, I don't think there is a problem with the eigenvalue computation:

mat = np.random.rand(5, 5)
Hxx, Hxy, Hyy = hessian_matrix(mat, sigma=5)

H = np.zeros((5, 5, 2, 2))
H[..., 0, 0] = Hxx
H[..., 0, 1] = Hxy
H[..., 1, 0] = Hxy
H[..., 1, 1] = Hyy

l1, l2 = hessian_matrix_eigvals(Hxx, Hxy, Hyy)

for i in range(5):
    for j in range(5):
        print np.linalg.det(H[i, j] - l1[i, j] * np.eye(2))
        print np.linalg.det(H[i, j] - l2[i, j] * np.eye(2))

Note, that hessian_matrix_det only computes the approximate determinant using integral images.

The problem of the hessian_matrix function at the moment is, that the generated kernel (second derivative of the Gaussian) is not correct for small sigmas due to discretization issues.

It would be highly appreciated, if that kernel function would be fixed. As a start, these references should be helpful: http://campar.in.tum.de/Chair/HaukeHeibelGaussianDerivatives, http://campar.in.tum.de/twiki/pub/Chair/HaukeHeibelGaussianDerivatives/gaussian2d.m

[stefanv]

This is a high priority; if you can make any progress on this I would be much obliged.

[mdmueller]

Sure, I can look into it. By any chance is there any sample data or some such where the Hessian is known, so that I can see where the error is for small sigma?

[mdmueller]

Sorry, I think I'm missing something--what is the error for small sigma? I notice that the kernel is extremely close to the identity, e.g.

array([[ -3.68287522e-42,   1.92874985e-22,  -3.68287522e-42],
       [ -1.90946235e-20,   1.00000000e+00,  -1.90946235e-20],
       [ -3.68287522e-42,   1.92874985e-22,  -3.68287522e-42]])

but if I'm not mistaken this is correct since sigma is small (0.1) and, e.g., np.exp(-1/(.1)**2) is 3.7e-44.

[stefanv]

@ahojnnes What are your thoughts around this PR?

[ahojnnes]

I am relatively certain that the Hessian matrix implementation is right now. Something is bogus with the hessian_matrix_det function though. E.g. the updated value in ll128-129 is not used at all...

[ahojnnes]

They seem to behave differently in terms of the return scale of values, e.g.:

import numpy as np
from skimage.feature import *
im = np.random.rand(100, 100)
Hxx, Hxy, Hyy = hessian_matrix(im, sigma=5)
det = hessian_matrix_det(im, sigma=5)
h0 = [[Hxx[50][50], Hxy[50][50]], [Hxy[50][50], Hyy[50][50]]]
det[50][50], h0[0][0]*h0[1][1] - h0[0][1]*h0[1][0]

im[49:51, 49:51] = 100
Hxx, Hxy, Hyy = hessian_matrix(im, sigma=5)
det = hessian_matrix_det(im, sigma=5)
h0 = [[Hxx[50][50], Hxy[50][50]], [Hxy[50][50], Hyy[50][50]]]
det[50][50], h0[0][0]*h0[1][1] - h0[0][1]*h0[1][0]