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