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