registration/clitkDemonsDeformableRegistrationGenericFilter.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://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 #ifndef __clitkDemonsDeformableRegistrationGenericFilter_txx
  19 #define __clitkDemonsDeformableRegistrationGenericFilter_txx
  20 #include "clitkDemonsDeformableRegistrationGenericFilter.h"
  21
  22
  23 namespace clitk
  24 {
  25
  26   //==============================================================================
  27   //Creating an observer class that allows us to change parameters at subsequent levels
  28   //==============================================================================
  29   template <typename MultiRegistrationFilterType>
  30   class CommandResolutionLevelUpdate : public itk::Command
  31   {
  32   public:
  33     typedef  CommandResolutionLevelUpdate   Self;
  34     typedef  itk::Command             Superclass;
  35     typedef  itk::SmartPointer<Self>  Pointer;
  36     itkNewMacro( Self );
  37
  38     //Typedefs
  39     itkStaticConstMacro(ImageDimension, unsigned int, MultiRegistrationFilterType::ImageDimension);
  40     typedef typename MultiRegistrationFilterType::FixedImageType FixedImageType;
  41     typedef typename MultiRegistrationFilterType::MovingImageType MovingImageType;
  42     typedef typename MultiRegistrationFilterType::DeformationFieldType DeformationFieldType;
  43     typedef typename MultiRegistrationFilterType::RegistrationType RegistrationType;
  44     typedef typename RegistrationType::FixedImageType InternalImageType;
  45     typedef itk::DiffeomorphicDemonsRegistrationFilter<InternalImageType, InternalImageType, DeformationFieldType> DiffFilterType;
  46     typedef itk::FastSymmetricForcesDemonsRegistrationFilter<InternalImageType, InternalImageType, DeformationFieldType> SymFilterType;
  47
  48   protected:
  49     CommandResolutionLevelUpdate()
  50     {
  51       m_CurrentLevel=0;
  52       m_MaxStep=2;
  53       m_ScaleStep=false;
  54       m_ScaleSD=false;
  55     }
  56
  57   public:
  58
  59     //Set
  60     void SetMaxRMSError(double* m){m_MaxRMSError=m;}
  61     void SetSD(double* m){m_SD=m;}
  62     void SetMaxStep(double m){m_MaxStep=m;}
  63     void SetScaleSD(bool m){m_ScaleSD=m;}
  64     void SetScaleStep(bool m){m_ScaleStep=m;}
  65     void SetRegistrationType(unsigned int i){m_RegistrationType=i;}
  66
  67     //Execute
  68     void Execute(const itk::Object * caller, const itk::EventObject & event )
  69     {
  70       std::cout << "Warning: The const Execute method shouldn't be called" << std::endl;
  71     }
  72     void Execute(itk::Object *caller, const itk::EventObject & event)
  73     {
  74
  75       // Check event
  76       if( !(itk::IterationEvent().CheckEvent( &event )) )
  77         {
  78           return;
  79         }
  80
  81       //Cast caller
  82       MultiRegistrationFilterType * filter = dynamic_cast<  MultiRegistrationFilterType * >( caller );
  83
  84       // Get the level
  85       m_CurrentLevel=filter->GetCurrentLevel();
  86       unsigned int numberOfLevels=filter->GetNumberOfLevels();
  87       unsigned int expandFactor=1<< (numberOfLevels-m_CurrentLevel-1);
  88       double * sd = m_SD;
  89       double maxStep=m_MaxStep;
  90
  91       // Scale the SD
  92       if (m_ScaleSD)
  93         {
  94           for (unsigned int i=0 ; i<ImageDimension ; i++) sd[i]*=(double)expandFactor;
  95           filter->GetRegistrationFilter()->SetStandardDeviations(sd);
  96           filter->GetRegistrationFilter()->SetUpdateFieldStandardDeviations(sd);
  97         }
  98
  99       // Scale max step: 2 & 3 only!
 100       if (m_ScaleStep)
 101         {
 102           maxStep= m_MaxStep*expandFactor;
 103           if (m_RegistrationType==2)
 104             {
 105               SymFilterType* symFilter=dynamic_cast< SymFilterType* >(filter->GetRegistrationFilter());
 106               symFilter->SetMaximumUpdateStepLength(maxStep);
 107             }
 108           else if (m_RegistrationType==3)
 109             {
 110               DiffFilterType* diffFilter=dynamic_cast< DiffFilterType* >(filter->GetRegistrationFilter());
 111               diffFilter->SetMaximumUpdateStepLength(maxStep);
 112             }
 113         }
 114
 115       // Set maxRMS
 116       filter->GetRegistrationFilter()->SetMaximumRMSError(m_MaxRMSError[m_CurrentLevel]);
 117
 118       //Print All
 119       std::cout << "--------------------------------------------------" << std::endl;
 120       std::cout << "Starting resolution level "<<m_CurrentLevel+1<<" of "<<filter->GetNumberOfLevels()<<"..."<<std::endl;
 121       std::cout << "Setting the standard deviations to [ "<<sd[0];
 122       for (unsigned int i=1; i<3; i++) std::cout<<", "<< sd[i];
 123       std::cout <<" ]..."<<std::endl;
 124       if ( (m_RegistrationType==2) || (m_RegistrationType==3) ) std::cout <<"Setting the maximum step size to "<< maxStep<<"..."<<std::endl;
 125       std::cout <<"Setting the maximum RMS field change to "<< m_MaxRMSError[m_CurrentLevel]<<"..."<<std::endl;
 126       std::cout << "--------------------------------------------------" << std::endl;
 127     }
 128
 129     double *m_MaxRMSError,*m_SD;
 130     double  m_MaxStep;
 131     unsigned int m_CurrentLevel, m_RegistrationType;
 132     bool m_ScaleSD, m_ScaleStep;
 133
 134   };
 135
 136
 137   //==============================================================================
 138   // Creating an observer class that allows output at each iteration
 139   //==============================================================================
 140   template <typename RegistrationFilterType, typename FixedImageType, typename MovingImageType>
 141   class CommandIterationUpdate : public itk::Command
 142   {
 143   public:
 144     typedef  CommandIterationUpdate   Self;
 145     typedef  itk::Command             Superclass;
 146     typedef  itk::SmartPointer<Self>  Pointer;
 147     typedef  itk::SmartPointer<const Self>  ConstPointer;
 148     itkNewMacro( Self );
 149
 150
 151     //find the multiresolution filter
 152     //     typedef typename  RegistrationFilterType::FixedImageType InternalImageType;
 153     //     typedef typename  RegistrationFilterType::MovingImageType MovingImageType;
 154 #if ITK_VERSION_MAJOR >= 4
 155     typedef typename  RegistrationFilterType::DisplacementFieldType DisplacementFieldType;
 156 #else
 157     typedef typename  RegistrationFilterType::DeformationFieldType DisplacementFieldType;
 158 #endif
 159     typedef clitk::MultiResolutionPDEDeformableRegistration<FixedImageType, MovingImageType, DisplacementFieldType> MultiResolutionRegistrationType;
 160     typedef CommandResolutionLevelUpdate<MultiResolutionRegistrationType> LevelObserver;
 161
 162   protected:
 163     CommandIterationUpdate(){};
 164
 165   public:
 166
 167     // Sets
 168     void SetStop( int* s){m_Stop=s;}
 169     void SetLevelObserver(LevelObserver* o ){m_LevelObserver=o;}
 170
 171
 172     //Execute
 173     void Execute(const itk::Object *, const itk::EventObject & )
 174     {
 175       std::cout << "Warning: The const Execute method shouldn't be called" << std::endl;
 176     }
 177
 178     void Execute(itk::Object *caller, const itk::EventObject & event)
 179     {
 180       if( !(itk::IterationEvent().CheckEvent( &event )) )
 181         {
 182           return;
 183         }
 184
 185       //Cast
 186       RegistrationFilterType * filter =  dynamic_cast<  RegistrationFilterType * >( caller );
 187
 188       if(filter)
 189         {
 190           // Get info
 191           m_Iteration=filter->GetElapsedIterations();
 192           m_Metric=filter->GetMetric();
 193
 194           // Output
 195           std::cout << m_Iteration<<"\t Field Tolerance= "<<filter->GetMaximumRMSError();
 196           std::cout <<"\t DVF Change= " << filter->GetRMSChange()<<"\t RMS= "<<m_Metric;
 197
 198           // Using stop criterion?
 199           if (m_Stop[m_LevelObserver->m_CurrentLevel]>=0)
 200             {
 201               // First iteration
 202               if(m_Iteration==1)
 203                 {
 204                   m_Minimum=m_Metric;
 205                   m_StopCounter=0;
 206                 }
 207
 208               // Less then minimum
 209               else if(m_Metric<m_Minimum)
 210                 {
 211                   m_StopCounter=0;
 212                   m_Minimum=m_Metric;
 213                 }
 214
 215               //Not less then minimum
 216               else
 217                 {
 218                   m_StopCounter++;
 219                   if (m_StopCounter>=m_Stop[m_LevelObserver->m_CurrentLevel])
 220                     filter->StopRegistration();
 221                 }
 222
 223               // Output
 224               std::cout <<"\t Stop= "<<m_StopCounter<<" / "<<m_Stop[m_LevelObserver->m_CurrentLevel]<<std::endl;
 225             }
 226
 227           // Not using stop criterion
 228           else std::cout <<std::endl;
 229         }
 230     }
 231
 232     double m_Minimum, m_Metric;
 233     int m_StopCounter, m_Iteration;
 234     int * m_Stop;
 235     LevelObserver* m_LevelObserver;
 236   };
 237
 238
 239   //==============================================================================
 240   // Update with the number of dimensions
 241   //==============================================================================
 242   template<unsigned int Dimension>
 243   void DemonsDeformableRegistrationGenericFilter::UpdateWithDim(std::string PixelType)
 244   {
 245     if(PixelType == "short"){
 246       if (m_Verbose) std::cout  << "Launching warp in "<< Dimension <<"D and signed short..." << std::endl;
 247       UpdateWithDimAndPixelType<Dimension, signed short>();
 248     }
 249     //    else if(PixelType == "unsigned_short"){
 250     //       if (m_Verbose) std::cout  << "Launching warp in "<< Dimension <<"D and unsigned_short..." << std::endl;
 251     //       UpdateWithDimAndPixelType<Dimension, unsigned short>();
 252     //     }
 253
 254     //     else if (PixelType == "unsigned_char"){
 255     //       if (m_Verbose) std::cout  << "Launching warp in "<< Dimension <<"D and unsigned_char..." << std::endl;
 256     //       UpdateWithDimAndPixelType<Dimension, unsigned char>();
 257     //     }
 258
 259     //     else if (PixelType == "char"){
 260     //       if (m_Verbose) std::cout  << "Launching warp in "<< Dimension <<"D and signed_char..." << std::endl;
 261     //       UpdateWithDimAndPixelType<Dimension, signed char>();
 262     //    }
 263     else {
 264       if (m_Verbose) std::cout  << "Launching warp in "<< Dimension <<"D and float..." << std::endl;
 265       UpdateWithDimAndPixelType<Dimension, float>();
 266     }
 267   }
 268
 269
 270
 271   //==============================================================================
 272   // Update with the number of dimensions and pixeltype
 273   //==============================================================================
 274   //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 275
 276   template<unsigned int ImageDimension, class PixelType>
 277   void DemonsDeformableRegistrationGenericFilter::UpdateWithDimAndPixelType()
 278   {
 279     //=======================================================
 280     // Run-time
 281     //=======================================================
 282     bool threadsGiven=m_ArgsInfo.threads_given;
 283     int threads=m_ArgsInfo.threads_arg;
 284
 285     typedef itk::Image< PixelType, ImageDimension >  FixedImageType;
 286     typedef itk::Image< PixelType, ImageDimension >  MovingImageType;
 287     typedef itk::Image< float, ImageDimension >  InternalImageType;
 288     typedef double TCoordRep;
 289
 290
 291     //=======================================================
 292     //Input
 293     //=======================================================
 294     typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
 295     typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
 296
 297     typename FixedImageReaderType::Pointer  fixedImageReader  = FixedImageReaderType::New();
 298     typename MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
 299
 300     fixedImageReader->SetFileName(  m_ArgsInfo.reference_arg );
 301     movingImageReader->SetFileName( m_ArgsInfo.target_arg );
 302     if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
 303     fixedImageReader->Update();
 304     movingImageReader->Update();
 305
 306     typename FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
 307     typename MovingImageType::Pointer movingImage =movingImageReader->GetOutput();
 308     typename FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
 309
 310
 311     //=======================================================
 312     //Output
 313     //=======================================================
 314     typedef itk::Vector< float, ImageDimension >    VectorPixelType;
 315     typedef itk::Image<  VectorPixelType, ImageDimension > DeformationFieldType;
 316
 317
 318     //=======================================================
 319     //Pyramids
 320     //=======================================================
 321     typedef itk::RecursiveMultiResolutionPyramidImageFilter< FixedImageType, InternalImageType>    FixedImagePyramidType;
 322     typedef itk::RecursiveMultiResolutionPyramidImageFilter< MovingImageType, InternalImageType>    MovingImagePyramidType;
 323     // typedef itk::MultiResolutionPyramidImageFilter< FixedImageType, FixedImageType>    FixedImagePyramidType;
 324     // typedef itk::MultiResolutionPyramidImageFilter< MovingImageType, MovingImageType>    MovingImagePyramidType;
 325     typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
 326     typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
 327
 328
 329     // =======================================================
 330     // Define the registation filters
 331     // =======================================================
 332     typedef itk::PDEDeformableRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType>   RegistrationFilterType;
 333     typename RegistrationFilterType::Pointer pdeFilter;
 334     typedef clitk::MultiResolutionPDEDeformableRegistration<FixedImageType,MovingImageType,DeformationFieldType >  MultiResolutionRegistrationFilterType;
 335     typename MultiResolutionRegistrationFilterType::Pointer multiResolutionFilter =  MultiResolutionRegistrationFilterType::New();
 336     typedef itk::ESMDemonsRegistrationFunction<InternalImageType, InternalImageType, DeformationFieldType> RegistrationFunctionType;
 337
 338
 339     // =======================================================
 340     // The multiresolution scheme
 341     // =======================================================
 342     if (threadsGiven) multiResolutionFilter->SetNumberOfThreads(threads);
 343     unsigned int nLevels=m_ArgsInfo.levels_arg;
 344     if (m_Verbose) std::cout<<"Setting the number of resolution levels to "<<nLevels<<"..."<<std::endl;
 345     multiResolutionFilter->SetFixedImage( fixedImage );
 346     multiResolutionFilter->SetMovingImage( movingImage );
 347     multiResolutionFilter->SetNumberOfLevels( nLevels );
 348     multiResolutionFilter->SetFixedImagePyramid( fixedImagePyramid );
 349     multiResolutionFilter->SetMovingImagePyramid( movingImagePyramid );
 350     if (threadsGiven) multiResolutionFilter->SetNumberOfThreads( threads );
 351
 352     //------------------------------------
 353     //Set the number of iterations
 354     //------------------------------------
 355     std::vector<unsigned int> nIterations(nLevels);
 356     for (unsigned int i=0 ; i<nLevels; i++)
 357       {
 358         if (m_ArgsInfo.maxIter_given==nLevels)
 359           {
 360             nIterations[i] =m_ArgsInfo.maxIter_arg[i];
 361           }
 362         else
 363           nIterations[i]=m_ArgsInfo.maxIter_arg[0];
 364       }
 365     multiResolutionFilter->SetNumberOfIterations( &(nIterations[0]) );
 366     if(m_Verbose) {
 367       std::cout<<"Setting the number of iterations to: "<<nIterations[0];
 368       for (unsigned int i=1; i<nLevels; i++)
 369         std::cout<<", "<<nIterations [i];
 370       std::cout<<std::endl;
 371     }
 372
 373     //------------------------------------
 374     //Set the max RMS error for the field update
 375     //------------------------------------
 376     std::vector<double> maxRMSError(nLevels);
 377     for (unsigned int i=0 ; i<nLevels; i++)
 378       {
 379         if (m_ArgsInfo.maxRMSError_given==nLevels)
 380           maxRMSError[i] =m_ArgsInfo.maxRMSError_arg[i];
 381         else
 382           maxRMSError[i]=m_ArgsInfo.maxRMSError_arg[0];
 383       }
 384     if(m_Verbose) {
 385       std::cout<<"Setting the max root mean squared error to: "<<maxRMSError[0];
 386       for (unsigned int i=1; i<nLevels; i++)
 387         std::cout<<", "<<maxRMSError[i];
 388       std::cout<<std::endl;
 389     }
 390
 391     //------------------------------------
 392     //Get the stop criterion
 393     //------------------------------------
 394     std::vector<int> stop(nLevels);
 395     for (unsigned int i=0; i<nLevels; i++)
 396       if (m_ArgsInfo.stop_given==nLevels)
 397         stop[i]=m_ArgsInfo.stop_arg[i];
 398       else
 399         stop[i]=m_ArgsInfo.stop_arg[0];
 400     if(m_Verbose) {
 401       std::cout<<"Setting the stop criterion to : "<<stop[0];
 402       for (unsigned int i=1; i<nLevels; i++)
 403         std::cout<<", "<<stop[i];
 404       std::cout<<std::endl;
 405     }
 406
 407     //------------------------------------
 408     //Grad type
 409     //------------------------------------
 410     typename RegistrationFunctionType::GradientType grad=RegistrationFunctionType::Symmetric;
 411     switch (m_ArgsInfo.gradType_arg)
 412       {
 413       case 0:
 414         grad= RegistrationFunctionType::Symmetric;
 415         break;
 416       case 1:
 417         grad= RegistrationFunctionType::Fixed;
 418         break;
 419       case 2:
 420         grad= RegistrationFunctionType::WarpedMoving;
 421         break;
 422       case 3:
 423         grad= RegistrationFunctionType::MappedMoving;
 424         break;
 425       }
 426
 427     //------------------------------------
 428     // Set smoothing standard deviations
 429     //------------------------------------
 430     double sd [ImageDimension];
 431     for (unsigned int i=0; i<ImageDimension; i++)
 432       if (m_ArgsInfo.sd_given==ImageDimension)
 433         sd[i]=m_ArgsInfo.sd_arg[i];
 434       else
 435         sd[i]=m_ArgsInfo.sd_arg[0];
 436     if(m_Verbose) {
 437       std::cout<<"Setting the final standard deviations for smoothing to: "<<sd[0];
 438       for (unsigned int i=1; i<ImageDimension; i++)
 439         std::cout<<", "<<sd[i];
 440       std::cout<<std::endl;
 441     }
 442
 443     // =======================================================
 444     // The level observer: adjust settings at each level
 445     // =======================================================
 446     typename CommandResolutionLevelUpdate<MultiResolutionRegistrationFilterType>::Pointer levelObserver =
 447       CommandResolutionLevelUpdate<MultiResolutionRegistrationFilterType>::New();
 448     multiResolutionFilter->AddObserver( itk::IterationEvent(), levelObserver );
 449     levelObserver->SetMaxRMSError(&(maxRMSError[0]));
 450     levelObserver->SetMaxStep(m_ArgsInfo.maxStep_arg);
 451     levelObserver->SetSD(sd);
 452     levelObserver->SetScaleStep(m_ArgsInfo.scaleStep_flag);
 453     levelObserver->SetScaleSD(m_ArgsInfo.scaleSD_flag);
 454     levelObserver->SetRegistrationType(m_ArgsInfo.demons_arg);
 455
 456
 457     // =======================================================
 458     // The type of filter
 459     // =======================================================
 460     switch (m_ArgsInfo.demons_arg){
 461
 462     case 0:
 463       {
 464         typedef itk::DemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
 465         typename DemonsFilterType::Pointer m =DemonsFilterType::New();
 466
 467         //Set Parameters for this filter
 468         m->SetIntensityDifferenceThreshold( m_ArgsInfo.intThreshold_arg);
 469         m->SetUseMovingImageGradient( m_ArgsInfo.movGrad_flag);
 470         typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
 471         observer->SetStop(&(stop[0]));
 472         observer->SetLevelObserver(levelObserver);
 473         m->AddObserver( itk::IterationEvent(), observer );
 474         if (m_Verbose)  std::cout<<"Using the demons registration filter..."<<std::endl;
 475           pdeFilter=m;
 476           break;
 477         }
 478
 479     case 1:
 480       {
 481         typedef itk::SymmetricForcesDemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
 482         typename DemonsFilterType::Pointer m =DemonsFilterType::New();
 483
 484         //Set Parameters for this filter
 485         m->SetIntensityDifferenceThreshold( m_ArgsInfo.intThreshold_arg);
 486         typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
 487         observer->SetStop(&(stop[0]));
 488         observer->SetLevelObserver(levelObserver);
 489         m->AddObserver( itk::IterationEvent(), observer );
 490         if (m_Verbose) std::cout<<"Using the symmetric forces demons registration filter..."<<std::endl;
 491         pdeFilter=m;
 492         break;
 493       }
 494
 495     case 2:
 496       {
 497         typedef itk::FastSymmetricForcesDemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
 498         typename DemonsFilterType::Pointer m = DemonsFilterType::New();
 499
 500         //Set Parameters for this filter
 501         m->SetIntensityDifferenceThreshold( m_ArgsInfo.intThreshold_arg);
 502         m->SetMaximumUpdateStepLength( m_ArgsInfo.maxStep_arg);
 503         m->SetUseGradientType(grad);
 504         typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
 505         observer->SetStop(&(stop[0]));
 506         observer->SetLevelObserver(levelObserver);
 507         m->AddObserver( itk::IterationEvent(), observer );
 508         if (m_Verbose) std::cout<<"Using the fast symmetric forces demons registration filter..."<<std::endl;
 509         pdeFilter=m;
 510         break;
 511       }
 512
 513     case 3:
 514       {
 515         typedef itk::DiffeomorphicDemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
 516         typename DemonsFilterType::Pointer m =  DemonsFilterType::New();
 517
 518         //Set Parameters for this filter
 519         m->SetMaximumUpdateStepLength( m_ArgsInfo.maxStep_arg);
 520         m->SetUseFirstOrderExp(m_ArgsInfo.firstOrder_flag);
 521         m->SetUseGradientType(grad);
 522         typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
 523         observer->SetStop(&(stop[0]));
 524         observer->SetLevelObserver(levelObserver);
 525         m->AddObserver( itk::IterationEvent(), observer );
 526         if (m_Verbose) std::cout<<"Using the diffeomorphic demons registration filter..."<<std::endl;
 527         pdeFilter=m;
 528         break;
 529       }
 530     }
 531
 532
 533
 534     //Set common options
 535     pdeFilter->SetStandardDeviations( sd );
 536     pdeFilter->SetUpdateFieldStandardDeviations( sd );
 537     //JV TODO
 538     // pdeFilter->SetMaximumError(m_ArgsInfo.maxError_arg);
 539     // pdeFilter->SetMaximumKernelWidth(m_ArgsInfo.maxError_arg);
 540 #if ITK_VERSION_MAJOR >= 4
 541     pdeFilter->SetSmoothDisplacementField(!m_ArgsInfo.fluid_flag);
 542 #else
 543     pdeFilter->SetSmoothDeformationField(!m_ArgsInfo.fluid_flag);
 544 #endif
 545     pdeFilter->SetSmoothUpdateField(m_ArgsInfo.fluid_flag);
 546     pdeFilter->SetUseImageSpacing( m_ArgsInfo.spacing_flag );
 547
 548     //Pass to the multi resolution scheme
 549     multiResolutionFilter->SetRegistrationFilter( pdeFilter );
 550
 551
 552     // =======================================================
 553     // The initial solution
 554     // =======================================================
 555     if (m_ArgsInfo.init_given)
 556       {
 557         typedef itk::ImageFileReader<DeformationFieldType> DeformationFieldReaderType;
 558         typename DeformationFieldReaderType::Pointer defReader=DeformationFieldReaderType::New();
 559         defReader->SetFileName(m_ArgsInfo.init_arg);
 560         defReader->Update();
 561         multiResolutionFilter->SetArbitraryInitialDeformationField(defReader->GetOutput());
 562       }
 563
 564
 565     // =======================================================
 566     // Execute
 567     // =======================================================
 568     try
 569       {
 570         multiResolutionFilter->Update();
 571       }
 572     catch( itk::ExceptionObject & excp )
 573       {
 574         std::cerr <<"Error executing the demons filter: "<< excp << std::endl;
 575         return;
 576       }
 577
 578
 579     //=======================================================
 580     // Get the output
 581     //=======================================================
 582     typename DeformationFieldType::Pointer deformationField=multiResolutionFilter->GetOutput();
 583
 584
 585     //=======================================================
 586     // Write the DVF
 587     //=======================================================
 588     typedef itk::ImageFileWriter< DeformationFieldType >  FieldWriterType;
 589     typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
 590     fieldWriter->SetInput( deformationField );
 591     fieldWriter->SetFileName( m_ArgsInfo.vf_arg );
 592     try
 593       {
 594         fieldWriter->Update();
 595       }
 596     catch( itk::ExceptionObject & excp )
 597       {
 598         std::cerr << "Exception thrown writing the DVF" << std::endl;
 599         std::cerr << excp << std::endl;
 600         return;
 601       }
 602
 603
 604     //=======================================================
 605     // Warp the moving image
 606     //=======================================================
 607     typedef itk::WarpImageFilter< MovingImageType, FixedImageType, DeformationFieldType >    WarpFilterType;
 608     typename WarpFilterType::Pointer warp = WarpFilterType::New();
 609
 610 #if ITK_VERSION_MAJOR >= 4
 611     warp->SetDisplacementField( deformationField );
 612 #else
 613     warp->SetDeformationField( deformationField );
 614 #endif
 615     warp->SetInput( movingImageReader->GetOutput() );
 616     warp->SetOutputOrigin(  fixedImage->GetOrigin() );
 617     warp->SetOutputSpacing( fixedImage->GetSpacing() );
 618     warp->SetOutputDirection( fixedImage->GetDirection() );
 619     warp->SetEdgePaddingValue( 0.0 );
 620     warp->Update();
 621
 622
 623     //=======================================================
 624     // Write the warped image
 625     //=======================================================
 626     typedef itk::ImageFileWriter< FixedImageType >  WriterType;
 627     typename WriterType::Pointer      writer =  WriterType::New();
 628     writer->SetFileName( m_ArgsInfo.output_arg );
 629     writer->SetInput( warp->GetOutput()    );
 630
 631     try
 632       {
 633         writer->Update();
 634       }
 635     catch( itk::ExceptionObject & err )
 636       {
 637         std::cerr << "ExceptionObject caught !" << std::endl;
 638         std::cerr << err << std::endl;
 639         return;
 640       }
 641
 642
 643     //=======================================================
 644     // Calculate the difference after the deformable transform
 645     //=======================================================
 646     typedef clitk::DifferenceImageFilter<  FixedImageType, FixedImageType> DifferenceFilterType;
 647     if (m_ArgsInfo.after_given)
 648       {
 649         typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
 650         difference->SetValidInput( fixedImage );
 651         difference->SetTestInput( warp->GetOutput() );
 652
 653         try
 654           {
 655             difference->Update();
 656           }
 657         catch( itk::ExceptionObject & err )
 658           {
 659             std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
 660             std::cerr << err << std::endl;
 661             return;
 662           }
 663
 664         typename WriterType::Pointer differenceWriter=WriterType::New();
 665         differenceWriter->SetInput(difference->GetOutput());
 666         differenceWriter->SetFileName(m_ArgsInfo.after_arg);
 667         differenceWriter->Update();
 668
 669       }
 670
 671
 672     //=======================================================
 673     // Calculate the difference before the deformable transform
 674     //=======================================================
 675     if( m_ArgsInfo.before_given )
 676       {
 677         typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
 678         difference->SetValidInput( fixedImage );
 679         difference->SetTestInput( movingImage );
 680
 681         try
 682           {
 683             difference->Update();
 684           }
 685         catch( itk::ExceptionObject & err )
 686           {
 687             std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
 688             std::cerr << err << std::endl;
 689             return;
 690           }
 691
 692         typename WriterType::Pointer differenceWriter=WriterType::New();
 693         writer->SetFileName( m_ArgsInfo.before_arg  );
 694         writer->SetInput( difference->GetOutput()  );
 695         writer->Update( );
 696       }
 697
 698     return;
 699   }
 700
 701 }//end namespace
 702
 703 #endif // __clitkDemonsDeformableRegistrationGenericFilter_txx