[Insight-users] extract filter problem

[michakuhn at gmx.ch]

Hi,

I've written the following program, which is supposed to create the gradient
image for a particular region of the input image.

----

#include <iostream>

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkGradientImageFilter.h"
#include "itkGradientRecursiveGaussianImageFilter.h"
#include "itkNumericTraits.h"
#include "itkImageRegionConstIterator.h"
#include "itkDiscreteGaussianImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkExtractImageFilter.h"

int main(int argc, char** argv)
{
	const unsigned int Dimension = 3;
	typedef short PixelType;
	typedef itk::Image<PixelType, Dimension> ImageType;

	typedef double RealType;

	typedef itk::Image<RealType, Dimension> IntermediateImageType;

	typedef itk::CovariantVector<RealType, Dimension> GradientPixelType;
	typedef itk::Image<GradientPixelType, Dimension> GradientImageType;
	
	
	typedef itk::ImageFileReader<ImageType> ReaderType;

	ReaderType::Pointer reader = ReaderType::New();

	if (argc < 2) {
		std::cout << "image in required" << std::endl;
		return -1;
	}
	reader->SetFileName(argv[1]);

	try {
		reader->Update();
	} catch (itk::ExceptionObject &e) {
		std::cout << e << std::endl;
		return -1;
	}

	typedef itk::CastImageFilter<ImageType, IntermediateImageType>
		CastFilterType;

	CastFilterType::Pointer castFilter = CastFilterType::New();

	castFilter->SetInput(reader->GetOutput());


	typedef itk::DiscreteGaussianImageFilter<IntermediateImageType,
		IntermediateImageType>	PreFilterType;

	PreFilterType::Pointer preFilter = PreFilterType::New();
	preFilter->SetInput(castFilter->GetOutput());
	preFilter->SetUseImageSpacing(false);
	preFilter->SetMaximumKernelWidth(8);
	std::cout << "default variance: " << preFilter->GetVariance() << std::endl;
	preFilter->SetVariance(1);

	typedef itk::GradientImageFilter<IntermediateImageType, RealType, 
		RealType> Filter1Type;

	Filter1Type::Pointer filter1 = Filter1Type::New();
	filter1->SetInput(preFilter->GetOutput());

	typedef itk::ExtractImageFilter<GradientImageType, GradientImageType>
		ExtractFilterType;
	ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();

	extractFilter->SetInput(filter1->GetOutput());

	typedef GradientImageType::RegionType RegionType;
	RegionType region;
	RegionType::SizeType size;
	RegionType::IndexType index;

	size[0] = 1;
	size[1] = 32;
	size[2] = 32;

	index[0] = 64;
	index[1] = 0;
	index[2] = 0;

	region.SetSize(size);
	region.SetIndex(index);

	extractFilter->SetExtractionRegion(region);

	try {
		extractFilter->Update();
	} catch (itk::ExceptionObject &e) {
		std::cout << e << std::endl;
		return -1;
	}

	GradientImageType::Pointer image1 = extractFilter->GetOutput();

	std::cout << "image1: " << image1 << std::endl;

	return 0;
}

----

This program crashes, if my input image has a size of 65x65x127 pixels. It
works fine, however, as soon as I change (in the source code) "index[0] =
64" to "index[0] = 63". It still works when increasing size[0] to 2 after
changing the index. Can somebody explain me, why I am not allowed to extract
this particular region I've set in the code? Or is there something else
wrong?

Thanks,

Michael

P.S.: Below the input image file (in.mhd):

----

ObjectType = Image
NDims = 3
BinaryData = True
BinaryDataByteOrderMSB = False
ElementSpacing = 0.127664 0.127664 0.127664
DimSize = 65 65 127
ElementType = MET_SHORT
ElementDataFile = in.raw

----

-- 
NEU: Bis zu 10 GB Speicher für e-mails & Dateien!
1 GB bereits bei GMX FreeMail http://www.gmx.net/de/go/mail