[Insight-users] problem with MI regsitration

[Luis Ibanez]

Hi Yannick,

Mutual Information is *very* noisy. The fact that you are seeing a
Metric plot that is *not monotonically decreasing* doesn't necessarily
mean that the optimizer is not *trying* to minimize the cost function.

Note also that every run of the Mutual Information metric will give
you different values because the points used for computing the metrics
are randomly selected. If you want to reproduce a run you have to make
sure that you initialize the seed of the random number generator.
Otherwise, every run will give you a different metric plot, even if you
set the exact same input parameters.

One way of getting around this is to use Evolutionary Optimizers that
better suited for noisy cost functions. You will find a use of the
OnePlusOneEvolutionary optimizer along with Mutual Information in the
files

        Insight/Examples/Registration/
                    ImageRegistration11.cxx
                    ImageRegistration14.cxx


Regards,


    Luis


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

> Hi Luis,
> 
> It is just that the optimizer was not "maximizing" the MI when I called the
> MaximizeOn() function of the optimizer but I did maximized it when I call the
> MinimizeOff(). I looked at the source code and the MinimizeOff() call the
> MaximizeOn()... I'm currently playing with the Learning rate of the optimizer
> to see if it is not the actual problem... maybe it was set a bit too high and
> therefore the registration was not converging properly. In fact I'm still
> observing some strange jump of the MI during optimization...
> 
> #iter  MI
> 
> 131   0.166983
> 132   0.174018
> 133   0.174946
> 134   0.183752
> 135   0.175778
> 136   0.178464
> 137   0.133462   [
> 138   0.18057
> 139   0.207989
> 140   0.206494
> 141   0.20218
> 142   0.184866
> 143   0.161572
> 144   0.182577
> 145   0.195277
> 
> 
> I'll let you know.
> 
> Thank you
> 
> Yannick
> 
> Selon Luis Ibanez <luis.ibanez at kitware.com>:
> 
> 
>>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();
>>>>>
>>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>
>>>
>>>
>>>
>>
>>
>>
>>
>>
> 
> 
> 
> 
>