[Insight-users] problem with MI regsitration
Koning, P.J.H. de (LKEB)
P.J.H.de_Koning at lumc.nl
Fri Nov 26 09:53:24 EST 2004
On Wed, 10 Nov 2004 18:09:18 -0500, Luis Ibanez <luis.ibanez at kitware.com> wrote:
You mentioned that you need to initialize the seed of the random number generator in order to reproduce a run. How do I do that? What kind of random number generator is used?
>
> 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();
>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>>
>>>>
>>>
>>>
>>>
>>>
>>>
>>
>>
>>
>>
>>
>
>
>
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>
--
Patrick de Koning
More information about the Insight-users
mailing list