[Insight-users] Registration of US and CT images

[Anders Almås]

Hi!

I am performing registration of US images and CT images in 3D. So far I
am only using Mattes mutual information metric. Now my ragistration is
only performed on a perfect aligned set of US and CT images. What I do
is to initially rotate and translate the CT image and try to registrate
it such that the translated/rotated image is aligned again. First i
used the Centered Affine transform, but that gave only results which
was sometimes ok for only rotation. Then when I only tested for
translations the translation transform worked quite good. At present
time I am testing the VersorRigid3DTransform. This one works ok at
translations and it seems that it manages rotations in one direction.
My question for VersorRigid3DTransform is how it should be initialized?
I send the code so that any one can comment it... Is it the
optimizerscales which are the problem. Or is it really that I perform
an affine transformation/rotation that makes it difficult to registrate
it back? I hope anyone can help me about this...
The size of the images is 256*256*176 with spacing 1.453,1.453 and 2

Here follows the code for the registration:

 std::cout<<"Here starts the registration process..."<<std::endl;

 
  //TranslationTransform
 // typedef itk::TranslationTransform<double, ImageDimension>
TranslationType;
  //Rigid3DTransform
  typedef itk::VersorRigid3DTransform<double> TranslationType;
 // typedef itk::CenteredAffineTransform< double, ImageDimension  > 
TranslationType;
  //Transform filters
  //Optimizer filters
   typedef itk::VersorRigid3DTransformOptimizer     OptimizerType;
  
  //  typedef itk::RegularStepGradientDescentOptimizer      
OptimizerType;

  //Image Function filters
  typedef itk::LinearInterpolateImageFunction<ImageType, double  >
InterpolatorType;
  
  
  //Metric filter
  typedef itk::MattesMutualInformationImageToImageMetric< ImageType,
ImageType >    MetricType;
  //typedef itk::MutualInformationImageToImageMetric<ImageType,ImageType
> MetricType;
  
  //RegistrationMethod
  //typedef itk::ImageRegistrationMethod< ImageType, ImageType    >
RegistrationType;
     typedef itk::MultiResolutionImageRegistrationMethod< ImageType,
ImageType >   RegistrationType;
   
	 typedef itk::MultiResolutionPyramidImageFilter<ImageType,ImageType >  
FixedImagePyramidType;
  typedef itk::MultiResolutionPyramidImageFilter<ImageType,ImageType >  
MovingImagePyramidType;
	 
	 typedef OptimizerType::ScalesType       OptimizerScalesType;

  TranslationType::Pointer transform2 = TranslationType::New();
  OptimizerType::Pointer      optimizer     = OptimizerType::New();
  InterpolatorType::Pointer   interpolator2  = InterpolatorType::New();
  RegistrationType::Pointer   registration  = RegistrationType::New();
  MetricType::Pointer         metric        = MetricType::New();


  FixedImagePyramidType::Pointer fixedImagePyramid =  
FixedImagePyramidType::New();
  MovingImagePyramidType::Pointer movingImagePyramid = 
MovingImagePyramidType::New();

 typedef itk::CenteredTransformInitializer< TranslationType, ImageType,
ImageType >  TransformInitializerType;
  TransformInitializerType::Pointer initializer =
TransformInitializerType::New();
  initializer->SetTransform(   transform2 );
 initializer->SetFixedImage(   range->GetOutput());
 //  initializer->SetFixedImage(  gaussian->GetOutput());
  initializer->SetMovingImage( filter->GetOutput() );
  //initializer->SetMovingImage( ctimage );
 initializer->MomentsOn();
  initializer->InitializeTransform();

  typedef TranslationType::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  );

  transform2->SetRotation( rotation );
  axis[0] = 1.0;
  axis[1] = 0.0;
  axis[2] = 0.0;

  rotation.Set(  axis, angle  );

  transform2->SetRotation( rotation );

  axis[0] = 0.0;
  axis[1] = 1.0;
  axis[2] = 0.0;

 // const double angle = 0;

  rotation.Set(  axis, angle  );

  transform2->SetRotation( rotation );
  //transform2->SetIdentity();
  registration->SetTransform(transform2);
  registration->SetOptimizer(     optimizer     );
  registration->SetInterpolator(  interpolator2 );
  registration->SetMetric( metric  );
  registration->SetFixedImagePyramid( fixedImagePyramid );
  registration->SetMovingImagePyramid( movingImagePyramid );
 
  OptimizerScalesType optimizerScales(
transform2->GetNumberOfParameters() );

  std::cout<<"Number of parametres:
"<<transform2->GetNumberOfParameters()<<std::endl;


  optimizerScales[0] = 1.0; 
  optimizerScales[1] = 1.0; 
  optimizerScales[2] = 1.0; 
  optimizerScales[3] = 1.0/1000; 

  optimizerScales[4] = 1.0/1000; 
  optimizerScales[5] = 1.0/1000 ;
  
  optimizer->SetScales( optimizerScales );
  optimizer->MaximizeOn();
  
 registration->SetFixedImage(    range->GetOutput()); //US IMAGE EDGE
MAP
   registration->SetMovingImage(  filter->GetOutput()  ); // CTimage
gaussian blurred translated and rotated...
  registration->SetFixedImageRegion(
range->GetOutput()->GetBufferedRegion() );

  typedef RegistrationType::ParametersType ParametersType;
  
  /*
  ParametersType initialParameters( transform->GetNumberOfParameters()
);

  initialParameters[0] = 0.0;  // Initial offset in mm along X
  initialParameters[1] = 0.0;  // Initial offset in mm along Y
  initialParameters[2] = 0.0;  // Initial offset in mm along Z
*/
  registration->SetInitialTransformParameters(
transform2->GetParameters() );

  metric->SetNumberOfHistogramBins( par[7] ); // 20
  metric->SetNumberOfSpatialSamples( par[8] ); //10000

  optimizer->SetNumberOfIterations( par[6] ); // Tested from 300 - 800
iterations

  CommandIterationUpdate::Pointer observer =
CommandIterationUpdate::New();
  optimizer->AddObserver( itk::IterationEvent(), observer );


  
  typedef RegistrationInterfaceCommand<RegistrationType> CommandType;
  CommandType::Pointer command = CommandType::New();
  registration->AddObserver( itk::IterationEvent(), command );

  registration->SetNumberOfLevels( par[9] ); //Tested from 3-5 levels..

  try 
    { 
    registration->StartRegistration(); 
    } 
  catch( itk::ExceptionObject & err ) 
    { 
    std::cout << "ExceptionObject caught !" << std::endl; 
    std::cout << err << std::endl; 
    
	return input;
    } 

   ParametersType finalParameters =
                  registration->GetLastTransformParameters();

  unsigned int numberOfIterations = optimizer->GetCurrentIteration();
  

double bestValue = optimizer->GetValue();

 
  // Print out results	
  //
  std::cout << "Result = " << std::endl;
  std::cout << " Iterations    = " << numberOfIterations << std::endl;
  std::cout << " Metric value  = " << bestValue          << std::endl;
  //std::cout << " Stop Condition  = " << optimizer->GetStopCondition()
<< std::endl;

I really hope somebody can help me with some advice. Maybe is it that my
parameters is set wrong??

Regards 

Anders