[Insight-users] bspline interpolation

[Luis Ibanez]

Hi Michael,

With BSplines in general you don't recover the same grid value when
you interpolate in a grid point position. That will only happen if
you using zero or first order BSplines, that are equivalent to
nearest-neighbor and linear interpolation respectively.

You may want to look at any introductory book on BSplines.
For example: Farin "Curves and Surfaces"

If you look at the diagrams of the recursive contruction of BSplines
in the Registration Techniques Course at RPI you will see why the
value interpolated at a grid point will not match the grid value.

     http://www.cs.rpi.edu/courses/spring04/imagereg/

Please look at lecture #26
http://www.cs.rpi.edu/courses/spring04/imagereg/lectureBSplines.ppt

Slides 10 to 24. Note that the slides have animations, so you
may want to see them in presentation/show mode.


   Regards,


     Luis


-------------------------
michakuhn at gmx.ch wrote:

> Hi,
> 
> I assumed, that the interpolated value at a grid position is the same as the
> pixel value at this position. However, I noticed that this is not always the
> case when bspline interpolation is used. I wonder whether this is the
> correct mathematical behaviour. I've written a little test program (see
> below) that counts the equal and none-equal values. The output for the
> BrainProtonDensitySliceBorder20.png looks as follows:
> 
> checking image...
> equal: 31811
> plus 1: 0
> minus 1: 24986
> 
>>+1: 0
> 
> < -1: 0
> 
> Can somebody explain me why there are grid positions where the interpolated
> values are not equal to the pixel values?
> 
> Thanks, Michael
> 
> P.S.: The program code that produced above results:
> 
> #include <iostream>
> 
> #include "itkImage.h"
> #include "itkImageFileReader.h"
> #include "itkBSplineInterpolateImageFunction.h"
> 
> 
> int main(int argc, char** argv)
> {
> 	const unsigned int Dimension = 2;
> 	typedef unsigned char PixelType;
> 	typedef itk::Image<PixelType, Dimension> ImageType;
> 
> 	typedef itk::ImageFileReader<ImageType> ReaderType;
> 
> 	ReaderType::Pointer reader = ReaderType::New();
> 
> 	if (argc < 2) {
> 		std::cout << "filename required." << std::endl;
> 	}
> 	reader->SetFileName(argv[1]);
> 
> 	try {
> 		reader->Update();
> 	} catch (itk::ExceptionObject &e) {
> 		std::cout << e << std::endl;
> 		return -1;
> 	}
> 
> 	ImageType::Pointer image = reader->GetOutput();
> 
> 	typedef itk::BSplineInterpolateImageFunction<ImageType, double, double>
> 		InterpolatorType;
> 
> 	InterpolatorType::Pointer interpolator = InterpolatorType::New();
> 	interpolator->SetSplineOrder(3);
> 
> 	interpolator->SetInputImage(image);
> 
> 	ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();
> 
> 
> 	std::cout << "checking image..." << std::endl;
> 
> 	int diffCount[5];
> 	for (int i = 0; i < 5; i++) {
> 		diffCount[i] = 0;
> 	}
> 
> 	ImageType::IndexType index;
> 
> 	for (int i = 0; i < size[0]; i++) {
> 		index[0] = i;
> 		for (int j = 0; j < size[1]; j++) {
> 			index[1] = j;
> 			int interpolated = interpolator->EvaluateAtContinuousIndex(index);
> 			int value = image->GetPixel(index);
> 			int diff = interpolated - value;
> 			if (diff > 2) {
> 				diff = 2;
> 			}
> 			if (diff < -2) {
> 				diff = -2;
> 			}
> 			diffCount[diff + 2]++;
> 		}
> 	}
> 
> 	std::cout << "equal: " << diffCount[2] << std::endl;
> 	std::cout << "plus 1: " << diffCount[3] << std::endl;
> 	std::cout << "minus 1: " << diffCount[1] << std::endl;
> 	std::cout << "> +1: " << diffCount[4] << std::endl;
> 	std::cout << "< -1: " << diffCount[0] << std::endl;
> 
> 	return 0;
> }
> 
>