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