[ITK Community] [Insight-users] 3D registration

[Massinissa Bandou]

Hi everyone,
can somebody tell me what's wrong with this code because I cant match the
source & moving image with the final matrix of transformation??

I would appreciate your help!

Massi

	template<typename pixelType1, typename pixelType2>
	vtkSmartPointer<vtkMatrix4x4> MutualInformation(const string source, const
string target,int iteration,pixelType1,pixelType2)
	{
		typedef itk::Image< pixelType1, 3 >  FixedImageType; 
		typedef itk::Image< pixelType2, 3 >  MovingImageType;

		typedef itk::VersorRigid3DTransform< double > TransformType;
		typedef itk::VersorRigid3DTransformOptimizer                                 
OptimizerType;
		typedef itk::MutualInformationImageToImageMetric< FixedImageType,
MovingImageType > MetricType;
		typedef itk::LinearInterpolateImageFunction< MovingImageType, double >       
InterpolatorType;
		typedef itk::ImageRegistrationMethod< FixedImageType, MovingImageType >      
RegistrationType;
	
		MetricType::Pointer         metric        = MetricType::New();
		OptimizerType::Pointer      optimizer     = OptimizerType::New();
		InterpolatorType::Pointer   interpolator  = InterpolatorType::New();
		TransformType::Pointer  transform = TransformType::New();

		RegistrationType::Pointer   registration  = RegistrationType::New();
		registration->SetMetric(        metric        );
		registration->SetOptimizer(     optimizer     );
		registration->SetInterpolator(  interpolator  );
		registration->SetTransform( transform );

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

		typedef itk::ImageFileReader<FixedImageType> FixedImageReaderType;
		typedef itk::ImageFileReader<MovingImageType> MovingImageReaderType;

		FixedImageReaderType::Pointer  fixedImageReader  =
FixedImageReaderType::New();
		MovingImageReaderType::Pointer movingImageReader =
MovingImageReaderType::New();
		fixedImageReader->SetFileName(source.c_str()); 
		movingImageReader->SetFileName(target.c_str());


		registration->SetFixedImage(fixedImageReader->GetOutput());
		registration->SetMovingImage(   movingImageReader->GetOutput()   );

		try{
			fixedImageReader->Update();
		}catch(itk::ExceptionObject & excp){
			std::cerr << excp << std::endl;
		}

	
registration->SetFixedImageRegion(fixedImageReader->GetOutput()->GetBufferedRegion()
);

		typedef
itk::CenteredTransformInitializer<TransformType,FixedImageType,MovingImageType> 
TransformInitializerType;
		TransformInitializerType::Pointer initializer =
TransformInitializerType::New();
	  
		initializer->SetTransform(transform);
		initializer->SetFixedImage(fixedImageReader->GetOutput());
		initializer->SetMovingImage(movingImageReader->GetOutput());
		initializer->MomentsOn();
		initializer->InitializeTransform();

		typedef TransformType::VersorType  VersorType;
		typedef VersorType::VectorType     VectorType;
		VersorType     rotation;
		VectorType     axis;
		axis[0] = 0.0;
		axis[1] = 0.0;
		axis[2] = 1.0;
		const double angle = 0;
		rotation.Set(  axis, angle  );
		transform->SetRotation( rotation );

		registration->SetInitialTransformParameters( transform->GetParameters() );

		typedef OptimizerType::ScalesType       OptimizerScalesType;
		OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );
		const double translationScale = 1.0 / 1000.0;
		optimizerScales[0] = 1.0;
		optimizerScales[1] = 1.0;
		optimizerScales[2] = 1.0;
		optimizerScales[3] = translationScale;
		optimizerScales[4] = translationScale;
		optimizerScales[5] = translationScale;
		optimizer->SetScales( optimizerScales );
		optimizer->SetMaximumStepLength( 0.2000  );
		optimizer->SetMinimumStepLength( 0.00000001 );
		optimizer->SetNumberOfIterations( iteration );
		/*monitor registration*/
		CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
		observer->GetTomo(this->in);
		optimizer->AddObserver(itk::IterationEvent(),observer);

		try{
			registration->Update();
			std::cout << "Optimizer stop condition: "<<
registration->GetOptimizer()->GetStopConditionDescription()<< std::endl;
		}
		catch( itk::ExceptionObject & err ){
			std::cerr << "ExceptionObject caught !" << std::endl;
			std::cerr << err << std::endl;
		}

		OptimizerType::ParametersType finalParameters =
registration->GetLastTransformParameters();
		const double versorX              = finalParameters[0];
		const double versorY              = finalParameters[1];
		const double versorZ              = finalParameters[2];
		const double finalTranslationX    = finalParameters[3];
		const double finalTranslationY    = finalParameters[4];
		const double finalTranslationZ    = finalParameters[5];
		const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
		const double bestValue = optimizer->GetValue();

		std::cout << std::endl << std::endl;
		std::cout << " Result = " << std::endl;
		std::cout << " versor X      = " << versorX  << std::endl;
		std::cout << " versor Y      = " << versorY  << std::endl;
		std::cout << " versor Z      = " << versorZ  << std::endl;
		std::cout << " Translation X = " << finalTranslationX  << std::endl;
		std::cout << " Translation Y = " << finalTranslationY  << std::endl;
		std::cout << " Translation Z = " << finalTranslationZ  << std::endl;
		std::cout << " Iterations    = " << numberOfIterations << std::endl;
		std::cout << " Metric value  = " << bestValue          << std::endl;

		transform->SetParameters( finalParameters );
		TransformType::MatrixType matrix = transform->GetMatrix();
		TransformType::OffsetType offset = transform->GetOffset();
		std::cout<<"Matrix = "<<std::endl<<matrix<<std::endl;
		std::cout<<"Offset =
"<<std::endl<<offset<<std::endl;

		vtkSmartPointer<vtkMatrix4x4> matrixOfTransformation =
vtkSmartPointer<vtkMatrix4x4>::New();
		for(unsigned int i=0;i<3;i++){
			for(unsigned int j=0;j<3;j++){
				matrixOfTransformation->SetElement(i,j,matrix(i,j));
			}
		}

		matrixOfTransformation->SetElement(0,3,finalTranslationX );
		matrixOfTransformation->SetElement(1,3,finalTranslationY );
		matrixOfTransformation->SetElement(2,3,finalTranslationZ );
		matrixOfTransformation->SetElement(3,0,0);
		matrixOfTransformation->SetElement(3,1,0);
		matrixOfTransformation->SetElement(3,2,0);
		matrixOfTransformation->SetElement(3,3,1);


		return matrixOfTransformation;
}




--
View this message in context: http://itk-users.7.n7.nabble.com/3D-registration-tp33441.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