[Insight-users] problem with MI regsitration

[Luis Ibanez]

Hi Yannick

Thanks for letting us know that you solved the problem.


Just to double-check:...

What you found is that calling

               optimizer->MaximizeOn();

is not equivalent to calling

              optimizer->MinimizeOff();    ??



If that's the case, that sounds like a bug in the optimizer...


Please let us know so we can trace this in case is a real bug.



    Thanks


       Luis


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

> Hi Luis,
> 
> Thank you for everything. It did work after the scaling but that was not the
> only problem. Somehow, even if I called the MaximizeOn function, the optimizer
> was still minimizing the MI, so I had to call the MinimizeOff function to solve
> that problem.
> 
> Have a nice day
> 
> Yannick
> 
> 
> Selon Luis Ibanez <luis.ibanez at kitware.com>:
> 
> 
>>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();
>>>
>>>
>>
>>
>>
>>
>>
>>
> 
> 
> 
> 
>