[ITK] [ITK-users] 3D Affine Registration - no rotation found

[Nicolas Gallego]

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>:

> 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 );
>
> ....
>
> // 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
>
> 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
> http://public.kitware.com/mailman/listinfo/community
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20140703/fc0aeaa0/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