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