<div>Hello,<br> <br>PROBLEM:<br> <br>I am trying to perform registration of two labelled volumes using ITK. I intend to get it to work with 3D medical images.<br> <br>STUDY:<br> <br>I went through the Image Registration Section of the ITK software guide and I came to know that the appropriate metric for my problem is the <br>
itk::MatchCardinalityImageToImageMetric. I also went through the corresponding example in itk --- "Examples\Registration\ImageRegistration10.cxx". <br> <br>I found that the MatchCardinalityImageToImageMetric does not provide analytical derivatives and simply counts the number of matches/mismatches. Hence the regular optimizers wouldnt work, and hence i tried to use the AmoebaOptimizer as shown in the example.<br>
<br>However i want to perform registration in the following order: Similarity Transform, Affine, Deformable<br> <br>To start with, i am just trying it with the Similarity3DTransform. <br> <br>itk::Similarity3DTransform derives from the itk::VersorRigid3DTransform and i read in the ITK software guide that the most appropriate optimizer for this transform is itk::VersorTransformOptimizer as it computes the versor derivatives as defined by hamilton. <br>
<br>QUESTION:<br> <br>Now Im confused as to which optimizer to use for my problem <br> <br>1. AmoebaOptimizer ( that is appropriate for MatchCardinalityImageToImageMetric)<br>2. VersorTransformOptimizer (that is appropriate for Similarity3DTransform)<br>
<br>Also, does the VersorTransformOptimizer require analytical derivates of the metric (which the MatchCardinalityImageToImageMetric does not provide).<br> <br>MY EXPERIMENTS:<br> <br>I tried to use the AmoebaOptimizer with the MatchCardinalityImageToImageMetric and Similarity3DTransform. <br>
<br>I was not clear how to assign the Initial Simplex for my transform parameters particularly those of the quaternion. So i let the optimizer determine it automatically.<br> <br>However, the registration method fails with an exception. </div>
<div> </div>
<div>My Input test data (i created "simple" toy images to test my code):</div>
<div> </div>
<div>> Fixed Image: Stack of two 2D slices with a rectangle positioned at the center of the image</div>
<div>> Moving Image: Stack of two 2D slices with an ellipse (of smaller size than a rectangle) position at the center of the image. The major axis of the ellipse is at 45 degrees about the image x-axis.</div>
<div> </div>
<div>I think it is failing when the optimizer is working with the scaling parameter ... I dont see a reason why it should do so .... </div>
<div> </div>
<div>Below is the command line output:<br> <br>****************************************************************************************<br> <br>Iter 1 : 0.1122 --- [0, 0, 0, 0, 0, 0, 1] <br>Iter 2 : 0.117374 --- [0.00025, 0, 0, 0, 0, 0, 1]<br>
Iter 3 : 0.114343 --- [0, 0.00025, 0, 0, 0, 0, 1]<br>Iter 4 : 0.114466 --- [0, 0, 0.00025, 0, 0, 0, 1]<br>Iter 5 : 0.113333 --- [0, 0, 0, 0.00025, 0, 0, 1]<br>Iter 6 : 0.113333 --- [0, 0, 0, 0, 0.00025, 0, 1]<br>Iter 7 : 0.1122 --- [0, 0, 0, 0, 0, 0.00025, 1]<br>
ERROR-REGISTRATION : <br>itk::ExceptionObject (013FFA90)<br>Location: "double __thiscall itk::MatchCardinalityImageToImageMetric<class itk::Image<unsigned char,3>,class itk::Image<unsigned char,3> >::GetNonconstValue(const class itk::Array<double> &)"<br>
File: c:\itk\insighttoolkit-3.10.2\code\algorithms\itkMatchCardinalityImageToImageMetric.txx<br>Line: 130<br>Description: itk::ERROR: MatchCardinalityImageToImageMetric(018DCEE0): All the points mapped to outside of the moving image<br>
<br>***************************************************************************************************<br> <br>At the end of this message is a section of the code that i wrote ... <br> <br>Any help would be greatly appreciated.<br>
<br>Thanks in advance.<br> <br>Regards,<br> <br>Deepak<br> <br> <br>****************************** CODE ******************************<br> <br>// Perform registration with Similarity3DTransform first <br>typedef itk::ImageRegistrationMethod<FixedImageType, MovingImageType> RegistrationType;<br>
RegistrationType::Pointer pRegistrationMethod = RegistrationType::New(); <br> // Metric<br> typedef itk::MatchCardinalityImageToImageMetric<FixedImageType,MovingImageType> MetricType;<br> MetricType::Pointer pMetric = MetricType::New();<br>
pMetric->MeasureMatchesOff();<br> pRegistrationMethod->SetMetric( pMetric );<br> <br> // Interpolator <br> typedef itk::NearestNeighborInterpolateImageFunction< MovingImageType, double> InterpolatorType;<br> InterpolatorType::Pointer pInterpolator = InterpolatorType::New();<br>
pRegistrationMethod->SetInterpolator( pInterpolator );<br> <br> // transform <br> typedef itk::Similarity3DTransform<double> SimilarityTransformType;<br> typedef itk::CenteredTransformInitializer<SimilarityTransformType,FixedImageType,MovingImageType> TransformInitializerType;<br>
<br> SimilarityTransformType::Pointer pTransform = SimilarityTransformType::New(); <br> TransformInitializerType::Pointer pTransformInitializer = TransformInitializerType::New();<br> <br> pTransformInitializer->SetFixedImage( pFixedImage ); <br>
pTransformInitializer->SetMovingImage( pMovingImage );<br> pTransformInitializer->SetTransform( pTransform );<br> pTransformInitializer->MomentsOn();<br> pTransformInitializer->InitializeTransform(); <br> <br>
// rotation<br> typedef SimilarityTransformType::VersorType VersorType;<br> typedef VersorType::VectorType VectorType;<br> VectorType rot_axis;<br> rot_axis[0] = 0.0; <br> rot_axis[1] = 0.0;<br> rot_axis[2] = 1.0;<br>
pTransform->SetRotation( rot_axis , 0.0 );<br> <br> // scaling <br> pTransform->SetScale( 1.0 );<br> pRegistrationMethod->SetInitialTransformParameters( pTransform->GetParameters() ); <br> pRegistrationMethod->SetTransform( pTransform ); <br>
// Optimizer <br> typedef itk::AmoebaOptimizer OptimizerType;<br> OptimizerType::Pointer pOptimizer = OptimizerType::New();<br> OptimizerType::ParametersType simplexDelta( pTransform->GetNumberOfParameters() );<br>// simplexDelta.Fill( 5.0 );<br>
// pOptimizer->AutomaticInitialSimplexOff();<br>// pOptimizer->SetInitialSimplexDelta( simplexDelta );<br> pOptimizer->SetParametersConvergenceTolerance( 0.25 ); // quarter pixel <br> pOptimizer->SetFunctionConvergenceTolerance( 0.001 ); // 0.1%<br>
pOptimizer->SetMaximumNumberOfIterations( 10 );<br> pRegistrationMethod->SetOptimizer( pOptimizer );<br> <br> // Observer <br> typedef CommandIterationUpdate ObserverType;<br> ObserverType::Pointer pObserver = ObserverType::New();<br>
pOptimizer->AddObserver( itk::IterationEvent() , pObserver);<br> <br>pRegistrationMethod->SetFixedImage( pFixedImage ); <br>pRegistrationMethod->SetMovingImage( pMovingImage );<br>pRegistrationMethod->SetFixedImageRegion( pFixedImage->GetLargestPossibleRegion() );<br>
<br>try <br>{<br> pRegistrationMethod->Initialize();<br> pRegistrationMethod->StartRegistration();<br>}<br>catch(itk::ExceptionObject & err)<br>{<br> std::cout << "\nERROR-REGISTRATION : " << err << std::endl;<br>
return EXIT_FAILURE; <br>}<br> <br>RegistrationType::ParametersType finalParameters;<br>finalParameters = pRegistrationMethod->GetLastTransformParameters();<br> <br>std::cout << "Final: " << pMetric->GetValue( finalParameters ) << " --- " << finalParameters << std::endl;<br>
<br>// Resampling -- Transform the moving image using the obtained transform parameters <br>typedef itk::ResampleImageFilter<MovingImageType,FixedImageType> ResampleFilterType;<br>ResampleFilterType::Pointer pResampleImageFilter = ResampleFilterType::New();<br>
// transform<br> SimilarityTransformType::Pointer pFinalTransform = SimilarityTransformType::New();<br> pFinalTransform->SetParameters( pRegistrationMethod->GetLastTransformParameters() );<br>pResampleImageFilter->SetTransform( pFinalTransform );<br>
pResampleImageFilter->SetInput( pMovingImage ); <br>pResampleImageFilter->SetSize( pFixedImage->GetLargestPossibleRegion().GetSize() );<br>pResampleImageFilter->SetOutputOrigin( pFixedImage->GetOrigin() );<br>
pResampleImageFilter->SetOutputSpacing( pFixedImage->GetSpacing() );<br>pResampleImageFilter->SetOutputDirection( pFixedImage->GetDirection() );<br>pResampleImageFilter->SetDefaultPixelValue( 0 );<br>pResampleImageFilter->SetInterpolator( pInterpolator );<br>
pResampleImageFilter->Update();<br> <br>************************************************************************************************<br clear="all"><br></div>