[Insight-developers] Re: Bug somewhere in ResampleImageFilter?
Mark Foskey
mark_foskey at unc . edu
Wed, 25 Jun 2003 22:22:38 -0400
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.
It seems to me that this is a lot to ask of somebody who may simply
want to, say, rotate and translate a mask before overlaying it onto
another image. Something as frequently used as a resampling filter
ought, if practical, to behave the way users would expect. Of course,
the "if practical" part can be a big issue.
>> Aylward:
>>>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?
You can recursively reduce an nD problem to an (n-1)D problem, but I
don't know how numerically nice that is.
For any given number of dimensions, you can rearrange the expression
analogous to
overlap * p1 + (1-overlap) * p2
so that you won't get rounding errors when all neighboring pixels are
identical. But I don't know how to code it in a general way.
>
> 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?
>>
>>
>>Hi Stephen,
>>
>>It seems that it is better to fix the current
>>linear interpolator rather than adding a new
>>class.
>>
>>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.
>>
>>
>> Luis
>>
>>
>>-------------------------
>>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 = 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 = 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 == 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
>>>>>
>>>>>
>>>>
>>>
>>>
>>
>>
>>_______________________________________________
>>Insight-developers mailing list
>>Insight-developers at itk . org
>>http://www . itk . org/mailman/listinfo/insight-developers
>
>
--
Mark Foskey (919) 843-5436 Computer-Aided Diagnosis and Display Lab
mark_foskey at unc . edu Department of Radiology, CB 7515, UNC
http://www . cs . unc . edu/~foskey Chapel Hill, NC 27599-7515