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