[ITK Community] [Insight-users] 3D registration
Matt McCormick
matt.mccormick at kitware.com
Sun Feb 23 14:50:36 EST 2014
Hi Massi,
Nothing pops out. It may help to plot the metric value and transform
parameters throughout the optimization process.
HTH,
Matt
On Sat, Feb 22, 2014 at 10:48 PM, Massinissa Bandou
<Massinissa.Bandou at usherbrooke.ca> wrote:
> 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
> _______________________________________________
> Community mailing list
> Community at itk.org
> http://public.kitware.com/cgi-bin/mailman/listinfo/community
_____________________________________
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
More information about the Community
mailing list