[Insight-users] mapping points from moving image to fixed image

John Drozd john.drozd at gmail.com
Thu Oct 15 13:01:21 EDT 2009


Hi Luis,

I made some small corrections in my email below *in red*.  Also this time I
remembered to carbon copy the insight users mailing list.

thanks,
john

On Thu, Oct 15, 2009 at 10:52 AM, John Drozd <john.drozd at gmail.com> wrote:

> Hello Luis,
>
> I used itkInverseDeformationFieldImageFilter.h to calculate an inverse
> deformation field from the moving subject to the fixed subject.   I tried to
> read in the deformation field in vtk format but reading the
> deformationField.vtk file took forever so I included calculating the inverse
> deformation field as the last part of the original deformable registration
> code that first calculated the deformation field (from the fixed subject to
> the moving subject). I was wondering about setting the index of the region.
> I arbitrarily set it to 0,0,0 as done in the testing code
> Insight/Testing/Code/BasicFilters/itkInverseDeformationFieldImageFilterTest.cxx
> Does the choice of the index make a difference?  Also it took almost 1 hour
> to calculate the inverse deformation field. I was wondering if there is a
> way to optimize the code so that it can run on the 8 processors on my
> computer. I noticed that the registration procedure used the 8 processors.
>
> To validate my code, I would like *to use the inverse deformation field* *that
> my code has calculated* to next translate physical coordinates of points
> from the moving subject to physical coordinates of points in the fixed
> subject and then back to the moving subject to see if I get back the
> original point (location).  Is there any example code that would show me how
> to do this? I am currently looking over ThinPlateSplineWarp.cxx which
> transforms some points.
>
> Eventually, I want to use the moving subject as an atlas and use the
> inverse deformation field to translate physical coordinates of points from
> ventricles in the moving subject atlas to physical coordinates of points of
> ventricles in the fixed subject, and then use these as seeds in the fixed
> subject in a region growing algorithm to segment the ventricles *and
> calculate the volumes* *of the ventricles* in the fixed subject. This way
> I can have an automated segmentation code.
>
> thank you
> john
>
> Below is my code
>
> #if defined(_MSC_VER)
> #pragma warning ( disable : 4786 )
> #endif
>
>
> /*
> ./DeformableRegistration8andinverse correctedsubject1.dcm
> correctedsubject2.dcm outputImage.mhd differenceOutput.mhd
> differenceBeforeRegistration.mhd deformationField.vtk
> inversedeformationField.vtk
> */
>
>
> // Software Guide : BeginLatex
> //
> // This example illustrates the use of the
> \doxygen{BSplineDeformableTransform}
> // class for performing registration of two $3D$ images and for the case of
> // multi-modality images. The image metric of choice in this case is the
> // \doxygen{MattesMutualInformationImageToImageMetric}.
> //
> // \index{itk::BSplineDeformableTransform}
> // \index{itk::BSplineDeformableTransform!DeformableRegistration}
> // \index{itk::LBFGSBOptimizer}
> //
> //
> // Software Guide : EndLatex
>
> #include "itkImageRegistrationMethod.h"
> #include "itkMattesMutualInformationImageToImageMetric.h"
> #include "itkLinearInterpolateImageFunction.h"
> #include "itkImage.h"
>
> #include "itkTimeProbesCollectorBase.h"
>
> //added per itkInverseDeformationFieldImageFilterTest.cxx
> #include "itkVector.h"
> #include "itkInverseDeformationFieldImageFilter.h"
> #include "itkImageRegionIteratorWithIndex.h"
> //#include "itkImageFileWriter.h"
> #include "itkFilterWatcher.h"
> //end of added code
>
> #ifdef ITK_USE_REVIEW
> #include "itkMemoryProbesCollectorBase.h"
> #define itkProbesCreate()  \
>   itk::TimeProbesCollectorBase chronometer; \
>   itk::MemoryProbesCollectorBase memorymeter
> #define itkProbesStart( text ) memorymeter.Start( text );
> chronometer.Start( text )
> #define itkProbesStop( text )  chronometer.Stop( text ); memorymeter.Stop(
> text  )
> #define itkProbesReport( stream )  chronometer.Report( stream );
> memorymeter.Report( stream  )
> #else
> #define itkProbesCreate()  \
>   itk::TimeProbesCollectorBase chronometer
> #define itkProbesStart( text ) chronometer.Start( text )
> #define itkProbesStop( text )  chronometer.Stop( text )
> #define itkProbesReport( stream )  chronometer.Report( stream )
> #endif
>
>
> //  Software Guide : BeginLatex
> //
> //  The following are the most relevant headers to this example.
> //
> //  \index{itk::BSplineDeformableTransform!header}
> //  \index{itk::LBFGSBOptimizer!header}
> //
> //  Software Guide : EndLatex
>
> // Software Guide : BeginCodeSnippet
> #include "itkBSplineDeformableTransform.h"
> #include "itkLBFGSBOptimizer.h"
> // Software Guide : EndCodeSnippet
>
> //  Software Guide : BeginLatex
> //
> //  The parameter space of the \code{BSplineDeformableTransform} is
> composed by
> //  the set of all the deformations associated with the nodes of the
> BSpline
> //  grid.  This large number of parameters makes possible to represent a
> wide
> //  variety of deformations, but it also has the price of requiring a
> //  significant amount of computation time.
> //
> //  \index{itk::BSplineDeformableTransform!header}
> //
> //  Software Guide : EndLatex
>
> #include "itkImageFileReader.h"
> #include "itkImageFileWriter.h"
>
> #include "itkResampleImageFilter.h"
> #include "itkCastImageFilter.h"
> #include "itkSquaredDifferenceImageFilter.h"
>
>
> #include "itkTransformFileReader.h"
>
> //added from ReadAtlasAndSubject.cxx
> #include "itkGDCMImageIO.h"
> //end of added code
>
> //  The following section of code implements a 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::LBFGSBOptimizer     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->GetCachedValue() << "   ";
>       std::cout << optimizer->GetInfinityNormOfProjectedGradient() <<
> 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 outputImagefile  ";
>     std::cerr << " [differenceOutputfile] [differenceBeforeRegistration] ";
>
>     //std::cerr << " [deformationField] ";
>     std::cerr << " [deformationField] " << " [inversedeformationField] ";
>
>     std::cerr << " [useExplicitPDFderivatives ] [useCachingBSplineWeights ]
> ";
>     std::cerr << " [filenameForFinalTransformParameters] ";
>     std::cerr << std::endl;
>     return EXIT_FAILURE;
>     }
>
>   const    unsigned int    ImageDimension = 3;
>   typedef  signed short    PixelType;
>
>   typedef itk::Image< PixelType, ImageDimension >  FixedImageType;
>   typedef itk::Image< PixelType, ImageDimension >  MovingImageType;
>
>
>   //  Software Guide : BeginLatex
>   //
>   //  We instantiate now the type of the \code{BSplineDeformableTransform}
> using
>   //  as template parameters the type for coordinates representation, the
>   //  dimension of the space, and the order of the BSpline.
>   //
>   //  \index{BSplineDeformableTransform!New}
>   //  \index{BSplineDeformableTransform!Instantiation}
>   //
>   //  Software Guide : EndLatex
>
>   // Software Guide : BeginCodeSnippet
>   const unsigned int SpaceDimension = ImageDimension;
>   const unsigned int SplineOrder = 3;
>   typedef double CoordinateRepType;
>
>   typedef itk::BSplineDeformableTransform<
>                             CoordinateRepType,
>                             SpaceDimension,
>                             SplineOrder >     TransformType;
>   // Software Guide : EndCodeSnippet
>
>
>   typedef itk::LBFGSBOptimizer       OptimizerType;
>
>
>   typedef itk::MattesMutualInformationImageToImageMetric<
>                                     FixedImageType,
>                                     MovingImageType >    MetricType;
>
>   typedef itk:: LinearInterpolateImageFunction<
>                                     MovingImageType,
>                                     double          >    InterpolatorType;
>
>   typedef itk::ImageRegistrationMethod<
>                                     FixedImageType,
>                                     MovingImageType >    RegistrationType;
>
>   MetricType::Pointer         metric        = MetricType::New();
>   OptimizerType::Pointer      optimizer     = OptimizerType::New();
>   InterpolatorType::Pointer   interpolator  = InterpolatorType::New();
>   RegistrationType::Pointer   registration  = RegistrationType::New();
>
>
>   registration->SetMetric(        metric        );
>   registration->SetOptimizer(     optimizer     );
>   registration->SetInterpolator(  interpolator  );
>
>
>   //  Software Guide : BeginLatex
>   //
>   //  The transform object is constructed below and passed to the
> registration
>   //  method.
>   //  \index{itk::RegistrationMethod!SetTransform()}
>   //
>   //  Software Guide : EndLatex
>
>   // Software Guide : BeginCodeSnippet
>   TransformType::Pointer  transform = TransformType::New();
>   registration->SetTransform( transform );
>   // Software Guide : EndCodeSnippet
>
>
>   //
>   //   In general, you must first solve an Affine registration between
>   //   the images before attempting to solve a deformable registration.
>   //   If you have solve an affine transform, it can be loaded into the
>   //   BSplineDeformableTransform as a "bulk" transform that will be
>   //   pre-composed with the deformation computed by the BSpline.
>   //   The following code loads one of such initial transforms if they
>   //   are available.
>   //
>   typedef itk::TransformFileReader        TransformReaderType;
>   typedef itk::AffineTransform<double, 3> AffineTransformType;
>
>   TransformReaderType::Pointer transformReader =
> TransformReaderType::New();
>
>   if( argc > 11 )
>     {
>     std::cout << "Loading Transform: " << argv[11] << std::endl;
>     transformReader->SetFileName( argv[11] );
>     transformReader->Update();
>
>     typedef TransformReaderType::TransformListType * TransformListType;
>     TransformListType transforms = transformReader->GetTransformList();
>     TransformReaderType::TransformListType::const_iterator tit =
> transforms->begin();
>     if( !strcmp((*tit)->GetNameOfClass(),"AffineTransform") )
>       {
>       AffineTransformType::Pointer affine_read =
>                   static_cast<AffineTransformType*>((*tit).GetPointer());
>       AffineTransformType::Pointer affine_transform =
>                   dynamic_cast< AffineTransformType * >(
> affine_read.GetPointer() );
>
>       if( affine_transform )
>         {
>         transform->SetBulkTransform( affine_transform );
>         }
>       }
>     else
>       {
>       std::cerr << "Bulk transform wasn't an affine transform." <<
> std::endl;
>       return EXIT_FAILURE;
>       }
>     }
>
>
>   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] );
>
>   //added from ReadAtlasAndSubject.cxx
>   typedef itk::GDCMImageIO           ImageIOTypefixed;
>   ImageIOTypefixed::Pointer gdcmImageIOfixed = ImageIOTypefixed::New();
>   fixedImageReader->SetImageIO( gdcmImageIOfixed );
>
>   typedef itk::GDCMImageIO           ImageIOTypemoving;
>   ImageIOTypemoving::Pointer gdcmImageIOmoving = ImageIOTypemoving::New();
>   movingImageReader->SetImageIO( gdcmImageIOmoving );
>   //end of added code
>
>   movingImageReader->Update();
>   fixedImageReader->Update();
>
>   FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();
>
>   //added
>   MovingImageType::ConstPointer movingImage =
> movingImageReader->GetOutput();
>
>   registration->SetFixedImage(  fixedImage   );
>   registration->SetMovingImage(   movingImageReader->GetOutput()   );
>
>   //fixedImageReader->Update();
>
>   FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
>
>  registration->SetFixedImageRegion( fixedRegion );
>
>   //  Software Guide : BeginLatex
>   //
>   //  Here we define the parameters of the BSplineDeformableTransform
> grid.  We
>   //  arbitrarily decide to use a grid with $5 \times 5$ nodes within the
> image.
>   //  The reader should note that the BSpline computation requires a
>   //  finite support region ( 1 grid node at the lower borders and 2
>   //  grid nodes at upper borders). Therefore in this example, we set
>   //  the grid size to be $8 \times 8$ and place the grid origin such that
>   //  grid node (1,1) coincides with the first pixel in the fixed image.
>   //
>   //  \index{BSplineDeformableTransform}
>   //
>   //  Software Guide : EndLatex
>
>   unsigned int numberOfGridNodesInOneDimension = 5;
>
>   if( argc > 10 )
>     {
>     numberOfGridNodesInOneDimension = atoi( argv[10] );
>     }
>
>   // Software Guide : BeginCodeSnippet
>   typedef TransformType::RegionType RegionType;
>   RegionType bsplineRegion;
>   RegionType::SizeType   gridSizeOnImage;
>   RegionType::SizeType   gridBorderSize;
>   RegionType::SizeType   totalGridSize;
>
>   gridSizeOnImage.Fill( numberOfGridNodesInOneDimension );
>   gridBorderSize.Fill( SplineOrder );    // Border for spline order = 3 ( 1
> lower, 2 upper )
>   totalGridSize = gridSizeOnImage + gridBorderSize;
>
>   bsplineRegion.SetSize( totalGridSize );
>   //  Software Guide : EndCodeSnippet
>
>
>   // Software Guide : BeginCodeSnippet
>   typedef TransformType::SpacingType SpacingType;
>   SpacingType spacing = fixedImage->GetSpacing();
>
>   typedef TransformType::OriginType OriginType;
>   OriginType origin = fixedImage->GetOrigin();;
>
>   FixedImageType::SizeType fixedImageSize = fixedRegion.GetSize();
>
>   for(unsigned int r=0; r<ImageDimension; r++)
>     {
>     spacing[r] *= static_cast<double>(fixedImageSize[r] - 1)  /
>                   static_cast<double>(gridSizeOnImage[r] - 1);
>     }
>
>   FixedImageType::DirectionType gridDirection = fixedImage->GetDirection();
>   SpacingType gridOriginOffset = gridDirection * spacing;
>
>   OriginType gridOrigin = origin - gridOriginOffset;
>
>   transform->SetGridSpacing( spacing );
>   transform->SetGridOrigin( gridOrigin );
>   transform->SetGridRegion( bsplineRegion );
>   transform->SetGridDirection( gridDirection );
>
>
>   typedef TransformType::ParametersType     ParametersType;
>
>   const unsigned int numberOfParameters =
>                transform->GetNumberOfParameters();
>
>   ParametersType parameters( numberOfParameters );
>
>   parameters.Fill( 0.0 );
>
>   transform->SetParameters( parameters );
>   //  Software Guide : EndCodeSnippet
>
>
>
>   //  Software Guide : BeginLatex
>   //
>   //  We now pass the parameters of the current transform as the initial
>   //  parameters to be used when the registration process starts.
>   //
>   //  Software Guide : EndLatex
>
>   // Software Guide : BeginCodeSnippet
>   registration->SetInitialTransformParameters( transform->GetParameters()
> );
>   // Software Guide : EndCodeSnippet
>
>
>   //  Software Guide : BeginLatex
>   //
>   //  Next we set the parameters of the LBFGSB Optimizer.
>   //
>   //  Software Guide : EndLatex
>
>
>   // Software Guide : BeginCodeSnippet
>   OptimizerType::BoundSelectionType boundSelect(
> transform->GetNumberOfParameters() );
>   OptimizerType::BoundValueType upperBound(
> transform->GetNumberOfParameters() );
>   OptimizerType::BoundValueType lowerBound(
> transform->GetNumberOfParameters() );
>
>   boundSelect.Fill( 0 );
>   upperBound.Fill( 0.0 );
>   lowerBound.Fill( 0.0 );
>
>   optimizer->SetBoundSelection( boundSelect );
>   optimizer->SetUpperBound( upperBound );
>   optimizer->SetLowerBound( lowerBound );
>
>   optimizer->SetCostFunctionConvergenceFactor( 1.e7 );
>   optimizer->SetProjectedGradientTolerance( 1e-6 );
>   optimizer->SetMaximumNumberOfIterations( 200 );
>   optimizer->SetMaximumNumberOfEvaluations( 30 );
>   optimizer->SetMaximumNumberOfCorrections( 5 );
>   // Software Guide : EndCodeSnippet
>
>   // Create the Command observer and register it with the optimizer.
>   //
>   CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
>   optimizer->AddObserver( itk::IterationEvent(), observer );
>
>
>   //  Software Guide : BeginLatex
>   //
>   //  Next we set the parameters of the Mattes Mutual Information Metric.
>   //
>   //  Software Guide : EndLatex
>
>   // Software Guide : BeginCodeSnippet
>   metric->SetNumberOfHistogramBins( 50 );
>
>   const unsigned int numberOfSamples =
>     static_cast<unsigned int>( fixedRegion.GetNumberOfPixels() * 20.0 /
> 100.0 );
>
>   metric->SetNumberOfSpatialSamples( numberOfSamples );
>   // Software Guide : EndCodeSnippet
>
>
>   //  Software Guide : BeginLatex
>   //
>   //  Given that the Mattes Mutual Information metric uses a random
> iterator in
>   //  order to collect the samples from the images, it is usually
> convenient to
>   //  initialize the seed of the random number generator.
>   //
>   //
> \index{itk::Mattes\-Mutual\-Information\-Image\-To\-Image\-Metric!ReinitializeSeed()}
>   //
>   //  Software Guide : EndLatex
>
>   // Software Guide : BeginCodeSnippet
>   metric->ReinitializeSeed( 76926294 );
>   // Software Guide : EndCodeSnippet
>
>   if( argc > 8 )
>     {
>     // Define whether to calculate the metric derivative by explicitly
>     // computing the derivatives of the joint PDF with respect to the
> Transform
>     // parameters, or doing it by progressively accumulating contributions
> from
>     // each bin in the joint PDF.
>     metric->SetUseExplicitPDFDerivatives( atoi( argv[8] ) );
>     }
>
>   if( argc > 9 )
>     {
>     // Define whether to cache the BSpline weights and indexes
> corresponding to
>     // each one of the samples used to compute the metric. Enabling caching
> will
>     // make the algorithm run faster but it will have a cost on the amount
> of memory
>     // that needs to be allocated. This option is only relevant when using
> the
>     // BSplineDeformableTransform.
>     metric->SetUseCachingOfBSplineWeights( atoi( argv[9] ) );
>     }
>
>
>   // Add time and memory probes
>   itkProbesCreate();
>
>   std::cout << std::endl << "Starting Registration" << std::endl;
>
>   try
>     {
>     itkProbesStart( "Registration" );
>     registration->StartRegistration();
>     itkProbesStop( "Registration" );
>     }
>   catch( itk::ExceptionObject & err )
>     {
>     std::cerr << "ExceptionObject caught !" << std::endl;
>     std::cerr << err << std::endl;
>     return EXIT_FAILURE;
>     }
>
>   OptimizerType::ParametersType finalParameters =
>                     registration->GetLastTransformParameters();
>
>
>   // Report the time and memory taken by the registration
>   itkProbesReport( std::cout );
>
>   // Software Guide : BeginCodeSnippet
>   transform->SetParameters( finalParameters );
>   // Software Guide : EndCodeSnippet
>
>
>   typedef itk::ResampleImageFilter<
>                             MovingImageType,
>                             FixedImageType >    ResampleFilterType;
>
>   ResampleFilterType::Pointer resample = ResampleFilterType::New();
>
>   resample->SetTransform( transform );
>   resample->SetInput( movingImageReader->GetOutput() );
>
>   resample->SetSize(    fixedImage->GetLargestPossibleRegion().GetSize() );
>   resample->SetOutputOrigin(  fixedImage->GetOrigin() );
>   resample->SetOutputSpacing( fixedImage->GetSpacing() );
>   resample->SetOutputDirection( fixedImage->GetDirection() );
>
>   // This value is set to zero in order to make easier to perform
>   // regression testing in this example. However, for didactic
>   // exercise it will be better to set it to a medium gray value
>   // such as 100 or 128.
>   //resample->SetDefaultPixelValue( 0 );
>   resample->SetDefaultPixelValue( 0 );
>
>   typedef  signed short  OutputPixelType;
>
>   typedef itk::Image< OutputPixelType, ImageDimension > OutputImageType;
>
>   typedef itk::CastImageFilter<
>                         FixedImageType,
>                         OutputImageType > CastFilterType;
>
>   typedef itk::ImageFileWriter< OutputImageType >  WriterType;
>
>
>   WriterType::Pointer      writer =  WriterType::New();
>   CastFilterType::Pointer  caster =  CastFilterType::New();
>
>   //added from ReadAtlasAndSubject.cxx
>   /*typedef itk::GDCMImageIO           ImageIOType;
>   ImageIOType::Pointer gdcmImageIO = ImageIOType::New();
>
>   writer->UseInputMetaDataDictionaryOff();
>   writer->SetImageIO( gdcmImageIO );*/
>   //end of added code
>
>   writer->SetFileName( argv[3] );
>
>   writer->UseInputMetaDataDictionaryOff();
>   //writer->SetImageIO( gdcmImageIOAtlas );
>
>   caster->SetInput( resample->GetOutput() );
>   writer->SetInput( caster->GetOutput()   );
>
>
>   try
>     {
>     writer->Update();
>     }
>   catch( itk::ExceptionObject & err )
>     {
>     std::cerr << "ExceptionObject caught !" << std::endl;
>     std::cerr << err << std::endl;
>     return EXIT_FAILURE;
>     }
>
>
>
>   typedef itk::SquaredDifferenceImageFilter<
>                                   FixedImageType,
>                                   FixedImageType,
>                                   OutputImageType > DifferenceFilterType;
>
>   DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
>
>   WriterType::Pointer writer2 = WriterType::New();
>
>   //added
>   writer2->UseInputMetaDataDictionaryOff();
>   //end of added code
>
>
>   writer2->SetInput( difference->GetOutput() );
>
>
>   // Compute the difference image between the
>   // fixed and resampled moving image.
>   if( argc > 4 )
>     {
>     difference->SetInput1( fixedImageReader->GetOutput() );
>     difference->SetInput2( resample->GetOutput() );
>     writer2->SetFileName( argv[4] );
>     try
>       {
>       writer2->Update();
>       }
>     catch( itk::ExceptionObject & err )
>       {
>       std::cerr << "ExceptionObject caught !" << std::endl;
>       std::cerr << err << std::endl;
>       return EXIT_FAILURE;
>       }
>     }
>
>
>   // Compute the difference image between the
>   // fixed and moving image before registration.
>   if( argc > 5 )
>     {
>     writer2->SetFileName( argv[5] );
>     difference->SetInput1( fixedImageReader->GetOutput() );
>     difference->SetInput2( movingImageReader->GetOutput() );
>     try
>       {
>       writer2->Update();
>       }
>     catch( itk::ExceptionObject & err )
>       {
>       std::cerr << "ExceptionObject caught !" << std::endl;
>       std::cerr << err << std::endl;
>       return EXIT_FAILURE;
>       }
>     }
>
>   // Generate the explicit deformation field and inverse deformation field
> resulting from
>   // the registration.
>   if( argc > 6 )
>     {
>
>     typedef itk::Vector< float, ImageDimension >  VectorType;
>     typedef itk::Image< VectorType, ImageDimension >  DeformationFieldType;
>
>     typedef itk::InverseDeformationFieldImageFilter<
>                                     DeformationFieldType,
>                                     DeformationFieldType
>                                              >  FilterType;
>
>     FilterType::Pointer filter = FilterType::New();
>
>     DeformationFieldType::Pointer field = DeformationFieldType::New();
>     field->SetRegions( fixedRegion );
>     field->SetOrigin( fixedImage->GetOrigin() );
>     field->SetSpacing( fixedImage->GetSpacing() );
>     field->SetDirection( fixedImage->GetDirection() );
>     field->Allocate();
>
>     typedef itk::ImageRegionIterator< DeformationFieldType > FieldIterator;
>     FieldIterator fi( field, fixedRegion );
>
>     fi.GoToBegin();
>
>     TransformType::InputPointType  fixedPoint;
>     TransformType::OutputPointType movingPoint;
>     DeformationFieldType::IndexType index;
>
>     VectorType displacement;
>
>     while( ! fi.IsAtEnd() )
>       {
>       index = fi.GetIndex();
>       field->TransformIndexToPhysicalPoint( index, fixedPoint );
>       movingPoint = transform->TransformPoint( fixedPoint );
>       displacement = movingPoint - fixedPoint;
>       fi.Set( displacement );
>       ++fi;
>       }
>
>     typedef itk::ImageFileWriter< DeformationFieldType >  FieldWriterType;
>     FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
>
>     fieldWriter->SetInput( field );
>
>     fieldWriter->SetFileName( argv[6] );
>     try
>       {
>       fieldWriter->Update();
>       }
>     catch( itk::ExceptionObject & excp )
>       {
>       std::cerr << "Exception thrown " << std::endl;
>       std::cerr << excp << std::endl;
>       return EXIT_FAILURE;
>       }
>
>
>   // Generate the explicit inverse deformation field resulting from
>   // the registration.
>
>   DeformationFieldType::SpacingType spacing;
>   //spacing.Fill( 1.0 );
>
>   DeformationFieldType::PointType origin;
>   //origin.Fill( 0.0 );
>
>   DeformationFieldType::RegionType     region;
>   DeformationFieldType::SizeType       size;
>   DeformationFieldType::IndexType      start;
>
>   /*size[0] = movingImage->GetLargestPossibleRegion().GetSize()[0];
>   size[1] = movingImage->GetLargestPossibleRegion().GetSize()[1];
>   size[2] = movingImage->GetLargestPossibleRegion().GetSize()[2];
>   */
>
>   size[0] = fixedImage->GetLargestPossibleRegion().GetSize()[0];
>   size[1] = fixedImage->GetLargestPossibleRegion().GetSize()[1];
>   size[2] = fixedImage->GetLargestPossibleRegion().GetSize()[2];
>
>
>   std::cout << "size[0] "<< size[0] << std::endl;
>   std::cout << "size[1] "<< size[1] << std::endl;
>   std::cout << "size[2] "<< size[2] << std::endl;
>
>   region.SetSize( size );
>
>   start[0] = 0;
>   start[1] = 0;
>   start[2] = 0;
>
>
>
>   region.SetIndex( start );
>
>   field->SetOrigin( field->GetOrigin() );
>   field->SetSpacing( field->GetSpacing() );
>
>   field->SetRegions( region );
>
>
>   field->Allocate();
>
>   VectorType pixelValue;
>
>   filter->SetOutputSpacing( field->GetSpacing() );
>
>   filter->SetOutputOrigin( field->GetOrigin() );
>
>   filter->SetSize( size );
>
>
>
>   filter->SetInput( field );
>
>
>   filter->SetSubsamplingFactor( 16 );
>
>   try
>     {
>     std::cout << "updating inverse deformation filter" << std::endl;
>     filter->UpdateLargestPossibleRegion();
>     }
>   catch( itk::ExceptionObject & excp )
>     {
>     std::cerr << "Exception thrown " << std::endl;
>     std::cerr << excp << std::endl;
>     }
>
>   std::cout << "finished updating inverse deformation filter" << std::endl;
>
>
>   // Write an image for regression testing
>   typedef itk::ImageFileWriter<  DeformationFieldType  > InvWriterType;
>
>   InvWriterType::Pointer writerinv = InvWriterType::New();
>
>   writerinv->SetInput (filter->GetOutput());
>   writerinv->SetFileName( argv[7] );
>
>   try
>     {
>     std::cout << "writing inverse deformation filter" << std::endl;
>     writerinv->Update();
>     }
>   catch( itk::ExceptionObject & excp )
>     {
>     std::cerr << "Exception thrown by writer" << std::endl;
>     std::cerr << excp << std::endl;
>     return EXIT_FAILURE;
>     }
>
>   std::cout << "finished writing inverse deformation filter" << std::endl;
>
>     }
>
>   // Optionally, save the transform parameters in a file
>   if( argc > 10 )
>     {
>     std::ofstream parametersFile;
>     parametersFile.open( argv[10] );
>     parametersFile << finalParameters << std::endl;
>     parametersFile.close();
>     }
>
>   return EXIT_SUCCESS;
> }
>
> #undef itkProbesCreate
> #undef itkProbesStart
> #undef itkProbesStop
> #undef itkProbesReport
>
>
>
>
> On Mon, Oct 12, 2009 at 12:29 PM, Luis Ibanez <luis.ibanez at kitware.com>wrote:
>
>> Hi John,
>>
>> 1) The Deformation Field resulting from a Deformable Registration
>>     process in ITK, contains the vectors that point from the physical
>>     coordinates of a point in the Fixed image reference system to
>>     its corresponding point in the Moving image reference system.
>>
>> 2)  You could combine that deformation field with the WarpImageFilter
>>     in order to map intensities from the Moving image reference system
>>     into the Fixed Image reference system. That is: to resample the
>>     Moving image into the image grid of the Fixed image.
>>
>> 3)  If you want to map points (not pixel values) from the Moving image
>>     to the Fixed image, then you need to invert the deformation field.
>>     Simply negating (multiplying by -1 ) the initial deformation field
>>     will NOT do the trick.
>>
>> 4)  You may want to look at the following two options for inverting
>>     deformation fields:
>>
>> 4.1)  Insight/Code/BasicFilters/itkInverseDeformationFieldImageFilter.h
>>
>> 4.2)
>>  Insight/Code/BasicFilters/itkIterativeInverseDeformationFieldImageFilter.h
>>
>>     The first method uses a KernelTransform ( such as ThinPlateSplines)
>>      as an intermediate helper for computing the inverse field.
>>
>>      The second method computes the inverse field iterative based
>>      on   F * R = I  (forward field composed with Reverse field should
>>      be equal to the Identity), by redistributing the residual errors of
>>      Error = F * R - I
>>
>>
>> 5) An alternative option can be found in
>>
>>     InsightApplications/InverseConsistentLandmarkRegistration
>>
>>     Based on the paper:
>>
>> Consistent Landmark and Intensity-Based Image Registration
>> by H. J. Johnson* and G. E. Christensen
>> IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 21, NO. 5, MAY 2002
>>
>>
>>
>> Please give it a try and let us know if you find
>> any problems.
>>
>>
>>
>>     Thanks
>>
>>
>>          Luis
>>
>>
>> --------------------------------------------------------
>> On Sat, Oct 10, 2009 at 4:12 AM, John Drozd <john.drozd at gmail.com> wrote:
>> > Hello,
>> >
>> > I used deformable registration using DeformableRegistration8.cxx to
>> register
>> > one brain subject from one person with another brain subject from
>> another
>> > person.  I now have a deformation field file that was outputted from
>> > deformableregistration8.cxx.  I was reading on the internet the
>> following
>> > link http://www.itk.org/pipermail/insight-users/2009-April/030004.html
>>  that
>> > describes how to map points from the fixed image to the moving image
>> using a
>> > deformation field.  If i want to map points from the moving image to the
>> > fixed image, can I use the same procedure but take the negative of the
>> > deformation field, and if so how can i negate the deformation field and
>> use
>> > it with ResampleImagefilter? My logic is that when you take the negative
>> of
>> > a vector, this reverses its direction.
>> >
>> > Thank you,
>> >
>> > john
>> >
>> >
>> >
>> >
>> > --
>> > John Drozd
>> > Postdoctoral Fellow
>> > Imaging Research Laboratories
>> > Robarts Research Institute
>> > Room 1256
>> > 100 Perth Drive
>> > London, Ontario, Canada
>> > N6A 5K8
>> > (519) 661-2111 ext. 25306
>> >
>>
>
>
>
> --
> John Drozd
> Postdoctoral Fellow
> Imaging Research Laboratories
> Robarts Research Institute
> Room 1256
> 100 Perth Drive
> London, Ontario, Canada
> N6A 5K8
> (519) 661-2111 ext. 25306
>



-- 
John Drozd
Postdoctoral Fellow
Imaging Research Laboratories
Robarts Research Institute
Room 1256
100 Perth Drive
London, Ontario, Canada
N6A 5K8
(519) 661-2111 ext. 25306
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20091015/972349de/attachment-0001.htm>


More information about the Insight-users mailing list