[Insight-users] Problem with Rigid3DTransform

[Serena Fabbri]

Hi All,

I am trying to register T1-MRI and T2-MRI images with

itkVersorRigid3DTransform
itkVersorRigid3DTransformOptimizer
itkMattesMutualInformationImageToImageMetric

fixed: T2
moving: T1

but i find that T1-reg has the last 4 slides corrupted.
the n. 160 is empty and the other have half full statistics.
I don't think that this task is difficult for MI because the 2 images look very similar and I did a test: i registered T1 with T1 and i found T1-reg corrupted.
I think I made a mistake but i don't recognize where it is.
maybe it is not possible to use  MattesMutualInformationImageToImageMetric with VersorRigid3DTransform, is it?



T1
dim 176 256 160
spacing 1 1 1
origin 0 0 0

T2
dim 392 512 160
spacing 0.5 0.5 1 
origin 0 0 0

they are mha format

this is my code.....any idea?

Thanks.

Serena.




const    unsigned int    Dimension = 3;
typedef  float           PixelType;

typedef itk::Image< PixelType, Dimension >  FixedImageType;
typedef itk::Image< PixelType, Dimension >  MovingImageType;



typedef itk::VersorRigid3DTransform< double > TransformType;




typedef itk::VersorRigid3DTransformOptimizer           OptimizerType;

typedef itk::MattesMutualInformationImageToImageMetric< 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();
RegistrationType::Pointer   registration  = RegistrationType::New();



registration->SetOptimizer(     optimizer     );
registration->SetInterpolator(  interpolator  );




registration->SetMetric(        metric        );
metric->SetNumberOfHistogramBins( 32 );
metric->SetNumberOfSpatialSamples( 20000 );

 	if( argc > 9 )
     {

 		metric->SetUseExplicitPDFDerivatives( atoi( argv[10] ) );
     }

 	if( argc > 10 )
     {

 		metric->SetUseCachingOfBSplineWeights( atoi( argv[11] ) );
     }



TransformType::Pointer  transform = TransformType::New();
registration->SetTransform( transform );



typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;

FixedImageReaderType::Pointer  fixedImageReader  = FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();

fixedImageReader->SetFileName(  argv[1] );  fixedImageReader->Update();
movingImageReader->SetFileName( argv[2] );  movingImageReader->Update();





typedef itk::RescaleIntensityImageFilter< FixedImageType, FixedImageType >   RescalerType;

RescalerType::Pointer intensityRescaler1 = RescalerType::New();

intensityRescaler1->SetInput( fixedImageReader->GetOutput() );
intensityRescaler1->SetOutputMinimum(   0 );
intensityRescaler1->SetOutputMaximum( 255 );
intensityRescaler1->Update();

RescalerType::Pointer intensityRescaler2 = RescalerType::New();

intensityRescaler2->SetInput( movingImageReader->GetOutput() );
intensityRescaler2->SetOutputMinimum(   0 );
intensityRescaler2->SetOutputMaximum( 255 );
intensityRescaler2->Update();

registration->SetFixedImage(    intensityRescaler1->GetOutput()    );
registration->SetMovingImage(   intensityRescaler2->GetOutput()   );






registration->SetFixedImageRegion( intensityRescaler1->GetOutput()->GetBufferedRegion() );



typedef itk::CenteredTransformInitializer< TransformType,
                                              FixedImageType,
                                              MovingImageType
                                                  >  TransformInitializerType;

TransformInitializerType::Pointer initializer =
                                           TransformInitializerType::New();

initializer->SetTransform(   transform );
initializer->SetFixedImage(  intensityRescaler1->GetOutput() );
initializer->SetMovingImage(  intensityRescaler2->GetOutput() );
initializer->GeometryOn();

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 / 2000.0;

   optimizerScales[0] = 1.0;//original
   optimizerScales[1] = 1.0;
   optimizerScales[2] = 1.0;
   optimizerScales[3] = translationScale;
   optimizerScales[4] = translationScale;
   optimizerScales[5] = translationScale;

   optimizer->SetScales( optimizerScales );

   optimizer->SetMaximumStepLength( 0.1000  );//original 0.2000
   optimizer->SetMinimumStepLength( 0.0001 );

   optimizer->SetNumberOfIterations( 500 );



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


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

   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;





   typedef itk::ResampleImageFilter<
                             MovingImageType,
                             FixedImageType >    ResampleFilterType;

   TransformType::Pointer finalTransform = TransformType::New();

   finalTransform->SetCenter( transform->GetCenter() );

   finalTransform->SetParameters( finalParameters );

   ResampleFilterType::Pointer resampler = ResampleFilterType::New();

   resampler->SetTransform( finalTransform );

   resampler->SetInput( intensityRescaler2->GetOutput() );


   FixedImageType::Pointer fixedImage = intensityRescaler1->GetOutput();

   resampler->SetSize(    fixedImage->GetLargestPossibleRegion().GetSize() );
   resampler->SetOutputOrigin(  fixedImage->GetOrigin() );
   resampler->SetOutputSpacing( fixedImage->GetSpacing() );
   resampler->SetOutputDirection( fixedImage->GetDirection() );
   resampler->SetDefaultPixelValue( 100 );

   typedef  unsigned char  OutputPixelType;

   typedef itk::Image< OutputPixelType, Dimension > OutputImageType;

   typedef itk::CastImageFilter<
                         FixedImageType,
                         OutputImageType > CastFilterType;

   typedef itk::ImageFileWriter< OutputImageType >  WriterType;


   WriterType::Pointer      writer =  WriterType::New();
   CastFilterType::Pointer  caster =  CastFilterType::New();


  writer->SetFileName( argv[3] );



   caster->SetInput( resampler->GetOutput() );
   writer->SetInput( caster->GetOutput()   );
   writer->Update();