registration/clitkShapedBLUTSpatioTemporalDIRGenericFilter.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://oncora1.lyon.fnclcc.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ======================================================================-====*/
  18 #ifndef clitkShapedBLUTSpatioTemporalDIRGenericFilter_txx
  19 #define clitkShapedBLUTSpatioTemporalDIRGenericFilter_txx
  20
  21 /* =================================================
  22  * @file   clitkShapedBLUTSpatioTemporalDIRGenericFilter.txx
  23  * @author
  24  * @date
  25  *
  26  * @brief
  27  *
  28  ===================================================*/
  29
  30 #include "clitkShapedBLUTSpatioTemporalDIRGenericFilter.h"
  31
  32 namespace clitk
  33 {
  34
  35   //==============================================================================
  36   // Creating an observer class that allows output at each iteration
  37   //==============================================================================
  38   template <typename TRegistration, class args_info_clitkShapedBLUTSpatioTemporalDIR>
  39   class CommandIterationUpdate : public itk::Command
  40   {
  41   public:
  42     typedef  CommandIterationUpdate   Self;
  43     typedef  itk::Command             Superclass;
  44     typedef  itk::SmartPointer<Self>  Pointer;
  45     itkNewMacro( Self );
  46
  47     // Registration
  48     typedef   TRegistration                              RegistrationType;
  49     typedef   RegistrationType *                         RegistrationPointer;
  50
  51     // Transform
  52     typedef typename RegistrationType::FixedImageType FixedImageType;
  53     typedef typename FixedImageType::RegionType RegionType;
  54     itkStaticConstMacro(ImageDimension, unsigned int,FixedImageType::ImageDimension);
  55     typedef clitk::ShapedBLUTSpatioTemporalDeformableTransform<double, ImageDimension, ImageDimension> TransformType;
  56     typedef clitk::ShapedBLUTSpatioTemporalDeformableTransformInitializer<TransformType, FixedImageType> InitializerType;
  57     typedef typename InitializerType::CoefficientImageType CoefficientImageType;
  58     typedef itk::CastImageFilter<CoefficientImageType, CoefficientImageType> CastImageFilterType;
  59     typedef typename TransformType::ParametersType ParametersType;
  60     typedef typename InitializerType::Pointer InitializerPointer;
  61
  62   protected:
  63     CommandIterationUpdate() { m_IterationCounter=0;}
  64   public:
  65     typedef   clitk::GenericOptimizer<args_info_clitkShapedBLUTSpatioTemporalDIR>     OptimizerType;
  66     typedef   const OptimizerType   *           OptimizerPointer;
  67
  68     // Execute
  69     void Execute(itk::Object *caller, const itk::EventObject & event)
  70     {
  71       Execute( (const itk::Object *)caller, event);
  72     }
  73
  74     void Execute(const itk::Object * object, const itk::EventObject & event)
  75     {
  76       if( !(itk::IterationEvent().CheckEvent( &event )) )
  77         {
  78           return;
  79         }
  80
  81       // Output
  82       m_Optimizer->OutputIterationInfo();
  83       m_IterationCounter++;
  84
  85       // Write intermediate result
  86       if (m_ArgsInfo.current_given && (m_IterationCounter> m_ArgsInfo.current_arg) )
  87         {
  88           // Write && Reset
  89           writeImage<CoefficientImageType>(m_Initializer->GetTransform()->GetCoefficientImage(), m_ArgsInfo.coeff_arg, m_ArgsInfo.verbose_flag );
  90           m_IterationCounter=0;
  91         }
  92     }
  93
  94    // Members
  95     void SetOptimizer(OptimizerPointer o){m_Optimizer=o;}
  96     OptimizerPointer m_Optimizer;
  97
  98     void SetInitializer(InitializerPointer i){m_Initializer=i;}
  99     InitializerPointer m_Initializer;
 100
 101     void SetArgsInfo(args_info_clitkShapedBLUTSpatioTemporalDIR a){m_ArgsInfo=a;}
 102     args_info_clitkShapedBLUTSpatioTemporalDIR m_ArgsInfo;
 103
 104     int m_IterationCounter;
 105   };
 106
 107
 108   //==============================================================================
 109   //Creating an observer class that allows us to change parameters at subsequent levels
 110   //==============================================================================
 111   template <typename TRegistration, class args_info_clitkShapedBLUTSpatioTemporalDIR>
 112   class RegistrationInterfaceCommand : public itk::Command
 113   {
 114   public:
 115     typedef RegistrationInterfaceCommand   Self;
 116     typedef itk::Command             Superclass;
 117     typedef itk::SmartPointer<Self>  Pointer;
 118     itkNewMacro( Self );
 119   protected:
 120     RegistrationInterfaceCommand() {};
 121   public:
 122
 123     // Registration
 124     typedef   TRegistration                              RegistrationType;
 125     typedef   RegistrationType *                         RegistrationPointer;
 126
 127     // Transform
 128     typedef typename RegistrationType::FixedImageType FixedImageType;
 129     typedef typename FixedImageType::RegionType RegionType;
 130     itkStaticConstMacro(ImageDimension, unsigned int,FixedImageType::ImageDimension);
 131     typedef clitk::ShapedBLUTSpatioTemporalDeformableTransform<double, ImageDimension, ImageDimension> TransformType;
 132     typedef clitk::ShapedBLUTSpatioTemporalDeformableTransformInitializer<TransformType, FixedImageType> InitializerType;
 133     typedef typename InitializerType::CoefficientImageType CoefficientImageType;
 134     typedef itk::CastImageFilter<CoefficientImageType, CoefficientImageType> CastImageFilterType;
 135     typedef typename TransformType::ParametersType ParametersType;
 136     typedef typename InitializerType::Pointer InitializerPointer;
 137     typedef  typename  CommandIterationUpdate< TRegistration, args_info_clitkShapedBLUTSpatioTemporalDIR>::Pointer CommandIterationUpdatePointer;
 138
 139     // Optimizer
 140     typedef clitk::GenericOptimizer<args_info_clitkShapedBLUTSpatioTemporalDIR> GenericOptimizerType;
 141     typedef typename GenericOptimizerType::Pointer GenericOptimizerPointer;
 142
 143     // Metric
 144     typedef typename RegistrationType::FixedImageType    InternalImageType;
 145     typedef clitk::GenericMetric<args_info_clitkShapedBLUTSpatioTemporalDIR, InternalImageType, InternalImageType> GenericMetricType;
 146     typedef typename GenericMetricType::Pointer GenericMetricPointer;
 147
 148     // Two arguments are passed to the Execute() method: the first
 149     // is the pointer to the object which invoked the event and the
 150     // second is the event that was invoked.
 151     void Execute(itk::Object * object, const itk::EventObject & event)
 152     {
 153       if( !(itk::IterationEvent().CheckEvent( &event )) )
 154         {
 155           return;
 156         }
 157
 158       // Get the levels
 159       RegistrationPointer registration = dynamic_cast<RegistrationPointer>( object );
 160       unsigned int numberOfLevels=registration->GetNumberOfLevels();
 161       unsigned int currentLevel=registration->GetCurrentLevel()+1;
 162
 163       // Output the levels
 164       std::cout<<std::endl;
 165       std::cout<<"========================================"<<std::endl;
 166       std::cout<<"Starting resolution level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
 167       std::cout<<"========================================"<<std::endl;
 168       std::cout<<std::endl;
 169
 170       // Higher level?
 171       if (currentLevel>1)
 172         {
 173           // fixed image region pyramid
 174           typedef clitk::MultiResolutionPyramidRegionFilter<InternalImageType> FixedImageRegionPyramidType;
 175           typename FixedImageRegionPyramidType::Pointer fixedImageRegionPyramid=FixedImageRegionPyramidType::New();
 176           fixedImageRegionPyramid->SetRegion(m_MetricRegion);
 177           fixedImageRegionPyramid->SetSchedule(registration->GetFixedImagePyramid()->GetSchedule());
 178
 179           // Reinitialize the metric (!= number of samples)
 180           m_GenericMetric= GenericMetricType::New();
 181           m_GenericMetric->SetArgsInfo(m_ArgsInfo);
 182           m_GenericMetric->SetFixedImage(registration->GetFixedImagePyramid()->GetOutput(registration->GetCurrentLevel()));
 183           if (m_ArgsInfo.referenceMask_given)  m_GenericMetric->SetFixedImageMask(registration->GetMetric()->GetFixedImageMask());
 184           m_GenericMetric->SetFixedImageRegion(fixedImageRegionPyramid->GetOutput(registration->GetCurrentLevel()));
 185           typedef itk::ImageToImageMetric< InternalImageType, InternalImageType >  MetricType;
 186           typename  MetricType::Pointer metric=m_GenericMetric->GetMetricPointer();
 187           registration->SetMetric(metric);
 188
 189           // Get the current coefficient image and make a COPY
 190           typename itk::ImageDuplicator<CoefficientImageType>::Pointer caster=itk::ImageDuplicator<CoefficientImageType>::New();
 191           caster->SetInputImage(m_Initializer->GetTransform()->GetCoefficientImage());
 192           caster->Update();
 193           typename CoefficientImageType::Pointer currentCoefficientImage=caster->GetOutput();
 194
 195           // Write the intermediate result?
 196           if (m_ArgsInfo.intermediate_given>=numberOfLevels)
 197             writeImage<CoefficientImageType>(currentCoefficientImage, m_ArgsInfo.intermediate_arg[currentLevel-2], m_ArgsInfo.verbose_flag);
 198
 199           // Set the new transform properties
 200           m_Initializer->SetImage(registration->GetFixedImagePyramid()->GetOutput(currentLevel-1));
 201           if( m_Initializer->m_ControlPointSpacingIsGiven)
 202             m_Initializer->SetControlPointSpacing(m_Initializer->m_ControlPointSpacingArray[registration->GetCurrentLevel()]);
 203           if( m_Initializer->m_NumberOfControlPointsIsGiven)
 204             m_Initializer->SetNumberOfControlPointsInsideTheImage(m_Initializer->m_NumberOfControlPointsInsideTheImageArray[registration->GetCurrentLevel()]);
 205
 206           // Reinitialize the transform
 207           if (m_ArgsInfo.verbose_flag) std::cout<<"Initializing transform for level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
 208           m_Initializer->InitializeTransform();
 209           ParametersType* newParameters= new typename TransformType::ParametersType(m_Initializer->GetTransform()->GetNumberOfParameters());
 210
 211           // Reinitialize an Optimizer (!= number of parameters)
 212           m_GenericOptimizer = GenericOptimizerType::New();
 213           m_GenericOptimizer->SetArgsInfo(m_ArgsInfo);
 214           m_GenericOptimizer->SetMaximize(m_Maximize);
 215           m_GenericOptimizer->SetNumberOfParameters(m_Initializer->GetTransform()->GetNumberOfParameters());
 216           typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
 217           OptimizerType::Pointer optimizer = m_GenericOptimizer->GetOptimizerPointer();
 218           optimizer->AddObserver( itk::IterationEvent(), m_CommandIterationUpdate);
 219           registration->SetOptimizer(optimizer);
 220           m_CommandIterationUpdate->SetOptimizer(m_GenericOptimizer);
 221
 222           // Set the previous transform parameters to the registration
 223           // if(m_Initializer->m_Parameters!=NULL )delete m_Initializer->m_Parameters;
 224           m_Initializer->SetInitialParameters(currentCoefficientImage,*newParameters);
 225           registration->SetInitialTransformParametersOfNextLevel(*newParameters);
 226         }
 227     }
 228
 229     void Execute(const itk::Object * , const itk::EventObject & )
 230     { return; }
 231
 232
 233     // Members
 234     void SetInitializer(InitializerPointer i){m_Initializer=i;}
 235     InitializerPointer m_Initializer;
 236
 237     void SetArgsInfo(args_info_clitkShapedBLUTSpatioTemporalDIR a){m_ArgsInfo=a;}
 238     args_info_clitkShapedBLUTSpatioTemporalDIR m_ArgsInfo;
 239
 240     void SetCommandIterationUpdate(CommandIterationUpdatePointer c){m_CommandIterationUpdate=c;};
 241     CommandIterationUpdatePointer m_CommandIterationUpdate;
 242
 243     GenericOptimizerPointer m_GenericOptimizer;
 244     void SetMaximize(bool b){m_Maximize=b;}
 245     bool m_Maximize;
 246
 247     GenericMetricPointer m_GenericMetric;
 248     void SetMetricRegion(RegionType i){m_MetricRegion=i;}
 249     RegionType m_MetricRegion;
 250
 251
 252   };
 253
 254
 255   //==============================================================================
 256   // Update with the number of dimensions
 257   //==============================================================================
 258   template<unsigned int Dimension>
 259   void
 260   ShapedBLUTSpatioTemporalDIRGenericFilter::UpdateWithDim(std::string PixelType)
 261   {
 262     if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
 263
 264     if(PixelType == "short"){
 265       if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
 266       UpdateWithDimAndPixelType<Dimension, signed short>();
 267     }
 268     //    else if(PixelType == "unsigned_short"){
 269     //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
 270     //       UpdateWithDimAndPixelType<Dimension, unsigned short>();
 271     //     }
 272
 273     //     else if (PixelType == "unsigned_char"){
 274     //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
 275     //       UpdateWithDimAndPixelType<Dimension, unsigned char>();
 276     //     }
 277
 278     //     else if (PixelType == "char"){
 279     //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
 280     //       UpdateWithDimAndPixelType<Dimension, signed char>();
 281     //    }
 282     else {
 283       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
 284       UpdateWithDimAndPixelType<Dimension, float>();
 285     }
 286   }
 287
 288
 289   //==============================================================================
 290   // Update with the number of dimensions and pixeltype
 291   //==============================================================================
 292   template<unsigned int ImageDimension, class PixelType>
 293   void
 294   ShapedBLUTSpatioTemporalDIRGenericFilter::UpdateWithDimAndPixelType()
 295   {
 296
 297
 298     //=============================================================================
 299     //Input
 300     //=============================================================================
 301     bool threadsGiven=m_ArgsInfo.threads_given;
 302     int threads=m_ArgsInfo.threads_arg;
 303
 304     typedef itk::Image< PixelType, ImageDimension >  FixedImageType;
 305     typedef itk::Image< PixelType, ImageDimension >  MovingImageType;
 306     const unsigned int SpaceDimension = ImageDimension-1;
 307     typedef double TCoordRep;
 308
 309
 310     //=======================================================
 311     //Input
 312     //=======================================================
 313     typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
 314     typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
 315
 316     typename FixedImageReaderType::Pointer  fixedImageReader  = FixedImageReaderType::New();
 317     typename MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
 318
 319     fixedImageReader->SetFileName(  m_ArgsInfo.reference_arg );
 320     movingImageReader->SetFileName( m_ArgsInfo.target_arg );
 321     if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
 322     fixedImageReader->Update();
 323     movingImageReader->Update();
 324
 325     typename FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
 326     typename MovingImageType::Pointer movingImage =movingImageReader->GetOutput();
 327     typename FixedImageType::Pointer croppedFixedImage=fixedImage;
 328
 329
 330     //=======================================================
 331     // Regions
 332     //=======================================================
 333
 334     // The original input region
 335     typename FixedImageType::RegionType fixedImageRegion = fixedImage->GetLargestPossibleRegion();
 336
 337     // The transform region with respect to the input region:
 338     // where should the transform be DEFINED (depends on mask)
 339     typename FixedImageType::RegionType transformRegion = fixedImage->GetLargestPossibleRegion();
 340     typename FixedImageType::RegionType::SizeType transformRegionSize=transformRegion.GetSize();
 341     typename FixedImageType::RegionType::IndexType transformRegionIndex=transformRegion.GetIndex();
 342     typename FixedImageType::PointType transformRegionOrigin=fixedImage->GetOrigin();
 343
 344     // The metric region with respect to the extracted transform region:
 345     // where should the metric be CALCULATED (depends on transform)
 346     typename FixedImageType::RegionType metricRegion = fixedImage->GetLargestPossibleRegion();
 347     typename FixedImageType::RegionType::SizeType metricRegionSize=metricRegion.GetSize();
 348     typename FixedImageType::RegionType::IndexType metricRegionIndex=metricRegion.GetIndex();
 349     typename FixedImageType::PointType metricRegionOrigin=fixedImage->GetOrigin();
 350
 351     //===========================================================================
 352     // If given, we connect a mask to reference or target
 353     //============================================================================
 354     typedef itk::ImageMaskSpatialObject<  ImageDimension >   MaskType;
 355     typename MaskType::Pointer  fixedMask=NULL;
 356     if (m_ArgsInfo.referenceMask_given)
 357       {
 358         fixedMask= MaskType::New();
 359         typedef itk::Image< unsigned char, ImageDimension >   ImageMaskType;
 360         typedef itk::ImageFileReader< ImageMaskType >    MaskReaderType;
 361         typename MaskReaderType::Pointer  maskReader = MaskReaderType::New();
 362         maskReader->SetFileName(m_ArgsInfo.referenceMask_arg);
 363         try
 364           {
 365             maskReader->Update();
 366           }
 367         catch( itk::ExceptionObject & err )
 368           {
 369             std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
 370             std::cerr << err << std::endl;
 371             return;
 372           }
 373         if (m_Verbose)std::cout <<"Reference image mask was read..." <<std::endl;
 374
 375         // Set the image to the spatialObject
 376         fixedMask->SetImage( maskReader->GetOutput() );
 377
 378         // Find the bounding box of the "inside" label
 379         typedef itk::LabelStatisticsImageFilter<ImageMaskType, ImageMaskType> StatisticsImageFilterType;
 380         typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
 381         statisticsImageFilter->SetInput(maskReader->GetOutput());
 382         statisticsImageFilter->SetLabelInput(maskReader->GetOutput());
 383         statisticsImageFilter->Update();
 384         typename StatisticsImageFilterType::BoundingBoxType boundingBox = statisticsImageFilter->GetBoundingBox(1);
 385
 386         // Limit the transform region to the mask
 387         for (unsigned int i=0; i<ImageDimension; i++)
 388           {
 389             transformRegionIndex[i]=boundingBox[2*i];
 390             transformRegionSize[i]=boundingBox[2*i+1]-boundingBox[2*i]+1;
 391           }
 392         transformRegion.SetSize(transformRegionSize);
 393         transformRegion.SetIndex(transformRegionIndex);
 394         fixedImage->TransformIndexToPhysicalPoint(transformRegion.GetIndex(), transformRegionOrigin);
 395
 396         // Crop the fixedImage to the bounding box to facilitate multi-resolution
 397         typedef itk::ExtractImageFilter<FixedImageType,FixedImageType> ExtractImageFilterType;
 398         typename ExtractImageFilterType::Pointer extractImageFilter=ExtractImageFilterType::New();
 399         extractImageFilter->SetInput(fixedImage);
 400         extractImageFilter->SetExtractionRegion(transformRegion);
 401         extractImageFilter->Update();
 402         croppedFixedImage=extractImageFilter->GetOutput();
 403
 404         // Update the metric region
 405         metricRegion = croppedFixedImage->GetLargestPossibleRegion();
 406         metricRegionIndex=metricRegion.GetIndex();
 407         metricRegionSize=metricRegion.GetSize();
 408         croppedFixedImage->TransformIndexToPhysicalPoint(metricRegionIndex, metricRegionOrigin);
 409
 410         // Set start index to zero (with respect to croppedFixedImage/transform region)
 411         metricRegionIndex.Fill(0);
 412         metricRegion.SetIndex(metricRegionIndex);
 413         croppedFixedImage->SetRegions(metricRegion);
 414         croppedFixedImage->SetOrigin(metricRegionOrigin);
 415
 416       }
 417
 418     typedef itk::ImageMaskSpatialObject<  ImageDimension >   MaskType;
 419     typename MaskType::Pointer  movingMask=NULL;
 420     if (m_ArgsInfo.targetMask_given)
 421       {
 422         movingMask= MaskType::New();
 423         typedef itk::Image< unsigned char, ImageDimension >   ImageMaskType;
 424         typedef itk::ImageFileReader< ImageMaskType >    MaskReaderType;
 425         typename MaskReaderType::Pointer  maskReader = MaskReaderType::New();
 426         maskReader->SetFileName(m_ArgsInfo.targetMask_arg);
 427         try
 428           {
 429             maskReader->Update();
 430           }
 431         catch( itk::ExceptionObject & err )
 432           {
 433             std::cerr << "ExceptionObject caught !" << std::endl;
 434             std::cerr << err << std::endl;
 435           }
 436         if (m_Verbose)std::cout <<"Target image mask was read..." <<std::endl;
 437
 438         movingMask->SetImage( maskReader->GetOutput() );
 439       }
 440
 441
 442     //=======================================================
 443     // Output Regions
 444     //=======================================================
 445
 446     if (m_Verbose)
 447       {
 448         // Fixed image region
 449         std::cout<<"The fixed image has its origin at "<<fixedImage->GetOrigin()<<std::endl
 450                  <<"The fixed image region starts at index "<<fixedImageRegion.GetIndex()<<std::endl
 451                  <<"The fixed image region has size "<< fixedImageRegion.GetSize()<<std::endl;
 452
 453         // Transform region
 454         std::cout<<"The transform has its origin at "<<transformRegionOrigin<<std::endl
 455                  <<"The transform region will start at index "<<transformRegion.GetIndex()<<std::endl
 456                  <<"The transform region has size "<< transformRegion.GetSize()<<std::endl;
 457
 458         // Metric region
 459         std::cout<<"The metric region has its origin at "<<metricRegionOrigin<<std::endl
 460                  <<"The metric region will start at index "<<metricRegion.GetIndex()<<std::endl
 461                  <<"The metric region has size "<< metricRegion.GetSize()<<std::endl;
 462
 463       }
 464
 465
 466     //=======================================================
 467     // Pyramids (update them for transform initializer)
 468     //=======================================================
 469     typedef clitk::RecursiveSpatioTemporalMultiResolutionPyramidImageFilter< FixedImageType, FixedImageType>    FixedImagePyramidType;
 470     typedef clitk::RecursiveSpatioTemporalMultiResolutionPyramidImageFilter< MovingImageType, MovingImageType>    MovingImagePyramidType;
 471     typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
 472     typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
 473     fixedImagePyramid->SetUseShrinkImageFilter(false);
 474     fixedImagePyramid->SetInput(croppedFixedImage);
 475     fixedImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
 476     movingImagePyramid->SetUseShrinkImageFilter(false);
 477     movingImagePyramid->SetInput(movingImage);
 478     movingImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
 479     if (m_Verbose) std::cout<<"Creating the image pyramid..."<<std::endl;
 480     fixedImagePyramid->Update();
 481     movingImagePyramid->Update();
 482     typedef clitk::MultiResolutionPyramidRegionFilter<FixedImageType> FixedImageRegionPyramidType;
 483     typename FixedImageRegionPyramidType::Pointer fixedImageRegionPyramid=FixedImageRegionPyramidType::New();
 484     fixedImageRegionPyramid->SetRegion(metricRegion);
 485     fixedImageRegionPyramid->SetSchedule(fixedImagePyramid->GetSchedule());
 486
 487
 488
 489
 490     //     //=======================================================
 491     //     // Rigid Transform
 492     //     //=======================================================
 493     //     typedef itk::Euler3DTransform <double> RigidTransformType;
 494     //     RigidTransformType::Pointer rigidTransform;
 495     //     if (m_ArgsInfo.rigid_given)
 496     //       {
 497     //  rigidTransform=RigidTransformType::New();
 498     //  itk::Matrix<double,4,4> rigidTransformMatrix=clitk::ReadMatrix3D(m_ArgsInfo.rigid_arg);
 499
 500     //  //Set the rotation
 501     //  itk::Matrix<double,3,3> finalRotation = clitk::GetRotationalPartMatrix3D(rigidTransformMatrix);
 502     //  rigidTransform->SetMatrix(finalRotation);
 503
 504     //  //Set the translation
 505     //  itk::Vector<double,3> finalTranslation = clitk::GetTranslationPartMatrix3D(rigidTransformMatrix);
 506     //  rigidTransform->SetTranslation(finalTranslation);
 507
 508     //       }
 509
 510
 511     //=======================================================
 512     // ShapedBSplineSpatioTemporal Transform
 513     //=======================================================
 514     typedef  clitk::ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep,ImageDimension, ImageDimension > TransformType;
 515     typename TransformType::Pointer transform= TransformType::New();
 516     transform->SetTransformShape(m_ArgsInfo.shape_arg);
 517     if (fixedMask) transform->SetMask( fixedMask );
 518     // if (rigidTransform) transform->SetBulkTransform( rigidTransform );
 519
 520     //-------------------------------------------------------------------------
 521     // The transform initializer
 522     //-------------------------------------------------------------------------
 523     typedef clitk::ShapedBLUTSpatioTemporalDeformableTransformInitializer< TransformType,FixedImageType> InitializerType;
 524     typename InitializerType::Pointer initializer = InitializerType::New();
 525     initializer->SetVerbose(m_Verbose);
 526     initializer->SetImage(fixedImagePyramid->GetOutput(0));
 527     initializer->SetTransform(transform);
 528
 529     //-------------------------------------------------------------------------
 530     // Order
 531     //-------------------------------------------------------------------------
 532     typename FixedImageType::RegionType::SizeType splineOrders ;
 533     splineOrders.Fill(3);
 534     if (m_ArgsInfo.order_given)
 535       for(unsigned int i=0; i<ImageDimension;i++)
 536         splineOrders[i]=m_ArgsInfo.order_arg[i];
 537     if (m_Verbose) std::cout<<"Setting the spline orders  to "<<splineOrders<<"..."<<std::endl;
 538     initializer->SetSplineOrders(splineOrders);
 539
 540     //-------------------------------------------------------------------------
 541     // Levels
 542     //-------------------------------------------------------------------------
 543
 544     // Spacing
 545     if (m_ArgsInfo.spacing_given)
 546       {
 547         initializer->m_ControlPointSpacingArray.resize(m_ArgsInfo.levels_arg);
 548         initializer->SetControlPointSpacing(m_ArgsInfo.spacing_arg);
 549         initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1]=initializer->m_ControlPointSpacing;
 550         if (m_Verbose) std::cout<<"Using a control point spacing of "<<initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1]
 551                                 <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
 552
 553         for (int i=1; i<m_ArgsInfo.levels_arg; i++ )
 554           {
 555             initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1-i]=initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-i]*2;
 556             initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1-i][ImageDimension-1]=2;
 557             if (m_Verbose) std::cout<<"Using a control point spacing of "<<initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1-i]
 558                                     <<" at level "<<m_ArgsInfo.levels_arg-i<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
 559           }
 560
 561       }
 562
 563     // Control
 564     if (m_ArgsInfo.control_given)
 565       {
 566         initializer->m_NumberOfControlPointsInsideTheImageArray.resize(m_ArgsInfo.levels_arg);
 567         initializer->SetNumberOfControlPointsInsideTheImage(m_ArgsInfo.control_arg);
 568         initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1]=initializer->m_NumberOfControlPointsInsideTheImage;
 569         if (m_Verbose) std::cout<<"Using "<< initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1]<<"control points inside the image"
 570                                 <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
 571
 572         for (int i=1; i<m_ArgsInfo.levels_arg; i++ )
 573           {
 574             for(unsigned int j=0;j<ImageDimension-1;j++)
 575               initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i][j]=ceil ((double)initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-i][j]/2.);
 576               //            initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i]=ceil ((double)initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-i]/2.);
 577             if (m_Verbose) std::cout<<"Using "<< initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i]<<"control points inside the image"
 578                                     <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
 579
 580           }
 581       }
 582
 583     // Inialize on the first level
 584     if (m_ArgsInfo.verbose_flag) std::cout<<"Initializing transform for level 1 of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
 585     if (m_ArgsInfo.spacing_given) initializer->SetControlPointSpacing(  initializer->m_ControlPointSpacingArray[0]);
 586     if (m_ArgsInfo.control_given) initializer->SetNumberOfControlPointsInsideTheImage(initializer->m_NumberOfControlPointsInsideTheImageArray[0]);
 587     if (m_ArgsInfo.samplingFactor_given) initializer->SetSamplingFactors(m_ArgsInfo.samplingFactor_arg);
 588
 589     // Initialize
 590     initializer->InitializeTransform();
 591
 592     //-------------------------------------------------------------------------
 593     // Initial parameters (passed by reference)
 594     //-------------------------------------------------------------------------
 595     typedef typename TransformType::ParametersType     ParametersType;
 596     const unsigned int numberOfParameters =    transform->GetNumberOfParameters();
 597     ParametersType parameters(numberOfParameters);
 598     parameters.Fill( 0.0 );
 599     transform->SetParameters( parameters );
 600     if (m_ArgsInfo.initCoeff_given) initializer->SetInitialParameters(m_ArgsInfo.initCoeff_arg, parameters);
 601
 602
 603     //=======================================================
 604     // Interpolator
 605     //=======================================================
 606     typedef clitk::GenericInterpolator<args_info_clitkShapedBLUTSpatioTemporalDIR, FixedImageType,TCoordRep > GenericInterpolatorType;
 607     typename   GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
 608     genericInterpolator->SetArgsInfo(m_ArgsInfo);
 609     typedef itk::InterpolateImageFunction< FixedImageType, TCoordRep >  InterpolatorType;
 610     typename  InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
 611
 612
 613     //=======================================================
 614     // Metric
 615     //=======================================================
 616     typedef clitk::GenericMetric< args_info_clitkShapedBLUTSpatioTemporalDIR, FixedImageType,MovingImageType > GenericMetricType;
 617     typename GenericMetricType::Pointer genericMetric=GenericMetricType::New();
 618     genericMetric->SetArgsInfo(m_ArgsInfo);
 619     genericMetric->SetFixedImage(fixedImagePyramid->GetOutput(0));
 620     if (fixedMask) genericMetric->SetFixedImageMask(fixedMask);
 621     genericMetric->SetFixedImageRegion(fixedImageRegionPyramid->GetOutput(0));
 622     typedef itk::ImageToImageMetric< FixedImageType, MovingImageType >  MetricType;
 623     typename  MetricType::Pointer metric=genericMetric->GetMetricPointer();
 624     if (movingMask) metric->SetMovingImageMask(movingMask);
 625
 626 #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
 627     if (threadsGiven) metric->SetNumberOfThreads( threads );
 628 #else
 629     if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
 630 #endif
 631
 632
 633     //=======================================================
 634     // Optimizer
 635     //=======================================================
 636     typedef clitk::GenericOptimizer<args_info_clitkShapedBLUTSpatioTemporalDIR> GenericOptimizerType;
 637     GenericOptimizerType::Pointer genericOptimizer = GenericOptimizerType::New();
 638     genericOptimizer->SetArgsInfo(m_ArgsInfo);
 639     genericOptimizer->SetMaximize(genericMetric->GetMaximize());
 640     genericOptimizer->SetNumberOfParameters(transform->GetNumberOfParameters());
 641     typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
 642     OptimizerType::Pointer optimizer = genericOptimizer->GetOptimizerPointer();
 643
 644
 645     //=======================================================
 646     // Registration
 647     //=======================================================
 648     typedef clitk::SpatioTemporalMultiResolutionImageRegistrationMethod<  FixedImageType, MovingImageType >    RegistrationType;
 649     typename RegistrationType::Pointer   registration  = RegistrationType::New();
 650     registration->SetMetric(        metric        );
 651     registration->SetOptimizer(     optimizer     );
 652     registration->SetInterpolator(  interpolator  );
 653     registration->SetTransform (transform);
 654     if(threadsGiven) registration->SetNumberOfThreads(threads);
 655     registration->SetFixedImage(  croppedFixedImage   );
 656     registration->SetMovingImage(  movingImage   );
 657     registration->SetFixedImageRegion( metricRegion );
 658     registration->SetFixedImagePyramid( fixedImagePyramid );
 659     registration->SetMovingImagePyramid( movingImagePyramid );
 660     registration->SetInitialTransformParameters( transform->GetParameters() );
 661     registration->SetNumberOfLevels( m_ArgsInfo.levels_arg );
 662     if (m_Verbose) std::cout<<"Setting the number of resolution levels to "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
 663
 664
 665     //================================================================================================
 666     // Observers
 667     //================================================================================================
 668     if (m_Verbose)
 669       {
 670         // Output iteration info
 671         typedef CommandIterationUpdate <RegistrationType,args_info_clitkShapedBLUTSpatioTemporalDIR> CommandIterationUpdateType;
 672         typename CommandIterationUpdateType::Pointer observer = CommandIterationUpdateType::New();
 673         observer->SetInitializer(initializer);
 674         observer->SetArgsInfo(m_ArgsInfo);
 675         observer->SetOptimizer(genericOptimizer);
 676         optimizer->AddObserver( itk::IterationEvent(), observer );
 677
 678         // Output level info
 679         typedef RegistrationInterfaceCommand<RegistrationType,args_info_clitkShapedBLUTSpatioTemporalDIR> CommandType;
 680         typename CommandType::Pointer command = CommandType::New();
 681         command->SetInitializer(initializer);
 682         command->SetArgsInfo(m_ArgsInfo);
 683         command->SetCommandIterationUpdate(observer);
 684         command->SetMaximize(genericMetric->GetMaximize());
 685         command->SetMetricRegion(metricRegion);
 686         registration->AddObserver( itk::IterationEvent(), command );
 687       }
 688
 689
 690     //=======================================================
 691     // Let's go
 692     //=======================================================
 693     if (m_Verbose) std::cout << std::endl << "Starting Registration" << std::endl;
 694
 695     try
 696       {
 697         registration->StartRegistration();
 698       }
 699     catch( itk::ExceptionObject & err )
 700       {
 701         std::cerr << "ExceptionObject caught while registering!" << std::endl;
 702         std::cerr << err << std::endl;
 703         return;
 704       }
 705
 706
 707     //=======================================================
 708     // Get the result
 709     //=======================================================
 710     OptimizerType::ParametersType finalParameters =  registration->GetLastTransformParameters();
 711     transform->SetParameters( finalParameters );
 712     if (m_Verbose)
 713       {
 714         std::cout<<"Stop condition description: "
 715                  <<registration->GetOptimizer()->GetStopConditionDescription()<<std::endl;
 716       }
 717
 718
 719     //=======================================================
 720     // Get the BSplineSpatioTemporal coefficient images and write them
 721     //=======================================================
 722     if (m_ArgsInfo.coeff_given)
 723       {
 724         typedef typename TransformType::CoefficientImageType CoefficientImageType;
 725         typename CoefficientImageType::Pointer coefficientImage =transform->GetCoefficientImage();
 726         typedef itk::ImageFileWriter<CoefficientImageType> CoeffWriterType;
 727         typename CoeffWriterType::Pointer coeffWriter=CoeffWriterType::New();
 728         coeffWriter->SetInput(coefficientImage);
 729         coeffWriter->SetFileName(m_ArgsInfo.coeff_arg);
 730         coeffWriter->Update();
 731
 732       }
 733     if (m_ArgsInfo.padCoeff_given)
 734       {
 735         typedef  itk::Image<itk::Vector<TCoordRep, SpaceDimension>, ImageDimension> CoefficientImageType;
 736         typename CoefficientImageType::Pointer coefficientImage =transform->GetPaddedCoefficientImage();
 737         typedef itk::ImageFileWriter<CoefficientImageType> CoeffWriterType;
 738         typename CoeffWriterType::Pointer coeffWriter=CoeffWriterType::New();
 739         coeffWriter->SetInput(coefficientImage);
 740         coeffWriter->SetFileName(m_ArgsInfo.padCoeff_arg);
 741         coeffWriter->Update();
 742       }
 743
 744
 745     //=======================================================
 746     // Generate the DVF (4D with 3 comps and 4D with 4 comps)
 747     //=======================================================
 748     typedef itk::Vector< float, SpaceDimension >  DisplacementType;
 749     typedef itk::Vector< float, ImageDimension >  Displacement4DType;
 750     typedef itk::Image< DisplacementType, ImageDimension >  DeformationFieldType;
 751     typedef itk::Image< Displacement4DType, ImageDimension >  DeformationField4DType;
 752
 753     typename DeformationFieldType::Pointer field = DeformationFieldType::New();
 754     typename DeformationField4DType::Pointer field4D = DeformationField4DType::New();
 755     field->SetRegions( fixedImageRegion );
 756     field->SetOrigin( fixedImage->GetOrigin() );
 757     field->SetSpacing( fixedImage->GetSpacing() );
 758     field->SetDirection( fixedImage->GetDirection() );
 759     field->Allocate();
 760     field4D->SetRegions( fixedImageRegion );
 761     field4D->SetOrigin( fixedImage->GetOrigin() );
 762     field4D->SetSpacing( fixedImage->GetSpacing() );
 763     field4D->SetDirection( fixedImage->GetDirection() );
 764     field4D->Allocate();
 765
 766     typedef itk::ImageRegionIteratorWithIndex< DeformationFieldType > FieldIterator;
 767     typedef itk::ImageRegionIteratorWithIndex< DeformationField4DType > Field4DIterator;
 768     FieldIterator fi( field, fixedImageRegion );
 769     Field4DIterator fi4D( field4D, fixedImageRegion );
 770     fi.GoToBegin();
 771     fi4D.GoToBegin();
 772
 773     typename TransformType::InputPointType  fixedPoint;
 774     typename TransformType::OutputPointType movingPoint;
 775     typename DeformationFieldType::IndexType index;
 776
 777     DisplacementType displacement;
 778     Displacement4DType displacement4D;
 779     displacement4D[ImageDimension-1]=0;
 780     while( ! fi.IsAtEnd() )
 781       {
 782         index = fi.GetIndex();
 783         field->TransformIndexToPhysicalPoint( index, fixedPoint );
 784         movingPoint = transform->TransformPoint( fixedPoint );
 785         for (unsigned int i=0; i<SpaceDimension; i++)
 786           {
 787             displacement[i] = movingPoint[i] - fixedPoint[i];
 788             displacement4D[i] = movingPoint[i] - fixedPoint[i];
 789           }
 790         fi.Set( displacement );
 791         fi4D.Set(displacement4D);
 792         ++fi;
 793         ++fi4D;
 794       }
 795
 796
 797     //=======================================================
 798     // Write the DVF
 799     //=======================================================
 800     typedef itk::ImageFileWriter< DeformationFieldType >  FieldWriterType;
 801     typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
 802     fieldWriter->SetFileName( m_ArgsInfo.vf_arg );
 803     fieldWriter->SetInput( field );
 804     try
 805       {
 806         fieldWriter->Update();
 807       }
 808     catch( itk::ExceptionObject & excp )
 809       {
 810         std::cerr << "Exception thrown writing the DVF" << std::endl;
 811         std::cerr << excp << std::endl;
 812         return;
 813       }
 814
 815
 816     //=======================================================
 817     // Resample the moving image
 818     //=======================================================
 819     typedef itk::WarpImageFilter< MovingImageType, FixedImageType, DeformationField4DType >    WarpFilterType;
 820     typename WarpFilterType::Pointer warp = WarpFilterType::New();
 821
 822     warp->SetDeformationField( field4D );
 823     warp->SetInput( movingImageReader->GetOutput() );
 824     warp->SetOutputOrigin(  fixedImage->GetOrigin() );
 825     warp->SetOutputSpacing( fixedImage->GetSpacing() );
 826     warp->SetOutputDirection( fixedImage->GetDirection() );
 827     warp->SetEdgePaddingValue( 0.0 );
 828     warp->Update();
 829
 830
 831     //=======================================================
 832     // Write the warped image
 833     //=======================================================
 834     typedef itk::ImageFileWriter< FixedImageType >  WriterType;
 835     typename WriterType::Pointer      writer =  WriterType::New();
 836     writer->SetFileName( m_ArgsInfo.output_arg );
 837     writer->SetInput( warp->GetOutput()    );
 838
 839     try
 840       {
 841         writer->Update();
 842       }
 843     catch( itk::ExceptionObject & err )
 844       {
 845         std::cerr << "ExceptionObject caught !" << std::endl;
 846         std::cerr << err << std::endl;
 847         return;
 848       }
 849
 850
 851     //=======================================================
 852     // Calculate the difference after the deformable transform
 853     //=======================================================
 854     typedef clitk::DifferenceImageFilter<  FixedImageType, FixedImageType> DifferenceFilterType;
 855     if (m_ArgsInfo.after_given)
 856       {
 857         typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
 858         difference->SetValidInput( fixedImage );
 859         difference->SetTestInput( warp->GetOutput() );
 860
 861         try
 862           {
 863             difference->Update();
 864           }
 865         catch( itk::ExceptionObject & err )
 866           {
 867             std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
 868             std::cerr << err << std::endl;
 869             return;
 870           }
 871
 872         typename WriterType::Pointer differenceWriter=WriterType::New();
 873         differenceWriter->SetInput(difference->GetOutput());
 874         differenceWriter->SetFileName(m_ArgsInfo.after_arg);
 875         differenceWriter->Update();
 876
 877       }
 878
 879
 880     //=======================================================
 881     // Calculate the difference before the deformable transform
 882     //=======================================================
 883     if( m_ArgsInfo.before_given )
 884       {
 885
 886         typename FixedImageType::Pointer moving=FixedImageType::New();
 887         if (m_ArgsInfo.rigid_given)
 888           {
 889             typedef itk::ResampleImageFilter<MovingImageType, FixedImageType> ResamplerType;
 890             typename ResamplerType::Pointer resampler=ResamplerType::New();
 891             resampler->SetInput(movingImage);
 892             resampler->SetOutputOrigin(fixedImage->GetOrigin());
 893             resampler->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
 894             resampler->SetOutputSpacing(fixedImage->GetSpacing());
 895             resampler->SetDefaultPixelValue( 0. );
 896             //resampler->SetTransform(rigidTransform);
 897             resampler->Update();
 898             moving=resampler->GetOutput();
 899           }
 900         else
 901           moving=movingImage;
 902
 903         typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
 904         difference->SetValidInput( fixedImage );
 905         difference->SetTestInput( moving );
 906
 907         try
 908           {
 909             difference->Update();
 910           }
 911         catch( itk::ExceptionObject & err )
 912           {
 913             std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
 914             std::cerr << err << std::endl;
 915             return;
 916           }
 917
 918         typename WriterType::Pointer differenceWriter=WriterType::New();
 919         writer->SetFileName( m_ArgsInfo.before_arg  );
 920         writer->SetInput( difference->GetOutput()  );
 921         writer->Update( );
 922       }
 923
 924     return;
 925
 926   }
 927 }//end clitk
 928
 929 #endif //#define clitkShapedBLUTSpatioTemporalDIRGenericFilter_txx