[Insight-users] Change in patient orientation after registration

[Diego Persano]

Hi all,

I'm trying to registered two MRI images in NifTi format (*.nii) using a
simple translation (the code is below). The metric converges and the final
translation seems to be OK, but the patient orientation of the transformed
image is changed as can be seen printing the image information with the
method image->Print(std::cout):

Fixed image:

(...)
  Dimension: 3
  Index: [0, 0, 0]
  Size: [512, 512, 132]
Spacing: [0.625, 0.625, 1.3]
Origin: [155.785, 162.614, -66.2875]
Direction:
        -1 0 0
        0 -1 0
        0 0 1
(...)


Moving image:

(...)
  Dimension: 3
  Index: [0, 0, 0]
  Size: [512, 512, 132]
Spacing: [0.625, 0.625, 1.3]
Origin: [155.785, 162.614, -66.2875]
Direction:
        -1 0 0
        0 -1 0
        0 0 1
(...)


Transformed image (output of the reg.):

(...)
  Dimension: 3
  Index: [0, 0, 0]
  Size: [512, 512, 132]
Spacing: [0.625, 0.625, 1.3]
Origin: [155.785, 162.614, -66.2875]
Direction:
        1 0 0
        0 1 0
        0 0 1
(...)


Any idea about what is happening?

If you need further information about the images or something else, please
let me know.

Thank you.
Diego


/**** Code ****/

#include "itkImageRegistrationMethod.h"
#include "itkTranslationTransform.h"
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkImage.h"

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"


//Command observer used to monitor the evolution of the registration
process
#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( typeid( event ) != typeid( itk::IterationEvent ) )
      {
        return;
      }
    std::cout << optimizer->GetCurrentIteration() << "   ";
    std::cout << optimizer->GetValue() << "   ";
    std::cout << optimizer->GetCurrentPosition() << std::endl;
  }
};


