registration/clitkBSplineDeformableRegistrationGenericFilter.txx

   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 __clitkBSplineDeformableRegistrationGenericFilter_txx
  19 #define __clitkBSplineDeformableRegistrationGenericFilter_txx
  20 #include "clitkBSplineDeformableRegistrationGenericFilter.h"
  21
  22
  23 namespace clitk
  24 {
  25
  26   //==============================================================================
  27   //Creating an observer class that allows us to change parameters at subsequent levels
  28   //==============================================================================
  29   template <typename TRegistration>
  30   class RegistrationInterfaceCommand : public itk::Command
  31   {
  32   public:
  33     typedef RegistrationInterfaceCommand   Self;
  34
  35
  36     typedef itk::Command             Superclass;
  37     typedef itk::SmartPointer<Self>  Pointer;
  38     itkNewMacro( Self );
  39   protected:
  40     RegistrationInterfaceCommand() {};
  41   public:
  42     typedef   TRegistration                              RegistrationType;
  43     typedef   RegistrationType *                         RegistrationPointer;
  44
  45     // Two arguments are passed to the Execute() method: the first
  46     // is the pointer to the object which invoked the event and the
  47     // second is the event that was invoked.
  48     void Execute(itk::Object * object, const itk::EventObject & event)
  49     {
  50       if( !(itk::IterationEvent().CheckEvent( &event )) )
  51         {
  52           return;
  53         }
  54       RegistrationPointer registration = dynamic_cast<RegistrationPointer>( object );
  55       unsigned int numberOfLevels=registration->GetNumberOfLevels();
  56       unsigned int currentLevel=registration->GetCurrentLevel()+1;
  57       std::cout<<std::endl;
  58       std::cout<<"========================================"<<std::endl;
  59       std::cout<<"Starting resolution level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
  60       std::cout<<"========================================"<<std::endl;
  61       std::cout<<std::endl;
  62     }
  63
  64     void Execute(const itk::Object * , const itk::EventObject & )
  65     { return; }
  66
  67   };
  68
  69
  70   //==============================================================================
  71   // Creating an observer class that allows output at each iteration
  72   //==============================================================================
  73   class CommandIterationUpdate : public itk::Command
  74   {
  75   public:
  76     typedef  CommandIterationUpdate   Self;
  77     typedef  itk::Command             Superclass;
  78     typedef  itk::SmartPointer<Self>  Pointer;
  79     itkNewMacro( Self );
  80   protected:
  81     CommandIterationUpdate() {};
  82   public:
  83     typedef   clitk::GenericOptimizer<args_info_clitkBSplineDeformableRegistration>     OptimizerType;
  84     typedef   const OptimizerType   *           OptimizerPointer;
  85
  86     // We set the generic optimizer
  87     void SetOptimizer(OptimizerPointer o){m_Optimizer=o;}
  88
  89     // Execute
  90     void Execute(itk::Object *caller, const itk::EventObject & event)
  91     {
  92       Execute( (const itk::Object *)caller, event);
  93     }
  94
  95     void Execute(const itk::Object * object, const itk::EventObject & event)
  96     {
  97       if( !(itk::IterationEvent().CheckEvent( &event )) )
  98         {
  99           return;
 100         }
 101
 102       m_Optimizer->OutputIterationInfo();
 103     }
 104
 105     OptimizerPointer m_Optimizer;
 106   };
 107
 108
 109   //==============================================================================
 110   // Update with the number of dimensions
 111   //==============================================================================
 112   template<unsigned int Dimension>
 113   void BSplineDeformableRegistrationGenericFilter::UpdateWithDim(std::string PixelType)
 114   {
 115
 116     if (m_Verbose) std::cout  << "Images were detected to be "<< Dimension << "D and " << PixelType << "..." << std::endl;
 117
 118     if(PixelType == "short"){
 119       if (m_Verbose) std::cout  << "Launching warp in "<< Dimension <<"D and signed short..." << std::endl;
 120       UpdateWithDimAndPixelType<Dimension, signed short>();
 121     }
 122     //    else if(PixelType == "unsigned_short"){
 123     //       if (m_Verbose) std::cout  << "Launching warp in "<< Dimension <<"D and unsigned_short..." << std::endl;
 124     //       UpdateWithDimAndPixelType<Dimension, unsigned short>();
 125     //     }
 126
 127     //     else if (PixelType == "unsigned_char"){
 128     //       if (m_Verbose) std::cout  << "Launching warp in "<< Dimension <<"D and unsigned_char..." << std::endl;
 129     //       UpdateWithDimAndPixelType<Dimension, unsigned char>();
 130     //     }
 131
 132     //     else if (PixelType == "char"){
 133     //       if (m_Verbose) std::cout  << "Launching warp in "<< Dimension <<"D and signed_char..." << std::endl;
 134     //       UpdateWithDimAndPixelType<Dimension, signed char>();
 135     //    }
 136     else {
 137       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
 138       UpdateWithDimAndPixelType<Dimension, float>();
 139     }
 140   }
 141
 142
 143
 144   //==============================================================================
 145   // Update with the number of dimensions and pixeltype
 146   //==============================================================================
 147   template<unsigned int ImageDimension, class PixelType>
 148   void BSplineDeformableRegistrationGenericFilter::UpdateWithDimAndPixelType()
 149   {
 150
 151     //=======================================================
 152     // Run-time
 153     //=======================================================
 154     bool threadsGiven=m_ArgsInfo.threads_given;
 155     int threads=m_ArgsInfo.threads_arg;
 156
 157     typedef itk::Image< PixelType, ImageDimension >  FixedImageType;
 158     typedef itk::Image< PixelType, ImageDimension >  MovingImageType;
 159     const unsigned int SpaceDimension = ImageDimension;
 160     typedef double TCoordRep;
 161
 162
 163     //=======================================================
 164     //Input
 165     //=======================================================
 166     typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
 167     typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
 168
 169     typename FixedImageReaderType::Pointer  fixedImageReader  = FixedImageReaderType::New();
 170     typename MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
 171
 172     fixedImageReader->SetFileName(  m_ArgsInfo.reference_arg );
 173     movingImageReader->SetFileName( m_ArgsInfo.target_arg );
 174     if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
 175     fixedImageReader->Update();
 176     movingImageReader->Update();
 177
 178     typename FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
 179     typename MovingImageType::Pointer movingImage =movingImageReader->GetOutput();
 180     typename FixedImageType::RegionType fixedImageRegion = fixedImage->GetLargestPossibleRegion();
 181
 182     // The metric region: where should the metric be CALCULATED (depends on mask)
 183     typename FixedImageType::RegionType metricRegion = fixedImage->GetLargestPossibleRegion();
 184     typename FixedImageType::RegionType::SizeType metricRegionSize=metricRegion.GetSize();
 185     typename FixedImageType::RegionType::IndexType metricRegionIndex=metricRegion.GetIndex();
 186     typename FixedImageType::PointType metricRegionOrigin=fixedImage->GetOrigin();
 187
 188     // The transform region: where should the transform be DEFINED (depends on mask)
 189     typename FixedImageType::RegionType transformRegion = fixedImage->GetLargestPossibleRegion();
 190     typename FixedImageType::RegionType::SizeType transformRegionSize=transformRegion.GetSize();
 191     typename FixedImageType::RegionType::IndexType transformRegionIndex=transformRegion.GetIndex();
 192     typename FixedImageType::PointType transformRegionOrigin=fixedImage->GetOrigin();
 193
 194
 195     //=======================================================
 196     // If given, we connect a mask to the fixed image
 197     //======================================================
 198     typedef itk::ImageMaskSpatialObject<  ImageDimension >   MaskType;
 199     typename MaskType::Pointer  spatialObjectMask=NULL;
 200
 201     if (m_ArgsInfo.mask_given)
 202       {
 203         typedef itk::Image< unsigned char, ImageDimension >   ImageMaskType;
 204         typedef itk::ImageFileReader< ImageMaskType >    MaskReaderType;
 205         typename MaskReaderType::Pointer  maskReader = MaskReaderType::New();
 206         maskReader->SetFileName(m_ArgsInfo.mask_arg);
 207
 208         try
 209           {
 210             maskReader->Update();
 211           }
 212         catch( itk::ExceptionObject & err )
 213           {
 214             std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
 215             std::cerr << err << std::endl;
 216             return;
 217           }
 218         if (m_Verbose)std::cout <<"Reference image mask was read..." <<std::endl;
 219
 220
 221         // Set the image to the spatialObject
 222         spatialObjectMask = MaskType::New();
 223         spatialObjectMask->SetImage( maskReader->GetOutput() );
 224
 225         // Find the bounding box of the "inside" label
 226         typedef itk::LabelStatisticsImageFilter<ImageMaskType, ImageMaskType> StatisticsImageFilterType;
 227         typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
 228         statisticsImageFilter->SetInput(maskReader->GetOutput());
 229         statisticsImageFilter->SetLabelInput(maskReader->GetOutput());
 230         statisticsImageFilter->Update();
 231         typename StatisticsImageFilterType::BoundingBoxType boundingBox = statisticsImageFilter->GetBoundingBox(1);
 232
 233         // Limit the transform region to the mask
 234         for (unsigned int i=0; i<ImageDimension; i++)
 235           {
 236             transformRegionIndex[i]=boundingBox[2*i];
 237             transformRegionSize[i]=boundingBox[2*i+1]-boundingBox[2*i]+1;
 238           }
 239         transformRegion.SetSize(transformRegionSize);
 240         transformRegion.SetIndex(transformRegionIndex);
 241         fixedImage->TransformIndexToPhysicalPoint(transformRegion.GetIndex(), transformRegionOrigin);
 242
 243         // Limit the metric region to the mask
 244         metricRegion=transformRegion;
 245         fixedImage->TransformIndexToPhysicalPoint(metricRegion.GetIndex(), metricRegionOrigin);
 246
 247      }
 248
 249
 250
 251     //=======================================================
 252     // Regions
 253     //=====================================================
 254     if (m_Verbose)
 255       {
 256         // Fixed image region
 257         std::cout<<"The fixed image has its origin at "<<fixedImage->GetOrigin()<<std::endl
 258                  <<"The fixed image region starts at index "<<fixedImageRegion.GetIndex()<<std::endl
 259                  <<"The fixed image region has size "<< fixedImageRegion.GetSize()<<std::endl;
 260
 261         // Transform region
 262         std::cout<<"The transform has its origin at "<<transformRegionOrigin<<std::endl
 263                  <<"The transform region will start at index "<<transformRegion.GetIndex()<<std::endl
 264                  <<"The transform region has size "<< transformRegion.GetSize()<<std::endl;
 265
 266         // Metric region
 267         std::cout<<"The metric region has its origin at "<<metricRegionOrigin<<std::endl
 268                  <<"The metric region will start at index "<<metricRegion.GetIndex()<<std::endl
 269                  <<"The metric region has size "<< metricRegion.GetSize()<<std::endl;
 270
 271       }
 272
 273
 274     //=======================================================
 275     //Pyramids
 276     //=======================================================
 277     typedef itk::RecursiveMultiResolutionPyramidImageFilter< FixedImageType, FixedImageType>    FixedImagePyramidType;
 278     typedef itk::RecursiveMultiResolutionPyramidImageFilter< MovingImageType, MovingImageType>    MovingImagePyramidType;
 279     typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
 280     typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
 281
 282
 283     //     //=======================================================
 284     //     // Rigid Transform
 285     //     //=======================================================
 286     //     typedef itk::Euler3DTransform <double> RigidTransformType;
 287     //     RigidTransformType::Pointer rigidTransform=RigidTransformType::New();
 288
 289     //     if (m_ArgsInfo.rigid_given)
 290     //       {
 291     //          itk::Matrix<double,4,4> rigidTransformMatrix=clitk::ReadMatrix3D(m_ArgsInfo.rigid_arg);
 292
 293     //          //Set the rotation
 294     //          itk::Matrix<double,3,3> finalRotation = clitk::GetRotationalPartMatrix3D(rigidTransformMatrix);
 295     //          rigidTransform->SetMatrix(finalRotation);
 296
 297     //          //Set the translation
 298     //          itk::Vector<double,3> finalTranslation = clitk::GetTranslationPartMatrix3D(rigidTransformMatrix);
 299     //          rigidTransform->SetTranslation(finalTranslation);
 300
 301     //       }
 302
 303
 304     //=======================================================
 305     // BSpline Transform
 306     //=======================================================
 307     typename FixedImageType::RegionType::SizeType splineOrders ;
 308
 309     //Default is cubic splines
 310     splineOrders.Fill(3);
 311     if (m_ArgsInfo.order_given)
 312       for(unsigned int i=0; i<ImageDimension;i++)
 313         splineOrders[i]=m_ArgsInfo.order_arg[i];
 314
 315     // BLUT or ITK FFD
 316     typedef itk::Transform<TCoordRep, ImageDimension, SpaceDimension> TransformType;
 317     typename TransformType::Pointer transform;
 318     typedef  itk::BSplineDeformableTransform<TCoordRep,SpaceDimension, 3> BSplineTransformType;
 319     typedef  BSplineTransformType* BSplineTransformPointer;
 320     typedef  clitk::BSplineDeformableTransform<TCoordRep,ImageDimension, SpaceDimension > BLUTTransformType;
 321     typedef  BLUTTransformType* BLUTTransformPointer;
 322
 323     // JV parameter array is passed by reference, create outside context so it exists afterwards!!!!!
 324     typedef typename TransformType::ParametersType     ParametersType;
 325     ParametersType parameters;
 326
 327
 328     // CLITK BLUT transform
 329     if(m_ArgsInfo.wlut_flag)
 330       {
 331         typename BLUTTransformType::Pointer  bsplineTransform = BLUTTransformType::New();
 332         if (m_Verbose) std::cout<<"Setting the spline orders  to "<<splineOrders<<"..."<<std::endl;
 333         bsplineTransform->SetSplineOrders(splineOrders);
 334
 335         //-------------------------------------------------------------------------
 336         // Define the region: Either the spacing or the number of CP should be given
 337         //-------------------------------------------------------------------------
 338
 339         // Region
 340         typedef typename BSplineTransformType::RegionType RegionType;
 341         RegionType bsplineRegion;
 342         typename RegionType::SizeType   gridSizeOnImage;
 343         typename RegionType::SizeType   gridBorderSize;
 344         typename RegionType::SizeType   totalGridSize;
 345
 346         // Spacing
 347         typedef typename BSplineTransformType::SpacingType SpacingType;
 348         SpacingType fixedImageSpacing, chosenSpacing, adaptedSpacing;
 349         fixedImageSpacing = fixedImage->GetSpacing();
 350
 351         // Only spacing given: adjust if necessary
 352         if (m_ArgsInfo.spacing_given && !m_ArgsInfo.control_given)
 353           {
 354             for(unsigned int r=0; r<ImageDimension; r++)
 355               {
 356                 chosenSpacing[r]= m_ArgsInfo.spacing_arg[r];
 357                 gridSizeOnImage[r] = ceil( (double) transformRegion.GetSize()[r] / ( itk::Math::Round<double>(chosenSpacing[r]/fixedImageSpacing[r]) ) );
 358                 adaptedSpacing[r]= ( itk::Math::Round<double>(chosenSpacing[r]/fixedImageSpacing[r]) *fixedImageSpacing[r] ) ;
 359               }
 360             if (m_Verbose) std::cout<<"The chosen control point spacing "<<chosenSpacing<<"..."<<std::endl;
 361             if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<adaptedSpacing<<"..."<<std::endl;
 362             if (m_Verbose) std::cout<<"The number of (internal) control points is "<<gridSizeOnImage<<"..."<<std::endl;
 363           }
 364
 365         // Only number of CP given: adjust if necessary
 366         else if (m_ArgsInfo.control_given && !m_ArgsInfo.spacing_given)
 367           {
 368             for(unsigned int r=0; r<ImageDimension; r++)
 369               {
 370                 gridSizeOnImage[r]= m_ArgsInfo.control_arg[r];
 371                 chosenSpacing[r]=fixedImageSpacing[r]*( (double)(transformRegion.GetSize()[r])  /
 372                                                         (double)(gridSizeOnImage[r]) );
 373                 adaptedSpacing[r]= fixedImageSpacing[r]* ceil( (double)(transformRegion.GetSize()[r] - 1)  /
 374                                                                (double)(gridSizeOnImage[r] - 1) );
 375               }
 376             if (m_Verbose) std::cout<<"The chosen control point spacing "<<chosenSpacing<<"..."<<std::endl;
 377             if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<adaptedSpacing<<"..."<<std::endl;
 378             if (m_Verbose) std::cout<<"The number of (internal) control points is "<<gridSizeOnImage<<"..."<<std::endl;
 379           }
 380
 381         // Spacing and number of CP given: no adjustment adjust, just warnings
 382         else if (m_ArgsInfo.control_given && m_ArgsInfo.spacing_given)
 383           {
 384             for(unsigned int r=0; r<ImageDimension; r++)
 385               {
 386                 adaptedSpacing[r]= m_ArgsInfo.spacing_arg[r];
 387                 gridSizeOnImage[r] =  m_ArgsInfo.control_arg[r];
 388                 if (gridSizeOnImage[r]*adaptedSpacing[r]< transformRegion.GetSize()[r]*fixedImageSpacing[r])
 389                   {
 390                     std::cout<<"WARNING: Specified control point region ("<<gridSizeOnImage[r]*adaptedSpacing[r]
 391                              <<"mm) does not cover the transform region ("<< transformRegion.GetSize()[r]*fixedImageSpacing[r]
 392                              <<"mm) for dimension "<<r<<"!" <<std::endl
 393                              <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
 394                   }
 395                 if (  fmod(adaptedSpacing[r], fixedImageSpacing[r]) )
 396                   {
 397                     std::cout<<"WARNING: Specified control point spacing for dimension "<<r
 398                              <<" does not allow exact representation of BLUT FFD!"<<std::endl
 399                              <<"Spacing ratio is non-integer: "<<adaptedSpacing[r]/ fixedImageSpacing[r]<<std::endl
 400                              <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
 401                   }
 402               }
 403             if (m_Verbose) std::cout<<"The control points spacing was set to "<<adaptedSpacing<<"..."<<std::endl;
 404             if (m_Verbose) std::cout<<"The number of (internal) control points spacing is "<<gridSizeOnImage<<"..."<<std::endl;
 405           }
 406
 407         //JV  border size should depend on spline order
 408         for(unsigned int r=0; r<ImageDimension; r++) gridBorderSize[r]=splineOrders[r]; // Border for spline order = 3 ( 1 lower, 2 upper )
 409         totalGridSize = gridSizeOnImage + gridBorderSize;
 410         bsplineRegion.SetSize( totalGridSize );
 411         if (m_Verbose) std::cout<<"The total control point grid size was set to "<<totalGridSize<<"..."<<std::endl;
 412
 413         // Direction
 414         typename FixedImageType::DirectionType gridDirection = fixedImage->GetDirection();
 415         SpacingType gridOriginOffset = gridDirection * adaptedSpacing;
 416
 417         // Origin: 1 CP border for spatial dimensions
 418         typedef typename BSplineTransformType::OriginType OriginType;
 419         OriginType gridOrigin = transformRegionOrigin - gridOriginOffset;
 420         if (m_Verbose) std::cout<<"The control point grid origin was set to "<<gridOrigin<<"..."<<std::endl;
 421
 422         // Set
 423         bsplineTransform->SetGridSpacing( adaptedSpacing );
 424         bsplineTransform->SetGridOrigin( gridOrigin );
 425         bsplineTransform->SetGridRegion( bsplineRegion );
 426         bsplineTransform->SetGridDirection( gridDirection );
 427
 428         //Bulk transform
 429         //if (m_Verbose) std::cout<<"Setting rigid transform..."<<std::endl;
 430         //bsplineTransform->SetBulkTransform( rigidTransform );
 431
 432         //Vector BSpline interpolator
 433         //bsplineTransform->SetOutputSpacing(fixedImage->GetSpacing());
 434         typename RegionType::SizeType samplingFactors;
 435         for (unsigned int i=0; i< ImageDimension; i++)
 436           {
 437             if (m_Verbose) std::cout<<"For dimension "<<i<<", the ideal sampling factor (if integer) is a multitude of "
 438                                     << (double)adaptedSpacing[i]/ (double) fixedImageSpacing[i]<<"..."<<std::endl;
 439             if (m_ArgsInfo.samplingFactor_given) samplingFactors[i]=m_ArgsInfo.samplingFactor_arg[i];
 440             else samplingFactors[i]=(int) ((double)adaptedSpacing[i]/ (double) movingImage->GetSpacing()[i]);
 441             if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<samplingFactors[i]<<"..."<<std::endl;
 442           }
 443         bsplineTransform->SetLUTSamplingFactors(samplingFactors);
 444
 445         //initial parameters
 446         if (m_ArgsInfo.init_given)
 447           {
 448             typedef itk::ImageFileReader<typename BLUTTransformType::CoefficientImageType> CoefficientReaderType;
 449             typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New();
 450             coeffReader->SetFileName(m_ArgsInfo.init_arg[0]);
 451             coeffReader->Update();
 452             bsplineTransform->SetCoefficientImage(coeffReader->GetOutput());
 453           }
 454         else
 455           {
 456             //typedef typename TransformType::ParametersType     ParametersType;
 457             const unsigned int numberOfParameters =    bsplineTransform->GetNumberOfParameters();
 458             parameters=ParametersType( numberOfParameters );
 459             parameters.Fill( 0.0 );
 460             bsplineTransform->SetParameters( parameters );
 461           }
 462
 463         // Mask
 464         if (spatialObjectMask) bsplineTransform->SetMask( spatialObjectMask );
 465
 466         // Pass
 467         transform=bsplineTransform;
 468       }
 469
 470     //ITK BSpline transform
 471     else
 472       {
 473         typename BSplineTransformType::Pointer  bsplineTransform = BSplineTransformType::New();
 474
 475         // Define the region
 476         typedef typename BSplineTransformType::RegionType RegionType;
 477         RegionType bsplineRegion;
 478         typename RegionType::SizeType   gridSizeOnImage;
 479         typename RegionType::SizeType   gridBorderSize;
 480         typename RegionType::SizeType   totalGridSize;
 481
 482 #if 0
 483         //Set the number of control points
 484         for(unsigned int r=0; r<ImageDimension; r++)  gridSizeOnImage[r]=m_ArgsInfo.control_arg[r];
 485         if (m_Verbose) std::cout<<"Setting the number of internal control points "<<gridSizeOnImage<<"..."<<std::endl;
 486         gridBorderSize.Fill( 3 );    // Border for spline order = 3 ( 1 lower, 2 upper )
 487         totalGridSize = gridSizeOnImage + gridBorderSize;
 488         bsplineRegion.SetSize( totalGridSize );
 489
 490         // Spacing
 491         typedef typename BSplineTransformType::SpacingType SpacingType;
 492         SpacingType spacing = fixedImage->GetSpacing();
 493         typename FixedImageType::SizeType fixedImageSize = fixedImageRegion.GetSize();
 494         if (m_ArgsInfo.spacing_given)
 495           {
 496
 497             for(unsigned int r=0; r<ImageDimension; r++)
 498               {
 499                 spacing[r] =m_ArgsInfo.spacing_arg[r];
 500               }
 501           }
 502         else
 503           {
 504             for(unsigned int r=0; r<ImageDimension; r++)
 505               {
 506                 spacing[r] *= static_cast<double>(fixedImageSize[r] - 1)  /
 507                   static_cast<double>(gridSizeOnImage[r] - 1);
 508               }
 509           }
 510         if (m_Verbose) std::cout<<"The control points spacing was set to "<<spacing<<"..."<<std::endl;
 511
 512         // Direction
 513         typename FixedImageType::DirectionType gridDirection = fixedImage->GetDirection();
 514         SpacingType gridOriginOffset = gridDirection * spacing;
 515
 516         // Origin
 517         typedef typename BSplineTransformType::OriginType OriginType;
 518         OriginType origin = fixedImage->GetOrigin();
 519         OriginType gridOrigin = origin - gridOriginOffset;
 520
 521         // Set
 522         bsplineTransform->SetGridSpacing( spacing );
 523         bsplineTransform->SetGridOrigin( gridOrigin );
 524         bsplineTransform->SetGridRegion( bsplineRegion );
 525         bsplineTransform->SetGridDirection( gridDirection );
 526
 527         // Bulk transform
 528         // if (m_Verbose) std::cout<<"Setting rigid transform..."<<std::endl;
 529         // bsplineTransform->SetBulkTransform( rigidTransform );
 530
 531         // Initial parameters
 532         if (m_ArgsInfo.init_given)
 533           {
 534             typedef itk::ImageFileReader<typename BSplineTransformType::ImageType> CoefficientReaderType;
 535             typename BSplineTransformType::ImageType::Pointer coeffImages[SpaceDimension];
 536             for(unsigned int i=0; i<SpaceDimension; i++)
 537               {
 538                 typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New();
 539                 coeffReader->SetFileName(m_ArgsInfo.init_arg[i]);
 540                 coeffReader->Update();
 541                 coeffImages[i]=coeffReader->GetOutput();
 542               }
 543             bsplineTransform->SetCoefficientImage(coeffImages);
 544           }
 545         else
 546           {
 547             const unsigned int numberOfParameters =    bsplineTransform->GetNumberOfParameters();
 548             parameters=ParametersType( numberOfParameters );
 549             parameters.Fill( 0.0 );
 550             bsplineTransform->SetParameters( parameters );
 551           }
 552
 553         // Pass
 554         transform=bsplineTransform;
 555 #else
 556         // Spacing
 557         typedef typename BSplineTransformType::SpacingType SpacingType;
 558         SpacingType fixedImageSpacing, chosenSpacing, adaptedSpacing;
 559         fixedImageSpacing = fixedImage->GetSpacing();
 560
 561         // Only spacing given: adjust if necessary
 562         if (m_ArgsInfo.spacing_given && !m_ArgsInfo.control_given)
 563           {
 564             for(unsigned int r=0; r<ImageDimension; r++)
 565               {
 566                 chosenSpacing[r]= m_ArgsInfo.spacing_arg[r];
 567                 gridSizeOnImage[r] = ceil( (double) transformRegion.GetSize()[r] / ( round(chosenSpacing[r]/fixedImageSpacing[r]) ) );
 568                 adaptedSpacing[r]= ( round(chosenSpacing[r]/fixedImageSpacing[r]) *fixedImageSpacing[r] ) ;
 569               }
 570             if (m_Verbose) std::cout<<"The chosen control point spacing "<<chosenSpacing<<"..."<<std::endl;
 571             if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<adaptedSpacing<<"..."<<std::endl;
 572             if (m_Verbose) std::cout<<"The number of (internal) control points is "<<gridSizeOnImage<<"..."<<std::endl;
 573           }
 574
 575         // Only number of CP given: adjust if necessary
 576         else if (m_ArgsInfo.control_given && !m_ArgsInfo.spacing_given)
 577           {
 578             for(unsigned int r=0; r<ImageDimension; r++)
 579               {
 580                 gridSizeOnImage[r]= m_ArgsInfo.control_arg[r];
 581                 chosenSpacing[r]=fixedImageSpacing[r]*( (double)(transformRegion.GetSize()[r])  /
 582                                                         (double)(gridSizeOnImage[r]) );
 583                 adaptedSpacing[r]= fixedImageSpacing[r]* ceil( (double)(transformRegion.GetSize()[r] - 1)  /
 584                                                                (double)(gridSizeOnImage[r] - 1) );
 585               }
 586             if (m_Verbose) std::cout<<"The chosen control point spacing "<<chosenSpacing<<"..."<<std::endl;
 587             if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<adaptedSpacing<<"..."<<std::endl;
 588             if (m_Verbose) std::cout<<"The number of (internal) control points is "<<gridSizeOnImage<<"..."<<std::endl;
 589           }
 590
 591         // Spacing and number of CP given: no adjustment adjust, just warnings
 592         else if (m_ArgsInfo.control_given && m_ArgsInfo.spacing_given)
 593           {
 594             for(unsigned int r=0; r<ImageDimension; r++)
 595               {
 596                 adaptedSpacing[r]= m_ArgsInfo.spacing_arg[r];
 597                 gridSizeOnImage[r] =  m_ArgsInfo.control_arg[r];
 598                 if (gridSizeOnImage[r]*adaptedSpacing[r]< transformRegion.GetSize()[r]*fixedImageSpacing[r])
 599                   {
 600                     std::cout<<"WARNING: Specified control point region ("<<gridSizeOnImage[r]*adaptedSpacing[r]
 601                              <<"mm) does not cover the transform region ("<< transformRegion.GetSize()[r]*fixedImageSpacing[r]
 602                              <<"mm) for dimension "<<r<<"!" <<std::endl
 603                              <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
 604                   }
 605                 if (  fmod(adaptedSpacing[r], fixedImageSpacing[r]) )
 606                   {
 607                     std::cout<<"WARNING: Specified control point spacing for dimension "<<r
 608                              <<" does not allow exact representation of BLUT FFD!"<<std::endl
 609                              <<"Spacing ratio is non-integer: "<<adaptedSpacing[r]/ fixedImageSpacing[r]<<std::endl
 610                              <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
 611                   }
 612               }
 613             if (m_Verbose) std::cout<<"The control points spacing was set to "<<adaptedSpacing<<"..."<<std::endl;
 614             if (m_Verbose) std::cout<<"The number of (internal) control points spacing is "<<gridSizeOnImage<<"..."<<std::endl;
 615           }
 616
 617         //JV  border size should depend on spline order
 618         for(unsigned int r=0; r<ImageDimension; r++) gridBorderSize[r]=splineOrders[r]; // Border for spline order = 3 ( 1 lower, 2 upper )
 619         totalGridSize = gridSizeOnImage + gridBorderSize;
 620         bsplineRegion.SetSize( totalGridSize );
 621         if (m_Verbose) std::cout<<"The total control point grid size was set to "<<totalGridSize<<"..."<<std::endl;
 622
 623         // Direction
 624         typename FixedImageType::DirectionType gridDirection = fixedImage->GetDirection();
 625         SpacingType gridOriginOffset = gridDirection * adaptedSpacing;
 626
 627         // Origin: 1 CP border for spatial dimensions
 628         typedef typename BSplineTransformType::OriginType OriginType;
 629         OriginType gridOrigin = transformRegionOrigin - gridOriginOffset;
 630         if (m_Verbose) std::cout<<"The control point grid origin was set to "<<gridOrigin<<"..."<<std::endl;
 631
 632         // Set
 633         bsplineTransform->SetGridSpacing( adaptedSpacing );
 634         bsplineTransform->SetGridOrigin( gridOrigin );
 635         bsplineTransform->SetGridRegion( bsplineRegion );
 636         bsplineTransform->SetGridDirection( gridDirection );
 637
 638         //Bulk transform
 639         //if (m_Verbose) std::cout<<"Setting rigid transform..."<<std::endl;
 640         //bsplineTransform->SetBulkTransform( rigidTransform );
 641
 642         //Vector BSpline interpolator
 643         //bsplineTransform->SetOutputSpacing(fixedImage->GetSpacing());
 644
 645         // Initial parameters
 646         if (m_ArgsInfo.init_given)
 647           {
 648             typedef itk::ImageFileReader<typename BSplineTransformType::ImageType> CoefficientReaderType;
 649             typename BSplineTransformType::ImageType::Pointer coeffImages[SpaceDimension];
 650             for(unsigned int i=0; i<SpaceDimension; i++)
 651               {
 652                 typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New();
 653                 coeffReader->SetFileName(m_ArgsInfo.init_arg[i]);
 654                 coeffReader->Update();
 655                 coeffImages[i]=coeffReader->GetOutput();
 656               }
 657             bsplineTransform->SetCoefficientImage(coeffImages);
 658           }
 659         else
 660           {
 661             const unsigned int numberOfParameters =    bsplineTransform->GetNumberOfParameters();
 662             parameters=ParametersType( numberOfParameters );
 663             parameters.Fill( 0.0 );
 664             bsplineTransform->SetParameters( parameters );
 665           }
 666
 667         // Pass
 668         transform=bsplineTransform;
 669 #endif
 670
 671       }
 672
 673
 674     //=======================================================
 675     // Interpolator
 676     //=======================================================
 677     typedef clitk::GenericInterpolator<args_info_clitkBSplineDeformableRegistration, FixedImageType, TCoordRep > GenericInterpolatorType;
 678     typename   GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
 679     genericInterpolator->SetArgsInfo(m_ArgsInfo);
 680     typedef itk::InterpolateImageFunction< FixedImageType, TCoordRep >  InterpolatorType;
 681     typename  InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
 682
 683
 684     //=======================================================
 685     // Metric
 686     //=======================================================
 687     typedef clitk::GenericMetric<args_info_clitkBSplineDeformableRegistration, FixedImageType,MovingImageType > GenericMetricType;
 688     typename GenericMetricType::Pointer genericMetric=GenericMetricType::New();
 689     genericMetric->SetArgsInfo(m_ArgsInfo);
 690     genericMetric->SetFixedImage(fixedImage);
 691     if (spatialObjectMask) genericMetric->SetFixedImageMask( spatialObjectMask );
 692     genericMetric->SetFixedImageRegion(metricRegion);
 693     typedef itk::ImageToImageMetric< FixedImageType, MovingImageType >  MetricType;
 694     typename  MetricType::Pointer metric=genericMetric->GetMetricPointer();
 695
 696 #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
 697     if (threadsGiven) metric->SetNumberOfThreads( threads );
 698 #else
 699     if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
 700 #endif
 701
 702
 703     //=======================================================
 704     // Optimizer
 705     //=======================================================
 706     typedef clitk::GenericOptimizer<args_info_clitkBSplineDeformableRegistration> GenericOptimizerType;
 707     GenericOptimizerType::Pointer genericOptimizer = GenericOptimizerType::New();
 708     genericOptimizer->SetArgsInfo(m_ArgsInfo);
 709     genericOptimizer->SetMaximize(genericMetric->GetMaximize());
 710     genericOptimizer->SetNumberOfParameters(transform->GetNumberOfParameters());
 711     typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
 712     OptimizerType::Pointer optimizer = genericOptimizer->GetOptimizerPointer();
 713
 714
 715     //=======================================================
 716     // Registration
 717     //=======================================================
 718     typedef itk::MultiResolutionImageRegistrationMethod<  FixedImageType, MovingImageType >    RegistrationType;
 719     typename RegistrationType::Pointer   registration  = RegistrationType::New();
 720     registration->SetMetric(        metric        );
 721     registration->SetOptimizer(     optimizer     );
 722     registration->SetInterpolator(  interpolator  );
 723     registration->SetTransform (transform);
 724     if(threadsGiven) registration->SetNumberOfThreads(threads);
 725     registration->SetFixedImage(  fixedImage   );
 726     registration->SetMovingImage(   movingImage   );
 727     registration->SetFixedImageRegion( metricRegion );
 728     registration->SetFixedImagePyramid( fixedImagePyramid );
 729     registration->SetMovingImagePyramid( movingImagePyramid );
 730     registration->SetInitialTransformParameters( transform->GetParameters() );
 731     registration->SetNumberOfLevels(m_ArgsInfo.levels_arg);
 732     if (m_Verbose) std::cout<<"Setting the number of resolution levels to "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
 733
 734
 735     //================================================================================================
 736     // Observers
 737     //================================================================================================
 738     if (m_Verbose)
 739       {
 740         // Output iteration info
 741         CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
 742         observer->SetOptimizer(genericOptimizer);
 743         optimizer->AddObserver( itk::IterationEvent(), observer );
 744
 745
 746         // Output level info
 747         typedef RegistrationInterfaceCommand<RegistrationType> CommandType;
 748         typename CommandType::Pointer command = CommandType::New();
 749         registration->AddObserver( itk::IterationEvent(), command );
 750       }
 751
 752
 753     //=======================================================
 754     // Let's go
 755     //=======================================================
 756     if (m_Verbose) std::cout << std::endl << "Starting Registration" << std::endl;
 757
 758     try
 759       {
 760         registration->StartRegistration();
 761       }
 762     catch( itk::ExceptionObject & err )
 763       {
 764         std::cerr << "ExceptionObject caught while registering!" << std::endl;
 765         std::cerr << err << std::endl;
 766         return;
 767       }
 768
 769
 770     //=======================================================
 771     // Get the result
 772     //=======================================================
 773     OptimizerType::ParametersType finalParameters =  registration->GetLastTransformParameters();
 774     transform->SetParameters( finalParameters );
 775
 776
 777     //=======================================================
 778     // Get the BSpline coefficient images and write them
 779     //=======================================================
 780     if (m_ArgsInfo.coeff_given)
 781       {
 782         if(m_ArgsInfo.wlut_flag)
 783           {
 784             BLUTTransformPointer bsplineTransform=dynamic_cast<BLUTTransformPointer>(registration->GetTransform());
 785             typedef  itk::Image<itk::Vector<TCoordRep, SpaceDimension>, ImageDimension> CoefficientImageType;
 786             typename CoefficientImageType::Pointer coefficientImage =bsplineTransform->GetCoefficientImage();
 787             typedef itk::ImageFileWriter<CoefficientImageType> CoeffWriterType;
 788             typename CoeffWriterType::Pointer coeffWriter=CoeffWriterType::New();
 789             coeffWriter->SetInput(coefficientImage);
 790             coeffWriter->SetFileName(m_ArgsInfo.coeff_arg[0]);
 791             coeffWriter->Update();
 792           }
 793         else
 794           {
 795             BSplineTransformPointer bsplineTransform=dynamic_cast<BSplineTransformPointer>(registration->GetTransform());
 796             typedef  itk::Image<TCoordRep, ImageDimension> CoefficientImageType;
 797 #if ITK_VERSION_MAJOR > 3
 798             typename BSplineTransformType::CoefficientImageArray coefficientImages = bsplineTransform->GetCoefficientImage();
 799 #else
 800             typename CoefficientImageType::Pointer *coefficientImages =bsplineTransform->GetCoefficientImage();
 801 #endif
 802             typedef itk::ImageFileWriter<CoefficientImageType> CoeffWriterType;
 803             for (unsigned int i=0;i<SpaceDimension; i ++)
 804               {
 805                 typename CoeffWriterType::Pointer coeffWriter=CoeffWriterType::New();
 806                 coeffWriter->SetInput(coefficientImages[i]);
 807                 coeffWriter->SetFileName(m_ArgsInfo.coeff_arg[i]);
 808                 coeffWriter->Update();
 809               }
 810           }
 811       }
 812
 813
 814     //=======================================================
 815     // Generate the DVF
 816     //=======================================================
 817     typedef itk::Vector< float, SpaceDimension >  DisplacementType;
 818     typedef itk::Image< DisplacementType, ImageDimension >  DeformationFieldType;
 819
 820     typename DeformationFieldType::Pointer field = DeformationFieldType::New();
 821     field->SetRegions( fixedImageRegion );
 822     field->SetOrigin( fixedImage->GetOrigin() );
 823     field->SetSpacing( fixedImage->GetSpacing() );
 824     field->SetDirection( fixedImage->GetDirection() );
 825     field->Allocate();
 826
 827     typedef itk::ImageRegionIteratorWithIndex< DeformationFieldType > FieldIterator;
 828     FieldIterator fi( field, fixedImageRegion );
 829     fi.GoToBegin();
 830
 831     typename TransformType::InputPointType  fixedPoint;
 832     typename TransformType::OutputPointType movingPoint;
 833     typename DeformationFieldType::IndexType index;
 834
 835     DisplacementType displacement;
 836     while( ! fi.IsAtEnd() )
 837       {
 838         index = fi.GetIndex();
 839         field->TransformIndexToPhysicalPoint( index, fixedPoint );
 840         movingPoint = transform->TransformPoint( fixedPoint );
 841         displacement = movingPoint - fixedPoint;
 842         fi.Set( displacement );
 843         ++fi;
 844       }
 845
 846
 847     //=======================================================
 848     // Write the DVF
 849     //=======================================================
 850     typedef itk::ImageFileWriter< DeformationFieldType >  FieldWriterType;
 851     typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
 852     fieldWriter->SetFileName( m_ArgsInfo.vf_arg );
 853     fieldWriter->SetInput( field );
 854     try
 855       {
 856         fieldWriter->Update();
 857       }
 858     catch( itk::ExceptionObject & excp )
 859       {
 860         std::cerr << "Exception thrown writing the DVF" << std::endl;
 861         std::cerr << excp << std::endl;
 862         return;
 863       }
 864
 865
 866     //=======================================================
 867     // Resample the moving image
 868     //=======================================================
 869     typedef itk::WarpImageFilter< MovingImageType, FixedImageType, DeformationFieldType >    WarpFilterType;
 870     typename WarpFilterType::Pointer warp = WarpFilterType::New();
 871
 872     warp->SetDeformationField( field );
 873     warp->SetInput( movingImageReader->GetOutput() );
 874     warp->SetOutputOrigin(  fixedImage->GetOrigin() );
 875     warp->SetOutputSpacing( fixedImage->GetSpacing() );
 876     warp->SetOutputDirection( fixedImage->GetDirection() );
 877     warp->SetEdgePaddingValue( 0.0 );
 878     warp->Update();
 879
 880
 881     //=======================================================
 882     // Write the warped image
 883     //=======================================================
 884     typedef itk::ImageFileWriter< FixedImageType >  WriterType;
 885     typename WriterType::Pointer      writer =  WriterType::New();
 886     writer->SetFileName( m_ArgsInfo.output_arg );
 887     writer->SetInput( warp->GetOutput()    );
 888
 889     try
 890       {
 891         writer->Update();
 892       }
 893     catch( itk::ExceptionObject & err )
 894       {
 895         std::cerr << "ExceptionObject caught !" << std::endl;
 896         std::cerr << err << std::endl;
 897         return;
 898       }
 899
 900
 901     //=======================================================
 902     // Calculate the difference after the deformable transform
 903     //=======================================================
 904     typedef clitk::DifferenceImageFilter<  FixedImageType, FixedImageType> DifferenceFilterType;
 905     if (m_ArgsInfo.after_given)
 906       {
 907         typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
 908         difference->SetValidInput( fixedImage );
 909         difference->SetTestInput( warp->GetOutput() );
 910
 911         try
 912           {
 913             difference->Update();
 914           }
 915         catch( itk::ExceptionObject & err )
 916           {
 917             std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
 918             std::cerr << err << std::endl;
 919             return;
 920           }
 921
 922         typename WriterType::Pointer differenceWriter=WriterType::New();
 923         differenceWriter->SetInput(difference->GetOutput());
 924         differenceWriter->SetFileName(m_ArgsInfo.after_arg);
 925         differenceWriter->Update();
 926
 927       }
 928
 929
 930     //=======================================================
 931     // Calculate the difference before the deformable transform
 932     //=======================================================
 933     if( m_ArgsInfo.before_given )
 934       {
 935
 936         typename FixedImageType::Pointer moving=FixedImageType::New();
 937         if (m_ArgsInfo.rigid_given)
 938           {
 939             typedef itk::ResampleImageFilter<MovingImageType, FixedImageType> ResamplerType;
 940             typename ResamplerType::Pointer resampler=ResamplerType::New();
 941             resampler->SetInput(movingImage);
 942             resampler->SetOutputOrigin(fixedImage->GetOrigin());
 943             resampler->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
 944             resampler->SetOutputSpacing(fixedImage->GetSpacing());
 945             resampler->SetDefaultPixelValue( 0. );
 946             //resampler->SetTransform(rigidTransform);
 947             resampler->Update();
 948             moving=resampler->GetOutput();
 949           }
 950         else
 951           moving=movingImage;
 952
 953         typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
 954         difference->SetValidInput( fixedImage );
 955         difference->SetTestInput( moving );
 956
 957         try
 958           {
 959             difference->Update();
 960           }
 961         catch( itk::ExceptionObject & err )
 962           {
 963             std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
 964             std::cerr << err << std::endl;
 965             return;
 966           }
 967
 968         typename WriterType::Pointer differenceWriter=WriterType::New();
 969         writer->SetFileName( m_ArgsInfo.before_arg  );
 970         writer->SetInput( difference->GetOutput()  );
 971         writer->Update( );
 972       }
 973
 974     return;
 975   }
 976 }
 977
 978 #endif // __clitkBSplineDeformableRegistrationGenericFilter_txx