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

[Stephen R. Aylward]

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
>>
>>



-- 
===========================================================
Dr. Stephen R. Aylward
Associate Professor of Radiology
Adjunct Associate Professor of Computer Science and Surgery
http://caddlab . rad . unc . edu
aylward at unc . edu
(919) 966-9695