[Insight-users] Change in patient orientation after registration
Diego Persano
persano at isi.it
Wed May 30 09:51:21 EDT 2007
Vincent,
Thank you very much for your reply. That solved the problem!
I was not aware of such a method because I used to use the Analyze image
format.
Thanks again.
Diego
On Tue, 29 May 2007, Vincent A. Magnotta wrote:
> 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
>
More information about the Insight-users
mailing list