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

[Gunther Sudra]

Hi,

at first: thank you Luis, for your prompt answers. Although I understand now the figure, I have still problems with my example:

All questions refer to my added example-code which is a smaller version of the ITKSoftwareGuide example Model-Based-Registration, using only one ellipse instead of three

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

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

3. Thank you for answering all these beginner-questions.



My code:

typedef itk::GroupSpatialObject< 2 > GroupType;
typedef itk::EllipseSpatialObject< 2 > EllipseType;

typedef itk::Image< float, 2 > ImageType;

EllipseType::Pointer ellipse1 = EllipseType::New();

double radius[2] = {10.0, 15.0};
ellipse1->SetRadius( radius );

EllipseType::TransformType::OffsetType offset;
offset[ 0 ] = 100.0;
offset[ 1 ] = 40.0;

ellipse1->GetObjectToParentTransform()->SetOffset(offset);
ellipse1->ComputeObjectToWorldTransform();

GroupType::Pointer group = GroupType::New();
group->AddSpatialObject( ellipse1 );


typedef itk::SpatialObjectToImageFilter< GroupType, ImageType > 
SpatialObjectToImageFilterType;

SpatialObjectToImageFilterType::Pointer imageFilter = 
SpatialObjectToImageFilterType::New();

imageFilter->SetInput( group );

ImageType::SizeType size;
size[ 0 ] = 200;
size[ 1 ] = 200;
imageFilter->SetSize( size );

imageFilter->Update();

typedef itk::DiscreteGaussianImageFilter< ImageType, ImageType > 
GaussianFilterType;
GaussianFilterType::Pointer gaussianFilter = GaussianFilterType::New();

gaussianFilter->SetInput( imageFilter->GetOutput() );

const double variance = 20;
gaussianFilter->SetVariance(variance);
gaussianFilter->Update();

typedef itk::ImageToSpatialObjectRegistrationMethod< ImageType, GroupType >
RegistrationType;
RegistrationType::Pointer registration = RegistrationType::New();

typedef SimpleImageToSpatialObjectMetric< ImageType, GroupType > MetricType;
MetricType::Pointer metric = MetricType::New();

typedef itk::LinearInterpolateImageFunction< ImageType, double > 
InterpolatorType;
InterpolatorType::Pointer interpolator = InterpolatorType::New();

typedef itk::OnePlusOneEvolutionaryOptimizer OptimizerType;
OptimizerType::Pointer optimizer = OptimizerType::New();

typedef itk::Euler2DTransform<> TransformType;
TransformType::Pointer transform = TransformType::New();

itk::Statistics::NormalVariateGenerator::Pointer generator 
= itk::Statistics::NormalVariateGenerator::New();

generator->Initialize(12345);

optimizer->SetNormalVariateGenerator( generator );

//100
optimizer->Initialize( 100 );

//10000
optimizer->SetMaximumIteration( 1000 );

TransformType::ParametersType parametersScale;
parametersScale.set_size(3);
parametersScale[0] = 1000; // angle scale

for( unsigned int i=1; i<3; i++ )
{
parametersScale[i] = 2; // offset scale
}
optimizer->SetScales( parametersScale );

typedef IterationCallback< OptimizerType > IterationCallbackType;
IterationCallbackType::Pointer callback = IterationCallbackType::New();
callback->SetOptimizer( optimizer );


//////////////////////////////////////////////////////////////////////////
// model with same radius which is translated by x and y
//////////////////////////////////////////////////////////////////////////
int x = 10;
int y = 10;


EllipseType::Pointer ellipse_model = EllipseType::New();

ellipse_model->SetRadius( radius );

offset[ 0 ] = 100.0 + x;
offset[ 1 ] = 40.0 + y;

ellipse_model->GetObjectToParentTransform()->SetOffset(offset);
ellipse_model->ComputeObjectToWorldTransform();


GroupType::Pointer group_model= GroupType::New();
group_model->AddSpatialObject( ellipse_model );


registration->SetFixedImage( gaussianFilter->GetOutput() );
registration->SetMovingSpatialObject( group_model );
registration->SetTransform( transform );
registration->SetInterpolator( interpolator );
registration->SetOptimizer( optimizer );
registration->SetMetric( metric );

TransformType::ParametersType initialParameters( 
transform->GetNumberOfParameters() );

initialParameters[0] = 0.0; // Angle
initialParameters[1] = 0.0; // Offset X
initialParameters[2] = 0.0; // Offset Y 
registration->SetInitialTransformParameters(initialParameters);

std::cout << "Initial Parameters : " << initialParameters << std::endl;

optimizer->MaximizeOn();

try 
{
registration->StartRegistration();
}
catch( itk::ExceptionObject & exp ) 
{
std::cerr << "Exception caught ! " << std::endl;
std::cerr << exp << std::endl;
}

RegistrationType::ParametersType finalParameters 
= registration->GetLastTransformParameters();

std::cout << "Final Solution is : " << finalParameters << std::endl;


itk::Array mPositionArray;
mPositionArray = registration->GetLastTransformParameters();

double angle = mPositionArray[0];
TransformType::MatrixType matrix;
matrix[0][0] = cos(angle); matrix[0][1] = -sin(angle);
matrix[1][0] = sin(angle); matrix[1][1] = cos(angle);

offset[ 0 ] = 100 - mPositionArray[1];
offset[ 1 ] = 40 - mPositionArray[2];
ellipse_model->GetObjectToParentTransform()->SetOffset(offset);
ellipse_model->GetObjectToParentTransform()->SetMatrix(matrix);
ellipse_model->ComputeObjectToWorldTransform();

imageFilter->SetInput( group_model );
imageFilter->Update();

______________________________________________________________________________
Ein Grund zum Feiern: Die PC Praxis ermittelt zwischen 10 grossen
Mailprovidern WEB.DE FreeMail als Testsieger http://f.web.de/?mc=021190