[Insight-users] Image Registration / Convergence of optimizer
Koning, P.J.H. de (LKEB)
P.J.H.de_Koning at lumc.nl
Tue Nov 9 03:01:04 EST 2004
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
>>
>>
>
>
>
>
--
Patrick de Koning
More information about the Insight-users
mailing list