[Insight-users] Buggy HessianRecursiveGaussian?

[Andreas Keil]

Dear all,

I'm pretty sure that the filter pipeline which is constructed in
HessianRecursiveGaussian does not what it's expected to do.

E.g. assume one wants to compute the Hessian of a 2D image, then
HessianRecursiveGaussian produces the following components:

Hxx = DxGx * DxGx * I
Hxy = DxGx * DyGy * I
Hyy = DyGy * DyGy * I

This is not what's expected since the smoothing is not always performed
once in x and y direction (at least for Hxx and Hyy). Instead, we'd like
to get

Hxx = DxDx * GxGy * I
Hxy = DxDy * GxGy * I
Hyy = DyDy * GxGy * I

and which could be computed in the following way, e.g.

Hxx = DxxGx * Gy * I
Hxy = DxGx * DyGy * I
Hyy = DyyGy * Gx * I

If I deduced this correct from walking through the code, I consider it a
serious bug which should be fixed immediately. The problem also manifests
in higher-dimensional images which affects the commonly used
Hessian3DToVesselnessMeasure, e.g. Here, we get
Hxx = DxGx * DxGx * Gy * I = Dxx * Gx*Gx*Gy * I
which also yields smoothing in the wrong direction.

A question on the side: What's the exact result of applying a
RecursiveGaussianImageFilter with first/second order? Is the derivative of
the Gaussian computed in a way that a division by scaling after applying
this filter is needed or not? If not, then the line

        ot.Set( it.Get() / factor );

in itkHessianRecursiveGaussian.txx would be superfluous.

Regards,
Andreas.