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