[Insight-users] ImageFilter and ImageFunction for Gaussian derivatives

[Matthias Schneider]

Hi.

I'm a little confused about the functional difference of 
itkDiscreteGaussianDerivativeImageFilter and 
itkDiscreteGaussianDerivativeImageFunction.

Looking at the ITK manual, I'd assume that the image function produces 
the same output as the image filter when evaluated at the corresponding 
(non-continuous) image indices (no interpolation).

Comparing the results of the two filters (with same settings) for 
3D-image, I obtain a weird sign shift, i.e., if the sum of the 
derivative orders is odd, the resulting filter values have an opposite 
sign. Do I miss something here?

Any help is appreciated very much.

Thanks,
Matthias



This is the code that I used for the comparison:

#include <stdlib.h>

#include <itkImageFileWriter.h>
#include <itkImageFileReader.h>
#include <itkDiscreteGaussianDerivativeImageFunction.h>
#include <itkDiscreteGaussianDerivativeImageFilter.h>

int main(int /*argc*/, char **argv) {
	const std::string fname = "/path/to/my/volume.nii";
	int x[3] = { 21, 16, 15 };	// image index to evaluate
	int o[3] = { 1, 0, 0 };		// derivative orders [dx dy dz]

	typedef itk::Image<unsigned char, 3> ImageType;
	typedef itk::Image<double, 3> DImageType;
	typedef itk::ImageFileReader<ImageType> ReaderType;

	itk::SmartPointer<ReaderType> reader = ReaderType::New();
	reader->SetFileName(fname.c_str());
	reader->Update();

	typedef itk::DiscreteGaussianDerivativeImageFilter<ImageType, 
DImageType> FilterType;
	typedef itk::DiscreteGaussianDerivativeImageFunction<ImageType, 
DImageType::PixelType> FunctorType;
	itk::SmartPointer<FilterType> filter = FilterType::New();
	itk::SmartPointer<FunctorType> func = FunctorType::New();

	FilterType::ArrayType variance, order, maxError;
	ImageType::RegionType roi;
	ImageType::IndexType idx;

	variance[0] = variance[1] = variance[2] = 20;
	maxError[0] = maxError[1] = maxError[2] = 1e-4;
	for (unsigned int i = 0; i < 3; ++i) {
		roi.SetSize(i, 1);
		roi.SetIndex(i, x[i]);
		idx.SetElement(i, x[i]);
		order[i] = o[i];
	}

	filter->SetInput(reader->GetOutput());
	filter->UseImageSpacingOff();
	filter->SetVariance(variance);
	filter->SetMaximumKernelWidth(48);
	filter->NormalizeAcrossScaleOn();
	filter->SetMaximumError(maxError);
	filter->SetOrder(order);
	filter->GetOutput()->SetRequestedRegion(roi);
	filter->Update();

	func->SetInputImage(reader->GetOutput());
	func->UseImageSpacingOff();
	func->SetVariance(variance);
	func->SetMaximumKernelWidth(48);
	func->NormalizeAcrossScaleOn();
	func->SetMaximumError(maxError[0]);
	func->SetOrder(order);
	func->Initialize();

     std::cout << "Filter: " << filter->GetOutput()->GetPixel(idx) << 
std::endl;
	std::cout << "Functor: " << func->EvaluateAtIndex(idx) << std::endl;

	return EXIT_SUCCESS;
}

The output looks like:
Filter: -2.12266
Functor: 2.12266