[Insight-users] Registration of Two Labeled Volumes !!!

[Deepak Roy Chittajallu]

Hello,

PROBLEM:

I am trying to perform registration of two labelled volumes using ITK.

STUDY:

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
itk::MatchCardinalityImageToImageMetric. I also went through the
corresponding example in itk ---
"Examples\Registration\ImageRegistration10.cxx".

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.

However i want to perform registration in the following order:  Similarity
Transform, Affine, Deformable

To start with, i am just trying it with the Similarity3DTransform.

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.

QUESTION:

Now Im confused as to which optimizer to use for my problem

1. AmoebaOptimizer ( that is appropriate for
MatchCardinalityImageToImageMetric)
2. VersorTransformOptimizer (that is appropriate for Similarity3DTransform)

Also, does the VersorTransformOptimizer require analytical derivates of the
metric (which the MatchCardinalityImageToImageMetric does not provide).

MY EXPERIMENTS:

I tried to use the AmoebaOptimizer with the
MatchCardinalityImageToImageMetric and Similarity3DTransform.

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.

However, the registration method fails with an exception. Below is the
command line output:

****************************************************************************************

Iter 1 : 0.1122 --- [0, 0, 0, 0, 0, 0, 1] Iter 2 : 0.117374 --- [0.00025, 0,
0, 0, 0, 0, 1]
Iter 3 : 0.114343 --- [0, 0.00025, 0, 0, 0, 0, 1]
Iter 4 : 0.114466 --- [0, 0, 0.00025, 0, 0, 0, 1]
Iter 5 : 0.113333 --- [0, 0, 0, 0.00025, 0, 0, 1]
Iter 6 : 0.113333 --- [0, 0, 0, 0, 0.00025, 0, 1]
Iter 7 : 0.1122 --- [0, 0, 0, 0, 0, 0.00025, 1]
ERROR-REGISTRATION : itk::ExceptionObject (013FFA90)
Location: "double __thiscall itk::MatchCardinalityImageToImageMetric<class
itk::Image<unsigned char,3>,class itk::Image<unsigned char,3>
>::GetNonconstValue(const class itk::Array<double> &)"
File:
c:\itk\insighttoolkit-3.10.2\code\algorithms\itkMatchCardinalityImageToImageMetric.txx
Line: 130
Description: itk::ERROR: MatchCardinalityImageToImageMetric(018DCEE0): All
the points mapped to outside of the moving image


***************************************************************************************************

At the end of this message is a section of the code that i wrote ...

Any help would be greatly appreciated.

Thanks in advance.

Regards,

Deepak


****************************** CODE ******************************

// Perform registration with Similarity3DTransform first typedef
itk::ImageRegistrationMethod<FixedImageType, MovingImageType>
RegistrationType;
RegistrationType::Pointer pRegistrationMethod =
RegistrationType::New();

 // Metric
 typedef
itk::MatchCardinalityImageToImageMetric<FixedImageType,MovingImageType>
MetricType;
 MetricType::Pointer pMetric = MetricType::New();
 pMetric->MeasureMatchesOff();
 pRegistrationMethod->SetMetric( pMetric );

 // Interpolator  typedef itk::NearestNeighborInterpolateImageFunction<
MovingImageType, double> InterpolatorType;
 InterpolatorType::Pointer pInterpolator = InterpolatorType::New();
 pRegistrationMethod->SetInterpolator( pInterpolator );

 // transform  typedef itk::Similarity3DTransform<double>
SimilarityTransformType;
 typedef
itk::CenteredTransformInitializer<SimilarityTransformType,FixedImageType,MovingImageType>
TransformInitializerType;

 SimilarityTransformType::Pointer pTransform =
SimilarityTransformType::New();   TransformInitializerType::Pointer
pTransformInitializer = TransformInitializerType::New();

 pTransformInitializer->SetFixedImage( pFixedImage );
pTransformInitializer->SetMovingImage(
pMovingImage );
 pTransformInitializer->SetTransform( pTransform );
 pTransformInitializer->MomentsOn();
 pTransformInitializer->InitializeTransform();


  // rotation
  typedef SimilarityTransformType::VersorType VersorType;
  typedef VersorType::VectorType VectorType;
  VectorType rot_axis;
  rot_axis[0] = 0.0;   rot_axis[1] = 0.0;
  rot_axis[2] = 1.0;
  pTransform->SetRotation( rot_axis , 0.0 );

 // scaling  pTransform->SetScale( 1.0 );
 pRegistrationMethod->SetInitialTransformParameters(
pTransform->GetParameters() );  pRegistrationMethod->SetTransform(
pTransform );

 // Optimizer  typedef itk::AmoebaOptimizer OptimizerType;
 OptimizerType::Pointer pOptimizer = OptimizerType::New();
 OptimizerType::ParametersType simplexDelta(
pTransform->GetNumberOfParameters() );

//   simplexDelta.Fill( 5.0 );
//   pOptimizer->AutomaticInitialSimplexOff();
//   pOptimizer->SetInitialSimplexDelta( simplexDelta );

 pOptimizer->SetParametersConvergenceTolerance( 0.25 ); // quarter
pixel  pOptimizer->SetFunctionConvergenceTolerance(
0.001 ); // 0.1%
 pOptimizer->SetMaximumNumberOfIterations( 10 );
 pRegistrationMethod->SetOptimizer( pOptimizer );

 // Observer  typedef CommandIterationUpdate ObserverType;
 ObserverType::Pointer pObserver = ObserverType::New();
 pOptimizer->AddObserver( itk::IterationEvent() , pObserver);


pRegistrationMethod->SetFixedImage( pFixedImage );
pRegistrationMethod->SetMovingImage(
pMovingImage );
pRegistrationMethod->SetFixedImageRegion(
pFixedImage->GetLargestPossibleRegion() );

try {
 pRegistrationMethod->Initialize();
 pRegistrationMethod->StartRegistration();
}
catch(itk::ExceptionObject & err)
{
 std::cout << "\nERROR-REGISTRATION : " << err << std::endl;
 return EXIT_FAILURE; }

RegistrationType::ParametersType finalParameters;
finalParameters = pRegistrationMethod->GetLastTransformParameters();

std::cout << "Final: " << pMetric->GetValue( finalParameters ) << " --- " <<
finalParameters << std::endl;

// Resampling -- Transform the moving image using the obtained transform
parameters typedef itk::ResampleImageFilter<MovingImageType,FixedImageType>
ResampleFilterType;
ResampleFilterType::Pointer pResampleImageFilter =
ResampleFilterType::New();

 // transform
 SimilarityTransformType::Pointer pFinalTransform =
SimilarityTransformType::New();
 pFinalTransform->SetParameters(
pRegistrationMethod->GetLastTransformParameters() );

pResampleImageFilter->SetTransform( pFinalTransform );
pResampleImageFilter->SetInput( pMovingImage ); pResampleImageFilter->SetSize(
pFixedImage->GetLargestPossibleRegion().GetSize() );
pResampleImageFilter->SetOutputOrigin(  pFixedImage->GetOrigin() );
pResampleImageFilter->SetOutputSpacing( pFixedImage->GetSpacing() );
pResampleImageFilter->SetOutputDirection( pFixedImage->GetDirection() );
pResampleImageFilter->SetDefaultPixelValue( 0 );
pResampleImageFilter->SetInterpolator( pInterpolator );
pResampleImageFilter->Update();

************************************************************************************************
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090312/52f2f5fc/attachment.htm>