[Insight-users] Registration using itk::OrientedImage
Luis Ibanez
luis.ibanez at kitware.com
Fri Sep 7 08:33:45 EDT 2007
Hi Robert,
As Dan pointed out, this is a known bug.
You will find its entry as Bug # 5081
http://public.kitware.com/Bug/view.php?id=5081
There are a couple of suggested solutions...
The challenge here is how to fix the computation for
the case of the OrientedImage without penalizing
the computation time for the standard itk::Image.
Suggestions are welcome.
Thanks
Luis
-------------------------------------------
Blezek, Daniel J (GE, Research) wrote:
> 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
> _______________________________________________
> 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