[Insight-developers] Re: Bug somewhere in ResampleImageFilter?

[Luis Ibanez]

Hi Lydia,

I agree with you in that interpolation in a binary image doesn't make
much sense. However the same rounding error will degrade the
interpolation of any image with integer pixel type (e.g. char, short
long, signed or unsigned).


Extending the trick to N-D should be possible, however the
code in itkLinearInterpolateImageFunction.txx will have
a very different look.

Currently in itkLinearInterpolateImageFunction.txx we
are computing the full contribution factor for each neighbor
of the continuous index to be evaluated. We end up computing
somehting like

interpolatedValue  =  0

for all neighbors
{
interpolatedValue +=  K1 * K2 * K3 * K4.... * KN  * neighborPixelValue
}


Where "N" is the image dimension and the terms Ki are
"(overlap)" or "(1-overlap)" along dimension "i".
The number of neighbors in the delaunay region used for linear
interpolation is 2^N.

The drawback of the current computation is that we are not
taking advantage of the potential factorization in the products
of the Ki terms,

So in an N-D cell we are computing

      N * (2 ^ N)         Multiplications
      N * (2^ (N-1) )   Sums

(e.g. in 3D == 24 mult and 12 sums )

when in principle the same task could be done with

    2^N-1             Multiplications
    2^(N+1) -2     Sums

(e.g. in 3D == 7 mult and 14 sums )

when the terms are factorized.

The implementation will require a bit of auxiliary
memory in which the interpolations will be done
along a single direction at a time, and each time
the dimension of the problem will be reduced by
one.


Since this is the most commonly used interpolator
in the registration framework, it is worth to give
it a try to implement the factorization.



    Luis


----------------------
Lydia Ng wrote:

>Hi Stephen and Luis,
>
>I am not totally convinced that we should classify the "rounding error"
>as a bug.
>
>I am not quite sure that using a linear interpolation for a binary image
>makes sense if what you really want back is a warp/transform binary
>image. One way to this more robustly is to create a signed distance map
>from the binary image, then warp/transform the map and threshold to get
>the warped binary image - ala
>S. P. Raya and J. K. Udupa, 
>"Shape-based interpolation of multidimensional objects",
>IEEE Trans. on Medical Imaging, vol. 9, no. 1, pp 32-42, 1990.
>
>
>  
>
>>>If instead you do
>>>v = pixel1 + overlap * (pixel2-pixel1)
>>>there are fewer multiplications (faster?) and the rounding error
>>>      
>>>
>will
>  
>
>>>not occur.
>>>      
>>>
>
>I agree this would help with rounding error. Does this trick extend to
>N-D?
>
>In any case, if we are to replace/modify the existing linear
>interpolator we really should do some performance testing so that we can
>systematically compare before and after.
>
>My 2 cents worth ....
>
>- Lydia
>
>
>
>  
>