int main( int argc, char *argv[] )
{

  //Verify the number of required parameters from the command line
  if( argc < 4 )
    {
      std::cerr << "Missing Parameters " << std::endl;
      std::cerr << "Usage: " << argv[0];
      std::cerr << " fixedImageFile movingImageFile ";
      std::cerr << "inputFile outputImagefile" << std::endl;
      return 1;
    }


  //Instantiation of the fixed and moving images
  const    unsigned int    Dimension = 3;
  typedef  signed short    PixelType;

  typedef itk::Image< PixelType, Dimension >  FixedImageType;
  typedef itk::Image< PixelType, Dimension >  MovingImageType;

  //Components required in the coregistration framework
  typedef itk::TranslationTransform< double, Dimension >                        
TransformType;
  typedef itk::RegularStepGradientDescentOptimizer                              
OptimizerType;
  typedef itk::LinearInterpolateImageFunction< MovingImageType, double >        
InterpolatorType;
  typedef itk::ImageRegistrationMethod< FixedImageType, MovingImageType >       
RegistrationType;
  typedef itk::MattesMutualInformationImageToImageMetric< FixedImageType,
MovingImageType > MetricType;

  //Creation of the coregistration components
  TransformType::Pointer      transform     = TransformType::New();
  OptimizerType::Pointer      optimizer     = OptimizerType::New();
  InterpolatorType::Pointer   interpolator  = InterpolatorType::New();
  RegistrationType::Pointer   registration  = RegistrationType::New();
  MetricType::Pointer         metric        = MetricType::New();

  registration->SetOptimizer(    optimizer    );
  registration->SetTransform(    transform    );
  registration->SetInterpolator( interpolator );
  registration->SetMetric(       metric       );


  //Read the fixed and moving images
  typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
  typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;

  FixedImageReaderType::Pointer  fixedImageReader  =
FixedImageReaderType::New();
  MovingImageReaderType::Pointer movingImageReader =
MovingImageReaderType::New();

  fixedImageReader->SetFileName(  argv[1] );
  movingImageReader->SetFileName( argv[2] );

  registration->SetFixedImage(  fixedImageReader->GetOutput() );
  registration->SetMovingImage( movingImageReader->GetOutput() );

  //Read the filter input file
  std::cout << std::endl;
  std::cout << "Reading filter input  " << argv[3] << std::endl;
  std::ifstream fin( argv[3] );

  if( !fin )
  {
    std::cout << "Cannot read" << argv[3] << std::endl;
    std::cout << std::endl;
    return 1;
  }

  //Read comments
  std::string buffer;
  getline( fin, buffer );
  while( buffer.at( buffer.find_first_not_of(" ") ) == '#' )
  {
    std::cout << buffer << std::endl;
    getline( fin, buffer );
  }

  //Put back the last read line: it's not any more a comment
  fin.seekg( fin.tellg() - std::istream::pos_type(buffer.size() + 1) );

  //Select a region of the fixed image as input to the metric computation

  //Coordinates of the region
  unsigned int startx, starty, startz;
  fin >> startx >> starty >> startz;
  std::cout << startx << " " << starty << " " << startz << std::endl;

  unsigned int sizex, sizey, sizez;
  fin >> sizex >> sizey >> sizez;
  std::cout << sizex << " " << sizey << " " << sizez << std::endl;

  FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
  fixedImageReader->Update();

  FixedImageType::IndexType start;

  start[0] = startx;   // first index on X
  start[1] = starty;   // first index on Y
  start[2] = startz;   // first index on Z

  FixedImageType::SizeType size;

  size[0] = sizex;   // size along X
  size[1] = sizey;   // size along Y
  size[2] = sizez;   // size along Z

  FixedImageType::RegionType fixedRegion;

  fixedRegion.SetSize(  size  );
  fixedRegion.SetIndex( start );

//   FixedImageType::RegionType fixedRegion =
fixedImage->GetBufferedRegion();

  registration->SetFixedImageRegion( fixedRegion );


  //Parameters required by the metric
  unsigned int numberOfBins;
  fin >> numberOfBins;
  std::cout << numberOfBins << std::endl;
  metric->SetNumberOfHistogramBins( numberOfBins );

  float spatialSamples;
  fin >> spatialSamples;
  std::cout << spatialSamples << std::endl;
  const unsigned int numberOfSamples = fixedRegion.GetNumberOfPixels() *
spatialSamples;
  metric->SetNumberOfSpatialSamples( numberOfSamples );

  //Initial parameters for the transform (x y z) translations
  typedef RegistrationType::ParametersType ParametersType;
  ParametersType initialParameters( transform->GetNumberOfParameters() );

  fin >> initialParameters[0] >> initialParameters[1] >>
initialParameters[2];
  std::cout << initialParameters[0] << " " << initialParameters[1] << " "
<< initialParameters[2] << std::endl;

  registration->SetInitialTransformParameters( initialParameters );


  //Parameters required by the optimizer
  float maxStepLength;
  fin >> maxStepLength;
  std::cout << maxStepLength << std::endl;
  optimizer->SetMaximumStepLength( maxStepLength );

  float minStepLength;
  fin >> minStepLength;
  std::cout << minStepLength << std::endl;
  optimizer->SetMinimumStepLength( minStepLength );

  unsigned int maxNumberOfIterations;
  fin >> maxNumberOfIterations;
  std::cout << maxNumberOfIterations << std::endl;
  optimizer->SetNumberOfIterations( maxNumberOfIterations );

  optimizer->SetGradientMagnitudeTolerance( 1e-7  );


  //Fix a voxel value for the region outside of the mapping
  unsigned int voxelValue;
  fin >> voxelValue;
  std::cout << voxelValue << std::endl;
  PixelType defaultVoxelValue = voxelValue;

  fin.close();
  std::cout << "Done" << std::endl;
  std::cout << std::endl;

  //Print size of the fixed image region
  FixedImageType::SizeType fixedImageSize = fixedRegion.GetSize();

  std::cout << "Size of fixed-image region in voxles: ";
  std::cout << fixedImageSize[0] << " x " << fixedImageSize[1] << " x " <<
fixedImageSize[2] << std::endl;

  std::cout << std::endl;
  std::cout << "Starting Registration" << std::endl;


  //Create the Command observer and register it with the optimizer
  CommandIterationUpdate::Pointer observer =
CommandIterationUpdate::New();
  optimizer->AddObserver( itk::IterationEvent(), observer );


  //Trigger the coregistration process
  try
    {
      registration->StartRegistration();
    }
  catch( itk::ExceptionObject & err )
    {
      std::cout << "ExceptionObject caught !" << std::endl;
      std::cout << err << std::endl;
      return -1;
    }


  //Get the last values of the coregistration process
  ParametersType finalParameters =
registration->GetLastTransformParameters();

  double TranslationAlongX = finalParameters[0];
  double TranslationAlongY = finalParameters[1];
  double TranslationAlongZ = finalParameters[2];

  unsigned int numberOfIterations = optimizer->GetCurrentIteration();

  double bestValue = optimizer->GetValue();


  //Print out results
  std::cout << std::endl;
  std::cout << "Result = " << std::endl;
  std::cout << " Translation X = " << TranslationAlongX  << std::endl;
  std::cout << " Translation Y = " << TranslationAlongY  << std::endl;
  std::cout << " Translation Z = " << TranslationAlongZ  << std::endl;
  std::cout << " Iterations    = " << numberOfIterations << std::endl;
  std::cout << " Metric value  = " << bestValue          << std::endl;
  std::cout << " Stop Condition  = " << optimizer->GetStopCondition() <<
std::endl;


  //Re-sampling of the moving image using the resulting transform
  typedef itk::ResampleImageFilter< MovingImageType, FixedImageType >
ResampleFilterType;

  TransformType::Pointer finalTransform = TransformType::New();

  finalTransform->SetParameters( finalParameters );

  ResampleFilterType::Pointer resample = ResampleFilterType::New();

  //Parameters required by the re-sampling filter
  resample->SetTransform(         finalTransform );
  resample->SetInput(             movingImageReader->GetOutput() );
  resample->SetSize(
fixedImage->GetLargestPossibleRegion().GetSize() );
  resample->SetOutputOrigin(      fixedImage->GetOrigin() );
  resample->SetOutputSpacing(     fixedImage->GetSpacing() );
  resample->SetDefaultPixelValue( defaultVoxelValue );


  //Instantiation of the output image and the cast and writer filter
  typedef signed short OutputPixelType;
  typedef itk::Image< OutputPixelType, Dimension >
OutputImageType;
  typedef itk::CastImageFilter< FixedImageType, OutputImageType >
CastFilterType;
  typedef itk::ImageFileWriter< OutputImageType >
WriterType;

  WriterType::Pointer      writer =  WriterType::New();
  CastFilterType::Pointer  caster =  CastFilterType::New();

  //Write output image file
  std::cout << std::endl;
  std::cout << "Writing output image  " << argv[4] << std::endl;

  writer->SetFileName( argv[4] );

  caster->SetInput( resample->GetOutput() );
  writer->SetInput( caster->GetOutput()   );
  writer->Update();

  std::cout << "Done" << std::endl;

  return 0;
}