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

[Bradley Lowekamp]

Luca,

I have not returned to resolve this issue.

I assume that you are looking at the normalization factor in the current RecursiveHessain?

According to my calculations the is a mistake with applying the pixel spacing incorrectly to the normalization. That is we are being normalized by (sigma/spacing) not the physical size of  the sigma. For isotropic voxels this has only a constant scalar factor for being off. If however the pixels are anisotropic, the results will be off it a more complex fashion.

One important principle behind this normalization is that the value of a feature should be independent of the size or scale of the feature. This implies that if you scale the pixel spacing and the sigma, but don't change the pixel, you should obtain the same value. This is easy to test, I believe that the Recursive Hessian will not pass this currently. This is an extremely important property because it allows the use of different resolution images, or image pyramids.

Luca did you have a chance to look at this test I wrote:
https://github.com/blowekamp/ITK/commit/fea94cfc99775abf8ee3f04169051bd8ee348b23

Brad

On Dec 17, 2010, at 3:34 AM, Luca Antiga wrote:

> 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:
>> 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/mailman/private/insight-developers/attachments/20101217/60fafade/attachment-0001.htm>