Add kernels for computing eigenvalues of hessian matrices in 2d and 3d images

[maarzt]

I added OpenCL kernels and tests for the computation of the hessian eigenvalues in a 2d and 3d images.

The code for the respective CLIJ ops still needs to be done.

[haesleinhuepf]

Hey Matthias @maarzt ,

fantastic! I was looking forward to this day. Thanks!

I just added the API-method to make the CLIJ-Ops thingy work. Would you please mind checking if the documentation is correct?

Furthermore, please run the demo macro I added and check if the results on the T1 head dataset appear reasonable.

Thanks again! I'm happy to see this evolving :-)

Cheers, Robert

[haesleinhuepf]

Also, do you think it would be easy / straightforward to implement Tubeness? All I know it has something to do with Hessian...

[maarzt]

Hi Robert,

You asked me to verify if the results look good for the T1 head example. This made me think a bit. I don't actually know what good look like. I'm sure that the eigenvalues are calculated correctly. But I'm a lot less certain about the calculation of the hessian matrix. I found this paper Schaar 2007, that does use slightly different derivative filters. The tubeness filter points to this paper Sato et al. 1997 which I should have a look at.

Do you know about a paper, that describes an algorithm for calculating the hessian of an image?

There are also this open points:

1. It is often suggested to calculate the hessian eigenvalues after smoothing the image with a gaussian blur. Should the CLIJ op have a build in gaussian blur?
2. The current implementation will give poor results for anisotropic images.

Bests, Matthias

[maarzt]

More specifically I think we still need to make the following decisions:

How do we calculate derivatives for the hessian matrix?

The general approach is to convolve the image with separable kernels:

But the question is which kernels should be used, all D₂ and D₁ is pretty clear, as they clearly calculate the first or second derivative. But the there a several options to add blurring kernels I₁ and I₂: a) Current solution: D₂= [1 -2 1], D₁=[-0.5 0 0.5], I₁=I₂=[1] b) Sobel filters: D₂= [1 -2 1], D₁=[-0.5 0 0.5], I₁=I₂=[0.25 0.5 0.25] c) Schaar 2007: D₂= [1 -2 1], D₁=[-0.5 0 0.5], I₁=[0.12026,0.75948,0.12026], I₂=[0.21478,0.57044,0.21478]

Should we consider integrating a Gaussian blur into the op?

a) Yes, add a parameter sigma. And run a kernel to blur the image if sigma is not zero. b) No.

Should we enable the plugin to deal with z-spacing correctly?

a) Always assume the image is isotropic. b) Add an parameter the specifies the z scaling. Scale I_zz, I_xz and I_yz accordingly. c) Add three parameters for the voxel size in x, y, z. Scale all derivatives accordingly.

[haesleinhuepf]

Do you know about a paper, that describes an algorithm for calculating the hessian of an image?

Hm. Not about computation of an image, but I think ilastik uses those images quite successfully. A quick look suggests they also use some kind of HessianOfGaussian filter. I dived a bit deeper and found they used the vigra library, but I can't find the correct implementation of the filter. I traced it until here. Maybe @k-dominik can give us a hint how they compute the Hessian images? (-:

How do we calculate derivatives for the hessian matrix?

I have no idea and zero experience. But I would be ok if we have multiple Hessian operations. One with good default parameters and one that takes a D1, D2, I1 and I2 or a matrix of those 4x3 numbers as input.

Should we consider integrating a Gaussian blur into the op?

a) Sure, why not! I like the idea of not blurring in case sigma=0

Should we enable the plugin to deal with z-spacing correctly?

a) If we mention that it needs an isotropic input in the docs, that should be ok. It wouldn't be the first operation which has this constraint. There is a function called makeIsotropic, I typically use in advance in such cases.

Thanks for your efforts!

Cheers, Robert

[maarzt]

Github search finds 200K matches for hessian eigenvalues.

It's a rabbit hole, let us stick with the current implementation. Here is a quick overview of what I found:

The fiji VIB-lib uses the same kernels I used, maybe we try to stay consistent with this implementation as the Tubeness filter is also based on it.

scikit-image has skimage.filters.hessian they recursively apply the first order finite central difference kernel [-0.5 0 0.5] to calculate the second order derivatives. Which certainly isn't optimal. But it's probably good enough.

imglib2-algorithm HessianMatrix does the same as scikit-image.

The ImageScience library samples the analytically calculated gaussian derivatives and uses this as a kernel for the calculation of the image derivatives. The problem here is that the error introduced by sampling the gaussian derivatives leads to poorly scaled derivatives.

Ilastik uses sampled gassian derivatives approach too. But the code compensates for the sampling error.

TrainableWekaSegmentation in 2d uses the sobel filter twice to calculate second order derivatives. TrainableWekaSegmentation in 3d uses the ImageScience library mentioned above.

And there is also the theoretical work by Hanno Scharr, mentioned int the wikipedia article on sobel filter.

[haesleinhuepf]

Alright. Thanks @maarzt !