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

Fri Dec 17 03:34:05 EST 2010

Hi all,
 just out of curiosity, have there been progress on testing Brad's fixes?
I took a look and apparently the normalization factor seems to be correctly following what we were writing.
I propose to run Iván's very good testing set-up on both implementation and eventually get it into ITK.
If somebody is already at it or is willing to do it, that's great. I can help if needed.
Best

Luca

On Dec 15, 2010, at 6:17 PM, Iván Macía wrote:

> Hi,
>  
> I just wanted to comment that the GaussianDerivativeOperator (in Review) is or should be consistent as well with this behavior. I remember that I struggled a bit with Lindeberg’s papers and normalization factors and took this into account when implementing the class but I am not 100% sure that the implementation is correct, especially regarding the spacing.
>  
> http://www.insight-journal.org/browse/publication/179
>  
> In GaussianDerivativeOperator::GenerateCoefficients()
>  
> CoefficientVector coeff;
>  
> // Use image spacing to modify variance
> m_Variance /= ( m_Spacing * m_Spacing );
>  
> // Calculate normalization factor for derivatives when necessary
> double norm = m_NormalizeAcrossScale && m_Order ? vcl_pow( m_Variance, m_Order/2.0 ) : 1.0;
>  
> // (sorry Luis about the evil ternary operator J this was long ago)
>  
> ……
>  
> sum *= norm / vcl_pow( m_Variance, static_cast<int>( m_Order ) );
>  
> Did also some pencil work at the time. I think it is correct but it would be nice if someone else could also have a look at the consistency of the calculations.
>  
> There are some other operators that are built using this one (see the paper), I think once this operator is correct the rest are consistent as they are combined correctly.
>  
> Thanks
>  
> Iván
>  
>  
> De: insight-developers-bounces at itk.org [mailto:insight-developers-bounces at itk.org] En nombre de Luca Antiga
> Enviado el: miércoles, 24 de noviembre de 2010 23:44
> Para: Luis Ibanez
> CC: Jim Miller; Insight Developers
> Asunto: Re: [Insight-developers] GaussianFilter::NormalizeAcrossScale not good for scale space analysis
>  
> 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/20101217/ca10d825/attachment.htm>