[Insight-developers] GaussianFilter::NormalizeAcrossScale not good for scale space analysis

[Luca Antiga]

Hi Brad, Luis, Jim,
 sorry I didn't catch the email earlier.

So it looks like my earlier fix to the Hessian did not go down to the root of the problem.

In relation to Brad's initial comment on why the normalization factor is scaled by the spacing, you're right, it really shouldn't be. In fact, since the recursive filter works in pixel units and the sigma at the exponent gets scaled by pixel units, the two scaling factors cancel out, and the integral still depends on the unscaled sigma.

As to normalization factors for the derivative operator, see for example ftp://ftp.nada.kth.se/CVAP/reports/Lin08-EncCompSci.pdf pages 7-9. The normalization factor is t^{\gamma/2} (=\sigma), where \gamma should be chosen appropriately. Clearly, choosing \gamma = 1 achieves scale invariance across scaling transformations, which is what is desirable most of the times and it's what Brad's code does now.

Being Brad fixes all appropriate and going back to the Hessian, the filter was invariant in scale-space, but just because errors were compensating (after my fix, earlier it wasn't working at all). 
In the end the final scaling for the second partial derivatives was correct (sigma^2, or t in Lindeberg's notation), and this was because at line 208 in itkHessianRecursiveGaussianImageFilter the two derivative filters were

(inline comments are contributions to scaling factors coming from the recursive Gaussian filter)

prior to Brad's fix:

 if ( dimb == dima )
        {
        m_DerivativeFilterA->SetOrder(DerivativeFilterAType::SecondOrder);   // sigma
        m_DerivativeFilterB->SetOrder(DerivativeFilterBType::ZeroOrder);         // sigma
        [...]
        }
else
        {
        m_DerivativeFilterA->SetOrder(DerivativeFilterAType::FirstOrder);         // sigma
        m_DerivativeFilterB->SetOrder(DerivativeFilterBType::FirstOrder);         // sigma
        [...]
        }

resulting in a total scaling of sigma^2.

after Brad's fix:

 if ( dimb == dima )
        {
        m_DerivativeFilterA->SetOrder(DerivativeFilterAType::SecondOrder);   // sigma^2
        m_DerivativeFilterB->SetOrder(DerivativeFilterBType::ZeroOrder);         // 1
        [...]
        }
else
        {
        m_DerivativeFilterA->SetOrder(DerivativeFilterAType::FirstOrder);         // sigma
        m_DerivativeFilterB->SetOrder(DerivativeFilterBType::FirstOrder);         // sigma
        [...]
        }

still resulting in a total scaling of sigma^2. The latter is the correct one.

Great that you dug this bug out and had a fix for it.

Thanks


Luca


PS: just an extra comment: Lindeberg suggests to choose \gamma according to the specific blob/ridge/etc detector at hand. Maybe the normalization factors could be exposed at the user level, should one need to follow Lindeberg and choose \gamma=1/2 for edge detection, or \gamma=3/4 for ridge detection, etc.



On Nov 23, 2010, at 6:24 PM, Luis Ibanez wrote:

> On Tue, Nov 23, 2010 at 11:21 AM, Bradley Lowekamp <blowekamp at mail.nih.gov> wrote:
> 
> My "n" refers to the order of the derivative not the dimension. This will perform exactly what is required for N-Dimensions of separated n-order derivatives.
> 
> 
>> 
>> So, by applying the current normalization, (if only
>> one Sigma per dimension) you end up with a 
>> Sigma^N normalization.
> 
> 
> We are talking about different Ns.
> 
> 
> 
> That clarifies my misunderstanding.
> 
> Can I suggest that we use "K" instead
> of "N" for the order for the derivative ?
> 
> 
> Scaling by:
> 
>                      Sigma^ K
> 
> where K is the order of the derivative,
> makes sense.
> 
> 
> I would suggest that we test this by 
> using an input image with an impulse
> 
> E.g. a image with all pixels set to zero,
> and a single pixel in the center, set to
> a nominal value (1000.0 for example ?).
> 
> In that way we should be able to verify
> the correct numerical behavior for every
> combination.
> 
> I'll take a closer look at your branch in 
> github.
> 
> 
>      Thanks
> 
> 
>           Luis
> 
> _______________________________________________
> 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://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-developers

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/mailman/private/insight-developers/attachments/20101124/3b0028dc/attachment.htm>