[Insight-users] Question about workings of ExponentialDeformationFieldImageFilter

[Tom Vercauteren]

Hi Anja,

Please see below.

On Wed, May 18, 2011 at 12:08, Anja Ende <anja.ende at googlemail.com> wrote:
> Hello all,
>
> I am trying to understand the inner workings of the
> ExponentialDeformationFieldImageFilter that is used in the
> diffeomorphic demons registration. I have a few questions as I try to
> understand the scaling and squaring method behind it.
>
> - The first question is regarding the automatic calculation of the
> number of iterations/squaring steps that is needed. In the default
> implementation of the code, it ensures that the initial first order
> approximation is diffeomorphic. Here is used the fact that
> exp(Phi/2^N) = Phi/2^N. However, shouldn't it be exp(Phi/2^N) =
> Identity + Phi/2^N assuming that we are doing the Taylor series
> expansion around 0?

Correct. Mathematically you have
  exp(Phi/2^N) \approx Identity + Phi/2^N
but in the code you always use a vector field which is a deviation
from the identity (displacement field) hence
you represent the spatial transformation Identity + Phi/2^N by
  Identity + Phi/2^N - Identity = Phi/2^N


> - Second thing that I am confused about is the actual calculation of N
> from the inequality max(norm(Phi)/2^N) < 0.5*pixelspacing. I tried to
> solve this inequality (I must confess I am a bit weak on Maths), but
> could not come to the expression that is used in the code which is:
>
> maxnorm2 /= vnl_math_sqr(minpixelspacing);
>
> InputPixelRealValueType numiterfloat = 2.0 + 0.5 *
> vcl_log(maxnorm2)/vnl_math::ln2;

  max(norm(Phi)/2^N) = max(norm(Phi))/2^N < 0.5*pixelspacing.
is equivalent to
  2^N > max(norm(Phi)) / (0.5*pixelspacing)
which is equivalent to
  N > log2(max(norm(Phi))/pixelspacing) - log2(0.5)
which is equivalent to
  N > 0.5*log2( (max(norm(Phi))/pixelspacing)^2 ) + 1