registration/clitkBLUTDIRGenericFilter.cxx

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