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://oncora1.lyon.fnclcc.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ======================================================================-====*/
  18 #ifndef __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] / ( round(chosenSpacing[r]/fixedImageSpacing[r]) ) );
 358                 adaptedSpacing[r]= ( round(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         //Set the number of control points
 483         for(unsigned int r=0; r<ImageDimension; r++)  gridSizeOnImage[r]=m_ArgsInfo.control_arg[r];
 484         if (m_Verbose) std::cout<<"Setting the number of internal control points "<<gridSizeOnImage<<"..."<<std::endl;
 485         gridBorderSize.Fill( 3 );    // Border for spline order = 3 ( 1 lower, 2 upper )
 486         totalGridSize = gridSizeOnImage + gridBorderSize;
 487         bsplineRegion.SetSize( totalGridSize );
 488
 489         // Spacing
 490         typedef typename BSplineTransformType::SpacingType SpacingType;
 491         SpacingType spacing = fixedImage->GetSpacing();
 492         typename FixedImageType::SizeType fixedImageSize = fixedImageRegion.GetSize();
 493         if (m_ArgsInfo.spacing_given)
 494           {
 495
 496             for(unsigned int r=0; r<ImageDimension; r++)
 497               {
 498                 spacing[r] =m_ArgsInfo.spacing_arg[r];
 499               }
 500           }
 501         else
 502           {
 503             for(unsigned int r=0; r<ImageDimension; r++)
 504               {
 505                 spacing[r] *= static_cast<double>(fixedImageSize[r] - 1)  /
 506                   static_cast<double>(gridSizeOnImage[r] - 1);
 507               }
 508           }
 509         if (m_Verbose) std::cout<<"The control points spacing was set to "<<spacing<<"..."<<std::endl;
 510
 511         // Direction
 512         typename FixedImageType::DirectionType gridDirection = fixedImage->GetDirection();
 513         SpacingType gridOriginOffset = gridDirection * spacing;
 514
 515         // Origin
 516         typedef typename BSplineTransformType::OriginType OriginType;
 517         OriginType origin = fixedImage->GetOrigin();
 518         OriginType gridOrigin = origin - gridOriginOffset;
 519
 520         // Set
 521         bsplineTransform->SetGridSpacing( spacing );
 522         bsplineTransform->SetGridOrigin( gridOrigin );
 523         bsplineTransform->SetGridRegion( bsplineRegion );
 524         bsplineTransform->SetGridDirection( gridDirection );
 525
 526         // Bulk transform
 527         // if (m_Verbose) std::cout<<"Setting rigid transform..."<<std::endl;
 528         // bsplineTransform->SetBulkTransform( rigidTransform );
 529
 530         // Initial parameters
 531         if (m_ArgsInfo.init_given)
 532           {
 533             typedef itk::ImageFileReader<typename BSplineTransformType::ImageType> CoefficientReaderType;
 534             typename BSplineTransformType::ImageType::Pointer coeffImages[SpaceDimension];
 535             for(unsigned int i=0; i<SpaceDimension; i++)
 536               {
 537                 typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New();
 538                 coeffReader->SetFileName(m_ArgsInfo.init_arg[i]);
 539                 coeffReader->Update();
 540                 coeffImages[i]=coeffReader->GetOutput();
 541               }
 542             bsplineTransform->SetCoefficientImage(coeffImages);
 543           }
 544         else
 545           {
 546             const unsigned int numberOfParameters =    bsplineTransform->GetNumberOfParameters();
 547             parameters=ParametersType( numberOfParameters );
 548             parameters.Fill( 0.0 );
 549             bsplineTransform->SetParameters( parameters );
 550           }
 551
 552         // Pass
 553         transform=bsplineTransform;
 554
 555       }
 556
 557
 558     //=======================================================
 559     // Interpolator
 560     //=======================================================
 561     typedef clitk::GenericInterpolator<args_info_clitkBSplineDeformableRegistration, FixedImageType, TCoordRep > GenericInterpolatorType;
 562     typename   GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
 563     genericInterpolator->SetArgsInfo(m_ArgsInfo);
 564     typedef itk::InterpolateImageFunction< FixedImageType, TCoordRep >  InterpolatorType;
 565     typename  InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
 566
 567
 568     //=======================================================
 569     // Metric
 570     //=======================================================
 571     typedef clitk::GenericMetric<args_info_clitkBSplineDeformableRegistration, FixedImageType,MovingImageType > GenericMetricType;
 572     typename GenericMetricType::Pointer genericMetric=GenericMetricType::New();
 573     genericMetric->SetArgsInfo(m_ArgsInfo);
 574     genericMetric->SetFixedImageRegion(metricRegion);
 575     typedef itk::ImageToImageMetric< FixedImageType, MovingImageType >  MetricType;
 576     typename  MetricType::Pointer metric=genericMetric->GetMetricPointer();
 577     if (spatialObjectMask) metric->SetFixedImageMask( spatialObjectMask );
 578
 579 #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
 580     if (threadsGiven) metric->SetNumberOfThreads( threads );
 581 #else
 582     if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
 583 #endif
 584
 585
 586     //=======================================================
 587     // Optimizer
 588     //=======================================================
 589     typedef clitk::GenericOptimizer<args_info_clitkBSplineDeformableRegistration> GenericOptimizerType;
 590     GenericOptimizerType::Pointer genericOptimizer = GenericOptimizerType::New();
 591     genericOptimizer->SetArgsInfo(m_ArgsInfo);
 592     genericOptimizer->SetMaximize(genericMetric->GetMaximize());
 593     genericOptimizer->SetNumberOfParameters(transform->GetNumberOfParameters());
 594     typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
 595     OptimizerType::Pointer optimizer = genericOptimizer->GetOptimizerPointer();
 596
 597
 598     //=======================================================
 599     // Registration
 600     //=======================================================
 601     typedef itk::MultiResolutionImageRegistrationMethod<  FixedImageType, MovingImageType >    RegistrationType;
 602     typename RegistrationType::Pointer   registration  = RegistrationType::New();
 603     registration->SetMetric(        metric        );
 604     registration->SetOptimizer(     optimizer     );
 605     registration->SetInterpolator(  interpolator  );
 606     registration->SetTransform (transform);
 607     if(threadsGiven) registration->SetNumberOfThreads(threads);
 608     registration->SetFixedImage(  fixedImage   );
 609     registration->SetMovingImage(   movingImage   );
 610     registration->SetFixedImageRegion( metricRegion );
 611     registration->SetFixedImagePyramid( fixedImagePyramid );
 612     registration->SetMovingImagePyramid( movingImagePyramid );
 613     registration->SetInitialTransformParameters( transform->GetParameters() );
 614     registration->SetNumberOfLevels(m_ArgsInfo.levels_arg);
 615     if (m_Verbose) std::cout<<"Setting the number of resolution levels to "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
 616
 617
 618     //================================================================================================
 619     // Observers
 620     //================================================================================================
 621     if (m_Verbose)
 622       {
 623         // Output iteration info
 624         CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
 625         observer->SetOptimizer(genericOptimizer);
 626         optimizer->AddObserver( itk::IterationEvent(), observer );
 627
 628
 629         // Output level info
 630         typedef RegistrationInterfaceCommand<RegistrationType> CommandType;
 631         typename CommandType::Pointer command = CommandType::New();
 632         registration->AddObserver( itk::IterationEvent(), command );
 633       }
 634
 635
 636     //=======================================================
 637     // Let's go
 638     //=======================================================
 639     if (m_Verbose) std::cout << std::endl << "Starting Registration" << std::endl;
 640
 641     try
 642       {
 643         registration->StartRegistration();
 644       }
 645     catch( itk::ExceptionObject & err )
 646       {
 647         std::cerr << "ExceptionObject caught while registering!" << std::endl;
 648         std::cerr << err << std::endl;
 649         return;
 650       }
 651
 652
 653     //=======================================================
 654     // Get the result
 655     //=======================================================
 656     OptimizerType::ParametersType finalParameters =  registration->GetLastTransformParameters();
 657     transform->SetParameters( finalParameters );
 658
 659
 660     //=======================================================
 661     // Get the BSpline coefficient images and write them
 662     //=======================================================
 663     if (m_ArgsInfo.coeff_given)
 664       {
 665         if(m_ArgsInfo.wlut_flag)
 666           {
 667             BLUTTransformPointer bsplineTransform=dynamic_cast<BLUTTransformPointer>(registration->GetTransform());
 668             typedef  itk::Image<itk::Vector<TCoordRep, SpaceDimension>, ImageDimension> CoefficientImageType;
 669             typename CoefficientImageType::Pointer coefficientImage =bsplineTransform->GetCoefficientImage();
 670             typedef itk::ImageFileWriter<CoefficientImageType> CoeffWriterType;
 671             typename CoeffWriterType::Pointer coeffWriter=CoeffWriterType::New();
 672             coeffWriter->SetInput(coefficientImage);
 673             coeffWriter->SetFileName(m_ArgsInfo.coeff_arg[0]);
 674             coeffWriter->Update();
 675           }
 676         else
 677           {
 678             BSplineTransformPointer bsplineTransform=dynamic_cast<BSplineTransformPointer>(registration->GetTransform());
 679             typedef  itk::Image<TCoordRep, ImageDimension> CoefficientImageType;
 680 #if ITK_VERSION_MAJOR > 3
 681             typename BSplineTransformType::CoefficientImageArray coefficientImages = bsplineTransform->GetCoefficientImage();
 682 #else
 683             typename CoefficientImageType::Pointer *coefficientImages =bsplineTransform->GetCoefficientImage();
 684 #endif
 685             typedef itk::ImageFileWriter<CoefficientImageType> CoeffWriterType;
 686             for (unsigned int i=0;i<SpaceDimension; i ++)
 687               {
 688                 typename CoeffWriterType::Pointer coeffWriter=CoeffWriterType::New();
 689                 coeffWriter->SetInput(coefficientImages[i]);
 690                 coeffWriter->SetFileName(m_ArgsInfo.coeff_arg[i]);
 691                 coeffWriter->Update();
 692               }
 693           }
 694       }
 695
 696
 697     //=======================================================
 698     // Generate the DVF
 699     //=======================================================
 700     typedef itk::Vector< float, SpaceDimension >  DisplacementType;
 701     typedef itk::Image< DisplacementType, ImageDimension >  DeformationFieldType;
 702
 703     typename DeformationFieldType::Pointer field = DeformationFieldType::New();
 704     field->SetRegions( fixedImageRegion );
 705     field->SetOrigin( fixedImage->GetOrigin() );
 706     field->SetSpacing( fixedImage->GetSpacing() );
 707     field->SetDirection( fixedImage->GetDirection() );
 708     field->Allocate();
 709
 710     typedef itk::ImageRegionIteratorWithIndex< DeformationFieldType > FieldIterator;
 711     FieldIterator fi( field, fixedImageRegion );
 712     fi.GoToBegin();
 713
 714     typename TransformType::InputPointType  fixedPoint;
 715     typename TransformType::OutputPointType movingPoint;
 716     typename DeformationFieldType::IndexType index;
 717
 718     DisplacementType displacement;
 719     while( ! fi.IsAtEnd() )
 720       {
 721         index = fi.GetIndex();
 722         field->TransformIndexToPhysicalPoint( index, fixedPoint );
 723         movingPoint = transform->TransformPoint( fixedPoint );
 724         displacement = movingPoint - fixedPoint;
 725         fi.Set( displacement );
 726         ++fi;
 727       }
 728
 729
 730     //=======================================================
 731     // Write the DVF
 732     //=======================================================
 733     typedef itk::ImageFileWriter< DeformationFieldType >  FieldWriterType;
 734     typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
 735     fieldWriter->SetFileName( m_ArgsInfo.vf_arg );
 736     fieldWriter->SetInput( field );
 737     try
 738       {
 739         fieldWriter->Update();
 740       }
 741     catch( itk::ExceptionObject & excp )
 742       {
 743         std::cerr << "Exception thrown writing the DVF" << std::endl;
 744         std::cerr << excp << std::endl;
 745         return;
 746       }
 747
 748
 749     //=======================================================
 750     // Resample the moving image
 751     //=======================================================
 752     typedef itk::WarpImageFilter< MovingImageType, FixedImageType, DeformationFieldType >    WarpFilterType;
 753     typename WarpFilterType::Pointer warp = WarpFilterType::New();
 754
 755     warp->SetDeformationField( field );
 756     warp->SetInput( movingImageReader->GetOutput() );
 757     warp->SetOutputOrigin(  fixedImage->GetOrigin() );
 758     warp->SetOutputSpacing( fixedImage->GetSpacing() );
 759     warp->SetOutputDirection( fixedImage->GetDirection() );
 760     warp->SetEdgePaddingValue( 0.0 );
 761     warp->Update();
 762
 763
 764     //=======================================================
 765     // Write the warped image
 766     //=======================================================
 767     typedef itk::ImageFileWriter< FixedImageType >  WriterType;
 768     typename WriterType::Pointer      writer =  WriterType::New();
 769     writer->SetFileName( m_ArgsInfo.output_arg );
 770     writer->SetInput( warp->GetOutput()    );
 771
 772     try
 773       {
 774         writer->Update();
 775       }
 776     catch( itk::ExceptionObject & err )
 777       {
 778         std::cerr << "ExceptionObject caught !" << std::endl;
 779         std::cerr << err << std::endl;
 780         return;
 781       }
 782
 783
 784     //=======================================================
 785     // Calculate the difference after the deformable transform
 786     //=======================================================
 787     typedef clitk::DifferenceImageFilter<  FixedImageType, FixedImageType> DifferenceFilterType;
 788     if (m_ArgsInfo.after_given)
 789       {
 790         typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
 791         difference->SetValidInput( fixedImage );
 792         difference->SetTestInput( warp->GetOutput() );
 793
 794         try
 795           {
 796             difference->Update();
 797           }
 798         catch( itk::ExceptionObject & err )
 799           {
 800             std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
 801             std::cerr << err << std::endl;
 802             return;
 803           }
 804
 805         typename WriterType::Pointer differenceWriter=WriterType::New();
 806         differenceWriter->SetInput(difference->GetOutput());
 807         differenceWriter->SetFileName(m_ArgsInfo.after_arg);
 808         differenceWriter->Update();
 809
 810       }
 811
 812
 813     //=======================================================
 814     // Calculate the difference before the deformable transform
 815     //=======================================================
 816     if( m_ArgsInfo.before_given )
 817       {
 818
 819         typename FixedImageType::Pointer moving=FixedImageType::New();
 820         if (m_ArgsInfo.rigid_given)
 821           {
 822             typedef itk::ResampleImageFilter<MovingImageType, FixedImageType> ResamplerType;
 823             typename ResamplerType::Pointer resampler=ResamplerType::New();
 824             resampler->SetInput(movingImage);
 825             resampler->SetOutputOrigin(fixedImage->GetOrigin());
 826             resampler->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
 827             resampler->SetOutputSpacing(fixedImage->GetSpacing());
 828             resampler->SetDefaultPixelValue( 0. );
 829             //resampler->SetTransform(rigidTransform);
 830             resampler->Update();
 831             moving=resampler->GetOutput();
 832           }
 833         else
 834           moving=movingImage;
 835
 836         typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
 837         difference->SetValidInput( fixedImage );
 838         difference->SetTestInput( moving );
 839
 840         try
 841           {
 842             difference->Update();
 843           }
 844         catch( itk::ExceptionObject & err )
 845           {
 846             std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
 847             std::cerr << err << std::endl;
 848             return;
 849           }
 850
 851         typename WriterType::Pointer differenceWriter=WriterType::New();
 852         writer->SetFileName( m_ArgsInfo.before_arg  );
 853         writer->SetInput( difference->GetOutput()  );
 854         writer->Update( );
 855       }
 856
 857     return;
 858   }
 859 }
 860
 861 #endif // __clitkBSplineDeformableRegistrationGenericFilter_txx