[Insight-users] Question about workings of ExponentialDeformationFieldImageFilter

[Tom Vercauteren]

Hi Anja,

Please see below.

On Fri, May 20, 2011 at 01:14, Anja Ende <anja.ende at googlemail.com> wrote:
>
> Hi Tom, itk-users
>
> One final bit of question regarding this filter:
>
> In the final bit of of the filter where you do the recursive squaring
> of the by composing the transformation by itself, why do you pass the
> result each iteration through an adder. From reading the paper, it
> seems that all one needs to do is do the scaling and recursive
> squaring and then use the result in the update step. The presence of
> the adder has me stumped. By looking at the code, it seems that at
> each iteration, we square the current displacement field and then add
> it to itself to compute the input to the warper for the next
> iteration. I am unable to understand why we need this addition step.
> Also, wouldn't this addition also might result in us leaving the space
> of diffeomorphisms? I am pretty sure I have a big hole in my
> understanding of the code here and I would be very grateful if you
> could comment on these queries.

This addition is need to compute the composition. Let u, v and w be
displacement fields, we want w such that
  Id+w = (Id + u) o (Id + v)
          = Id + v + u_warped_by_v

meaning
  w = v + u_warped_by_v


> Also, just to clarify,
> - the input to the filter is a displacement field and NOT a position
> field, am I correct? I have managed to confuse myself regarding this
> as well.

correct

>
> Many thanks and I would be really grateful for any reply you can give me.
>
> Thanks,
>
> Anja
>
> On 19 May 2011 11:24, Anja Ende <anja.ende at googlemail.com> wrote:
> > Hi Tom,
> >
> > I think there is a problem in the last line of the code for this calculation:
> >
> >>  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
> >
> > Well, first thing is that the change of base (as in code) can be done
> > without the last line:
> >
> > So, N can be simply
> > N > log(max(norm(Phi))/pixelspacing)/vnl_log_2 + 1
> >
> > In your code, you do:
> > maxnorm2 /= vnl_math_sqr(minpixelspacing);

maxnorm2 is the square of the norm. Hence
  maxnorm2 /= vnl_math_sqr(minpixelspacing);
provides
  (max(norm(Phi))/pixelspacing)^2


> > However, according to the math it should be:
> > maxnorm2 /= minpixelspacing;
> > maxnorm2 = vnl_math_sqr(maxnorm2);
> >
> > I think the code segment should be as follows:
> >
> > maxnorm2 /= minpixelspacing;
> > maxnorm2 = vnl_math_sqr(maxnorm2);
> >
> > InputPixelRealValueType numiterfloat = 2.0 + 0.5 *
> > vcl_log(maxnorm2)/vnl_math::ln2;
> >
> > or simply:
> > maxnorm2 /= minpixelspacing;
> > InputPixelRealValueType numiterfloat = 2.0 + vcl_log(maxnorm2)/vnl_math::ln2;
> >
> > Maybe I am mistaken and have messed something up?? I would be grateful
> > for your comment on this.

Your computations look wrong. Please see above.

> I have a small related question. I was reading your SPIE paper from
> last year where you use NMI forces to adapt the registration for
> inter-modal registrations. In this paper, you use a conjugate gradient
> based optimisation scheme where you use the NMI gradient as the
> driving force. I wanted to ask whether it is theoretically sound to
> use the conjugate gradient scheme as there is a line in the mentioned
> paper about maybe this approach being on a bit shaky mathematical
> foundation. I was wondering if you would be kind enough to comment on
> that as I am quite keen to try and implement this in the ITK
> framework. I was looking for a subsequent publication regarding this
> but there seems to be very little work on this. I would like to try
> and code this up in the ITK framework.

I haven't put myself back into the NMI forces. The mathematical
justifications are thus still unclear to me.

Hope this helps,
Tom