[Insight-users] Registration using itk::OrientedImage
Rupert Brooks
rupe.brooks at gmail.com
Thu Sep 6 21:31:45 EDT 2007
Hi,
I seem to be having a problem using the itk::OrientedImage class in
the registration framework. I don't know if this is a bug, or if I'm
using it wrong.
The problem is that image derivatives, seem to be calculated with
respect to the pixel orientations. That is, an image derivative of
[1, -1] means that the image increases as you move along the first
pixel axis, and decreases as you move along the second. However,
for the image registration update step, you need [dI/dx, dI/dy] - that
is the derivative with respect to spatial X and Y. If the
transformation matrix (direction cosines) are near the identity, the
difference is unimportant. But if the transformation matrix is
relatively far from the identity, this is a big problem.
I'm attaching an example program that shows this problem. When run,
it expects two arguments
OrientedImageExample rotate 2Dimagefilename
rotate is either 0 or 1
0 - in which case it leaves the direction cosine of the image alone
1 - in which case it sets the direction cosines to [0 1;-1 0] (90
degree rotation)
It then does a simple registration of the image to itself from a small
distance away, from the same starting position, transformed
accordingly. Before registering, it prints out the starting value and
derivative of the metric
Not surprisingly, the starting metric value is the same in each case,
but the derivative is also the same, (except with a sign change) which
seems like a problem to me - it should also be transformed i think.
The registration fails in the rotated case.
Rupert B.
Output, when run on the typical example file BrainProtonDensitySliceBorder20
rupert at marita ~/development/registration/testing
$ $BUILD/Release/OrientedImageExample 0 $W/BrainProtonDensitySliceBorder20.png
I'm leaving the images as read
Probing metric at initial position:
Value: 2149.26 derivative: [207.904, -211.837]
0 = 2149.26 : [0.198196, -2.1452]
1 = 806.405 : [-0.392464, 1.81095]
2 = 655.604 : [0.344304, -0.0483956]
3 = 47.2975 : [-0.65047, 0.0537061]
4 = 163.683 : [-0.151342, 0.0241907]
5 = 9.38078 : [0.345085, -0.0354765]
6 = 47.0304 : [0.0957721, -0.0169548]
7 = 3.81136 : [-0.152001, 0.0163413]
8 = 9.26022 : [-0.0273984, 0.00637869]
9 = 0.324928 : [0.095644, -0.0156568]
10 = 3.77619 : [0.0336195, -0.00796226]
11 = 0.489949 : [-0.0278694, 0.0032345]
12 = 0.315358 : [0.00326033, 0.000495335]
13 = 0.00439 : [-0.0122514, -0.00138299]
Registration done !
Number of iterations = 15
Parameters: [-0.0122514, -0.00138299]
Optimal metric value = 0.0607644
rupert at marita ~/development/registration/testing
$ $BUILD/Release/OrientedImageExample 1 $W/BrainProtonDensitySliceBorder20.png
I'm rotating the images 90 degrees
Probing metric at initial position:
Value: 2149.26 derivative: [-207.904, 211.837]
0 = 2149.26 : [7.8018, 0.145198]
1 = 2448.22 : [8.56937, -3.78047]
2 = 3004.69 : [5.37094, -6.18257]
3 = 2862.28 : [2.4309, -8.89479]
[...]
198 = 4354.61 : [-14.1694, -10.946]
199 = 4340.54 : [-14.1899, -10.9224]
Registration done !
Number of iterations = 200
Parameters: [-14.1899, -10.9224]
Optimal metric value = 4340.54
Code:
/*=====================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#include "itkImageRegistrationMethod.h"
#include "itkTranslationTransform.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkSubtractImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkOrientedImage.h"
#include "itkChangeInformationImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkCommand.h"
class CommandIterationUpdate : public itk::Command
{
public:
typedef CommandIterationUpdate Self;
typedef itk::Command Superclass;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro( Self );
protected:
CommandIterationUpdate() {};
public:
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
typedef const OptimizerType *OptimizerPointer;
void Execute(itk::Object *caller, const itk::EventObject & event)
{
Execute( (const itk::Object *)caller, event);
}
void Execute(const itk::Object * object, const itk::EventObject & event)
{
OptimizerPointer optimizer =
dynamic_cast< OptimizerPointer >( object );
if( ! itk::IterationEvent().CheckEvent( &event ) )
{
return;
}
std::cout << optimizer->GetCurrentIteration() << " = ";
std::cout << optimizer->GetValue() << " : ";
std::cout << optimizer->GetCurrentPosition() << std::endl;
}
};
int main( int argc, char *argv[] )
{
if( argc < 3 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " flip testImageFile ";
std::cerr << "outputImagefile [differenceImage]" << std::endl;
return 1;
}
bool flip=true;
const unsigned int Dimension = 2;
typedef unsigned short PixelType;
typedef itk::OrientedImage< PixelType, Dimension > FixedImageType;
typedef itk::OrientedImage< PixelType, Dimension > MovingImageType;
typedef itk::ChangeInformationImageFilter<FixedImageType>
FixedChangeFilterType;
typedef itk::ChangeInformationImageFilter<MovingImageType>
MovingChangeFilterType;
typedef itk::TranslationTransform< double, Dimension > TransformType;
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
typedef itk::LinearInterpolateImageFunction<
MovingImageType,
double > InterpolatorType;
typedef itk::ImageRegistrationMethod<
FixedImageType,
MovingImageType > RegistrationType;
typedef itk::MeanSquaresImageToImageMetric<
FixedImageType,
MovingImageType > MetricType;
TransformType::Pointer transform = TransformType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
InterpolatorType::Pointer interpolator = InterpolatorType::New();
RegistrationType::Pointer registration = RegistrationType::New();
registration->SetOptimizer( optimizer );
MetricType::Pointer metric = MetricType::New();
registration->SetMetric( metric );
typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
FixedImageReaderType::Pointer fixedImageReader =
FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader =
MovingImageReaderType::New();
flip=atoi(argv[1]);
if(flip) {
std::cout<<"I'm rotating the images 90 degrees"<<std::endl;
} else {
std::cout<<"I'm leaving the images as read"<<std::endl;
}
fixedImageReader->SetFileName( argv[2] );
movingImageReader->SetFileName( argv[2] );
FixedChangeFilterType::Pointer fixedChanger=FixedChangeFilterType::New();
MovingChangeFilterType::Pointer movingChanger=MovingChangeFilterType::New();
movingChanger->SetInput(movingImageReader->GetOutput());
fixedChanger->SetInput(fixedImageReader->GetOutput());
FixedChangeFilterType::DirectionType directions;
directions[0][0]=0;
directions[0][1]=1;
directions[1][0]=-1;
directions[1][1]=0;
movingChanger->SetOutputDirection(directions);
fixedChanger->SetOutputDirection(directions);
if(flip) {
fixedChanger->ChangeDirectionOn();
movingChanger->ChangeDirectionOn();
}
fixedChanger->Update(); // This is needed to make the BufferedRegion
below valid.
typedef RegistrationType::ParametersType ParametersType;
ParametersType initialParameters=transform->GetParameters();
if(flip) {
initialParameters[0]=5;
initialParameters[1]=3;
} else {
initialParameters[0]=3;
initialParameters[1]=-5;
}
metric->SetFixedImage( fixedChanger->GetOutput() );
metric->SetMovingImage( movingChanger->GetOutput() );
metric->SetFixedImageRegion(
fixedChanger->GetOutput()->GetBufferedRegion() );
metric->SetTransform( transform );
metric->SetInterpolator( interpolator );
MetricType::DerivativeType derivative(transform->GetNumberOfParameters());
double value;
try{
metric->Initialize();
metric->GetValueAndDerivative(initialParameters,value,derivative);
}
catch( itk::ExceptionObject & err )
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return -1;
}
std::cout<<"Probing metric at initial position:"<<std::endl;
std::cout<<"Value: "<<value<<" derivative: "<<derivative<<std::endl;
registration->SetInitialTransformParameters( initialParameters );
registration->SetFixedImage( fixedChanger->GetOutput() );
registration->SetMovingImage( movingChanger->GetOutput() );
registration->SetTransform( transform );
registration->SetInterpolator( interpolator );
registration->SetFixedImageRegion(
fixedImageReader->GetOutput()->GetBufferedRegion() );
optimizer->SetMaximumStepLength( 4.00 );
optimizer->SetMinimumStepLength( 0.01 );
optimizer->SetNumberOfIterations( 200 );
optimizer->MaximizeOff();
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver( itk::IterationEvent(), observer );
try
{
registration->StartRegistration();
}
catch( itk::ExceptionObject & err )
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return -1;
}
ParametersType finalParameters = registration->GetLastTransformParameters();
const double TranslationAlongX = finalParameters[0];
const double TranslationAlongY = finalParameters[1];
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
const double bestValue = optimizer->GetValue();
std::cout << "Registration done !" << std::endl;
std::cout << "Number of iterations = " << numberOfIterations << std::endl;
std::cout << "Parameters: "<<finalParameters<< std::endl;
std::cout << "Optimal metric value = " << bestValue << std::endl;
return 0;
}
--
--------------------------------------------------------------
Rupert Brooks
McGill Centre for Intelligent Machines (www.cim.mcgill.ca)
Ph.D Student, Electrical and Computer Engineering
http://www.cyberus.ca/~rbrooks
More information about the Insight-users
mailing list