[Insight-users] problem with MI regsitration

[Luis Ibanez]

Hi Yannick,


Thanks for posting your code.


You are not setting the scaling parameters for the transform.


Those values are fundamental when you use the AffineTransform
because the dynamic range of the coefficients in the rotational
part of the matrix is *much smaller* than the dynamic range of
the coefficients in the translational part.



Please read the Chapter on Image Registration from the ITK
Software Guide

      http://www.itk.org/ItkSoftwareGuide.pdf


and pay particular attention to sections:

    8.5.1 "Rigid Registration in 2D"
           pdf-page 263, and its code in pdf-page 266.


                   AND


   8.6.2 "Parameter Tunning" in pdf-page 291.


If you don't to set the parameter scaling, the optimizer treats
angular units at the same range of translations, so you can easily
end up using rotations of 5 radians  !!!


Note also that in general you should prefer the use of the
CenteredAffineTransform over the standard AffineTransform.



and...

as a standard issue you should always connect an observer to your
optimizer, so you can track its progress on the parametric space.




    Regards,



         Luis



---------------------------
Yannick Allard wrote:

> Hi Luis,
> 
> here the code portion that I'm using to regsiter the output of the warper with
> the fixed image. With that, all the samples points map outside the moving
> image... the output of the warper is perfect and the two image now have the
> same dimension and spacing, so I really dont know what is the problem... I'll
> add the observer to see what wrong in the iteration :
> 
>     //Mutual information maximisation for coregistration refinement
> 
>     typedef itk::AffineTransform< double, Dimension > TransformType;
>     typedef itk::GradientDescentOptimizer OptimizerType;
>     typedef itk::LinearInterpolateImageFunction< FixedImageType, double >
> InterpolatorTypeMI;
>     typedef itk::ImageRegistrationMethod< FixedImageType, FixedImageType >
> RegistrationType;
> 
>     TransformType::Pointer transform = TransformType::New();
>     OptimizerType::Pointer optimizer = OptimizerType::New();
>     InterpolatorTypeMI::Pointer interpolatorMI = InterpolatorTypeMI::New();
>     RegistrationType::Pointer registration = RegistrationType::New();
>     registration->SetOptimizer( optimizer );
>     registration->SetTransform( transform );
>     registration->SetInterpolator( interpolatorMI );
> 
>     typedef itk::MutualInformationImageToImageMetric< FixedImageType,
> FixedImageType > MetricType;
>     MetricType::Pointer metric = MetricType::New();
>     registration->SetMetric( metric );
> 
>     metric->SetFixedImageStandardDeviation( 0.4 );
>     metric->SetMovingImageStandardDeviation( 0.4 );
>     metric->SetNumberOfSpatialSamples( 100 );
> 
>     typedef itk::NormalizeImageFilter< FixedImageType, FixedImageType >
> FixedNormalizeFilterType;
>     typedef itk::NormalizeImageFilter< MovingImageType, FixedImageType
> 
>>MovingNormalizeFilterType;
> 
>     FixedNormalizeFilterType::Pointer fixedNormalizer =
> FixedNormalizeFilterType::New();
>     MovingNormalizeFilterType::Pointer movingNormalizer =
> MovingNormalizeFilterType::New();
> 
>     typedef itk::DiscreteGaussianImageFilter< FixedImageType, FixedImageType >
> GaussianFilterType;
>     GaussianFilterType::Pointer fixedSmoother = GaussianFilterType::New();
>     GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();
> 
>     fixedSmoother->SetVariance( 2.0 );
>     movingSmoother->SetVariance( 2.0 );
> 
>     //Modify the following twolines of code to get the data from the databuffers
>     fixedNormalizer->SetInput( l_FixedImage );
>    movingNormalizer->SetInput( warper->GetOutput() );
>     movingNormalizer->Update();
>     fixedNormalizer->Update();
> 
>     fixedSmoother->SetInput( fixedNormalizer->GetOutput() );
>     movingSmoother->SetInput( movingNormalizer->GetOutput() );
>     fixedSmoother->Update();
>     movingSmoother->Update();
> 
>     registration->SetFixedImage( fixedSmoother->GetOutput() );
>     registration->SetMovingImage( movingSmoother->GetOutput() );
> 
>     registration->SetFixedImageRegion(
> fixedNormalizer->GetOutput()->GetBufferedRegion() );
>     typedef RegistrationType::ParametersType ParametersType;
> 
>     transform->SetIdentity();
> 
> 
>     registration->SetInitialTransformParameters( transform->GetParameters() );
> 
>     optimizer->SetLearningRate( 100.0 );
>     optimizer->SetNumberOfIterations( 200 );
>     optimizer->MaximizeOn();
> 
>     try
>     {
>       registration->StartRegistration();
>     }
>     catch( itk::ExceptionObject & err )
>     {
>       std::cerr << "ExceptionObject caught !" << std::endl;
>       std::cerr << err << std::endl;
>       return;
>     }
> 
>     ParametersType finalParameters = registration->GetLastTransformParameters();
>     typedef itk::ResampleImageFilter< MovingImageType, FixedImageType >   
> ResampleFilterType;
> 
>     TransformType::Pointer finalTransform = TransformType::New();
> 
>     // Print out results
>     for( int Index = 0; Index < transform->GetNumberOfParameters(); Index++ )
>     {
>       std::cerr<<"Tranform parameters : "<<finalParameters[ Index ]<<std::endl;
>     }
> 
>     finalTransform->SetParameters( finalParameters );
> 
>     ResampleFilterType::Pointer resample = ResampleFilterType::New();
> 
>     resample->SetTransform( finalTransform );
>     resample->SetInput( warper->GetOutput() );
>     resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
>     resample->SetOutputOrigin(  fixedImage->GetOrigin() );
>     resample->SetOutputSpacing( fixedImage->GetSpacing() );
>     resample->SetDefaultPixelValue( 0.0 );
>     resample->Update();
> 
>