[Insight-users] Question about workings of ExponentialDeformationFieldImageFilter

[Anja Ende]

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