[Insight-users] bspline interpolation

[michakuhn at gmx.ch]

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


-- 
NEU: WLAN-Router für 0,- EUR* - auch für DSL-Wechsler!
GMX DSL = supergünstig & kabellos http://www.gmx.net/de/go/dsl