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

[Lydia Ng]

> I tried it under Linux with gcc3.3 and verified
> your observation:
>=20
>    The time is not significantly different.

That's interesting.
On XP/.NET I get the lerp version being 4 times slower!
(e.g. 16 sec versus 4 sec both in Release mode ).

- Lydia

> -----Original Message-----
> From: Luis Ibanez [mailto:luis . ibanez at kitware . com]
> Sent: Wednesday, June 25, 2003 10:45 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?
>=20
> Hi Lydia,
>=20
> Thanks for the implementation of the Lerp,
> and the performance test.
>=20
> I tried it under Linux with gcc3.3 and verified
> your observation:
>=20
>    The time is not significantly different.
>=20
>=20
>=20
> It seems that the pixel access and neighbor
> computations takes much longer than the
> actual evaluation of the interpolation
> shape function.
>=20
> The timing were basically the same when
> using the current LinearInterpolator and
> when using your implementation of Lerp.
>=20
> 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.
>=20
> 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.
>=20
> Interpolators can easily be a bottleneck in
> the process of evaluating ImageMetrics.
>=20
>=20
>=20
>          Luis
>=20
>=20
>=20
> ----------------------
> Lydia Ng wrote:
>=20
> >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  =3D  0
> >>
> >>for all neighbors
> >>{
> >>interpolatedValue +=3D  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 =3D=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 =3D=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 =3D 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
> >>>
> >>>
> >>>
> >>>
> >>>
> >>>
> >>>
> >
> >
> >
> >
>=20
>=20
>=20