[Insight-users] Registration of US and CT images
Luis Ibanez
luis.ibanez@kitware.com
Tue May 11 22:49:54 EDT 2004
Hi Anders
The easy way to initialize the VersorRigid3DTransform
is to first call SetIdentity() on it, and then use
the CenteredTransformInitializer in order to set the
initial translation.
The Versor is not limited to manage rotations in one
direction. It will represent a rotation around an
arbitrary axis in 3D space.
Note that the VersorRigid3D transform can only represent
Rigid transforms not generic Affine transforms. If you
applied an Affine transform to your dataset and this
Affine transform included anisotropic scaling or shearing,
the VersorRigid3DTransform *will not* be able to correct
for those deformations. You are free to use an Affine
transform for initializing the dataset as long as you
only set rotations and translation in that transform.
Actually, if you want to configure a test, the easy way
to go is to miss-register the image on purpose using
a VersorRigid3DTransform. Then you perform the registration
and verify if the solution is the inverse of the transform
that you manually applied at the beginning.
---
The optimzerScales are the usual suspect in any
problem related to Transforms combining translations,
rotations and scaling. Your values seems to be ok for
the scale of the image though.
The next usual suspect is the step length of the optimizer.
I couldn't find in your code the setup for this parameters.
Did I missed it ?
or is that you omitted to set this parameter ?
The step length is critical in the optimization process.
You should start with very small values, for example
0.01, and only increase them if you see that the transform
parameters move too slowly.
The best way for you to get a grip on the appropriate
parameters for the optimizer is to manually subsample the
images to the lowest resolution level of the pyramid, and
feed those images in a single-resolution level program.
Fine tune the parameters of the optimizer for that level
until you get convergence.
Then you can move to the next resolution level, and so on...
Note that *there is no reason* to expect that the optimization
parameters that were useful at one resolution level will
also be appropriate for another resolution level. They are
a good starting point though.
Once you figure out the parameters for every resolution level
you can come back to the multi-resolution level code and
make the Command/Observer set up the appropriate parameters
for each level as part of its response to the iteration (the
change of resolution level).
Regards,
Luis
-----------------
Anders Almås wrote:
> Hi!
>
> I am performing registration of US images and CT images in 3D. So far I
> am only using Mattes mutual information metric. Now my ragistration is
> only performed on a perfect aligned set of US and CT images. What I do
> is to initially rotate and translate the CT image and try to registrate
> it such that the translated/rotated image is aligned again. First i
> used the Centered Affine transform, but that gave only results which
> was sometimes ok for only rotation. Then when I only tested for
> translations the translation transform worked quite good. At present
> time I am testing the VersorRigid3DTransform. This one works ok at
> translations and it seems that it manages rotations in one direction.
> My question for VersorRigid3DTransform is how it should be initialized?
> I send the code so that any one can comment it... Is it the
> optimizerscales which are the problem. Or is it really that I perform
> an affine transformation/rotation that makes it difficult to registrate
> it back? I hope anyone can help me about this...
> The size of the images is 256*256*176 with spacing 1.453,1.453 and 2
>
> Here follows the code for the registration:
>
> std::cout<<"Here starts the registration process..."<<std::endl;
>
>
> //TranslationTransform
> // typedef itk::TranslationTransform<double, ImageDimension>
> TranslationType;
> //Rigid3DTransform
> typedef itk::VersorRigid3DTransform<double> TranslationType;
> // typedef itk::CenteredAffineTransform< double, ImageDimension >
> TranslationType;
> //Transform filters
> //Optimizer filters
> typedef itk::VersorRigid3DTransformOptimizer OptimizerType;
>
> // typedef itk::RegularStepGradientDescentOptimizer
> OptimizerType;
>
> //Image Function filters
> typedef itk::LinearInterpolateImageFunction<ImageType, double >
> InterpolatorType;
>
>
> //Metric filter
> typedef itk::MattesMutualInformationImageToImageMetric< ImageType,
> ImageType > MetricType;
> //typedef itk::MutualInformationImageToImageMetric<ImageType,ImageType
>
>>MetricType;
>
>
> //RegistrationMethod
> //typedef itk::ImageRegistrationMethod< ImageType, ImageType >
> RegistrationType;
> typedef itk::MultiResolutionImageRegistrationMethod< ImageType,
> ImageType > RegistrationType;
>
> typedef itk::MultiResolutionPyramidImageFilter<ImageType,ImageType >
> FixedImagePyramidType;
> typedef itk::MultiResolutionPyramidImageFilter<ImageType,ImageType >
> MovingImagePyramidType;
>
> typedef OptimizerType::ScalesType OptimizerScalesType;
>
> TranslationType::Pointer transform2 = TranslationType::New();
> OptimizerType::Pointer optimizer = OptimizerType::New();
> InterpolatorType::Pointer interpolator2 = InterpolatorType::New();
> RegistrationType::Pointer registration = RegistrationType::New();
> MetricType::Pointer metric = MetricType::New();
>
>
> FixedImagePyramidType::Pointer fixedImagePyramid =
> FixedImagePyramidType::New();
> MovingImagePyramidType::Pointer movingImagePyramid =
> MovingImagePyramidType::New();
>
> typedef itk::CenteredTransformInitializer< TranslationType, ImageType,
> ImageType > TransformInitializerType;
> TransformInitializerType::Pointer initializer =
> TransformInitializerType::New();
> initializer->SetTransform( transform2 );
> initializer->SetFixedImage( range->GetOutput());
> // initializer->SetFixedImage( gaussian->GetOutput());
> initializer->SetMovingImage( filter->GetOutput() );
> //initializer->SetMovingImage( ctimage );
> initializer->MomentsOn();
> initializer->InitializeTransform();
>
> typedef TranslationType::VersorType VersorType;
> typedef VersorType::VectorType VectorType;
>
> VersorType rotation;
> VectorType axis;
>
> axis[0] = 0.0;
> axis[1] = 0.0;
> axis[2] = 1.0;
>
> const double angle = 0;
>
> rotation.Set( axis, angle );
>
> transform2->SetRotation( rotation );
> axis[0] = 1.0;
> axis[1] = 0.0;
> axis[2] = 0.0;
>
> rotation.Set( axis, angle );
>
> transform2->SetRotation( rotation );
>
> axis[0] = 0.0;
> axis[1] = 1.0;
> axis[2] = 0.0;
>
> // const double angle = 0;
>
> rotation.Set( axis, angle );
>
> transform2->SetRotation( rotation );
> //transform2->SetIdentity();
> registration->SetTransform(transform2);
> registration->SetOptimizer( optimizer );
> registration->SetInterpolator( interpolator2 );
> registration->SetMetric( metric );
> registration->SetFixedImagePyramid( fixedImagePyramid );
> registration->SetMovingImagePyramid( movingImagePyramid );
>
> OptimizerScalesType optimizerScales(
> transform2->GetNumberOfParameters() );
>
> std::cout<<"Number of parametres:
> "<<transform2->GetNumberOfParameters()<<std::endl;
>
>
> optimizerScales[0] = 1.0;
> optimizerScales[1] = 1.0;
> optimizerScales[2] = 1.0;
> optimizerScales[3] = 1.0/1000;
>
> optimizerScales[4] = 1.0/1000;
> optimizerScales[5] = 1.0/1000 ;
>
> optimizer->SetScales( optimizerScales );
> optimizer->MaximizeOn();
>
> registration->SetFixedImage( range->GetOutput()); //US IMAGE EDGE
> MAP
> registration->SetMovingImage( filter->GetOutput() ); // CTimage
> gaussian blurred translated and rotated...
> registration->SetFixedImageRegion(
> range->GetOutput()->GetBufferedRegion() );
>
> typedef RegistrationType::ParametersType ParametersType;
>
> /*
> ParametersType initialParameters( transform->GetNumberOfParameters()
> );
>
> initialParameters[0] = 0.0; // Initial offset in mm along X
> initialParameters[1] = 0.0; // Initial offset in mm along Y
> initialParameters[2] = 0.0; // Initial offset in mm along Z
> */
> registration->SetInitialTransformParameters(
> transform2->GetParameters() );
>
> metric->SetNumberOfHistogramBins( par[7] ); // 20
> metric->SetNumberOfSpatialSamples( par[8] ); //10000
>
> optimizer->SetNumberOfIterations( par[6] ); // Tested from 300 - 800
> iterations
>
> CommandIterationUpdate::Pointer observer =
> CommandIterationUpdate::New();
> optimizer->AddObserver( itk::IterationEvent(), observer );
>
>
>
> typedef RegistrationInterfaceCommand<RegistrationType> CommandType;
> CommandType::Pointer command = CommandType::New();
> registration->AddObserver( itk::IterationEvent(), command );
>
> registration->SetNumberOfLevels( par[9] ); //Tested from 3-5 levels..
>
> try
> {
> registration->StartRegistration();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cout << "ExceptionObject caught !" << std::endl;
> std::cout << err << std::endl;
>
> return input;
> }
>
> ParametersType finalParameters =
> registration->GetLastTransformParameters();
>
> unsigned int numberOfIterations = optimizer->GetCurrentIteration();
>
>
> double bestValue = optimizer->GetValue();
>
>
> // Print out results
> //
> std::cout << "Result = " << std::endl;
> std::cout << " Iterations = " << numberOfIterations << std::endl;
> std::cout << " Metric value = " << bestValue << std::endl;
> //std::cout << " Stop Condition = " << optimizer->GetStopCondition()
> << std::endl;
>
> I really hope somebody can help me with some advice. Maybe is it that my
> parameters is set wrong??
>
> Regards
>
> Anders
>
> _______________________________________________
> Insight-users mailing list
> Insight-users@itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>
More information about the Insight-users
mailing list