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

[Jerome Plumat]

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