[Insight-users] Registration using itk::OrientedImage
Blezek, Daniel J (GE, Research)
blezek at crd.ge.com
Fri Sep 7 08:09:01 EDT 2007
Hi Rupert,
This is a known, but perhaps not documented, problem with
OrientedImage and registration. Since OrientedImage was added after
much of ITK was written, almost none of the filters take
orientation/directional cosines into account. Most importantly, the
image derivatives assume that the images are axial, i.e. direction
matrix is the identity. This is a real bummer when you are trying to
register coronal MR to axial CT.
The work around is to not use OrientedImages, but initialize your
transform using the fixed orientation matrix multiplied by the inverse
of the moving orientation matrix (off the top of my head, it may be the
other way around). ITK will treat the Images is axial, but the initial
transform (in the "ITK coordinate system") will properly reflect the RAS
orientation of the images.
This is a pain! But unless someone refactors all the filters, you are
stuck here. It is an admitted design mistake to add orientation so late
in the development of ITK.
-dan
-----Original Message-----
From: insight-users-bounces+blezek=crd.ge.com at itk.org
[mailto:insight-users-bounces+blezek=crd.ge.com at itk.org] On Behalf Of
Rupert Brooks
Sent: Thursday, September 06, 2007 9:32 PM
To: insight-users at itk.org
Subject: [Insight-users] Registration using itk::OrientedImage
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
_______________________________________________
Insight-users mailing list
Insight-users at itk.org
http://www.itk.org/mailman/listinfo/insight-users
More information about the Insight-users
mailing list