[Insight-users] Hessian format for SymmetricEigenAnalysisImageFilter

[Luis Ibanez]

Hi Alex,

If I understand correctly, you are attempting to compare

                      Gxx     versus      Gx( Gx )

If so, part of the problem is that when we do Gxx, the Gaussian
smoothing is applied only once,  while, when using Gx(Gx) the
Gaussian Smooting is applied twice.

And... depending on how you are doing it,... it may not be the
equivalent to Gxx, since the Gxx that we do requires:

                                 Dxx( G )

which means that along X we do the equivalent of convolving
with the second derivative of the Gaussian, while at the same
time we still do the equivalent of convolve with a Gaussian along
Y and along Z.


The Scale normalization by the Gaussian will also apply
twice on the Gx(Gx) version.

All that said,
It will be very healthy to run this on a synthetic image for
which the theoretical values can be computed, so we have
a reference value.

It may well be that you might have uncovered a bug...


   Regards,


       Luis

---------------------------------------------------------------------------
On Mon, Dec 14, 2009 at 7:30 PM, Oleksandr Dzyubak <adzyubak at gmail.com> wrote:
>
> Hi Luis,
>
> Looks like I am terribly missing something here...
>
> OK. For my testings I built objects with Gaussian profiles.
> These are simple: Gaussian lines with various sigmas extruded
> into Z-direction thus I got "Gaussian planes".
> Objects: 3D Gaussian planes with Sigma=1,2,3,4,5,10 and intensities I=100.
>
> Now I use RecursiveGaussians and HessianRecursiveGaussian from the ITK stock
> library
> and calculate the second order derivative across the Gaussian profiles:
>  Hxx -- Hessian x-components, Gxx - Gaussian second order derivative,
> and GxGx - Gaussian second order derivative taken sequentially as Gx(Gx)
> using Gaussian first order derivative.
>
> Of course, all that was done with "NormalizeAcrossScale" ON and OFF.
> Afterwards the maximum of absolute values were measured.
> The results I summarized in the table below.
>
> Sigma    ITK_Hxx_OFF    ITK_Gxx_OFF    ITK_Hxx_ON    ITK_Gxx_ON
> ITK_GxGx_OFF    ITK_GxGx_ON
>
> 1                 35.77                 35.77                   35.77
>     35.77                   19.47               19.47
> 2                  8.92                   8.92                    71.33
>     71.33                     4.87               19.47
> 3                   3.96                  3.96                  106.95
>   106.95                    2.16                19.47
> 4                   2.23                  2.23                  142.56
>   142.56                    1.22                19.47
> 5                   1.43                  1.43                  178.18
>     178.18                    0.78                19.47
> 10                 0.36                  0.36                   356.29
>     356.29                    0.19                19.47
>
>
> As you can see from the table, Hxx is always equal to Gxx for both
> "NormalizeAcrossScale" ON and OFF.
> Ok. No surprises here. So far so good.
> Should not be the values Scale indipendent though?
>
> Well, now lets take the second order derivative sequentially using Gx(Gx).
> Now there is a dramatic difference here.
> a) GxGx values do not match Gxx and Hxx anymore.
> b) "ITK_GxGx_ON" version apparently demonstrates Scale Space invariance as
> it should be.
>
>
> Luis, I would highly appreciate any comments on that.
>
> Best regards,
>
> Alex
>
>
> On Fri, Dec 11, 2009 at 7:39 PM, Luis Ibanez <luis.ibanez at kitware.com>
> wrote:
>>
>> Hi Alex,
>>
>>
>> The order for the Hessian follows the same
>> pattern of the class
>>
>>   Insight/Code/Common/
>>          itkSymmetricSecondRankTensor.h
>>
>> The elements store only the upper right triangle
>> of the actual matrix.
>>
>> For example, in 3D you get
>> the following pattern
>>
>>                0  1  2
>>                    3  4
>>                        5
>>
>> All that to say:
>>
>>   Yes,
>>   the order that you are using is correct:
>>
>>                Dxx  Dxy    Dxz
>>                         Dyy   Dyz
>>                                  Dzz
>>
>>
>> BTW:
>> You may want to use (or at least look) the filter:
>>
>>   Insight/Code/BasicFilters/
>>        itkHessianRecursiveGaussianImageFilter.h
>>
>>
>>
>>     Regards,
>>
>>
>>             Luis
>>
>>
>> ---------------------------------------------------------------
>> On Fri, Dec 11, 2009 at 12:33 PM, Oleksandr Dzyubak <adzyubak at gmail.com>
>> wrote:
>> >
>> > Dear users,
>> >
>> > After I built the individual Hessian components,
>> > in what order should I fill them up to further feed
>> > into the SymmetricEigenAnalysisImageFilter?
>> >
>> > At the moment I am using the order below.
>> >
>> > H[0] = Dxx
>> > H[1] = Dxy
>> > H[2] = Dxz
>> > H[3] = Dyy
>> > H[4] = Dyz
>> > H[5] = Dzz
>> >
>> > Is this correct order?
>> >
>> > Thanks,
>> >
>> > Alex
>> >
>> >
>> > _____________________________________
>> > Powered by www.kitware.com
>> >
>> > Visit other Kitware open-source projects at
>> > http://www.kitware.com/opensource/opensource.html
>> >
>> > Kitware offers ITK Training Courses, for more information visit:
>> > http://www.kitware.com/products/protraining.html
>> >
>> > Please keep messages on-topic and check the ITK FAQ at:
>> > http://www.itk.org/Wiki/ITK_FAQ
>> >
>> > Follow this link to subscribe/unsubscribe:
>> > http://www.itk.org/mailman/listinfo/insight-users
>> >
>
>