[Insight-users] Translation transform isn't initializing?

[David Welch]

Hello all,

I'm working on a simple (I thought) registration problem, very similar
to the example ImageRegistration13.cxx.

I have two tiff test images as the fixed image and moving images. Both
images are in RGB format, but one is in grayscale (fixed) while the
second is in color. They are 100x100 pixels and I set the spacing of
the images to an isotropic 0.33 mm.

Casting them to scalar pixels was done and the resulting images were
registered. There are a number of output statements and an observer
I've used to track the registration, but no matter what I try, I can't
get the moving image to translate remotely close to the fixed image. I
must be missing something essential, but I just can't see it. When I
output the moving image before registration, there is rotation, but no
translation. Using MomentsOn() for the TransformInitializer gives me a
value of [0.118895, -0.0253759], whereas GeometryOn() gives the expected
[0, 0].

The pertinent code is below. Any help is appreciated!

#include "itkImageRegistrationMethod.h"
#include "itkCenteredRigid2DTransform.h"
#include "itkCenteredTransformInitializer.h"
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkVectorResampleImageFilter.h"
...
typedef itk::ImageRegistrationMethod< IntensityImageType,
IntensityImageType > RegistrationType;
typedef itk::CenteredRigid2DTransform< double > TransformType;
float rotation2D = -30.0; // in degrees
int xOffset = 0;
int yOffset = 0;
typedef itk::CenteredTransformInitializer< TransformType,
IntensityImageType, IntensityImageType > InitialTransformType;
typedef itk::MattesMutualInformationImageToImageMetric<
IntensityImageType, IntensityImageType > MetricType;
const unsigned int numberOfBins = 256;
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
int numberOfIterations = 100;
typedef itk::LinearInterpolateImageFunction< IntensityImageType, double
 > InterpType;
...
MetricType::Pointer setupMetric( IntensityImageType::RegionType
fixedRegion )
{
/*############################################################################# 


Set up the metric
############################################################################*/ 

const unsigned int numberOfPixels = fixedRegion.GetNumberOfPixels();
const unsigned int numberOfSamples = static_cast< unsigned int >(
numberOfPixels * 0.5 );

MetricType::Pointer metric = MetricType::New();
metric->SetUseExplicitPDFDerivatives( 1 );
metric->SetNumberOfSpatialSamples( numberOfSamples );
metric->SetNumberOfHistogramBins( numberOfBins );

return metric;
}

TransformType::Pointer setupTransform( IntensityImageType::Pointer
fixed, IntensityImageType::Pointer moving )
{
/*############################################################################# 


Set up the transform
############################################################################*/ 

TransformType::Pointer transform = TransformType::New();
InitialTransformType::Pointer initializer =
InitialTransformType::New();
initializer->SetFixedImage( fixed );
initializer->SetMovingImage( moving );
initializer->SetTransform( transform );
initializer->MomentsOn();
initializer->InitializeTransform();
transform->SetAngleInDegrees( rotation2D );

return transform;
}

OptimizerType::Pointer setupOptimizer( TransformType::Pointer transform )
{
/*############################################################################# 


Set up optimizer and observer
############################################################################*/ 

OptimizerType::Pointer optimizer = OptimizerType::New();
optimizer->SetMaximumStepLength( 1.00 );
optimizer->SetMinimumStepLength( 0.001 );
optimizer->SetNumberOfIterations( numberOfIterations );
optimizer->SetRelaxationFactor( 0.8 );
optimizer->MinimizeOn();

typedef OptimizerType::ScalesType OptimizerScalesType;
OptimizerScalesType optimizerScales(
transform->GetNumberOfParameters() );

const double translationScale = 1.0 / 128.0;
const double centerScale = 1.0 / 100.0;
optimizerScales[0] = 1.0;
optimizerScales[1] = centerScale;
optimizerScales[2] = centerScale;
optimizerScales[3] = translationScale;
optimizerScales[4] = translationScale;

optimizer->SetScales( optimizerScales );

return optimizer;
}

int main( int argc, char * argv[] )
{
...
IntensityImageType::RegionType fixedImageRegion =
fixedImage->GetBufferedRegion();

// Set up the registration components
TransformType::Pointer transform = setupTransform(
fixedImage, movingImage );
OptimizerType::Pointer optimizer = setupOptimizer(
transform );
InterpType::Pointer interpolator = InterpType::New();

MetricType::Pointer metric = setupMetric(
fixedImageRegion );
metric->SetFixedImage( fixedImage );
metric->SetMovingImage( movingImage );
metric->SetInterpolator( interpolator );
metric->SetTransform( transform );

CommandIterationUpdate::Pointer observer =
CommandIterationUpdate::New();
optimizer->AddObserver( itk::IterationEvent(), observer );

RegistrationType::Pointer registration =
RegistrationType::New();
registration->SetMetric( metric );
registration->SetTransform( transform );
registration->SetInitialTransformParameters(
transform->GetParameters() );
registration->SetOptimizer( optimizer );
registration->SetInterpolator( interpolator);
registration->SetFixedImage( fixedImage );
registration->SetFixedImageRegion( fixedImageRegion );
registration->SetMovingImage( movingImage );

registration->Update();
registration->StartRegistration();
...
TransformType::Pointer finalTransform = TransformType::New();

finalTransform->SetParameters( finalParameters );
finalTransform->SetFixedParameters( transform->GetFixedParameters() );

RGBPixelType defaultRGB;
int rgb = 255;
defaultRGB[0] = rgb;
defaultRGB[1] = rgb;
defaultRGB[2] = rgb;

ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform( finalTransform );
resample->SetInput( cameraImage );
resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
resample->SetOutputOrigin( fixedImage->GetOrigin() );
resample->SetOutputSpacing( fixedImage->GetSpacing() );
resample->SetOutputDirection( fixedImage->GetDirection() );
resample->SetDefaultPixelValue( defaultRGB );

resample->Update();
WriterType::Pointer writerRGB = WriterType::New();
writerRGB->SetFileName( "registered.tif" ); //outputVolume );
writerRGB->SetInput( resample->GetOutput() );

writerRGB->Update();
...
return EXIT_SUCCESS;
}


Thank you,

Dave Welch
University of Iowa