[Insight-users] Analyze Image Orientation (Again)
Christian Marshall Rieck
rieck at stud.ntnu.no
Thu Mar 22 04:57:21 EST 2007
> 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?)
with the output file specified as *.hdr. (I use IBSR data where they are
called .hdr). Changing this to *.mhd makes it work.
Second: Same as above, but I rescale the intensities to 0..255 when I read
them from disk. The registration works, but writing the file does not.
Comment out the rescaling and it works.
Hopefully some of this may be of help.
Christian.
>
> 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/
>
> _______________________________________________
> 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