[Insight-users] VersorRigid3DTransform and VersorTransformOptimizer
Luis Ibanez
luis . ibanez at kitware . com
Mon, 19 May 2003 09:51:38 -0400
Hi Peter,
It looks like there is some confusion here.
The itk::VersorTransformOptimizer is not intended to be used with the
itk::VersorRigid3DTransform but with the itk::VersortTransform.
VersorTransform --> performs pure 3D rotations
VersorRigid3DTransform --> performs 3D rotations plus translations
The VersorTransform optimizer was introduced in order to cope
with the fact that the space of Quaternions is not a vector space.
Typical gradient descent optimizers assume that the parametric space
is vectorial. (i.e you can modify any parameters independently of
the others).
It is actually surprising that the combination runs, since the
VectorRigid3DTransform produces an array of parameters with 6
components (3 for the versor and 3 for the translation).
The VersorTransformOptimizer is hard-coded to work on 3 parameters
since the VersorTransform only needs 3 parameters to encode a Versor
(a unit quaternion). So probably the 3 translation parameters are
never touched by the optimizer.
We are missing an optimizer specific for the VersorRigid3DTransform.
The difficulty in implementing such optimizer is that it will have to
treat 3 parameters as non vectorial and the other 3 parameters as
vectorial. This is not impossible, but may result a bit tricky to
get it working right.
This new optimizer could be a sister class of the
VersorTransfromOptimizer
http://www.itk.org/Insight/Doxygen/html/classitk_1_1VersorTransformOptimizer.html
or even derive from it.
Note that there is only one method in the implementation of the
VersorTransformOptimizer. "StepAlongGradient()", this method is
implemented following the definition of derivatives in quaternion
spaces (which is different from a vector space).
The rest of the optimizer behavior is in the base class, the
RegularStepGradientDescentBaseOptimizer
http://www.itk.org/Insight/Doxygen/html/classitk_1_1RegularStepGradientDescentBaseOptimizer.html
Implementing a VersorRigid3DTransform will require to implement
an appropriate version of the "StepAlongGradient()" method.
Regards,
Luis
---------------------------------
Peter Lorenzen wrote:
> 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
> _______________________________________________
> Insight-users mailing list
> Insight-users@public.kitware.com
> http://public.kitware.com/mailman/listinfo/insight-users
>