registration/clitkAffineRegistrationGenericFilter.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://oncora1.lyon.fnclcc.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ======================================================================-====*/
  18
  19 #ifndef CLITKAFFINEREGISTRATIONGENERICFILTER_TXX
  20 #define CLITKAFFINEREGISTRATIONGENERICFILTER_TXX
  21
  22 namespace clitk
  23 {
  24
  25 //==============================================================================
  26 //Creating an observer class that allows us to monitor the registration
  27 //================================================================================
  28 class CommandIterationUpdate : public itk::Command
  29 {
  30 public:
  31   typedef  CommandIterationUpdate   Self;
  32   typedef  itk::Command             Superclass;
  33   typedef  itk::SmartPointer<Self>  Pointer;
  34   itkNewMacro( Self );
  35 protected:
  36   CommandIterationUpdate() {};
  37 public:
  38   typedef   clitk::GenericOptimizer<args_info_clitkAffineRegistration>     OptimizerType;
  39   typedef   const OptimizerType   *           OptimizerPointer;
  40
  41   // Set the generic optimizer
  42   void SetOptimizer(OptimizerPointer o) {
  43     m_Optimizer=o;
  44   }
  45
  46   // Execute
  47   void Execute(itk::Object *caller, const itk::EventObject & event) {
  48     Execute( (const itk::Object *)caller, event);
  49   }
  50
  51   void Execute(const itk::Object * object, const itk::EventObject & event) {
  52     if( !(itk::IterationEvent().CheckEvent( &event )) ) {
  53       return;
  54     }
  55
  56     m_Optimizer->OutputIterationInfo();
  57   }
  58
  59   OptimizerPointer m_Optimizer;
  60 };
  61
  62
  63 //==============================================================================
  64 //Creating an observer class that allows us to change parameters at subsequent levels
  65 //==============================================================================
  66 template <typename TRegistration, class args_info_clitkAffineRegistration>
  67 class RegistrationInterfaceCommand : public itk::Command
  68 {
  69 public:
  70   typedef  RegistrationInterfaceCommand   Self;
  71   typedef  itk::Command             Superclass;
  72   typedef itk::SmartPointer<Self>  Pointer;
  73   itkNewMacro( Self );
  74 protected:
  75   RegistrationInterfaceCommand() {};
  76 public:
  77
  78   // Registration
  79   typedef   TRegistration                              RegistrationType;
  80   typedef   RegistrationType *                         RegistrationPointer;
  81
  82   // Metric
  83   typedef typename RegistrationType::FixedImageType    InternalImageType;
  84   typedef clitk::GenericMetric<args_info_clitkAffineRegistration, InternalImageType, InternalImageType> GenericMetricType;
  85
  86   // Two arguments are passed to the Execute() method: the first
  87   // is the pointer to the object which invoked the event and the
  88   // second is the event that was invoked.
  89   void Execute(itk::Object * object, const itk::EventObject & event) {
  90     if( !(itk::IterationEvent().CheckEvent( &event )) ) {
  91       return;
  92     }
  93
  94     // Get the levels
  95     RegistrationPointer registration = dynamic_cast<RegistrationPointer>( object );
  96     unsigned int numberOfLevels=registration->GetNumberOfLevels();
  97     unsigned int currentLevel=registration->GetCurrentLevel()+1;
  98
  99     // Output the levels
 100     std::cout<<std::endl;
 101     std::cout<<"========================================"<<std::endl;
 102     std::cout<<"Starting resolution level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
 103     std::cout<<"========================================"<<std::endl;
 104     std::cout<<std::endl;
 105
 106     // Higher level?
 107     if (currentLevel>1) {
 108       // Reinitialize the metric (!= number of samples)
 109       typename GenericMetricType::Pointer genericMetric= GenericMetricType::New();
 110       genericMetric->SetArgsInfo(m_ArgsInfo);
 111       genericMetric->SetFixedImage(registration->GetFixedImagePyramid()->GetOutput(registration->GetCurrentLevel()));
 112       if (m_ArgsInfo.referenceMask_given)  genericMetric->SetFixedImageMask(registration->GetMetric()->GetFixedImageMask());
 113       typedef itk::ImageToImageMetric< InternalImageType, InternalImageType >  MetricType;
 114       typename  MetricType::Pointer metric=genericMetric->GetMetricPointer();
 115       registration->SetMetric(metric);
 116     }
 117   }
 118
 119   void Execute(const itk::Object * , const itk::EventObject & ) {
 120     return;
 121   }
 122
 123   void SetArgsInfo(args_info_clitkAffineRegistration a) {
 124     m_ArgsInfo=a;
 125   }
 126   args_info_clitkAffineRegistration m_ArgsInfo;
 127 };
 128
 129
 130 //==============================================================================
 131 // Update with the number of dimensions
 132 //==============================================================================
 133 template <unsigned int Dimension>
 134 void AffineRegistrationGenericFilter::UpdateWithDim(std::string PixelType)
 135 {
 136   if (m_Verbose) std::cout  << "Images were detected to be "<< Dimension << "D and " << PixelType << "..." << std::endl;
 137   if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
 138   UpdateWithDimAndPixelType<Dimension, float>();
 139   // }
 140
 141 }
 142
 143
 144 //==============================================================================
 145 // Update with the number of dimensions and pixeltype
 146 //==============================================================================
 147 template <unsigned int Dimension, class PixelType>
 148 void
 149 AffineRegistrationGenericFilter::UpdateWithDimAndPixelType()
 150 {
 151   //=============================================================================
 152   //Input
 153   //=============================================================================
 154   bool threadsGiven=m_ArgsInfo.threads_given;
 155   int threads=m_ArgsInfo.threads_arg;
 156
 157   //Coordinate Representation
 158   typedef double TCoordRep;
 159
 160   //The pixeltype of the fixed image will be used for output
 161   typedef itk::Image< PixelType, Dimension > FixedImageType;
 162
 163   //Whatever the pixel type, internally we work with an image represented in float
 164   typedef   PixelType     InternalPixelType;
 165   typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
 166
 167   //Read in the reference/fixed image
 168   typedef itk::ImageFileReader< InternalImageType > ReaderType;
 169   typename ReaderType::Pointer  fixedImageReader  = ReaderType::New();
 170   fixedImageReader->SetFileName( m_ArgsInfo.reference_arg);
 171
 172
 173   //Read in the object/moving image
 174   typename ReaderType::Pointer movingImageReader = ReaderType::New();
 175   movingImageReader->SetFileName( m_ArgsInfo.target_arg );
 176   if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
 177   fixedImageReader->Update();
 178   movingImageReader->Update();
 179
 180   if (m_Verbose) std::cout  << "Reading images... " << std::endl;
 181
 182   //we connect pointers to these internal images
 183   typename InternalImageType::Pointer fixedImage= fixedImageReader->GetOutput();
 184   typename InternalImageType::Pointer movingImage= movingImageReader->GetOutput();
 185
 186   //We keep the images used for input for possible output
 187   typename InternalImageType::Pointer inputFixedImage= fixedImageReader->GetOutput();
 188   typename InternalImageType::Pointer inputMovingImage= movingImageReader->GetOutput();
 189
 190
 191   //============================================================================
 192   // Preprocessing
 193   //============================================================================
 194
 195   //If given, the intensities of both images are first normalized to a zero mean and SD of 1
 196   // (usefull for MI, not necessary for Mattes' MI but performed anyway for the ouput)
 197   if ( m_ArgsInfo.normalize_flag ) {
 198     typedef itk::NormalizeImageFilter< InternalImageType,InternalImageType >  NormalizeFilterType;
 199
 200     typename  NormalizeFilterType::Pointer  fixedNormalizeFilter = NormalizeFilterType::New();
 201     typename  NormalizeFilterType::Pointer  movingNormalizeFilter = NormalizeFilterType::New();
 202
 203     fixedNormalizeFilter->SetInput( fixedImage );
 204     movingNormalizeFilter->SetInput( movingImage );
 205
 206     fixedNormalizeFilter->Update();
 207     movingNormalizeFilter->Update();
 208
 209     //We keep the images used for input for possible output
 210     inputFixedImage= fixedNormalizeFilter->GetOutput();
 211     inputMovingImage= movingNormalizeFilter->GetOutput();
 212
 213     //the pointers are reconnected for further output
 214     fixedImage=fixedNormalizeFilter->GetOutput();
 215     movingImage=movingNormalizeFilter->GetOutput();
 216
 217     if (m_Verbose)  std::cout <<  "Normalizing image intensities to zero mean and SD of 1..." << std::endl;
 218   }
 219
 220
 221   //If given, the images are blurred before processing
 222   if ( m_ArgsInfo.blur_arg!= 0) {
 223     typedef itk::DiscreteGaussianImageFilter<InternalImageType,InternalImageType> GaussianFilterType;
 224     typename GaussianFilterType::Pointer fixedSmoother = GaussianFilterType::New();
 225     typename GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();
 226     fixedSmoother->SetVariance( m_ArgsInfo.blur_arg );
 227     movingSmoother->SetVariance(m_ArgsInfo.blur_arg );
 228
 229     fixedSmoother->SetInput( fixedImage );
 230     movingSmoother->SetInput( movingImage );
 231
 232     fixedSmoother->Update();
 233     movingSmoother->Update();
 234
 235     fixedImage=fixedSmoother->GetOutput();
 236     movingImage=movingSmoother->GetOutput();
 237
 238     if (m_Verbose)  std::cout <<  "Blurring images with a Gaussian with standard deviation of " << m_ArgsInfo.blur_arg <<"..." << std::endl;
 239   }
 240
 241
 242   //============================================================================
 243   // Setting up the moving image in a reference system
 244   //============================================================================
 245   const itk::Vector<double, Dimension> movingResolution = movingImage->GetSpacing();
 246   typename InternalImageType::RegionType movingRegion = movingImage->GetLargestPossibleRegion();
 247   typename InternalImageType::RegionType::SizeType  movingSize = movingRegion.GetSize();
 248
 249   // Print the parameters of the moving image
 250   if (m_Verbose) {
 251     std::cout << "Object or Moving image:"<<std::endl;
 252     std::cout << "Size: " << movingSize[0] << ", " << movingSize[1];
 253     if (Dimension==3) std::cout<<", " << movingSize[2];
 254     std::cout << std::endl;
 255
 256     std::cout<< "Resolution: "<< movingResolution[0] << ", " << movingResolution[1];
 257     if (Dimension==3) std::cout<< ", " << movingResolution[2];
 258     std::cout << std::endl;
 259   }
 260
 261
 262   //============================================================================
 263   // Setting up the fixed image in a reference system
 264   //============================================================================
 265   const itk::Vector<double, Dimension> fixedResolution = fixedImage->GetSpacing();
 266   typename InternalImageType::RegionType fixedRegion = fixedImage->GetLargestPossibleRegion();
 267   typename InternalImageType::RegionType::SizeType fixedSize = fixedRegion.GetSize();
 268
 269   // Print the parameters of the moving image and the transform
 270   if (m_Verbose) {
 271     std::cout << "Target or Moving image:"<<std::endl;
 272     std::cout << "Size: " << fixedSize[0] << ", " << fixedSize[1];
 273     if (Dimension==3) std::cout<<", " << fixedSize[2];
 274     std::cout << std::endl;
 275
 276     std::cout<< "Resolution: "<< fixedResolution[0] << ", " << fixedResolution[1];
 277     if (Dimension==3) std::cout<< ", " << fixedResolution[2];
 278     std::cout << std::endl;
 279   }
 280
 281
 282
 283   //===========================================================================
 284   // If given, we connect a mask to reference or target
 285   //============================================================================
 286   typedef itk::ImageMaskSpatialObject<  Dimension >   MaskType;
 287   typename MaskType::Pointer  fixedMask=NULL;
 288   if (m_ArgsInfo.referenceMask_given) {
 289     fixedMask= MaskType::New();
 290     typedef itk::Image< unsigned char, Dimension >   ImageMaskType;
 291     typedef itk::ImageFileReader< ImageMaskType >    MaskReaderType;
 292     typename MaskReaderType::Pointer  maskReader = MaskReaderType::New();
 293     maskReader->SetFileName(m_ArgsInfo.referenceMask_arg);
 294     try {
 295       maskReader->Update();
 296     } catch( itk::ExceptionObject & err ) {
 297       std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
 298       std::cerr << err << std::endl;
 299       return;
 300     }
 301     if (m_Verbose)std::cout <<"Reference image mask was read..." <<std::endl;
 302
 303     // Set the image to the spatialObject
 304     fixedMask->SetImage( maskReader->GetOutput() );
 305   }
 306
 307   typedef itk::ImageMaskSpatialObject<  Dimension >   MaskType;
 308   typename MaskType::Pointer  movingMask=NULL;
 309   if (m_ArgsInfo.targetMask_given) {
 310     movingMask= MaskType::New();
 311     typedef itk::Image< unsigned char, Dimension >   ImageMaskType;
 312     typedef itk::ImageFileReader< ImageMaskType >    MaskReaderType;
 313     typename MaskReaderType::Pointer  maskReader = MaskReaderType::New();
 314     maskReader->SetFileName(m_ArgsInfo.targetMask_arg);
 315     try {
 316       maskReader->Update();
 317     } catch( itk::ExceptionObject & err ) {
 318       std::cerr << "ExceptionObject caught !" << std::endl;
 319       std::cerr << err << std::endl;
 320     }
 321     if (m_Verbose)std::cout <<"Target image mask was read..." <<std::endl;
 322
 323     movingMask->SetImage( maskReader->GetOutput() );
 324   }
 325
 326
 327   //============================================================================
 328   // The image pyramids
 329   //============================================================================
 330   typedef itk::RecursiveMultiResolutionPyramidImageFilter<InternalImageType,InternalImageType >  FixedImagePyramidType;
 331   typedef itk::RecursiveMultiResolutionPyramidImageFilter<InternalImageType,InternalImageType >  MovingImagePyramidType;
 332   typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
 333   typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
 334   fixedImagePyramid->SetUseShrinkImageFilter(false);
 335   fixedImagePyramid->SetInput(fixedImage);
 336   fixedImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
 337   movingImagePyramid->SetUseShrinkImageFilter(false);
 338   movingImagePyramid->SetInput(movingImage);
 339   movingImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
 340   if (m_Verbose) std::cout<<"Creating the image pyramid..."<<std::endl;
 341   fixedImagePyramid->Update();
 342   movingImagePyramid->Update();
 343
 344
 345
 346   //============================================================================
 347   // We retrieve the type of metric from the command line
 348   //============================================================================
 349   typedef clitk::GenericMetric<args_info_clitkAffineRegistration, InternalImageType, InternalImageType> GenericMetricType;
 350   typename GenericMetricType::Pointer genericMetric=GenericMetricType::New();
 351   genericMetric->SetArgsInfo(m_ArgsInfo);
 352   genericMetric->SetFixedImage(fixedImagePyramid->GetOutput(0));
 353   if (fixedMask) genericMetric->SetFixedImageMask(fixedMask);
 354   typedef itk::ImageToImageMetric< InternalImageType, InternalImageType >  MetricType;
 355   typename  MetricType::Pointer metric=genericMetric->GetMetricPointer();
 356   if (movingMask) metric->SetMovingImageMask(movingMask);
 357
 358 #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
 359   if (threadsGiven) metric->SetNumberOfThreads( threads );
 360 #else
 361   if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
 362 #endif
 363
 364   //============================================================================
 365   // Initialize using image moments.
 366   //============================================================================
 367   if (m_ArgsInfo.moment_flag) {
 368     typedef itk::ImageMomentsCalculator< InternalImageType > CalculatorType;
 369     typename CalculatorType::Pointer fixedCalculator= CalculatorType::New();
 370
 371     typename InternalImageType::Pointer fixedThresh;
 372     if (m_ArgsInfo.intThreshold_given) {
 373       typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
 374       typename ThresholdImageFilterType::Pointer thresholder = ThresholdImageFilterType::New();
 375       thresholder->SetInput(fixedImage);
 376       thresholder->SetLower(m_ArgsInfo.intThreshold_arg);
 377       thresholder->Update();
 378       fixedThresh=thresholder->GetOutput();
 379     } else fixedThresh=fixedImage;
 380
 381     fixedCalculator->SetImage(fixedThresh);
 382     fixedCalculator->Compute();
 383     Vector<double, Dimension> fixedCenter=fixedCalculator->GetCenterOfGravity();
 384     if(m_Verbose)std::cout<<"The fixed center of gravity is "<<fixedCenter<<"..."<<std::endl;
 385
 386     typedef itk::ImageMomentsCalculator< InternalImageType > CalculatorType;
 387     typename CalculatorType::Pointer movingCalculator= CalculatorType::New();
 388
 389     typename InternalImageType::Pointer movingThresh;
 390     if (m_ArgsInfo.intThreshold_given) {
 391       typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
 392       typename ThresholdImageFilterType::Pointer thresholder = ThresholdImageFilterType::New();
 393       thresholder->SetInput(movingImage);
 394       thresholder->SetLower(m_ArgsInfo.intThreshold_arg);
 395       thresholder->Update();
 396       movingThresh=thresholder->GetOutput();
 397     } else movingThresh=movingImage;
 398
 399     movingCalculator->SetImage(movingThresh);
 400     movingCalculator->Compute();
 401     Vector<double, Dimension> movingCenter=movingCalculator->GetCenterOfGravity();
 402     if(m_Verbose)std::cout<<"The moving center of gravity is "<<movingCenter<<"..."<<std::endl;
 403
 404     Vector<double, Dimension> shift= movingCenter-fixedCenter;
 405     if(m_Verbose)std::cout<<"The initial shift applied is "<<shift<<"..."<<std::endl;
 406
 407     m_ArgsInfo.transX_arg= shift [0];
 408     m_ArgsInfo.transY_arg= shift [1];
 409     if (Dimension==3) m_ArgsInfo.transZ_arg=shift [2];
 410   }
 411
 412   //============================================================================
 413   // Transform
 414   //============================================================================
 415   typedef clitk::GenericAffineTransform<args_info_clitkAffineRegistration, TCoordRep, Dimension > GenericAffineTransformType;
 416   typename GenericAffineTransformType::Pointer genericAffineTransform = GenericAffineTransformType::New();
 417   genericAffineTransform->SetArgsInfo(m_ArgsInfo);
 418   typedef itk::Transform< double, Dimension, Dimension > TransformType;
 419   typename TransformType::Pointer transform = genericAffineTransform->GetTransform();
 420
 421
 422   //=======================================================
 423   // Interpolator
 424   //=======================================================
 425   typedef clitk::GenericInterpolator<args_info_clitkAffineRegistration, InternalImageType,TCoordRep > GenericInterpolatorType;
 426   typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
 427   genericInterpolator->SetArgsInfo(m_ArgsInfo);
 428   typedef itk::InterpolateImageFunction< InternalImageType, TCoordRep >  InterpolatorType;
 429   typename  InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
 430
 431
 432   //============================================================================
 433   // Optimizer
 434   //============================================================================
 435   typedef clitk::GenericOptimizer<args_info_clitkAffineRegistration> GenericOptimizerType;
 436   unsigned int nParam = transform->GetNumberOfParameters();
 437   GenericOptimizerType::Pointer genericOptimizer=GenericOptimizerType::New();
 438   genericOptimizer->SetArgsInfo(m_ArgsInfo);
 439   genericOptimizer->SetOutputIteration(m_Verbose);
 440   genericOptimizer->SetOutputPosition(m_Verbose);
 441   genericOptimizer->SetOutputValue(m_Verbose);
 442   genericOptimizer->SetOutputGradient(m_ArgsInfo.gradient_flag);
 443   genericOptimizer->SetMaximize(genericMetric->GetMaximize());
 444   genericOptimizer->SetNumberOfParameters(nParam);
 445   typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
 446   OptimizerType::Pointer optimizer=genericOptimizer->GetOptimizerPointer();
 447
 448   // Scales
 449   itk::Optimizer::ScalesType scales( nParam );
 450   for (unsigned int i=nParam-Dimension; i<nParam; i++) //Translations
 451     scales[i] = m_ArgsInfo.tWeight_arg;
 452   for (unsigned int i=0; i<nParam-Dimension; i++)      //Rest
 453     scales[i] = m_ArgsInfo.rWeight_arg*180./M_PI;
 454   optimizer->SetScales(scales);
 455
 456   //============================================================================
 457   // Multiresolution registration
 458   //============================================================================
 459   typedef itk::MultiResolutionImageRegistrationMethod< InternalImageType,InternalImageType >  RegistrationType;
 460   typename  RegistrationType::Pointer registration = RegistrationType::New();
 461   registration->SetFixedImage( fixedImage  );
 462   registration->SetFixedImageRegion(fixedImage->GetLargestPossibleRegion());
 463   registration->SetMovingImage(  movingImage );
 464   registration->SetFixedImagePyramid( fixedImagePyramid );
 465   registration->SetMovingImagePyramid( movingImagePyramid );
 466   registration->SetTransform( transform );
 467   registration->SetInitialTransformParameters( transform->GetParameters() );
 468   registration->SetInterpolator( interpolator );
 469   registration->SetMetric(metric);
 470   registration->SetOptimizer(optimizer);
 471   registration->SetNumberOfLevels( m_ArgsInfo.levels_arg );
 472   if (m_Verbose) std::cout << "Setting "<< m_ArgsInfo.levels_arg <<" resolution levels..." << std::endl;
 473   if (m_Verbose) std::cout << "Initial Transform: "<< registration->GetInitialTransformParameters()<<std::endl;
 474
 475   //============================================================================
 476   // Connecting the commander to the registration to monitor it
 477   //============================================================================
 478   if (m_Verbose) {
 479
 480     // Output iteration info
 481     CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
 482     observer->SetOptimizer(genericOptimizer);
 483     optimizer->AddObserver( itk::IterationEvent(), observer );
 484
 485
 486     // Output level info
 487     typedef RegistrationInterfaceCommand<RegistrationType, args_info_clitkAffineRegistration> CommandType;
 488     typename CommandType::Pointer command = CommandType::New();
 489     command->SetArgsInfo(m_ArgsInfo);
 490     registration->AddObserver( itk::IterationEvent(), command );
 491
 492   }
 493
 494
 495   //============================================================================
 496   // Finally we can start the registration with the given amount of multiresolution levels
 497   //============================================================================
 498   if (m_Verbose) std::cout << "Starting the registration now..." << std::endl;
 499
 500   try {
 501     registration->StartRegistration();
 502   } catch( itk::ExceptionObject & err ) {
 503     std::cerr << "ExceptionObject caught !" << std::endl;
 504     std::cerr << err << std::endl;
 505   }
 506
 507
 508   //============================================================================
 509   // Processing the result of the registration
 510   //============================================================================
 511   OptimizerType::ParametersType finalParameters = registration->GetLastTransformParameters();
 512   std::cout<< "Result : " <<std::setprecision(12)<<std::endl;
 513
 514   for (unsigned int i=nParam-Dimension; i<nParam; i++) //Translations
 515     std::cout << " Translation " << i << " = " << finalParameters[i];
 516   for (unsigned int i=0; i<nParam-Dimension; i++)      //Rest
 517     std::cout << " Other parameter " << i << " = " << finalParameters[i];
 518
 519
 520   itk::Matrix<double,Dimension+1,Dimension+1> matrix;
 521   if(m_ArgsInfo.transform_arg == 3) {
 522     for(unsigned int i=0; i<Dimension; i++) {
 523       matrix[i][3] =  finalParameters[nParam-Dimension+i];
 524       for(unsigned int j=0; j<Dimension; j++) {
 525         matrix[i][j] = finalParameters[i*3+j];
 526       }
 527       matrix[3][3] = 1.0;
 528     }
 529   } else {
 530 //    matrix = clitk::GetBackwardAffineMatrix3D(finalParameters);
 531   }
 532
 533   std::cout << " Affine transform matrix =" << std::endl;
 534   std::cout << matrix <<std::setprecision(6)<< std::endl;
 535
 536   // Write matrix to a file
 537   if(m_ArgsInfo.matrix_given) {
 538     std::ofstream mFile;
 539     mFile.open(m_ArgsInfo.matrix_arg);
 540     mFile<<std::setprecision(12)<<matrix<< std::setprecision(6)<<std::endl;
 541     mFile.close();
 542   }
 543
 544   //============================================================================
 545   // Prepare the resampling filter in order to transform the moving image.
 546   //============================================================================
 547   if (m_ArgsInfo.output_given || m_ArgsInfo.checker_after_given || m_ArgsInfo.after_given ) {
 548     transform->SetParameters( finalParameters );
 549     typedef itk::ResampleImageFilter< InternalImageType,InternalImageType >    ResampleFilterType;
 550     typename    ResampleFilterType::Pointer resampler = ResampleFilterType::New();
 551
 552     resampler->SetTransform( transform );
 553     resampler->SetInput( movingImageReader->GetOutput() );
 554     resampler->SetSize( fixedImageReader ->GetOutput()->GetLargestPossibleRegion().GetSize() );
 555     resampler->SetOutputOrigin(  fixedImageReader ->GetOutput()->GetOrigin() );
 556     resampler->SetOutputSpacing( fixedImageReader ->GetOutput()->GetSpacing() );
 557     resampler->SetDefaultPixelValue( 0 );
 558
 559     //Output?
 560     if (m_ArgsInfo.output_given) {
 561       //We write an output in the same pixeltype then the input
 562       typedef itk::ImageFileWriter< FixedImageType >  WriterType;
 563       typename WriterType::Pointer outputWriter =  WriterType::New();
 564       outputWriter->SetFileName(m_ArgsInfo.output_arg );
 565       outputWriter->SetInput( resampler->GetOutput()   );
 566       outputWriter->Update();
 567     }
 568
 569
 570     //============================================================================
 571     // Checker after?
 572     //============================================================================
 573     if (m_ArgsInfo.checker_after_given) {
 574       //To display correctly the checkerboard image, the intensities must lie in the same range (normalized)
 575       //We write the image in the internal image type
 576       typedef itk::ResampleImageFilter< InternalImageType,InternalImageType >    ResampleFilterType;
 577       typename    ResampleFilterType::Pointer internalResampler = ResampleFilterType::New();
 578       internalResampler->SetTransform( transform );
 579       internalResampler->SetInput( inputMovingImage );
 580       internalResampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
 581       internalResampler->SetOutputOrigin(  fixedImage->GetOrigin() );
 582       internalResampler->SetOutputSpacing( fixedImage->GetSpacing() );
 583       internalResampler->SetDefaultPixelValue( 0 );
 584
 585       //We pass the normalized images to the checker filter
 586       typedef itk::CheckerBoardImageFilter< InternalImageType > CheckerBoardFilterType;
 587       typename CheckerBoardFilterType::Pointer checkerFilter= CheckerBoardFilterType::New();
 588
 589       checkerFilter->SetInput1(inputFixedImage);
 590       checkerFilter->SetInput2(internalResampler->GetOutput());
 591       typedef itk::ImageFileWriter< InternalImageType >  InternalWriterType;
 592       typename  InternalWriterType::Pointer checkerWriter =  InternalWriterType::New();
 593       checkerWriter->SetFileName(m_ArgsInfo.checker_after_arg);
 594       checkerWriter->SetInput( checkerFilter->GetOutput() );
 595       checkerWriter->Update();
 596     }
 597
 598
 599     //============================================================================
 600     // Checker before?
 601     //============================================================================
 602     if (m_ArgsInfo.checker_before_given) {
 603       //To display correctly the checkerboard image, the intensities must lie in the same range (normalized)
 604       //We write the image in the internal image type
 605       //We pass the normalized images to the checker filter
 606       typedef itk::CheckerBoardImageFilter< InternalImageType > CheckerBoardFilterType;
 607       typename CheckerBoardFilterType::Pointer checkerFilter= CheckerBoardFilterType::New();
 608
 609       checkerFilter->SetInput1(inputFixedImage);
 610       checkerFilter->SetInput2(inputMovingImage);
 611       typedef itk::ImageFileWriter< InternalImageType >  InternalWriterType;
 612       typename  InternalWriterType::Pointer checkerWriter =  InternalWriterType::New();
 613       checkerWriter->SetFileName(m_ArgsInfo.checker_before_arg);
 614       checkerWriter->SetInput( checkerFilter->GetOutput() );
 615       checkerWriter->Update();
 616     }
 617
 618
 619     //============================================================================
 620     // Difference After?
 621     //============================================================================
 622     if (m_ArgsInfo.after_given) {
 623       typedef itk::SubtractImageFilter< InternalImageType, FixedImageType,FixedImageType > DifferenceImageFilterType;
 624       typename DifferenceImageFilterType::Pointer differenceAfterFilter= DifferenceImageFilterType::New();
 625
 626       differenceAfterFilter->SetInput1(fixedImageReader ->GetOutput());
 627       differenceAfterFilter->SetInput2(resampler->GetOutput());
 628
 629       // Prepare a writer to write the difference image
 630       typedef itk::ImageFileWriter< FixedImageType >  WriterType;
 631       typename   WriterType::Pointer     differenceAfterWriter =  WriterType::New();
 632       differenceAfterWriter->SetFileName(m_ArgsInfo.after_arg );
 633       differenceAfterWriter->SetInput( differenceAfterFilter->GetOutput()   );
 634       differenceAfterWriter->Update();
 635     }
 636   }
 637
 638   //============================================================================
 639   // Difference Before?
 640   //============================================================================
 641   if (m_ArgsInfo.before_given) {
 642     typedef itk::CastImageFilter< InternalImageType,FixedImageType > CastFilterType;
 643     typename    CastFilterType::Pointer  caster =  CastFilterType::New();
 644     caster->SetInput( movingImageReader->GetOutput() );
 645
 646     typedef itk::SubtractImageFilter< InternalImageType, FixedImageType, FixedImageType > DifferenceImageFilterType;
 647     typename DifferenceImageFilterType::Pointer differenceBeforeFilter= DifferenceImageFilterType::New();
 648
 649
 650     differenceBeforeFilter->SetInput1(fixedImageReader ->GetOutput());
 651     differenceBeforeFilter->SetInput2(caster->GetOutput());
 652
 653     // Prepare a writer to write the difference image
 654     typedef itk::ImageFileWriter< FixedImageType >  WriterType;
 655     typename WriterType::Pointer     differenceBeforeWriter =  WriterType::New();
 656     differenceBeforeWriter->SetFileName(m_ArgsInfo.before_arg);
 657     differenceBeforeWriter->SetInput( differenceBeforeFilter->GetOutput()   );
 658     differenceBeforeWriter->Update();
 659   }
 660
 661 }
 662 }
 663 #endif  //#define  CLITKAFFINEREGISTRATIONGENERICFILTER_TXX
 664