[Insight-users] Strange Behaviour of LinearInterpolateImageFunction

[Matthias Kirschner]

Hi,

I had some access violations when using LinearInterpolateImageFunction. 
I first thought that they occured because I evaluated the function "out 
of range" of the Image. Later I found out that the function usually 
returns values even if the index is out of range (I suppose this is 
desired behaviour). My second idea was that it might be caused by very 
small values and numerical instability. After trying to reproduce the 
error with small test code, the behaviour seems to me more and more obscure.

The attached code crashes with an access violation, unless I execute it 
in debug mode! Clearly, (214, 214) is out of image range, but the code 
works well if I change the coordinates to (213, 213), which is out of 
range, too. With other test code, my code also crashes in debug mode. If 
it does, the access violation always happens in 
LinearInterpolateImageFunction::EvaluateAtContinuousIndex(...) at the 
following position

   // get neighbor value only if overlap is not zero
    if( overlap )
      {
      value += overlap * static_cast<RealType>( 
this->GetInputImage()->GetPixel( neighIndex ) );
      totalOverlap += overlap;
      }

The reason for the access violation is that neighIndex is out of range, 
and I suppose this is a bug and should be avoided in the code of 
LinearInterpolateImageFunction.

Though my test code to reproduce the error is quite artificial, the 
problem occured in a meaningful program in the first place.  In case you 
can't reproduce the error: I'm using Visual Studio 2005.

Regards,
Matthias

Here is the code:

typedef itk::Image< unsigned char, 2> ImageType;
ImageType::Pointer image = ImageType::New();
ImageType::IndexType startindex;
startindex[0] = 0;
startindex[1] = 0;
ImageType::SizeType size;
size[0] = 100;
size[1] = 100;
ImageType::RegionType region;
region.SetSize( size );
region.SetIndex( startindex );
image->SetRegions( region );
image->Allocate();
image->FillBuffer(255);

typedef itk::LinearInterpolateImageFunction<ImageType, double> 
InterpolatorType;
typedef InterpolatorType::ContinuousIndexType                 IndexType;
       
InterpolatorType::Pointer interp = InterpolatorType::New();
interp->SetInputImage(image);

IndexType index;
index[0] = 214.0;
index[1] = 214.0;
interp->EvaluateAtContinuousIndex(index);