[Insight-users] Buggy HessianRecursiveGaussian?

[Andreas Keil]

Wow Luis,

regarding A, you're just too fast :-)   Thank you really for this very
fast reaction to my report! I wanted to suggest another solution for this
bug which should reduce the lines of code and be easier readable:

What about setting up a pipeline of filters

   Input -> F_0 -> F_2 -> ... -> F_N-1 -> Output,

where N is the input image's dimension and we do

   for (i = 0; i < N; ++i) F_i.SetDirection(i)

After this general setup (staying the same for all entries of the
Hessian), we already obtain the correct smoothing. Inside the loop

   for ( unsigned int dima=0; dima < ImageDimension; dima++ )
     for ( unsigned int dimb=dima; dimb < ImageDimension; dimb++ )

we'd then only need to set the order of derivatives by specifying

   for (i = 0; i < N; ++i) {
     if ((i == dima) && (i == dimb)) {
       F_i->SetOrder(GaussianFilterType::SecondOrder);
     } else if ((i == dima && i != dimb) || (i != dima && i == dimb)) {
       F_i->SetOrder(GaussianFilterType::FirstOrder);
     } else if ((i != dima) && (i != dimb)) {
       F_i->SetOrder(GaussianFilterType::ZeroOrder);
     } else {
       // some error routine since this case can't happen
     }
   }

No more complicated while loops inside the for loops now. What do you
think?

Regarding B: Since Sigma has to be specified in "world coordinates"
(comment to Set/GetSigma), I assumed that the recursive implementation of
Gaussian smoothing already handles the spacing correctly. And building the
derivative of a continuous Gaussian doesn't need any division by scaling.
If the derivative of the recursive Gaussian implementation would be
implemented by a simple difference, I'd understand the division, but
(without knowing any details of R. Deriche's recursive approximation) I
assumed that the derivative is built differently...

Best regards,
Andreas.



-----Original Message-----
From: Luis Ibanez [mailto:luis.ibanez at kitware.com] 
Sent: Thursday, April 26, 2007 00:19
To: Andreas Keil
Cc: insight-users at itk.org
Subject: Re: [Insight-users] Buggy HessianRecursiveGaussian?


Hi Andreas,


A) Regarding the computation of the Hessian diagonal elements:

    You seem to be right in your analysis of the code.


         This indeed seems to be a bug.

    This issue has now been logged in the bug tracker as Bug # 4934:
    http://public.kitware.com/Bug/bug.php?op=show&bugid=4934&pos=0


    We are working on a fix for it, that most likely
    will involve to have a different configuration
    of the minipipeline for computing the diagonal
    elements of the Hessian.


B) Regarding the factor correction, we still have to double
    check the result in an experiment, but, my recollection
    is that the intensity values of the gradients are computed
    in a pixel-spacing independent way. We did a similar
    division in the GradientRecursiveGaussianImageFilter.
    Of course, that doesn't means that the operation is correct.

    We will run a test that will tell us for sure.



  Thanks


     Luis


-----------------------
Andreas Keil wrote:
> 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.