[Insight-users] Question about workings of ExponentialDeformationFieldImageFilter

[Anja Ende]

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.

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.

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);
>
> 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.
>
> Many thanks,
>
> Anja
>