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