[ITK Community] [Insight-users] registration using Mutual Information

[Massinissa Bandou]

Hello everyone,

With Qt+vtk+itk, I tried to perform a mutual information registration
between 2 vtkImageData (source: dicom series, target .raw file) as you can
see in the picture. but I got this error:
*Qt has caught an exception thrown from an event handler. Throwing
exceptions from an event handler is not supported in Qt. You must
reimplement QApplication::notify() and catch all exceptions there.**

<http://itk-users.7.n7.nabble.com/file/n32972/test.png> 

and my code is:


	typedef itk::VTKImageToImageFilter<ImageType> VTKToITKConnector;
        VTKToITKConnector::Pointer fixed = VTKToITKConnector::New();
	VTKToITKConnector::Pointer moving = VTKToITKConnector::New();
	fixed->SetInput(source);  //source vtkImageData
	moving->SetInput(target);  //target vtkImageData

	ImageType::Pointer  fixedImage = fixed->GetOutput();
	ImageType::Pointer  movingImage = moving->GetOutput();

	typedef float InternalPixelType;

	typedef itk::Image<float,3> InternalImageType;
	
	typedef itk::NormalizeImageFilter<ImageType, InternalImageType>
NormalizeFilterType;
	
	NormalizeFilterType::Pointer fixedNormalizer = NormalizeFilterType::New();
	NormalizeFilterType::Pointer movingNormalizer = NormalizeFilterType::New();
	
	fixedNormalizer->SetInput(fixedImage);
	movingNormalizer->SetInput(movingImage);

	  // Smooth the normalized images
	typedef
itk::DiscreteGaussianImageFilter<InternalImageType,InternalImageType>
GaussianFilterType;
	GaussianFilterType::Pointer fixedSmoother  = GaussianFilterType::New();
	GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();
	
	fixedSmoother->SetVariance(2.0);
	movingSmoother->SetVariance(2.0);
	
	fixedSmoother->SetInput(fixedNormalizer->GetOutput());
	movingSmoother->SetInput(movingNormalizer->GetOutput());

	typedef itk::TranslationTransform<double,3> TransformType;
	typedef itk::GradientDescentOptimizer OptimizerType;
	typedef itk::LinearInterpolateImageFunction<InternalImageType,double>
InterpolatorType;
	typedef itk::ImageRegistrationMethod<InternalImageType,InternalImageType>
RegistrationType;
	typedef
itk::MutualInformationImageToImageMetric<InternalImageType,InternalImageType>
MetricType;
 
	TransformType::Pointer transform = TransformType::New();
	OptimizerType::Pointer optimizer = OptimizerType::New();
	InterpolatorType::Pointer interpolator = InterpolatorType::New();
	RegistrationType::Pointer registration = RegistrationType::New();
 
	registration->SetOptimizer(optimizer);
	registration->SetTransform(transform);
	registration->SetInterpolator(interpolator);
 
	MetricType::Pointer metric = MetricType::New();
	registration->SetMetric(metric);


	metric->SetFixedImageStandardDeviation(0.4);
	metric->SetMovingImageStandardDeviation(0.4);

	registration->SetFixedImage(fixedSmoother->GetOutput());
	registration->SetMovingImage(movingSmoother->GetOutput());
	fixedNormalizer->Update();
	 
	ImageType::RegionType fixedImageRegion =
fixedNormalizer->GetOutput()->GetBufferedRegion();
	registration->SetFixedImageRegion(fixedImageRegion);
 
	typedef RegistrationType::ParametersType ParametersType;
	ParametersType initialParameters(transform->GetNumberOfParameters());

	initialParameters[0] = 0.0;  // Initial offset along X
	initialParameters[1] = 0.0;  // Initial offset along Y
	initialParameters[2] = 0.0;  // Initial offset along Z

	registration->SetInitialTransformParameters(initialParameters);


	const unsigned int numberOfPixels = fixedImageRegion.GetNumberOfPixels();
	const unsigned int numberOfSamples = static_cast<unsigned
int>(numberOfPixels * 0.01);
	metric->SetNumberOfSpatialSamples(numberOfSamples);


	optimizer->SetLearningRate(15.0);
	optimizer->SetNumberOfIterations(200);
	optimizer->MaximizeOn();


		registration->Update();
		std::cout<<"Optimizer stop condition:
"<<registration->GetOptimizer()->GetStopConditionDescription()<<std::endl;


	ParametersType finalParameters =
registration->GetLastTransformParameters();
	double TranslationAlongX = finalParameters[0];
	double TranslationAlongY = finalParameters[1];
	double TranslationAlongZ = finalParameters[2];

	unsigned int numberOfIterations = optimizer->GetCurrentIteration();
	double bestValue = optimizer->GetValue();
	// Print out results
	std::cout<<std::endl;
	std::cout<<"Result = "<<std::endl;
	std::cout<<" Translation X =
"<<TranslationAlongX<<std::endl;
	std::cout<<" Translation Y =
"<<TranslationAlongY<<std::endl;
	std::cout<<" Translation Z =
"<<TranslationAlongZ<<std::endl;
	std::cout<<" Iterations    =
"<<numberOfIterations<<std::endl;
	std::cout<<" Metric value  =
"<<bestValue<<std::endl;
	std::cout<<" Numb. Samples =
"<<numberOfSamples<<std::endl;

	typedef itk::ResampleImageFilter<ImageType,ImageType>
ResampleFilterType;
	
	TransformType::Pointer finalTransform = TransformType::New();
	finalTransform->SetParameters(finalParameters);
	finalTransform->SetFixedParameters(transform->GetFixedParameters());
	
	ResampleFilterType::Pointer resample = ResampleFilterType::New();
	resample->SetTransform(finalTransform);
	resample->SetInput(movingImage);

	resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
	resample->SetOutputOrigin(fixedImage->GetOrigin());
	resample->SetOutputSpacing(fixedImage->GetSpacing());
	resample->SetOutputDirection(fixedImage->GetDirection());
	resample->SetDefaultPixelValue(100);

	// connect to VTK
	typedef itk::ImageToVTKImageFilter<ImageType> ConnectorType;
	ConnectorType::Pointer connector = ConnectorType::New();
	connector->SetInput(resample->GetOutput());

	return connector->GetOutput();


Any help please!

Massi



--
View this message in context: http://itk-users.7.n7.nabble.com/registration-using-Mutual-Information-tp32972.html
Sent from the ITK - Users mailing list archive at Nabble.com.
_____________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://www.itk.org/mailman/listinfo/insight-users