[ITK] [ITK-users] 3D Affine Registration - no rotation found
Jerome Plumat
j.plumat at auckland.ac.nz
Thu Jul 3 19:23:08 EDT 2014
Hi Nicolás,
Thank for the answer and advices. Indeed you are right: the value 1/1E7
is too small. I set it to 1/1E2 and the registration is better
performed. I also increased the number of samples but the scale value
has a huge influence on the results.
However, I don't understand why this value has such an impact on the
rotation optimization. If I understand well the user guide section 8.7.2
top of page 359 the first optimizerScales[0] to
optimizerScales[Dimension*Dimension-1] are relative to the moments of
the affine transform while optimizerScales[Dimension*Dimension] to
optimizerScales[(Dimensions+1)*Dimension-1] should influence the
translations.
Any light?
Cheers,
Jerome Plumat
-------
School of Medical Sciences
University of Auckland
If I am not for myself, who will be for me? And if I am only for myself, then what an I? And if not now when? – Hillel HaZaken
On 03/07/14 20:12, Nicolas Gallego wrote:
> Hi Jérôme,
>
> Have you tried changing the optimizer scales?
> That parameter changes the sensitivity of the optimization process
> between the translation and rotation of the transform, it seams to me
> that the value 1.0/1e7 is probably too small and the optimizer is not
> sensitive enough to rotations.
>
> For test purposes you can also set the mutual information metric to
> compute histograms on the whole image, that may give you a much
> detailed metric map for the optimizer as well.
>
> metric->SetNumberOfSpatialSamples(<some value>);
>
> Hope that helps,
>
> Nicolás Gallego-Ortiz
> Université catholique de Louvain, Belgium
>
>
> 2014-07-03 1:29 GMT+02:00 Jerome Plumat <j.plumat at auckland.ac.nz
> <mailto:j.plumat at auckland.ac.nz>>:
>
> Hi everyone,
> I'm facing to a problem with 3D affine multi resolution registration.
> The pipeline registers two 3D MRI head and necks volumes
> (600x400x40). These are scans of the same patient at two times.
> Due to the protocol, volumes are similar expect for MRI artifacts.
> In fine the algorithm would be used to registered an atlas defined
> on the fixed image on moving image.
>
> The main elements of the pipeline are listed below. It's
> implemented with itk4 and is mostly based on
> MutliResImageRegistration2.cxx, I use the Mutual Information
> metric, linear interpolator, the optimizer is
> RegularStepGradientDescentOptimizer and CenteredTransformInitializer.
>
> The registration is performed without any error and the
> translation is correctly found but the algorithm can't provide any
> rotation and I can't figure why.
> The final results are
> versor X = [1.00007, 0.000485777, 1.5912e-05]
> versor Y = [-0.000126483, 0.999922, 1.34996e-05]
> versor Z = [2.35376e-05, -8.28481e-05, 0.999986]
> Translation X = 1.35204
> Translation Y = 0.647584
> Translation Z = -0.132708
> Iterations = 96
> Metric value = -0.738412
>
> Following
> http://www.cmake.org/pipermail/insight-users/2006-August/019025.html
> the final rotation is 0.0003 along axis [-0.155456, -0.0122989,
> -0.987766].
> Moreover, visual validation clearly highlight the fact the
> rotation is not sufficient.
>
> Hope any one has advices.
> Thanks in advances.
>
> Jerome.
>
>
> Part of code:
>
> ....
>
> // transform
> typedef itk::AffineTransform< double, Dimension> TransformType;
> typedef itk::CenteredTransformInitializer< TransformType,
> FixedImageType, MovingImageType > TransformInitializerType;
> // optimizer
> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
>
> // metric
> typedef itk::MattesMutualInformationImageToImageMetric<
> MovingImageType,
> MovingImageType > MetricType;
> // interpolator
> typedef itk:: LinearInterpolateImageFunction<
> MovingImageType,
> double > InterpolatorType;
>
> // registration
> typedef itk::MultiResolutionImageRegistrationMethod<
> FixedImageType,
> MovingImageType > RegistrationType;
>
>
> typedef itk::MultiResolutionPyramidImageFilter<
> FixedImageType,
> FixedImageType > FixedImagePyramidType;
> typedef itk::MultiResolutionPyramidImageFilter<
> MovingImageType,
> MovingImageType > MovingImagePyramidType;
>
>
> ...
>
> metric->SetNumberOfHistogramBins( 256 );
> metric->SetNumberOfSpatialSamples( 50000 );
> metric->ReinitializeSeed( 76926294 <tel:%28%2076926294> );
>
> ....
>
> // connect the inputs
> registration->SetFixedImage( fixedImage );
> registration->SetMovingImage( movingImage );
>
>
> // set the regions
> registration->SetFixedImageRegion(
> fixedImage->GetLargestPossibleRegion() );
>
> // configure the initializer
> initializer->SetTransform( transform );
> initializer->SetFixedImage( fixedImage );
> initializer->SetMovingImage( movingImage );
> initializer->MomentsOn();
> initializer->InitializeTransform();
>
> registration->SetInitialTransformParameters(
> transform->GetParameters() );
>
> // optimizer scales
> typedef OptimizerType::ScalesType OptimizerScalesType;
> OptimizerScalesType optimizerScales(
> transform->GetNumberOfParameters() );
> const double momentScale = 1;
> const double translationScale = 1.0 / 1e7;
>
> //number of dimensions: N*(N+1)
> for(unsigned int i=0; i<(Dimension*Dimension); i++ )
> {
> optimizerScales[i] = momentScale;
> }
>
> for( unsigned int i=(Dimension*Dimension);
> i<((Dimension+1)*Dimension); i++ )
> {
> optimizerScales[i] = translationScale;
> }
> optimizer->SetScales( optimizerScales );
>
>
> // performed registration
> ...
>
> //
> // deformed the moving image and save it
> //
> typedef itk::ResampleImageFilter<
> MovingImageType,
> FixedImageType > ResampleFilterType;
> typedef itk::CastImageFilter<
> MovingImageType,
> OutputImageType > CastFilterType;
> typedef itk::ImageFileWriter< OutputImageType > WriterType;
>
> ResampleFilterType::Pointer resampler = ResampleFilterType::New();
> WriterType::Pointer writer = WriterType::New();
> CastFilterType::Pointer caster = CastFilterType::New();
>
> TransformType::Pointer finalTransform = TransformType::New();
>
> finalTransform->SetParameters(
> registration->GetLastTransformParameters() );
> finalTransform->SetFixedParameters(
> transform->GetFixedParameters() );
>
> // set transform and initial moving image
> resampler->SetTransform( finalTransform );
> resampler->SetInput( movingImage );
>
> resampler->SetSize(
> fixedImage->GetLargestPossibleRegion().GetSize() );
> resampler->SetOutputOrigin( fixedImage->GetOrigin() );
> resampler->SetOutputSpacing( fixedImage->GetSpacing() );
> resampler->SetOutputDirection( fixedImage->GetDirection() );
> resampler->SetDefaultPixelValue( defaultOutputPixel );
>
> writer->SetFileName( argv[argvOutputImageFile] );
> caster->SetInput( resampler->GetOutput() );
> writer->SetInput( caster->GetOutput() );
>
> ...
>
> --
>
> Jerome
> -------
> School of Medical Sciences
> University of Auckland
>
> If I am not for myself, who will be for me? And if I am only for
> myself, then what an I? And if not now when? – Hillel HaZaken
>
> _____________________________________
> Powered by www.kitware.com <http://www.kitware.com>
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.php
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://public.kitware.com/mailman/listinfo/insight-users
> _______________________________________________
> Community mailing list
> Community at itk.org <mailto:Community at itk.org>
> http://public.kitware.com/mailman/listinfo/community
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20140704/98556415/attachment-0002.html>
-------------- next part --------------
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php
Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ
Follow this link to subscribe/unsubscribe:
http://public.kitware.com/mailman/listinfo/insight-users
More information about the Community
mailing list