[Insight-users] Hessian format for SymmetricEigenAnalysisImageFilter
Oleksandr Dzyubak
adzyubak at gmail.com
Tue Dec 22 10:25:54 EST 2009
Hi Luis,
I have already done all possible tests using both formulas
and synthetic images in 1D, 2D, and 3D space.
1) Symbolic simulations (using Maxima).
2) Numerical simulations (using Matlab).
3) Test on synthetic images (Matlab, ITK)
All the tests confirmed what I said before and summarized in the table
below.
I would dare to conclude that it is a bug which affects all the algorithms
that use high order Gaussian derivatives including Hessians, of course.
When you form Hessians, you have to use second order and cross derivatives
thus you get incorrect scaling anyway.
Since the incorrect scaling, the results from
SymmetricEigenAnalysisImageFilter
are also affected introducing the errors into consecutive algorithms.
I hope we can solve this problem since we know what to do.
Thanks for your help and support and
Merry Christmas!
Alex
Luis Ibanez wrote:
> 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
>>>>
>>>>
>>
More information about the Insight-users
mailing list