[Insight-users] Change in patient orientation after registration
Vincent A. Magnotta
vincent-magnotta at uiowa.edu
Tue May 29 09:59:55 EDT 2007
Diego,
You must specify the output directions for the image in the
ResampleImageFilter. This can be done with the SetOutputDirection()
method.
Vince
On Tue, 2007-05-29 at 11:48 +0200, Diego Persano wrote:
> 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;
> }
>
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
--
Assistant Professor
Department of Radiology
0453-D JCP
200 Hawkins Drive
Iowa City, IA 52242
E-mail: vincent-magnotta at uiowa.edu
Phone: 319-356-8255
Fax: 319-353-6275
Website: http://www.radiology.uiowa.edu
More information about the Insight-users
mailing list