[Insight-users] Analyze Image Orientation (Again)

Geoff Topping g_topping at hotmail.com
Fri Mar 23 18:43:13 EST 2007


>Date: Thu, 22 Mar 2007 10:57:21 +0100 (MET)
>From: Christian Marshall Rieck <rieck at stud.ntnu.no>
>Subject: Re: [Insight-users] Analyze Image Orientation (Again)
>To: Geoff Topping <g_topping at hotmail.com>
>Cc: insight-users at itk.org
>Message-ID: ><Pine.LNX.4.63.0703221050020.5008 at gaupe.stud.ntnu.no>
>Content-Type: TEXT/PLAIN; charset=US-ASCII
>
>>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?
>
>I have come across this myself, in two different ways:
>Running the code below (it apperas to be examples/imageRegistration8.cxx?)

It was the Hello World registration example, ImageRegistration1.cxx, with 
most of the Latex comments removed, the image dimensions changed to 3D from 
2D, and some added debug output and a bit of other code attempting to work 
around the problem by the method suggested in previous discussions on this 
topic (setting the image orientation).

>with the output file specified as *.hdr. (I use IBSR data where they are 
>called .hdr). Changing this to *.mhd makes it work.

After updating to ITK 3.2.0, I've tried this, but get the same results as 
with other filename extensions:


WARNING: In 
/home/gtopping/InsightToolkit-3.2.0/Code/IO/itkAnalyzeImageIO.cxx, line 1182
AnalyzeImageIO (0xb1f640): ERROR: Analyze 7.5 File Format Only Allows RPI, 
PIR, and RIP Orientation

terminate called after throwing an instance of 'itk::ExceptionObject'
  what():  
/home/gtopping/InsightToolkit-3.2.0/Code/IO/itkAnalyzeImageIO.cxx:1372:
Error opening image data file for writing../201_all_on_177_all_as_AVW.mhd
Aborted


I get the same results regardless of whether I use an OrientImageFilter on 
the output image before writing, or if I used the output of the resampler 
directly.


Any other suggestions?

Thanks,
Geoff

My latest 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 "itkAnalyzeImageIO.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;
  typedef itk::AnalyzeImageIO ImageIOType;


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

  writer->SetImageIO(analyzeIO);

  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 /*resampler->GetOutput()*/ );

  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;

  return 0;
}

_________________________________________________________________
Your Space. Your Friends. Your Stories. Share your world with Windows Live 
Spaces. http://spaces.live.com/?mkt=en-ca



More information about the Insight-users mailing list