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