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