[Insight-users] Model-Based-Registration results with one ellipse (code included)

[Julien Jomier]

Hi Gunther,

Thanks for providing the code. I just tried and here are some comments.
The problem using only one ellipse results in an ambiguity between the
translation and the rotation (you can get the same results by applying a
large rotation and small translation or the inverse), assuming the =
center of
rotation is at [0,0].

To answer your first question:

> 1. Question: I found out that I get very good results when I=20
> feed the optimizer with the same synthetic image and=20
> incorrect small initial values, e.g. [0,10,10]. But when I=20
> move the ellipse in the synthetic image by 10 in x-direction=20
> (variable x), and 10 in y-direction (variable y) and set the=20
> initial values to [0,0,0] I get bad results. What am I doing=20
> wrong? What is the difference?

If you set the initial parameters to [0,10,10] you should offset the =
image
by [-10,-10] to get the same results. In your case, you offset the image =
by
[10,10] and I think the optimizer is trying to optimize the rotation =
more
than the translation and that's why you get different results.

and the second one:

> 2. Question: A strong increase of the optimizer parameters=20
> (optimizer->Initialize( 10000 ) and=20
> optimizer->SetMaximumIteration( 10000 ) seem to have almost=20
> no effect on=20
> optimizer->the quality of the result.

You should initialize the optimizer with the maximum displacement you =
want
to optimize for. This means if you expect an offset of 10 you should
Initialize(10).
The 1+1 evolutionary optimizer is setting the scale of the parameters
automatically during the optimization, that's why even if you set a =
large
scale for rotation, it will optimize it anyway. On the other hand, this
optimizer converges quickly so increasing the number of iterations won't
really help.

You should try to solve the ambiguity with the rotation by setting the
center of rotation as the center of the ellipse. This should be done in =
the
metric where the points are being transformed.

hope this helps,

Julien


> -----Original Message-----
> From: insight-users-admin at itk.org=20
> [mailto:insight-users-admin at itk.org] On Behalf Of Gunther Sudra
> Sent: Thursday, February 19, 2004 4:33 PM
> To: insight-users at itk.org
> Subject: [Insight-users] Model-Based-Registration results=20
> with one ellipse (code included)
>=20
>=20
> Hi,
>=20
> at first: thank you Luis, for your prompt answers. Although I=20
> understand now the figure, I have still problems with my example:
>=20
> All questions refer to my added example-code which is a=20
> smaller version of the ITKSoftwareGuide example=20
> Model-Based-Registration, using only one ellipse instead of three
>=20
> 1. Question: I found out that I get very good results when I=20
> feed the optimizer with the same synthetic image and=20
> incorrect small initial values, e.g. [0,10,10]. But when I=20
> move the ellipse in the synthetic image by 10 in x-direction=20
> (variable x), and 10 in y-direction (variable y) and set the=20
> initial values to [0,0,0] I get bad results. What am I doing=20
> wrong? What is the difference?
>=20
> 2. Question: A strong increase of the optimizer parameters=20
> (optimizer->Initialize( 10000 ) and=20
> optimizer->SetMaximumIteration( 10000 ) seem to have almost=20
> no effect on=20
> optimizer->the quality of the result.
>=20
> 3. Thank you for answering all these beginner-questions.
>=20
>=20
>=20
> My code:
>=20
> typedef itk::GroupSpatialObject< 2 > GroupType;
> typedef itk::EllipseSpatialObject< 2 > EllipseType;
>=20
> typedef itk::Image< float, 2 > ImageType;
>=20
> EllipseType::Pointer ellipse1 =3D EllipseType::New();
>=20
> double radius[2] =3D {10.0, 15.0};
> ellipse1->SetRadius( radius );
>=20
> EllipseType::TransformType::OffsetType offset;
> offset[ 0 ] =3D 100.0;
> offset[ 1 ] =3D 40.0;
>=20
> ellipse1->GetObjectToParentTransform()->SetOffset(offset);
> ellipse1->ComputeObjectToWorldTransform();
>=20
> GroupType::Pointer group =3D GroupType::New();
> group->AddSpatialObject( ellipse1 );
>=20
>=20
> typedef itk::SpatialObjectToImageFilter< GroupType, ImageType >=20
> SpatialObjectToImageFilterType;
>=20
> SpatialObjectToImageFilterType::Pointer imageFilter =3D=20
> SpatialObjectToImageFilterType::New();
>=20
> imageFilter->SetInput( group );
>=20
> ImageType::SizeType size;
> size[ 0 ] =3D 200;
> size[ 1 ] =3D 200;
> imageFilter->SetSize( size );
>=20
> imageFilter->Update();
>=20
> typedef itk::DiscreteGaussianImageFilter< ImageType, ImageType >=20
> GaussianFilterType;
> GaussianFilterType::Pointer gaussianFilter =3D=20
> GaussianFilterType::New();
>=20
> gaussianFilter->SetInput( imageFilter->GetOutput() );
>=20
> const double variance =3D 20;
> gaussianFilter->SetVariance(variance);
> gaussianFilter->Update();
>=20
> typedef itk::ImageToSpatialObjectRegistrationMethod<=20
> ImageType, GroupType > RegistrationType;=20
> RegistrationType::Pointer registration =3D RegistrationType::New();
>=20
> typedef SimpleImageToSpatialObjectMetric< ImageType,=20
> GroupType > MetricType; MetricType::Pointer metric =3D=20
> MetricType::New();
>=20
> typedef itk::LinearInterpolateImageFunction< ImageType, double >=20
> InterpolatorType;
> InterpolatorType::Pointer interpolator =3D InterpolatorType::New();
>=20
> typedef itk::OnePlusOneEvolutionaryOptimizer OptimizerType;=20
> OptimizerType::Pointer optimizer =3D OptimizerType::New();
>=20
> typedef itk::Euler2DTransform<> TransformType;=20
> TransformType::Pointer transform =3D TransformType::New();
>=20
> itk::Statistics::NormalVariateGenerator::Pointer generator=20
> =3D itk::Statistics::NormalVariateGenerator::New();
>=20
> generator->Initialize(12345);
>=20
> optimizer->SetNormalVariateGenerator( generator );
>=20
> //100
> optimizer->Initialize( 100 );
>=20
> //10000
> optimizer->SetMaximumIteration( 1000 );
>=20
> TransformType::ParametersType parametersScale;=20
> parametersScale.set_size(3); parametersScale[0] =3D 1000; // angle =
scale
>=20
> for( unsigned int i=3D1; i<3; i++ )
> {
> parametersScale[i] =3D 2; // offset scale
> }
> optimizer->SetScales( parametersScale );
>=20
> typedef IterationCallback< OptimizerType >=20
> IterationCallbackType; IterationCallbackType::Pointer=20
> callback =3D IterationCallbackType::New();
> callback->SetOptimizer( optimizer );
>=20
>=20
> //////////////////////////////////////////////////////////////
> ////////////
> // model with same radius which is translated by x and y=20
> //////////////////////////////////////////////////////////////
> ////////////
> int x =3D 10;
> int y =3D 10;
>=20
>=20
> EllipseType::Pointer ellipse_model =3D EllipseType::New();
>=20
> ellipse_model->SetRadius( radius );
>=20
> offset[ 0 ] =3D 100.0 + x;
> offset[ 1 ] =3D 40.0 + y;
>=20
> ellipse_model->GetObjectToParentTransform()->SetOffset(offset);
> ellipse_model->ComputeObjectToWorldTransform();
>=20
>=20
> GroupType::Pointer group_model=3D GroupType::New();=20
> group_model->AddSpatialObject( ellipse_model );
>=20
>=20
> registration->SetFixedImage( gaussianFilter->GetOutput() );=20
> registration->SetMovingSpatialObject( group_model ); SetTransform(=20
> registration->transform ); SetInterpolator( interpolator );
> registration->SetOptimizer( optimizer );
> registration->SetMetric( metric );
>=20
> TransformType::ParametersType initialParameters(=20
> transform->GetNumberOfParameters() );
>=20
> initialParameters[0] =3D 0.0; // Angle
> initialParameters[1] =3D 0.0; // Offset X
> initialParameters[2] =3D 0.0; // Offset Y=20
> registration->SetInitialTransformParameters(initialParameters);
>=20
> std::cout << "Initial Parameters : " << initialParameters <<=20
> std::endl;
>=20
> optimizer->MaximizeOn();
>=20
> try=20
> {
> registration->StartRegistration();
> }
> catch( itk::ExceptionObject & exp )=20
> {
> std::cerr << "Exception caught ! " << std::endl;
> std::cerr << exp << std::endl;
> }
>=20
> RegistrationType::ParametersType finalParameters=20
> =3D registration->GetLastTransformParameters();
>=20
> std::cout << "Final Solution is : " << finalParameters << std::endl;
>=20
>=20
> itk::Array mPositionArray;
> mPositionArray =3D registration->GetLastTransformParameters();
>=20
> double angle =3D mPositionArray[0];
> TransformType::MatrixType matrix;
> matrix[0][0] =3D cos(angle); matrix[0][1] =3D -sin(angle);=20
> matrix[1][0] =3D sin(angle); matrix[1][1] =3D cos(angle);
>=20
> offset[ 0 ] =3D 100 - mPositionArray[1];
> offset[ 1 ] =3D 40 - mPositionArray[2];=20
> ellipse_model->GetObjectToParentTransform()->SetOffset(offset);
> ellipse_model->GetObjectToParentTransform()->SetMatrix(matrix);
> ellipse_model->ComputeObjectToWorldTransform();
>=20
> imageFilter->SetInput( group_model );
> imageFilter->Update();
>=20
> ______________________________________________________________
> ________________
> Ein Grund zum Feiern: Die PC Praxis ermittelt zwischen 10=20
> grossen Mailprovidern WEB.DE FreeMail als Testsieger=20
http://f.web.de/?mc=3D021190

_______________________________________________
Insight-users mailing list
Insight-users at itk.org http://www.itk.org/mailman/listinfo/insight-users