[Insight-users] VersorRigid3DTransform and VersorTransformOptimizer
Peter Lorenzen
lorenzen at cs . unc . edu
Sat, 17 May 2003 11:25:07 -0400
Hi All,
I am trying my hand at Rigid3D registration with ITK. Based on the
itkSoftwareGuide and several of the examples in the Insight Examples
directory I decided to use the transform class
itkVersorRigid3DTransform with the optimizer class
itkVersorTransformOptimizer.
Unfortunately, I have not been able to set the optimizer parameters in
a such a manner as to get a simple registration (e.g. registration
between an image and a rotated (by 2 degrees about a principal axis)
image) to work. The optimization appears to oscillate with cost
function value. The images I am using are simple cubes of uniform
intensity.
Additionally, the rotational component, in matrix form, of the final
transform is not always a rotation matrix (this may be a result of the
above observation but I am not sure).
Any suggestions on how to get the itkVersorTransformOptimizer working
properly with the itkVersorRigid3DTransform (if that is, indeed, the
correct thing to do in this case)?
Thanks,
- Peter
---------------------------- Code Snippet --------------------------
////////////////////////////////////
// Define registration components //
////////////////////////////////////
typedef itk::VersorRigid3DTransform<double> TransformType;
typedef itk::VersorTransformOptimizer OptimizerType;
typedef itk::MeanSquaresImageToImageMetric< ImageType, ImageType >
MetricType;
typedef itk::LinearInterpolateImageFunction< ImageType, double >
InterpolatorType;
typedef itk::ImageRegistrationMethod< ImageType, ImageType >
RegistrationType;
////////////////////////////////////
// Create registration components //
////////////////////////////////////
MetricType::Pointer metric = MetricType::New();
TransformType::Pointer transform = TransformType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
InterpolatorType::Pointer interpolator = InterpolatorType::New();
RegistrationType::Pointer registration = RegistrationType::New();
////////////////////////////////////////////////////////////////
// Connect the components to the instance of the registration //
////////////////////////////////////////////////////////////////
registration->SetMetric( metric );
registration->SetOptimizer( optimizer );
registration->SetTransform( transform );
registration->SetInterpolator( interpolator );
////////////////
// Set images //
////////////////
registration->SetFixedImage( fixedImage ); // Of type ImageType
registration->SetMovingImage( movingImage ); // Of type ImageType
registration->SetFixedImageRegion( fixedImage->GetBufferedRegion() );
//////////////////////////////////
// Translational initialization //
//////////////////////////////////
typedef itk::CenteredTransformInitializer< TransformType, ImageType,
ImageType > TransformInitializerType;
TransformInitializerType::Pointer initializer =
TransformInitializerType::New();
initializer->SetTransform( transform );
initializer->SetFixedImage( fixedImage );
initializer->SetMovingImage( movingImage );
initializer->MomentsOn();
initializer->InitializeTransform();
///////////////////////////////
// Rotational initialization //
///////////////////////////////
typedef TransformType::VersorType VersorType;
typedef VersorType::VectorType VectorType;
VectorType axis;
axis[0] = 0.0;
axis[1] = 0.0;
axis[2] = 1.0;
const double angle = 0;
VersorType rotation;
rotation.Set( axis, angle );
transform->SetRotation( rotation );
registration->SetInitialTransformParameters(
transform->GetParameters() );
///////////////////////////////
// Fine tuning the optimizer //
///////////////////////////////
optimizer->SetMaximumStepLength( 0.1 );
optimizer->SetMinimumStepLength( 0.001 );
optimizer->SetNumberOfIterations( maxIter );
///////////////////////////////////////////////////////////////////
// Set up a command observers and register it with the optimizer //
///////////////////////////////////////////////////////////////////
CommandIterationUpdate::Pointer observer =
CommandIterationUpdate::New();
optimizer->AddObserver( itk::IterationEvent(), observer );
//////////////////////////////
// Perform the registration //
//////////////////////////////
try {
std::cout << "Running registration..." << std::endl;
registration->StartRegistration();
}
catch ( itk::ExceptionObject& err ) {
std::cout << err << std::endl;
return( 1 );
}
-------- Sample Oscillate Output from attempting to register on image
with a rotated (by 2 degrees about the z-axis) version of itself ------
+---------------------------------------------------+
| Image Estimate Affine |
+---------------------------------------------------+
Running registration...
Current rotation = [ 0, 0, 0, 1 ]
factor = -0.10262
axis.Norm = 0.938404
Gradient rotaion = [ 0.000457015, 0.000499151, -0.0481264, 0.998841 ]
New rotaion = [ 0.000457015, 0.000499151, -0.0481264, 0.998841 ]
0 208.763 [0.000457015, 0.000499151, -0.0481264]
Current rotation = [ 0.000457015, 0.000499151, -0.0481264, 0.998841 ]
factor = -2.10831e-06
axis.Norm = 47429.9
Gradient rotaion = [ -0.0356385, -0.0348612, -0.0035152, 0.99875 ]
New rotaion = [ -0.0368202, -0.0326055, -0.0515755, 0.997457 ]
1 3157.85 [-0.0368202, -0.0326055, -0.0515755]
Current rotation = [ -0.0368202, -0.0326055, -0.0515755, 0.997457 ]
factor = -6.17673e-06
axis.Norm = 16188.7
Gradient rotaion = [ -0.00454395, -0.0496817, 0.00294461, 0.99875 ]
New rotaion = [ -0.043965, -0.0817774, -0.0468928, 0.994576 ]
2 1655.68 [-0.043965, -0.0817774, -0.0468928]
Current rotation = [ -0.043965, -0.0817774, -0.0468928, 0.994576 ]
factor = -2.70495e-06
axis.Norm = 36968.6
Gradient rotaion = [ -0.0367289, 0.000337154, -0.0338926, 0.99875 ]
New rotaion = [ -0.0776522, -0.0811077, -0.0835613, 0.990156 ]
3 2960.37 [-0.0776522, -0.0811077, -0.0835613]
Current rotation = [ -0.0776522, -0.0811077, -0.0835613, 0.990156 ]
factor = -1.3797e-05
axis.Norm = 7247.86
Gradient rotaion = [ -0.00506945, -0.0484109, -0.0113382, 0.99875 ]
New rotaion = [ -0.0857004, -0.129397, -0.0913354, 0.983651 ]
4 2467.91 [-0.0857004, -0.129397, -0.0913354]
Current rotation = [ -0.0857004, -0.129397, -0.0913354, 0.983651 ]
factor = -3.55664e-06
axis.Norm = 28115.7
Gradient rotaion = [ -0.0405131, 0.00119833, -0.0292413, 0.99875 ]
New rotaion = [ -0.121551, -0.126863, -0.12533, 0.976434 ]
5 4019.68 [-0.121551, -0.126863, -0.12533]
Current rotation = [ -0.121551, -0.126863, -0.12533, 0.976434 ]
factor = -2.40684e-05
axis.Norm = 4154.79
Gradient rotaion = [ -0.00669484, -0.0476132, -0.013639, 0.99875 ]
New rotaion = [ -0.132173, -0.174014, -0.133552, 0.96665 ]
6 3533.24 [-0.132173, -0.174014, -0.133552]
Current rotation = [ -0.132173, -0.174014, -0.133552, 0.96665 ]
factor = -5.30093e-06
axis.Norm = 18863.7
Gradient rotaion = [ -0.042266, 0.0122401, -0.0236948, 0.99875 ]
New rotaion = [ -0.167106, -0.159452, -0.165263, 0.958822 ]
7 4816.81 [-0.167106, -0.159452, -0.165263]
Current rotation = [ -0.167106, -0.159452, -0.165263, 0.958822 ]
factor = -6.62802e-06
axis.Norm = 7543.66
Gradient rotaion = [ 0.0232129, -0.00631942, 0.00678872, 0.999688 ]
New rotaion = [ -0.146924, -0.168163, -0.153945, 0.962515 ]
8 4335.18 [-0.146924, -0.168163, -0.153945]
Current rotation = [ -0.146924, -0.168163, -0.153945, 0.962515 ]
factor = -2.18005e-06
axis.Norm = 11467.3
Gradient rotaion = [ -0.0112489, -0.00290099, -0.00461309, 0.999922 ]
New rotaion = [ -0.157411, -0.169888, -0.159838, 0.959589 ]
9 4358.47 [-0.157411, -0.169888, -0.159838]
+---------------------------------------------------+
| Image Estimate Affine...Done! |
+---------------------------------------------------+
Matrix =
0.891179 0.360242 -0.275725
-0.253274 0.899347 0.356408
0.376366 -0.24779 0.89272
Offset =
-0.149975 0.0447058 -0.0572304