[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