registration/clitkAffineRegistrationGenericFilter.cxx

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