[Insight-users] Analyze Image Orientation (Again)

[Geoff Topping]

Hello all,

As discussed here,

http://public.kitware.com/pipermail/insight-users/2007-February/021144.html

I am having problems with Analyze image orientation when trying to write 
images.  I'm trying to get the registration Hello World example working, but 
despite adding code similar to that suggested in the above-linked post, I 
continue to get the same error message:

terminate called after throwing an instance of 'itk::ExceptionObject'
  what():  
/home/gtopping/InsightToolkit-3.0.1/Code/IO/itkAnalyzeImageIO.cxx:1180:
itk::ERROR: AnalyzeImageIO(0x98a1690): ERROR: Analyze 7.5 File Format Only 
Allows RPI, PIR, and RIP Orientation
Aborted

The errors occur when writer->Update() is called, as verified by adding 
debug output to the example code (see below).

Can anyone suggests how to get this working?

Thanks,
Geoff

My current code:

// Software Guide : BeginCodeSnippet
#include "itkImageRegistrationMethod.h"
#include "itkTranslationTransform.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkImage.h"
// Software Guide : EndCodeSnippet


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

#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkSubtractImageFilter.h"
#include "itkOrientImageFilter.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 < 4 )
    {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " fixedImageFile  movingImageFile ";
    std::cerr << "outputImagefile [differenceImageAfter]";
    std::cerr << "[differenceImageBefore]" << std::endl;
    return 1;
    }


  const    unsigned int    Dimension = 3;
  typedef  float           PixelType;


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


  typedef itk::TranslationTransform< double, Dimension > TransformType;

  typedef itk::RegularStepGradientDescentOptimizer       OptimizerType;

  typedef itk::MeanSquaresImageToImageMetric<
                                    FixedImageType,
                                    MovingImageType >    MetricType;

  typedef itk:: LinearInterpolateImageFunction<
                                    MovingImageType,
                                    double          >    InterpolatorType;

  typedef itk::ImageRegistrationMethod<
                                    FixedImageType,
                                    MovingImageType >    RegistrationType;


  MetricType::Pointer         metric        = MetricType::New();
  TransformType::Pointer      transform     = TransformType::New();
  OptimizerType::Pointer      optimizer     = OptimizerType::New();
  InterpolatorType::Pointer   interpolator  = InterpolatorType::New();
  RegistrationType::Pointer   registration  = RegistrationType::New();


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


  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] );


  // Set Orientation of Images.  
http://public.kitware.com/pipermail/insight-users/2007-February/021144.html
  // fixed image...
  typedef itk::OrientImageFilter<FixedImageType, FixedImageType> 
OrientFixedImageFilterType;
  OrientFixedImageFilterType::Pointer fixedImageOrienter = 
OrientFixedImageFilterType::New();
  fixedImageOrienter->SetInput(fixedImageReader->GetOutput());
  fixedImageOrienter->SetUseImageDirection(true);
  
fixedImageOrienter->SetDesiredCoordinateOrientation(itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RPI);

  // moving image
  typedef itk::OrientImageFilter<MovingImageType, MovingImageType> 
OrientMovingImageFilterType;
  OrientMovingImageFilterType::Pointer movingImageOrienter = 
OrientMovingImageFilterType::New();
  movingImageOrienter->SetInput(movingImageReader->GetOutput());
  movingImageOrienter->SetUseImageDirection(true);
  
