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

[Lydia Ng]

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,=20
"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



> -----Original Message-----
> From: Luis Ibanez [mailto:luis . ibanez at kitware . com]
> Sent: Wednesday, June 25, 2003 1:27 PM
> To: Stephen R. Aylward
> Cc: 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
>=20
> Hi Stephen,
>=20
> It seems that it is better to fix the current
> linear interpolator rather than adding a new
> class.
>=20
> Since your proposed operation is mathematically
> equivalent to the one being used now in the
> interpolator, and has the advantage of not being
> vulnerable to the float representation limitations,
> this is actually like a bug fix.
>=20
>=20
>    Luis
>=20
>=20
> -------------------------
> Stephen R. Aylward wrote:
> > Hi,
> >
> > There is a bit of a problem with the linear interpolator with
rounding
> > errors.   See Mark's email below this one.
> >
> > Current implementation is doing:
> > v =3D overlap * pixel1 + (1-overlap) * pixel2
> > to get the value v that occurs overlap distance between pixels
pixel1
> > and pixel2.   The rounding occurs because of the two multiplications
by
> > overlap and 1-overlap.   Specifically, even if pixel1 and pixel2
have
> > the same value, you'll sometimes get a value less than them when
> > interpolating!   For example, if you have a binary image encoded in
an
> > unsigned char as 0s and 1s, linear interpolation can result in a
"false"
> > being interpolated between two "trues."
> >
> > If instead you do
> > v =3D pixel1 + overlap * (pixel2-pixel1)
> > there are fewer multiplications (faster?) and the rounding error
will
> > not occur.
> >
> > Any suggestions/comments?   I propose to check in a second filter
that
> > implements this, and I am open to naming suggestions!
> >
> > Thanks,
> > Stephen
> >
> >
> >
> > Mark Foskey wrote:
> >
> >> It's rounding.  Along with truncation.  I'm testing an image that
is
> >> constant at 2.  I got the ResampleImageFilter to print out all
pixel
> >> values not equal to 2.0 (before casting to the output pixel type).
I
> >> get 1.9999999999999998, but not for every pixel -- about the same
> >> number as the holes I see when I view it as a byte image.
> >>
> >> If you cast 1.9999999999999998 to unsigned char, you of course get
1.
> >> But if you cast it to float, you get 2.0.  And that's the behavior
I
> >> see when I do the test with a float image.
> >>
> >> If you use NearestNeighborInterpolationFunction, you get no holes,
as
> >> expected since the pixel values aren't modified.  But that has more
> >> problems with aliasing.
> >>
> >> BSplineInterpolationFunction, the other choice, has more roundoff
> >> problems -- that is, more pixels come out to be slightly different
> >> than 2.  This makes sense since there are probably more numerical
ops
> >> involved.  BSpline may have some advantages from a signal
processing
> >> point of view, but I don't know how important they would be for our
> >> purposes.
> >>
> >> You do get the errors with translation and scaling, as long as you
> >> feed it suitably bad numerical values.  Scaling by 0.5 is fine, but
by
> >> something like 1.471 leads to roundoff errors.
> >>
> >> Derek and I talked and he's looking at working with float images as
> >> the output type of the resample image filter.
> >>
> >>>
> >>> If you try a different interpolation scheme, what happens?
> >>>
> >>> I've looked over the linear interpolator - my only concern is that
> >>> they have an early termination criterion (totalOverlap =3D=3D =
1.0).
The
> >>> cost of this at every pixel isn't worth the occasional time that
it
> >>> saves a loop iteration or two (if ever).  I really don't see any
> >>> "real" errors.
> >>>
> >>> Is it really the interpolator?   What happens if you translate a
> >>> really really small amount (again to reveal a rounding error)?
> >>>
> >>> s
> >>>
> >>>
> >
> >
> >
>=20
>=20
>=20
> _______________________________________________
> Insight-developers mailing list
> Insight-developers at itk . org
> http://www . itk . org/mailman/listinfo/insight-developers