[Insight-users] Fwd: Registration Speed using Normalized Mutual Information

[Alexia Rodríguez Ruano]

>
>Hi
>I'm modifying the example ImageRegistration8.cxx to reguister 3D 
>images using Normalized Mutual Information.
>The time it takes to calculate each iteration is too slow (1 minute 
>each one), and I am working with release version.

The images are identical but one of them is translated and 
rotated  and with a size of 181x271x181
The number of histogram bins is 32
>How can I solve that??

The NormalizedMutualInformationHistogramImageToImageMetric is 
implemented directly from the histogram or  follows any specific 
method (Viola, Mattes, etc)


>-----------------------------------------
>         ...
>         ....
>   const    unsigned int    Dimension = 3;
>   typedef  unsigned char   PixelType;
>
>   typedef itk::Image< PixelType, Dimension >  FixedImageType;
>   typedef itk::Image< PixelType, Dimension >  MovingImageType;
>   typedef itk::VersorRigid3DTransform< double > TransformType;
>   typedef itk::VersorRigid3DTransformOptimizer   OptimizerType;
>   typedef itk::NormalizedMutualInformationHistogramImageToImageMetric<
>                                           FixedImageType,
>                                           MovingImageType >    MetricType;
>   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->SetMetric(        metric        );
>   registration->SetOptimizer(     optimizer     );
>   registration->SetInterpolator(  interpolator  );
>
>   unsigned int numberOfHistogramBins =  *argv6;
>   MetricType::HistogramType::SizeType histogramSize;
>   histogramSize[0] = numberOfHistogramBins;
>   histogramSize[1] = numberOfHistogramBins;
>   metric->SetHistogramSize( histogramSize );
>
>   TransformType::Pointer  transform = TransformType::New();
>   registration->SetTransform( transform );
>
>         const unsigned int numberOfParameters = 
> transform->GetNumberOfParameters();
>         typedef MetricType::ScalesType ScalesType;
>         ScalesType scales( numberOfParameters );
>         const double transscale = 1.0 / 1000.0;
>         scales.Fill( 1.0 );
>         scales[0] = 1.0;
>         scales[1] = 1.0;
>         scales[2] = 1.0;
>         scales[3] = transscale;
>         scales[4] = transscale;
>         scales[5] = transscale;
>         metric->SetDerivativeStepLengthScales(scales);
>
>         typedef itk::ImportImageFilter< PixelType, 3  > 
> FixedImageImportFilterType;
>         typedef itk::ImportImageFilter< PixelType, 3  > 
> MovingImageImportFilterType;
>
>         FixedImageImportFilterType::Pointer  FixedImageImport  = 
> FixedImageImportFilterType::New();
>         MovingImageImportFilterType::Pointer  MovingImageImport  = 
> MovingImageImportFilterType::New();
>
>         FixedImageImportFilterType::SizeType  size;
>         size[0]  = argv4[0];  // size along X
>         size[1]  = argv4[1];  // size along Y
>         size[2] = argv4[2]; // size along Z
>         FixedImageImportFilterType::IndexType start;
>         start.Fill( 0 );
>         FixedImageImportFilterType::RegionType region;
>         region.SetIndex( start );
>         region.SetSize(  size  );
>         FixedImageImport->SetRegion( region );
>         const unsigned int numberOfPixels =  size[0]*size[1]*size[2];
>
>         MovingImageImportFilterType::SizeType  sizem;
>         sizem[0]  = argv5[0];  // size along X
>         sizem[1]  = argv5[1];  // size along Y
>         sizem[2]  = argv5[2]; // size along Z
>         MovingImageImportFilterType::IndexType startm;
>         startm.Fill( 0 );
>         MovingImageImportFilterType::RegionType regionm;
>         regionm.SetIndex( startm );
>         regionm.SetSize(  sizem  );
>         MovingImageImport->SetRegion( regionm );
>         const unsigned int numberOfPixelsm =  sizem[0]*sizem[1]*sizem[2];
>
>         const bool importImageFilterWillOwnTheBuffer = false;
>         FixedImageImport->SetImportPointer( argv1, 
> numberOfPixels,importImageFilterWillOwnTheBuffer );
>         MovingImageImport->SetImportPointer( argv2, 
> numberOfPixelsm,importImageFilterWillOwnTheBuffer );
>
>   registration->SetFixedImage(    FixedImageImport->GetOutput()    );
>   registration->SetMovingImage(   MovingImageImport->GetOutput()   );
>
>         FixedImageImport->Update();
>   registration->SetFixedImageRegion( 
> FixedImageImport->GetOutput()->GetBufferedRegion() );
>
>   typedef itk::CenteredTransformInitializer< TransformType,
>                                              FixedImageType,
>                                              MovingImageType
>                                                  >  TransformInitializerType;
>
>   TransformInitializerType::Pointer initializer =
>                                           TransformInitializerType::New();
>
>   initializer->SetTransform(   transform );
>   initializer->SetFixedImage(  FixedImageImport->GetOutput() );
>   initializer->SetMovingImage( MovingImageImport->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  );
>         std::cout<< "transform->SetRotation( rotation )" << std::endl;
>   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.2  );
>   optimizer->SetMinimumStepLength( 0.001 );
>   optimizer->SetRelaxationFactor(  0.90 );//Alexia
>   optimizer->SetNumberOfIterations( 100 );
>   optimizer->MaximizeOn();
>
>   CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
>
>   optimizer->AddObserver( itk::IterationEvent(), observer );
>
>   registration->ReleaseDataFlagOn();
>
>   try
>     {
>     registration->StartRegistration();
>     }
>   catch( itk::ExceptionObject & err )
>     {
>     std::cerr << "ExceptionObject caught !" << std::endl;
>     std::cerr << err << std::endl;
>     return -1;
>     }
>...
>...
>...