[Insight-users] [ITK Community] 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


More information about the Insight-users mailing list