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

[Luis Ibanez]

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