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     typedef typename  RegistrationFilterType::DisplacementFieldType DisplacementFieldType;
 155     typedef clitk::MultiResolutionPDEDeformableRegistration<FixedImageType, MovingImageType, DisplacementFieldType> MultiResolutionRegistrationType;
 156     typedef CommandResolutionLevelUpdate<MultiResolutionRegistrationType> LevelObserver;
 157
 158   protected:
 159     CommandIterationUpdate(){};
 160
 161   public:
 162
 163     // Sets
 164     void SetStop( int* s){m_Stop=s;}
 165     void SetLevelObserver(LevelObserver* o ){m_LevelObserver=o;}
 166
 167
 168     //Execute
 169     void Execute(const itk::Object *, const itk::EventObject & )
 170     {
 171       std::cout << "Warning: The const Execute method shouldn't be called" << std::endl;
 172     }
 173
 174     void Execute(itk::Object *caller, const itk::EventObject & event)
 175     {
 176       if( !(itk::IterationEvent().CheckEvent( &event )) )
 177         {
 178           return;
 179         }
 180
 181       //Cast
 182       RegistrationFilterType * filter =  dynamic_cast<  RegistrationFilterType * >( caller );
 183
 184       if(filter)
 185         {
 186           // Get info
 187           m_Iteration=filter->GetElapsedIterations();
 188           m_Metric=filter->GetMetric();
 189
 190           // Output
 191           std::cout << m_Iteration<<"\t Field Tolerance= "<<filter->GetMaximumRMSError();
 192           std::cout <<"\t DVF Change= " << filter->GetRMSChange()<<"\t RMS= "<<m_Metric;
 193
 194           // Using stop criterion?
 195           if (m_Stop[m_LevelObserver->m_CurrentLevel]>=0)
 196             {
 197               // First iteration
 198               if(m_Iteration==1)
 199                 {
 200                   m_Minimum=m_Metric;
 201                   m_StopCounter=0;
 202                 }
 203
 204               // Less then minimum
 205               else if(m_Metric<m_Minimum)
 206                 {
 207                   m_StopCounter=0;
 208                   m_Minimum=m_Metric;
 209                 }
 210
 211               //Not less then minimum
 212               else
 213                 {
 214                   m_StopCounter++;
 215                   if (m_StopCounter>=m_Stop[m_LevelObserver->m_CurrentLevel])
 216                     filter->StopRegistration();
 217                 }
 218
 219               // Output
 220               std::cout <<"\t Stop= "<<m_StopCounter<<" / "<<m_Stop[m_LevelObserver->m_CurrentLevel]<<std::endl;
 221             }
 222
 223           // Not using stop criterion
 224           else std::cout <<std::endl;
 225         }
 226     }
 227
 228     double m_Minimum, m_Metric;
 229     int m_StopCounter, m_Iteration;
 230     int * m_Stop;
 231     LevelObserver* m_LevelObserver;
 232   };
 233
 234
 235   //==============================================================================
 236   // Update with the number of dimensions
 237   //==============================================================================
 238   template<unsigned int Dimension>
 239   void DemonsDeformableRegistrationGenericFilter::UpdateWithDim(std::string PixelType)
 240   {
 241     if(PixelType == "short"){
 242       if (m_Verbose) std::cout  << "Launching warp in "<< Dimension <<"D and signed short..." << std::endl;
 243       UpdateWithDimAndPixelType<Dimension, signed short>();
 244     }
 245     //    else if(PixelType == "unsigned_short"){
 246     //       if (m_Verbose) std::cout  << "Launching warp in "<< Dimension <<"D and unsigned_short..." << std::endl;
 247     //       UpdateWithDimAndPixelType<Dimension, unsigned short>();
 248     //     }
 249
 250     //     else if (PixelType == "unsigned_char"){
 251     //       if (m_Verbose) std::cout  << "Launching warp in "<< Dimension <<"D and unsigned_char..." << std::endl;
 252     //       UpdateWithDimAndPixelType<Dimension, unsigned char>();
 253     //     }
 254
 255     //     else if (PixelType == "char"){
 256     //       if (m_Verbose) std::cout  << "Launching warp in "<< Dimension <<"D and signed_char..." << std::endl;
 257     //       UpdateWithDimAndPixelType<Dimension, signed char>();
 258     //    }
 259     else {
 260       if (m_Verbose) std::cout  << "Launching warp in "<< Dimension <<"D and float..." << std::endl;
 261       UpdateWithDimAndPixelType<Dimension, float>();
 262     }
 263   }
 264
 265
 266
 267   //==============================================================================
 268   // Update with the number of dimensions and pixeltype
 269   //==============================================================================
 270   //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 271
 272   template<unsigned int ImageDimension, class PixelType>
 273   void DemonsDeformableRegistrationGenericFilter::UpdateWithDimAndPixelType()
 274   {
 275     //=======================================================
 276     // Run-time
 277     //=======================================================
 278     bool threadsGiven=m_ArgsInfo.threads_given;
 279     int threads=m_ArgsInfo.threads_arg;
 280
 281     typedef itk::Image< PixelType, ImageDimension >  FixedImageType;
 282     typedef itk::Image< PixelType, ImageDimension >  MovingImageType;
 283     typedef itk::Image< float, ImageDimension >  InternalImageType;
 284     typedef double TCoordRep;
 285
 286
 287     //=======================================================
 288     //Input
 289     //=======================================================
 290     typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
 291     typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
 292
 293     typename FixedImageReaderType::Pointer  fixedImageReader  = FixedImageReaderType::New();
 294     typename MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
 295
 296     fixedImageReader->SetFileName(  m_ArgsInfo.reference_arg );
 297     movingImageReader->SetFileName( m_ArgsInfo.target_arg );
 298     if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
 299     fixedImageReader->Update();
 300     movingImageReader->Update();
 301
 302     typename FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
 303     typename MovingImageType::Pointer movingImage =movingImageReader->GetOutput();
 304     typename FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
 305
 306
 307     //=======================================================
 308     //Output
 309     //=======================================================
 310     typedef itk::Vector< float, ImageDimension >    VectorPixelType;
 311     typedef itk::Image<  VectorPixelType, ImageDimension > DeformationFieldType;
 312
 313
 314     //=======================================================
 315     //Pyramids
 316     //=======================================================
 317     typedef itk::RecursiveMultiResolutionPyramidImageFilter< FixedImageType, InternalImageType>    FixedImagePyramidType;
 318     typedef itk::RecursiveMultiResolutionPyramidImageFilter< MovingImageType, InternalImageType>    MovingImagePyramidType;
 319     // typedef itk::MultiResolutionPyramidImageFilter< FixedImageType, FixedImageType>    FixedImagePyramidType;
 320     // typedef itk::MultiResolutionPyramidImageFilter< MovingImageType, MovingImageType>    MovingImagePyramidType;
 321     typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
 322     typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
 323
 324
 325     // =======================================================
 326     // Define the registation filters
 327     // =======================================================
 328     typedef itk::PDEDeformableRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType>   RegistrationFilterType;
 329     typename RegistrationFilterType::Pointer pdeFilter;
 330     typedef clitk::MultiResolutionPDEDeformableRegistration<FixedImageType,MovingImageType,DeformationFieldType >  MultiResolutionRegistrationFilterType;
 331     typename MultiResolutionRegistrationFilterType::Pointer multiResolutionFilter =  MultiResolutionRegistrationFilterType::New();
 332     typedef itk::ESMDemonsRegistrationFunction<InternalImageType, InternalImageType, DeformationFieldType> RegistrationFunctionType;
 333
 334
 335     // =======================================================
 336     // The multiresolution scheme
 337     // =======================================================
 338     if (threadsGiven) {
 339 #if ITK_VERSION_MAJOR <= 4
 340       multiResolutionFilter->SetNumberOfThreads(threads);
 341 #else
 342       multiResolutionFilter->SetNumberOfWorkUnits(threads);
 343 #endif
 344     }
 345     unsigned int nLevels=m_ArgsInfo.levels_arg;
 346     if (m_Verbose) std::cout<<"Setting the number of resolution levels to "<<nLevels<<"..."<<std::endl;
 347     multiResolutionFilter->SetFixedImage( fixedImage );
 348     multiResolutionFilter->SetMovingImage( movingImage );
 349     multiResolutionFilter->SetNumberOfLevels( nLevels );
 350     multiResolutionFilter->SetFixedImagePyramid( fixedImagePyramid );
 351     multiResolutionFilter->SetMovingImagePyramid( movingImagePyramid );
 352     if (threadsGiven) {
 353 #if ITK_VERSION_MAJOR <= 4
 354       multiResolutionFilter->SetNumberOfThreads( threads );
 355 #else
 356       multiResolutionFilter->SetNumberOfWorkUnits( threads );
 357 #endif
 358     }
 359
 360     //------------------------------------
 361     //Set the number of iterations
 362     //------------------------------------
 363     std::vector<unsigned int> nIterations(nLevels);
 364     for (unsigned int i=0 ; i<nLevels; i++)
 365       {
 366         if (m_ArgsInfo.maxIter_given==nLevels)
 367           {
 368             nIterations[i] =m_ArgsInfo.maxIter_arg[i];
 369           }
 370         else
 371           nIterations[i]=m_ArgsInfo.maxIter_arg[0];
 372       }
 373     multiResolutionFilter->SetNumberOfIterations( &(nIterations[0]) );
 374     if(m_Verbose) {
 375       std::cout<<"Setting the number of iterations to: "<<nIterations[0];
 376       for (unsigned int i=1; i<nLevels; i++)
 377         std::cout<<", "<<nIterations [i];
 378       std::cout<<std::endl;
 379     }
 380
 381     //------------------------------------
 382     //Set the max RMS error for the field update
 383     //------------------------------------
 384     std::vector<double> maxRMSError(nLevels);
 385     for (unsigned int i=0 ; i<nLevels; i++)
 386       {
 387         if (m_ArgsInfo.maxRMSError_given==nLevels)
 388           maxRMSError[i] =m_ArgsInfo.maxRMSError_arg[i];
 389         else
 390           maxRMSError[i]=m_ArgsInfo.maxRMSError_arg[0];
 391       }
 392     if(m_Verbose) {
 393       std::cout<<"Setting the max root mean squared error to: "<<maxRMSError[0];
 394       for (unsigned int i=1; i<nLevels; i++)
 395         std::cout<<", "<<maxRMSError[i];
 396       std::cout<<std::endl;
 397     }
 398
 399     //------------------------------------
 400     //Get the stop criterion
 401     //------------------------------------
 402     std::vector<int> stop(nLevels);
 403     for (unsigned int i=0; i<nLevels; i++)
 404       if (m_ArgsInfo.stop_given==nLevels)
 405         stop[i]=m_ArgsInfo.stop_arg[i];
 406       else
 407         stop[i]=m_ArgsInfo.stop_arg[0];
 408     if(m_Verbose) {
 409       std::cout<<"Setting the stop criterion to : "<<stop[0];
 410       for (unsigned int i=1; i<nLevels; i++)
 411         std::cout<<", "<<stop[i];
 412       std::cout<<std::endl;
 413     }
 414
 415     //------------------------------------
 416     //Grad type
 417     //------------------------------------
 418     typename RegistrationFunctionType::GradientType grad=RegistrationFunctionType::Symmetric;
 419     switch (m_ArgsInfo.gradType_arg)
 420       {
 421       case 0:
 422         grad= RegistrationFunctionType::Symmetric;
 423         break;
 424       case 1:
 425         grad= RegistrationFunctionType::Fixed;
 426         break;
 427       case 2:
 428         grad= RegistrationFunctionType::WarpedMoving;
 429         break;
 430       case 3:
 431         grad= RegistrationFunctionType::MappedMoving;
 432         break;
 433       }
 434
 435     //------------------------------------
 436     // Set smoothing standard deviations
 437     //------------------------------------
 438     double sd [ImageDimension];
 439     for (unsigned int i=0; i<ImageDimension; i++)
 440       if (m_ArgsInfo.sd_given==ImageDimension)
 441         sd[i]=m_ArgsInfo.sd_arg[i];
 442       else
 443         sd[i]=m_ArgsInfo.sd_arg[0];
 444     if(m_Verbose) {
 445       std::cout<<"Setting the final standard deviations for smoothing to: "<<sd[0];
 446       for (unsigned int i=1; i<ImageDimension; i++)
 447         std::cout<<", "<<sd[i];
 448       std::cout<<std::endl;
 449     }
 450
 451     // =======================================================
 452     // The level observer: adjust settings at each level
 453     // =======================================================
 454     typename CommandResolutionLevelUpdate<MultiResolutionRegistrationFilterType>::Pointer levelObserver =
 455       CommandResolutionLevelUpdate<MultiResolutionRegistrationFilterType>::New();
 456     multiResolutionFilter->AddObserver( itk::IterationEvent(), levelObserver );
 457     levelObserver->SetMaxRMSError(&(maxRMSError[0]));
 458     levelObserver->SetMaxStep(m_ArgsInfo.maxStep_arg);
 459     levelObserver->SetSD(sd);
 460     levelObserver->SetScaleStep(m_ArgsInfo.scaleStep_flag);
 461     levelObserver->SetScaleSD(m_ArgsInfo.scaleSD_flag);
 462     levelObserver->SetRegistrationType(m_ArgsInfo.demons_arg);
 463
 464
 465     // =======================================================
 466     // The type of filter
 467     // =======================================================
 468     switch (m_ArgsInfo.demons_arg){
 469
 470     case 0:
 471       {
 472         typedef itk::DemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
 473         typename DemonsFilterType::Pointer m =DemonsFilterType::New();
 474
 475         //Set Parameters for this filter
 476         m->SetIntensityDifferenceThreshold( m_ArgsInfo.intThreshold_arg);
 477         m->SetUseMovingImageGradient( m_ArgsInfo.movGrad_flag);
 478         typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
 479         observer->SetStop(&(stop[0]));
 480         observer->SetLevelObserver(levelObserver);
 481         m->AddObserver( itk::IterationEvent(), observer );
 482         if (m_Verbose)  std::cout<<"Using the demons registration filter..."<<std::endl;
 483           pdeFilter=m;
 484           break;
 485         }
 486
 487     case 1:
 488       {
 489         typedef itk::SymmetricForcesDemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
 490         typename DemonsFilterType::Pointer m =DemonsFilterType::New();
 491
 492         //Set Parameters for this filter
 493         m->SetIntensityDifferenceThreshold( m_ArgsInfo.intThreshold_arg);
 494         typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
 495         observer->SetStop(&(stop[0]));
 496         observer->SetLevelObserver(levelObserver);
 497         m->AddObserver( itk::IterationEvent(), observer );
 498         if (m_Verbose) std::cout<<"Using the symmetric forces demons registration filter..."<<std::endl;
 499         pdeFilter=m;
 500         break;
 501       }
 502
 503     case 2:
 504       {
 505         typedef itk::FastSymmetricForcesDemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
 506         typename DemonsFilterType::Pointer m = DemonsFilterType::New();
 507
 508         //Set Parameters for this filter
 509         m->SetIntensityDifferenceThreshold( m_ArgsInfo.intThreshold_arg);
 510         m->SetMaximumUpdateStepLength( m_ArgsInfo.maxStep_arg);
 511         m->SetUseGradientType(grad);
 512         typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
 513         observer->SetStop(&(stop[0]));
 514         observer->SetLevelObserver(levelObserver);
 515         m->AddObserver( itk::IterationEvent(), observer );
 516         if (m_Verbose) std::cout<<"Using the fast symmetric forces demons registration filter..."<<std::endl;
 517         pdeFilter=m;
 518         break;
 519       }
 520
 521     case 3:
 522       {
 523         typedef itk::DiffeomorphicDemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
 524         typename DemonsFilterType::Pointer m =  DemonsFilterType::New();
 525
 526         //Set Parameters for this filter
 527         m->SetMaximumUpdateStepLength( m_ArgsInfo.maxStep_arg);
 528         m->SetUseFirstOrderExp(m_ArgsInfo.firstOrder_flag);
 529         m->SetUseGradientType(grad);
 530         typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
 531         observer->SetStop(&(stop[0]));
 532         observer->SetLevelObserver(levelObserver);
 533         m->AddObserver( itk::IterationEvent(), observer );
 534         if (m_Verbose) std::cout<<"Using the diffeomorphic demons registration filter..."<<std::endl;
 535         pdeFilter=m;
 536         break;
 537       }
 538     }
 539
 540
 541
 542     //Set common options
 543     pdeFilter->SetStandardDeviations( sd );
 544     pdeFilter->SetUpdateFieldStandardDeviations( sd );
 545     //JV TODO
 546     // pdeFilter->SetMaximumError(m_ArgsInfo.maxError_arg);
 547     // pdeFilter->SetMaximumKernelWidth(m_ArgsInfo.maxError_arg);
 548     pdeFilter->SetSmoothDisplacementField(!m_ArgsInfo.fluid_flag);
 549     pdeFilter->SetSmoothUpdateField(m_ArgsInfo.fluid_flag);
 550     pdeFilter->SetUseImageSpacing( m_ArgsInfo.spacing_flag );
 551
 552     //Pass to the multi resolution scheme
 553     multiResolutionFilter->SetRegistrationFilter( pdeFilter );
 554
 555
 556     // =======================================================
 557     // The initial solution
 558     // =======================================================
 559     if (m_ArgsInfo.init_given)
 560       {
 561         typedef itk::ImageFileReader<DeformationFieldType> DeformationFieldReaderType;
 562         typename DeformationFieldReaderType::Pointer defReader=DeformationFieldReaderType::New();
 563         defReader->SetFileName(m_ArgsInfo.init_arg);
 564         defReader->Update();
 565         multiResolutionFilter->SetArbitraryInitialDeformationField(defReader->GetOutput());
 566       }
 567
 568
 569     // =======================================================
 570     // Execute
 571     // =======================================================
 572     try
 573       {
 574         multiResolutionFilter->Update();
 575       }
 576     catch( itk::ExceptionObject & excp )
 577       {
 578         std::cerr <<"Error executing the demons filter: "<< excp << std::endl;
 579         return;
 580       }
 581
 582
 583     //=======================================================
 584     // Get the output
 585     //=======================================================
 586     typename DeformationFieldType::Pointer deformationField=multiResolutionFilter->GetOutput();
 587
 588
 589     //=======================================================
 590     // Write the DVF
 591     //=======================================================
 592     typedef itk::ImageFileWriter< DeformationFieldType >  FieldWriterType;
 593     typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
 594     fieldWriter->SetInput( deformationField );
 595     fieldWriter->SetFileName( m_ArgsInfo.vf_arg );
 596     try
 597       {
 598         fieldWriter->Update();
 599       }
 600     catch( itk::ExceptionObject & excp )
 601       {
 602         std::cerr << "Exception thrown writing the DVF" << std::endl;
 603         std::cerr << excp << std::endl;
 604         return;
 605       }
 606
 607
 608     //=======================================================
 609     // Warp the moving image
 610     //=======================================================
 611     typedef itk::WarpImageFilter< MovingImageType, FixedImageType, DeformationFieldType >    WarpFilterType;
 612     typename WarpFilterType::Pointer warp = WarpFilterType::New();
 613
 614     warp->SetDisplacementField( deformationField );
 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