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

Mon Dec 20 09:54:10 EST 2010

Dear all,

I had some recent conversation with Brad about these details. I explain some
of our current conclusions. We are quite sure that the normalization when
using spacing is not correct due to some weird results when spacing differs
from 1.0 as compared with the recursive implementation.

When I was implementing itk::GaussianDerivativeOperator it was a
modification of the existing itk::GaussianOperator. However, the latter
didn’t have the option to use the image spacing (and still does not). Then I
looked at the implementation of the RecursiveGaussian which divided the
sigma by the spacing when using image spacing. This is where the following
line comes from:

m_Variance /= ( m_Spacing * m_Spacing ); 

since variance is the square of sigma. On the other hand, when one
calculates the derivatives of a Gaussian it results in a polynomial of N
order which coincides with the order of the derivative. I found a rule to
calculate the coefficients of this polynomial which is implemented there.
When replacing x in the polynomial with the actual distance from the origin
I multiplied it by the spacing. 

       for( j=1, k = (m_Order-1)/2; j<polyCoeffs.size(); j+=2, k-- )

          {

          sum += polyCoeffs[j] * vcl_pow( m_Variance,(double)k ) * vcl_pow(
i*m_Spacing,(double)j );

          }

When spacing is 1.0, in the position x = 2 in the kernel the distance is 2
mm. but when the spacing is 1.5 mm then the position should be 3.0 mm. This
seemed obvious, but we found that this multiplication is causing a large
deviation with inversion of signs in the eigenvalues with respect to the
recursive version. Maybe this is also inconsistent with dividing the
variance by the spacing (shouldn’t it be multiplied instead??? or x divided
by the spacing???).

If this multiplication is removed, as Brad has confirmed, results seem to
differ only in magnitude. We should also confirm that the normalization in
terms of multiplication factors etc. is the same. This could account for
some differences in magnitude in both implementations (recursive and
discrete). For example the discrete Gaussian is implemented such that it
approximates a unit integral.

Best regards

Iván

De: Luca Antiga [mailto:luca.antiga at gmail.com] 
Enviado el: viernes, 17 de diciembre de 2010 18:27
Para: Bradley Lowekamp
CC: Iván Macía; 'Luis Ibanez'; 'Jim Miller'; 'Insight Developers'
Asunto: Re: [Insight-developers] GaussianFilter::NormalizeAcrossScale not
good for scale space analysis

Hi Brad,

On Dec 17, 2010, at 3:51 PM, Bradley Lowekamp wrote:

Luca,

I have not returned to resolve this issue.

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

Correct.

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.

I actually agree with the fix (see my previous email) and I think it should
be committed after we run a few cross-checks.

I was now checking the consistency with Iván's code. The normalization in
itkGaussianDerivativeOperator is

vcl_pow( m_Variance, m_Order/2.0 )

which is consistent with your fix in terms of exponents.

However, the sigma is being divided by the spacing before all operations are
performed

m_Variance /= ( m_Spacing * m_Spacing ); 

but then later on variance and spacing are going back in the neighborhood
operator coefficients

       for( j=1, k = (m_Order-1)/2; j<polyCoeffs.size(); j+=2, k-- )
          {
          sum += polyCoeffs[j] * vcl_pow( m_Variance,(double)k ) * vcl_pow(
i*m_Spacing,(double)j );
          }

Iván, it would be great if you provided a quick hint on what the
normalization ends up being in the end, I'll eventually do the maths myself
when I find some time.

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/fea94cfc99775abf8ee3f04169051bd8ee34
8b23

I didn't have a chance to run it yet, but I'm now collecting the code, fixes
and classes to make a sense of all this.

Thanks a lot

Luca

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/20101220/61b07a1b/attachment.htm>