[Insight-users] Image Registration / Convergence of optimizer
Luis Ibanez
luis.ibanez at kitware.com
Tue Nov 9 16:22:33 EST 2004
Hi Patrick,
Thanks for posting your code and additional information.
The reason why you see a different number of iteration in the
output of the Command/Observer and the final value returned
from the optimizer is that the optimizer that you are using
send "IterationEvents" in its "AdvanceOneStep()" method in
line 117 of
Insight/Code/Numerics/
itkRegularStepradientDescentBaseOptimizer.cxx
and after that,
the iteration counter is increased in line 119 of the same file.
The method AdvanceOneStep() is the one that checks for
termination (convergence) conditions, so it may return
*before* invoking the event for its current iteration value.
(that is, you will miss the signal for iteration #28), and
when returning from this function, the value of the iteration
will be incremented by 1 in line 119.
It is arguable that we should move the increment of the iteration
counter to line 245, just *after* invoking the event. In that way the
user-perception of the iteration
counter will be consistent with the number of IterationEvents
sent.
You may want to give it a try to this code change and
let us know if that works for you. In which case we will
enter this as a bug report and proceed to fix the code
in the CVS repository.
Please let us know what you find,
Thanks,
Luis
-------------------------------
Koning, P.J.H. de (LKEB) wrote:
> On Mon, 08 Nov 2004 15:21:01 -0500, Luis Ibanez
> <luis.ibanez at kitware.com> wrote:
>
>>
>> Hi Patrick,
>>
>> We may need more information from you
>> in order to help you with your question.
>>
>>
>> 1) Are you using any of the Examples in
>> Insight/Examples/Registration ?
>>
>> If not, please post to the list your code where
>> you are defining the Command/Observer for this
>> optimization.
>>
> Hi Luis,
>
> This is the code I use for the Observer and optimization
>
> class CommandIterationUpdate : public itk::Command
> {
> public:
> typedef CommandIterationUpdate Self;
> typedef itk::Command Superclass;
> typedef itk::SmartPointer<Self> Pointer;
> itkNewMacro( Self );
>
> protected:
> CommandIterationUpdate() {};
>
> public:
> typedef itk::VersorRigid3DTransformOptimizer OptimizerType;
> typedef const OptimizerType *OptimizerPointer;
> void Execute(itk::Object *caller, const itk::EventObject &
> event){Execute((const itk::Object *)caller, event);}
>
> void Execute(const itk::Object * object, const itk::EventObject &
> event)
> {
> OptimizerPointer optimizer =
> dynamic_cast<OptimizerPointer>(object);
>
> if( typeid( event ) != typeid( itk::IterationEvent ) ) {return;}
>
> std::cerr << optimizer->GetCurrentIteration() << " = ";
> std::cerr << optimizer->GetValue() << " : \t";
> std::cerr << optimizer->GetCurrentPosition() << std::endl;
> }
> };
>
> template <typename PixelType>
> void RigidMIRegistration<PixelType>::Execute()
> {
> //////////////////////////////////////////////////////////////////////////
>
> typename NormalizerType::Pointer normalizer;
> typename CentererType::Pointer centerer;
> typename TransformType::Pointer postTransform, preTransform, transform;
> typename TransformType::OffsetType offset;
>
> postTransform = TransformType::New();
> preTransform = TransformType::New();
>
> // Process the fixed image
> centerer = CentererType::New();
> centerer->CenterImageOn();
> centerer->SetInput(m_FixedImage);
> centerer->Update();
>
> preTransform->SetIdentity();
> for(unsigned int j = 0; j < ImageDimension; j++)
> offset[j] = centerer->GetOutput()->GetOrigin()[j];
> preTransform->SetOffset(offset);
>
> normalizer = NormalizerType::New();
> normalizer->SetInput(centerer->GetOutput());
> normalizer->Update();
> m_OutputFixedImage = normalizer->GetOutput();
>
> // Process the moving image
> centerer = CentererType::New();
> centerer->CenterImageOn();
> centerer->SetInput(m_MovingImage);
> centerer->Update();
>
> postTransform->SetIdentity();
> for(unsigned int j = 0; j < ImageDimension; j++)
> offset[j] = -1.0 * centerer->GetOutput()->GetOrigin()[j];
> postTransform->SetOffset(offset);
>
> normalizer = NormalizerType::New();
> normalizer->SetInput(centerer->GetOutput());
> normalizer->Update();
> m_OutputMovingImage = normalizer->GetOutput();
>
> //////////////////////////////////////////////////////////////////////////
>
> // Set up internal registrator with default components
> transform = TransformType::New();
> m_RigidOptimizer = RigidOptimizerType::New();
> m_Metric = MetricType::New();
> m_Interpolator = InterpolatorType::New();
> m_FixedImagePyramid = FixedImagePyramidType::New();
> m_MovingImagePyramid = MovingImagePyramidType::New();
> m_Registration = RegistrationType::New();
>
> m_Registration->SetTransform(transform);
> m_Registration->SetOptimizer(m_RigidOptimizer);
> m_Registration->SetMetric(m_Metric);
> m_Registration->SetInterpolator(m_Interpolator);
> m_Registration->SetFixedImagePyramid(m_FixedImagePyramid);
> m_Registration->SetMovingImagePyramid(m_MovingImagePyramid);
>
> // Setup an registration observer
> typedef itk::SimpleMemberCommand<RigidMIRegistration> CommandType;
> typename CommandType::Pointer command = CommandType::New();
> command->SetCallbackFunction(this,
> &RigidMIRegistration::StartNewLevel);
> m_Registration->AddObserver(itk::IterationEvent(), command);
>
> CommandIterationUpdate::Pointer observer =
> CommandIterationUpdate::New();
> m_RigidOptimizer->AddObserver( itk::IterationEvent(), observer );
>
> ParametersType initialParameters =
> ParametersType(transform->GetNumberOfParameters());
> initialParameters.Fill(0.0);
>
> // Setup the optimizer
> typename OptimizerType::ScalesType
> scales(transform->GetNumberOfParameters());
> scales.Fill(1.0);
>
> for (int j = 3; j < 6; j++)
> scales[j] = m_TranslationScale[0];
>
> m_RigidOptimizer->SetScales(scales);
> m_RigidOptimizer->MaximizeOn();
>
> // Setup the metric
> m_Metric->SetMovingImageStandardDeviation(m_MovingImageStandardDeviation);
>
> m_Metric->SetFixedImageStandardDeviation(m_FixedImageStandardDeviation);
>
> m_Metric->SetNumberOfSpatialSamples(m_NumberOfSpatialSamples);
>
> // Setup the image pyramids
> m_FixedImagePyramid->SetNumberOfLevels(m_NumberOfLevels);
> m_FixedImagePyramid->SetStartingShrinkFactors(m_FixedImageShrinkFactors.GetDataPointer());
>
>
> m_MovingImagePyramid->SetNumberOfLevels(m_NumberOfLevels);
> m_MovingImagePyramid->SetStartingShrinkFactors(m_MovingImageShrinkFactors.GetDataPointer());
>
>
> // Setup the registrator
> m_Registration->SetFixedImage(m_OutputFixedImage);
> m_Registration->SetMovingImage(m_OutputMovingImage);
> m_Registration->SetNumberOfLevels(m_NumberOfLevels);
> m_Registration->SetInitialTransformParameters(initialParameters);
> m_Registration->SetFixedImageRegion(m_FixedImage->GetBufferedRegion());
>
> try
> {
> m_Registration->StartRegistration();
> }
> catch(itk::ExceptionObject & err)
> {
> std::cerr << "ExceptionObject caught! " << err << std::endl;
> throw err;
> }
>
> std::cerr << m_RigidOptimizer->GetCurrentIteration() << " = ";
> std::cerr << m_RigidOptimizer->GetValue() << " : \t";
> std::cerr << m_RigidOptimizer->GetCurrentPosition() << std::endl;
>
> m_TotalTransform = TransformType::New();
> m_TotalTransform->SetIdentity();
> m_TotalTransform->Compose(postTransform, true);
> m_TotalTransform->Compose(transform, true);
> m_TotalTransform->Compose(preTransform, true);
> }
>
>
>>
>> 2) Did you set up the Scaling parameters for the
>> optimizer ?
>>
>> If so, please post to the list the values that
>> you used for the scaling parameters.
>
>
> These are the scaling parameters:
> scaling[0] = scaling[1] = scaling[2] = 1;
> scaling[3] = scaling[4] = scaling[5] = 1e-5
>
>>
>>
>> 3) It is common for Mutual Information to return
>> *very* noisy values. You may want to look at
>> the ITK Software Guide
>>
>> http://www.itk.org/ItkSoftwareGuide.pdf
>>
>> to see the typical Metric plots resulting from
>> Mutual Information. This noisy behavior is related
>> to the fact that the metric is computed from a
>> small number of randomly selected samples, this
>> samples are different at each iteration of the
>> optimizer.
>>
>> Look for example at Figure 8.11 in pdf-page 260.
>>
>
> Thanks, I missed that figure
>
>>
>>
>> 4) Note that the Mattes implementation of Mutual
>> Information produces much smoother results,
>> as you can see from Figure 8.13 in pdf-page 264.
>>
>> When you use Mutual Information along with
>> gradient descent optimizer, you probably want
>> to use the Mattes implementation instead of
>> the Viola-Wells implementation.
>>
>
> I will give that a try today, thanks for the advice.
>
>>
>>
>>
>> Regards,
>>
>>
>> Luis
>>
>>
>>
>>
>> --------------------------------
>> Koning, P.J.H. de (LKEB) wrote:
>>
>>> Hello,
>>>
>>> I am trying to register two 3D brain images. I am using the
>>> MutualInformationImageToImageMetric combined with a
>>> VersorRigid3DTransform and the VersorRigid3DTransformOptimizer. I am
>>> getting some good results, but when I look at the GetValue() of the
>>> optimizer during the iterations, I can't see it convergence (although
>>> the parameters of the rigid transform do show a convergence).
>>>
>>> Is this value not supposed to converge to a maximum value (I set
>>> MaximizeOn() for the optimzer)?
>>>
>>> Secondly, when I look at the values of the optimizer after the
>>> optimization has finished, I see that the current iteration has
>>> increased by 2, the transfor parameters stay constant and the GetValue()
>>> of the optimizer has changed ?!? How did this happen?
>>>
>>> A hope you can shed some light on this (for me) strange behaviour.
>>>
>>>
>>> Patrick de Koning
>>>
>>>
>>>
>>>
>>> P.S. Below you can see the output of the optimizer at each iteration:
>>>
>>> (Iteration number = optimizer value : [transform])
>>> 0 = 0.119745 : [-1.38545e-005 -2.53032e-006 ,7.4684e-006 ,-0.186788
>>> ,0.179258 ,4.99329
>>> 1 = 0.0503912 : [-1.71608e-005, -8.88087e-007, 9.23982e-006,
>>> -0.191543, 0.331362, 2.49793]
>>> 2 = 0.55825 : [-1.74154e-005, -3.9016e-007, 9.85893e-006, -0.204029,
>>> 0.326522, 3.74786]
>>> 3 = 0.395624 : [-1.77931e-005, 1.00222e-006, 1.00231e-005,
>>> -0.193226, 0.306919, 3.12326]
>>> 4 = 0.3617 : [-1.71888e-005, -1.22213e-006, 7.93444e-006, -0.200906,
>>> 0.316677, 2.49838]
>>> 5 = 0.335809 : [-1.81236e-005, -2.72138e-006, 6.68041e-006,
>>> -0.16859, 0.298884, 2.8087]
>>> 6 = 0.414256 : [-1.81705e-005, -2.66581e-006, 7.14423e-006,
>>> -0.176486, 0.29259, 2.65277]
>>> 7 = -0.0230083 : [-1.8131e-005, -2.49584e-006, 6.81443e-006,
>>> -0.179083, 0.28925, 2.49658]
>>> 8 = 0.328098 : [-1.82165e-005, -2.45698e-006, 7.11142e-006,
>>> -0.174855, 0.283865, 2.5744]
>>> 9 = -0.0404819 : [-1.81304e-005, -2.5165e-006, 7.35304e-006,
>>> -0.177024, 0.286958, 2.53553]
>>> 10 = 0.0277824 : [-1.81428e-005, -2.51756e-006, 7.50534e-006,
>>> -0.175715, 0.28891, 2.49653]
>>> 11 = 0.241717 : [-1.81675e-005, -2.4732e-006, 7.59313e-006,
>>> -0.174454, 0.290408, 2.45752]
>>> 12 = 0.198569 : [-1.81956e-005, -2.50076e-006, 7.62124e-006,
>>> -0.174316, 0.289997, 2.47705]
>>> 13 = 0.211859 : [-1.81711e-005, -2.48456e-006, 7.67883e-006,
>>> -0.175431, 0.28911, 2.46739]
>>> 14 = 0.271563 : [-1.81712e-005, -2.50097e-006, 7.66833e-006,
>>> -0.175269, 0.288967, 2.45762]
>>> 15 = 0.22593 : [-1.81791e-005, -2.48652e-006, 7.63361e-006,
>>> -0.175413, 0.288978, 2.44786]
>>> 16 = 0.954839 : [-1.81924e-005, -2.48619e-006, 7.62184e-006,
>>> -0.175887, 0.289044, 2.4381]
>>> 17 = 0.245999 : [-1.81852e-005, -2.46701e-006, 7.60605e-006,
>>> -0.175916, 0.289297, 2.44298]
>>> 18 = 0.163771 : [-1.81861e-005, -2.46174e-006, 7.60433e-006,
>>> -0.175958, 0.289195, 2.44054]
>>> 19 = 0.472102 : [-1.82019e-005, -2.47041e-006, 7.59436e-006,
>>> -0.175931, 0.289398, 2.43811]
>>> 20 = 0.306671 : [-1.8204e-005, -2.47435e-006, 7.58475e-006,
>>> -0.175902, 0.289338, 2.43567]
>>> 21 = 0.362819 : [-1.8203e-005, -2.47463e-006, 7.5833e-006,
>>> -0.175974, 0.289261, 2.43688]
>>> 22 = 0.33775 : [-1.82027e-005, -2.47465e-006, 7.58338e-006,
>>> -0.175984, 0.289256, 2.43627]
>>> 23 = 0.062375 : [-1.82026e-005, -2.47405e-006, 7.58279e-006,
>>> -0.175974, 0.289248, 2.43566]
>>> 24 = 0.294642 : [-1.82027e-005, -2.47339e-006, 7.58462e-006,
>>> -0.17596, 0.289238, 2.43505]
>>> 25 = 0.375321 : [-1.82031e-005, -2.47345e-006, 7.58648e-006,
>>> -0.175964, 0.289261, 2.43536]
>>> 26 = 0.295915 : [-1.82031e-005, -2.47429e-006, 7.58632e-006,
>>> -0.175959, 0.289254, 2.43521]
>>> 27 = 0.488161 : [-1.82028e-005, -2.47451e-006, 7.58712e-006,
>>> -0.175985, 0.289253, 2.43506]
>>>
>>> [End of optimization]
>>>
>>> 29 = 0.366311 : [-1.82028e-005, -2.47451e-006, 7.58712e-006,
>>> -0.175985, 0.289253, 2.43506]
>>> _______________________________________________
>>> Insight-users mailing list
>>> Insight-users at itk.org
>>> http://www.itk.org/mailman/listinfo/insight-users
>>>
>>>
>>
>>
>>
>>
>
>
>
More information about the Insight-users
mailing list