movingImageOrienter->SetDesiredCoordinateOrientation(itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RPI);



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


  fixedImageReader->Update();
  registration->SetFixedImageRegion(
                    fixedImageReader->GetOutput()->GetBufferedRegion() );


  typedef RegistrationType::ParametersType ParametersType;
  ParametersType initialParameters( transform->GetNumberOfParameters() );

  initialParameters[0] = 0.0;  // Initial offset in mm along X
  initialParameters[1] = 0.0;  // Initial offset in mm along Y

  registration->SetInitialTransformParameters( initialParameters );


  optimizer->SetMaximumStepLength( 4.00 );
  optimizer->SetMinimumStepLength( 0.01 );


  optimizer->SetNumberOfIterations( 2 );


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



  try
    {
    registration->Update();
    }
  catch( itk::ExceptionObject & err )
    {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << 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 << "Result = " << std::endl;
  std::cout << " Translation X = " << TranslationAlongX  << std::endl;
  std::cout << " Translation Y = " << TranslationAlongY  << std::endl;
  std::cout << " Iterations    = " << numberOfIterations << std::endl;
  std::cout << " Metric value  = " << bestValue          << std::endl;




  typedef itk::ResampleImageFilter<
                            MovingImageType,
                            FixedImageType >    ResampleFilterType;


  ResampleFilterType::Pointer resampler = ResampleFilterType::New();
  resampler->SetInput( movingImageReader->GetOutput() );


  resampler->SetTransform( registration->GetOutput()->Get() );

  std::cout << "set resampler transform" << std::endl;


  FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
  resampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
  resampler->SetOutputOrigin(  fixedImage->GetOrigin() );
  resampler->SetOutputSpacing( fixedImage->GetSpacing() );
  resampler->SetDefaultPixelValue( 100 );

  std::cout << "set resampler size, output origin, output spacing and 
default pixel value" << std::endl;


  typedef itk::Image< PixelType, Dimension > OutputImageType;
  typedef itk::CastImageFilter<
                        FixedImageType,
                        OutputImageType > CastFilterType;


  typedef itk::OrientImageFilter<OutputImageType, OutputImageType> 
OrientOutputImageFilterType;
  OrientOutputImageFilterType::Pointer outputImageOrienter = 
OrientOutputImageFilterType::New();

  typedef itk::ImageFileWriter< OutputImageType >  WriterType;


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

  writer->SetFileName( argv[3] );

  std::cout << "set writer file name" << std::endl;


  OutputImageType* outputImage = resampler->GetOutput();

  std::cout << "got output image from resampler" << std::endl;

  outputImageOrienter->SetInput(outputImage);
  outputImageOrienter->SetUseImageDirection(true);
  
outputImageOrienter->SetDesiredCoordinateOrientation(itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RPI);

  std::cout << "set output image orienter input, image direction and 
orientation" << std::endl;

  caster->SetInput(outputImage);

  std::cout << "set caster input" << std::endl;

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

  std::cout << "set writer input" << std::endl;

  writer->Update();

  std::cout << "updated writer" << std::endl;


  typedef itk::SubtractImageFilter<
                                  FixedImageType,
                                  FixedImageType,
                                  FixedImageType > DifferenceFilterType;

  DifferenceFilterType::Pointer difference = DifferenceFilterType::New();

  difference->SetInput1( fixedImageReader->GetOutput() );
  difference->SetInput2( resampler->GetOutput() );


  typedef itk::RescaleIntensityImageFilter<
                                  FixedImageType,
                                  OutputImageType >   RescalerType;

  RescalerType::Pointer intensityRescaler = RescalerType::New();

  intensityRescaler->SetInput( difference->GetOutput() );
  intensityRescaler->SetOutputMinimum(   0 );
  intensityRescaler->SetOutputMaximum( 255 );

  resampler->SetDefaultPixelValue( 1 );


  WriterType::Pointer writer2 = WriterType::New();
  writer2->SetInput( intensityRescaler->GetOutput() );


  if( argc > 4 )
    {
    writer2->SetFileName( argv[4] );
    writer2->Update();
    }


  TransformType::Pointer identityTransform = TransformType::New();
  identityTransform->SetIdentity();
  resampler->SetTransform( identityTransform );


  if( argc > 5 )
    {
    writer2->SetFileName( argv[5] );
    writer2->Update();
    }


  return 0;
}

_________________________________________________________________
Win a trip for four to a concert anywhere in the world! 
http://www.mobilelivetour.ca/