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

[Luis Ibanez]

Hi Lydia,

Thanks for the implementation of the Lerp,
and the performance test.

I tried it under Linux with gcc3.3 and verified
your observation:

   The time is not significantly different.



It seems that the pixel access and neighbor
computations takes much longer than the
actual evaluation of the interpolation
shape function.

The timing were basically the same when
using the current LinearInterpolator and
when using your implementation of Lerp.

It was actually surprising that the timings
didn't change when compiling for debugging
or for optimization (with -O3).  That leads
to think that we have a strong overhead in
neighbor trafficking and overlap computation.

Whenever we got into performance optimization
for the registration framework, it seems to be
worth to revisit this issue.  It may even be
interesting to create explicit especializations
of the interpolator per image dimension.

Interpolators can easily be a bottleneck in
the process of evaluating ImageMetrics.


 
         Luis



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

>Hi Mark and Luis,
>
>A couple of years back I did compare the current version with a version
>which explicitly reduces it to 1D lerp (hard coded for 3D) and as I
>recall I didn't get much of a speedup.
>
>I did play around with the interpolator this afternoon - using both
>reduction to 1D lerp and Stephen's trick. It seems to be actually slower
>- but it could be my implementation. I have attached the code for those
>who would like to play with it :-)
>
>Cheers,
>Lydia
>
>  
>
>>-----Original Message-----
>>From: Luis Ibanez [mailto:luis . ibanez at kitware . com]
>>Sent: Wednesday, June 25, 2003 8:42 PM
>>To: Lydia Ng
>>Cc: Stephen R. Aylward; Mark Foskey; Derek W Cool;
>>    
>>
>jisung at email . unc . edu;
>  
>
>>piloo at email . unc . edu; Insight-Developers (E-mail)
>>Subject: Re: [Insight-developers] Re: Bug somewhere in
>>ResampleImageFilter?
>>
>>
>>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
>>>
>>>
>>>
>>>
>>>
>>>      
>>>
>
>
>  
>