registration/clitkShapedBLUTSpatioTemporalDeformableTransformInitializer.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 __clitkShapedBLUTSpatioTemporalDeformableTransformInitializer_txx
  19 #define __clitkShapedBLUTSpatioTemporalDeformableTransformInitializer_txx
  20 #include "clitkShapedBLUTSpatioTemporalDeformableTransformInitializer.h"
  21
  22 namespace clitk
  23 {
  24
  25
  26   template < class TTransform, class TImage >
  27   ShapedBLUTSpatioTemporalDeformableTransformInitializer<TTransform, TImage >
  28   ::ShapedBLUTSpatioTemporalDeformableTransformInitializer()
  29   {
  30     this->m_NumberOfControlPointsInsideTheImage.Fill( 5 );
  31     this->m_ControlPointSpacing.Fill(64.);
  32     m_ControlPointSpacing[InputDimension-1]=2;
  33     this->m_ChosenSpacing.Fill(64.);
  34     m_ChosenSpacing[InputDimension-1]=2;
  35     m_ControlPointSpacing[InputDimension-1]=2;
  36     m_NumberOfControlPointsIsGiven=false;
  37     m_ControlPointSpacingIsGiven=false;
  38     m_TransformRegionIsGiven=false;
  39     m_SamplingFactorIsGiven=false;
  40     m_SplineOrders.Fill(3);
  41     m_BC1=true;
  42     m_BC2=true;
  43     m_TrajectoryShape=0;
  44   }
  45
  46   template < class TTransform, class TImage >
  47   void
  48   ShapedBLUTSpatioTemporalDeformableTransformInitializer<TTransform, TImage >
  49   ::InitializeTransform()
  50   {
  51     // Sanity check:
  52     // The image is required for the region and spacing
  53     if( ! this->m_Image )
  54       {
  55         itkExceptionMacro( "Reference Image has not been set" );
  56         return;
  57       }
  58
  59     // A pointer to the transform is required
  60     if( ! this->m_Transform )
  61       {
  62         itkExceptionMacro( "Transform has not been set" );
  63         return;
  64       }
  65
  66     // If the image come from a filter, then update that filter.
  67     if( this->m_Image->GetSource() )
  68       {
  69         this->m_Image->GetSource()->Update();
  70       }
  71
  72
  73     //--------------------------------------
  74     // Spacing & Size on image
  75     //--------------------------------------
  76     SpacingType fixedImageSpacing = m_Image->GetSpacing();
  77     SizeType fixedImageSize=m_Image->GetLargestPossibleRegion().GetSize();
  78     typename RegionType::SizeType   gridBorderSize;
  79     typename RegionType::SizeType   totalGridSize;
  80
  81     // Only spacing given: adjust if necessary
  82     if (m_ControlPointSpacingIsGiven  && !m_NumberOfControlPointsIsGiven)
  83       {
  84         for(unsigned int r=0; r<InputDimension; r++)
  85           {
  86             // JV
  87             m_ChosenSpacing[r]= m_ControlPointSpacing[r];
  88             m_ControlPointSpacing[r]= ( itk::Math::Round<double>(m_ChosenSpacing[r]/fixedImageSpacing[r]) *fixedImageSpacing[r] ) ;
  89             m_NumberOfControlPointsInsideTheImage[r] = ceil( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] );
  90             if ( (  ( ceil( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] ) )
  91                  == ( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] ) )
  92                  &&  (r!= InputDimension-1)  )
  93               {
  94                 m_NumberOfControlPointsInsideTheImage[r] +=1;
  95                 DD("Adding control point");
  96               }
  97           }
  98       }
  99
 100     // Only number of CP given: adjust if necessary
 101     // JV -1 ?
 102     else if (!m_ControlPointSpacingIsGiven  && m_NumberOfControlPointsIsGiven)
 103       {
 104         for(unsigned int r=0; r<InputDimension; r++)
 105           {
 106             m_ChosenSpacing[r]=fixedImageSpacing[r]*( (double)(fixedImageSize[r])  /
 107                                                       (double)(m_NumberOfControlPointsInsideTheImage[r]) );
 108             m_ControlPointSpacing[r]= fixedImageSpacing[r]* ceil( (double)(fixedImageSize[r] - 1)  /
 109                                                                   (double)(m_NumberOfControlPointsInsideTheImage[r] - 1) );
 110           }
 111       }
 112
 113     // Both or none of Spacing and number of CP given: no adjustment adjust, just warnings
 114     else
 115       {
 116         for(unsigned int r=0; r<InputDimension; r++)
 117           {
 118             m_ChosenSpacing[r]= m_ControlPointSpacing[r];
 119             if (m_NumberOfControlPointsInsideTheImage[r]*m_ControlPointSpacing[r]< fixedImageSize[r]*fixedImageSpacing[r])
 120               {
 121                 std::cout<<"WARNING: Specified control point region ("<<m_NumberOfControlPointsInsideTheImage[r]*m_ControlPointSpacing[r]
 122                          <<"mm) does not cover the fixed image region ("<< fixedImageSize[r]*fixedImageSpacing[r]
 123                          <<"mm) for dimension "<<r<<"!" <<std::endl
 124                          <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
 125               }
 126             if (  fmod(m_ControlPointSpacing[r], fixedImageSpacing[r]) )
 127               {
 128                 std::cout<<"WARNING: Specified control point spacing for dimension "<<r
 129                          <<" does not allow exact representation of BLUT FFD!"<<std::endl
 130                          <<"Spacing ratio is non-integer: "<<m_ControlPointSpacing[r]/ fixedImageSpacing[r]<<std::endl
 131                          <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
 132               }
 133           }
 134       }
 135
 136     //--------------------------------------
 137     // LUT sampling factors
 138     //--------------------------------------
 139     for (unsigned int i=0; i< InputDimension; i++)
 140       {
 141         if (m_Verbose) std::cout<<"For dimension "<<i<<", the ideal sampling factor (if integer) is a multitude of "
 142                                 << m_ControlPointSpacing[i]/ fixedImageSpacing[i]<<"..."<<std::endl;
 143         if ( !m_SamplingFactorIsGiven) m_SamplingFactors[i]= (int) ( m_ControlPointSpacing[i]/ fixedImageSpacing[i]);
 144       }
 145
 146     //--------------------------------------
 147     // Set the transform properties
 148     //--------------------------------------
 149     RegionType gridRegion;
 150     OriginType gridOrigin;
 151     DirectionType gridDirection;
 152
 153     //--------------------------------------
 154     // Border size
 155     //--------------------------------------
 156     for(unsigned int r=0; r<InputDimension; r++) gridBorderSize[r]=m_SplineOrders[r]; // Border for spline order = 3 ( 1 lower, 2 upper )
 157     totalGridSize = m_NumberOfControlPointsInsideTheImage + gridBorderSize;
 158
 159     //--------------------------------------
 160     // Origin
 161     //--------------------------------------
 162     OriginType fixedImageOrigin;
 163
 164     // In case of non-zero index
 165     m_Image->TransformIndexToPhysicalPoint(m_Image->GetLargestPossibleRegion().GetIndex(),fixedImageOrigin);
 166     gridDirection = m_Image->GetDirection();
 167     itk::FixedArray<int, InputDimension> orderShift;
 168
 169     // Spacing is 2.5 : manually modify the props
 170     if( !(m_Transform->GetTransformShape()%2) )
 171       {
 172         m_ControlPointSpacing[InputDimension-1]=2.5;
 173         m_SamplingFactors[InputDimension-1]=5;
 174         totalGridSize[InputDimension-1]-=1;
 175         gridRegion.SetSize(totalGridSize);
 176       }
 177
 178
 179     switch(m_Transform->GetTransformShape()){
 180
 181       // The egg shape
 182     case 0:
 183     case 1:
 184       {
 185         if (m_Verbose) std::cout<<"Using the egg shape..."<<std::endl;
 186
 187         //--------------------------------------
 188         // Border size
 189         //--------------------------------------
 190         // Boundary condition1: cyclic
 191         if (m_BC1)
 192           totalGridSize[InputDimension-1]-=3;
 193
 194         // Boundary condition1: reference is on trajectory
 195         if (m_BC2)
 196           totalGridSize[InputDimension-1]-=1;
 197
 198         gridRegion.SetSize(totalGridSize);
 199
 200         //--------------------------------------
 201         // Origin
 202         //--------------------------------------
 203         // Shift depends on order
 204         for(unsigned int r=0; r<InputDimension; r++)
 205           {
 206             orderShift[r] = m_SplineOrders[r] / 2;
 207             if (r==InputDimension-1)
 208               {
 209                 // Boundary condition1: cyclic
 210                 if (m_BC1)
 211                   orderShift[r]-=1;
 212
 213                 // Boundary condition1: reference is on trajectory
 214                 if (m_BC2)
 215                   orderShift[r]-=1;
 216               }
 217
 218             gridOrigin[r] = fixedImageOrigin[r]- (double) orderShift[r]* m_ControlPointSpacing[r];
 219           }
 220         break;
 221       }
 222
 223       // The rabit: a diamond with fixed extreme values
 224     case 2:
 225     case 3:
 226       {
 227
 228         if (m_Verbose) std::cout<<"Using the rabbit shape..."<<std::endl;
 229
 230         //--------------------------------------
 231         // Border size
 232         //--------------------------------------
 233         // Boundary condition1: cyclic
 234         if (m_BC1)
 235           totalGridSize[InputDimension-1]-=1;
 236
 237         // Boundary condition1: reference is on trajectory
 238         if (m_BC2)
 239           totalGridSize[InputDimension-1]-=1;
 240
 241         // Extreme values are constrained
 242         totalGridSize[InputDimension-1]-=2;
 243
 244         gridRegion.SetSize(totalGridSize);
 245
 246         //--------------------------------------
 247         // Origin
 248         //--------------------------------------
 249         // Shift depends on order
 250         for(unsigned int r=0; r<InputDimension; r++)
 251           {
 252             orderShift[r] = m_SplineOrders[r] / 2;
 253             if (r==InputDimension-1)
 254               {
 255                 // Boundary condition1: cyclic
 256                 orderShift[r]=0;
 257               }
 258             gridOrigin[r] = fixedImageOrigin[r]- (double) orderShift[r]* m_ControlPointSpacing[r];
 259           }
 260
 261         break;
 262       }
 263
 264       // The sputnik: a diamond with one indepent extreme value and 4DOF
 265     case 4:
 266     case 5:
 267       {
 268         if (m_Verbose) std::cout<<"Using the sputnik shape..."<<std::endl;
 269
 270         //--------------------------------------
 271         // Border size
 272         //--------------------------------------
 273         // Boundary condition1: cyclic
 274         if (m_BC1)
 275           totalGridSize[InputDimension-1]-=1;
 276
 277         // Boundary condition1: reference is on trajectory
 278         if (m_BC2)
 279           totalGridSize[InputDimension-1]-=1;
 280
 281         // One extreme value is fixed
 282         totalGridSize[InputDimension-1]-=1;
 283         gridRegion.SetSize(totalGridSize);
 284
 285         //--------------------------------------
 286         // Origin
 287         //--------------------------------------
 288         // Shift depends on order
 289         for(unsigned int r=0; r<InputDimension; r++)
 290           {
 291             orderShift[r] = m_SplineOrders[r] / 2;
 292             if (r==InputDimension-1)
 293               {
 294                 // Boundary condition1: cyclic
 295                 orderShift[r]=0;
 296               }
 297             gridOrigin[r] = fixedImageOrigin[r]- (double) orderShift[r]* m_ControlPointSpacing[r];
 298           }
 299         break;
 300       }
 301
 302       // The diamond
 303     case 6:
 304     case 7:
 305       {
 306         if (m_Verbose) std::cout<<"Using the diamond shape..."<<std::endl;
 307
 308         //--------------------------------------
 309         // Border size
 310         //--------------------------------------
 311         // Boundary condition1: cyclic
 312         if (m_BC1)
 313           totalGridSize[InputDimension-1]-=1;
 314
 315         // Boundary condition1: reference is on trajectory
 316         if (m_BC2)
 317           totalGridSize[InputDimension-1]-=1;
 318
 319         gridRegion.SetSize(totalGridSize);
 320
 321         //--------------------------------------
 322         // Origin
 323         //--------------------------------------
 324         // Shift depends on order
 325         for(unsigned int r=0; r<InputDimension; r++)
 326           {
 327             orderShift[r] = m_SplineOrders[r] / 2;
 328             gridOrigin[r] = fixedImageOrigin[r]- (double) orderShift[r]* m_ControlPointSpacing[r];
 329           }
 330
 331         break;
 332
 333       } //end case 7
 334
 335    // The sputnik: a diamond with one indepent extreme value and 4DOF
 336     case 8:
 337     case 9:
 338       {
 339         if (m_Verbose) std::cout<<"Using the sputnik shape..."<<std::endl;
 340
 341         //--------------------------------------
 342         // Border size
 343         //--------------------------------------
 344         // Boundary condition1: cyclic
 345         if (m_BC1)
 346           totalGridSize[InputDimension-1]-=1;
 347
 348         // Boundary condition1: reference is on trajectory
 349         if (m_BC2)
 350           totalGridSize[InputDimension-1]-=1;
 351
 352         // One extreme value is fixed
 353         totalGridSize[InputDimension-1]-=1;
 354         gridRegion.SetSize(totalGridSize);
 355
 356         //--------------------------------------
 357         // Origin
 358         //--------------------------------------
 359         // Shift depends on order
 360         for(unsigned int r=0; r<InputDimension; r++)
 361           {
 362             orderShift[r] = m_SplineOrders[r] / 2;
 363             if (r==InputDimension-1)
 364               {
 365                 // Boundary condition1: cyclic
 366                 orderShift[r]=0;
 367               }
 368             gridOrigin[r] = fixedImageOrigin[r]- (double) orderShift[r]* m_ControlPointSpacing[r];
 369           }
 370         break;
 371       }
 372
 373     default:
 374       {
 375         DD("Shape not available");
 376       }
 377
 378     }
 379
 380
 381
 382     if (m_Verbose) std::cout<<"The chosen control point spacing "<<m_ChosenSpacing<<"..."<<std::endl;
 383     if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<m_ControlPointSpacing<<"..."<<std::endl;
 384     if (m_Verbose) std::cout<<"The number of (internal) control points is "<<m_NumberOfControlPointsInsideTheImage<<"..."<<std::endl;
 385     if (m_Verbose) std::cout<<"Setting sampling factors to "<<m_SamplingFactors<<"..."<<std::endl;
 386     if (m_Verbose) std::cout<<"The total control point grid size was set to "<<totalGridSize<<"..."<<std::endl;
 387     if (m_Verbose) std::cout<<"The control point grid origin was set to "<<gridOrigin<<"..."<<std::endl;
 388
 389     //--------------------------------------
 390     // Set
 391     //--------------------------------------
 392     this->m_Transform->SetSplineOrders(m_SplineOrders);
 393     this->m_Transform->SetGridSpacing( m_ControlPointSpacing );
 394     this->m_Transform->SetGridOrigin( gridOrigin );
 395     this->m_Transform->SetGridDirection( gridDirection );
 396     this->m_Transform->SetGridRegion( gridRegion );
 397     this->m_Transform->SetLUTSamplingFactors(m_SamplingFactors);
 398
 399   }
 400
 401   template < class TTransform, class TImage >
 402   void
 403   ShapedBLUTSpatioTemporalDeformableTransformInitializer<TTransform, TImage >
 404   ::SetInitialParameters( const std::string& s, ParametersType& params)
 405   {
 406     //---------------------------------------
 407     // Read Initial parameters
 408     //---------------------------------------
 409     typedef itk::ImageFileReader<CoefficientImageType> CoefficientReaderType;
 410     typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New();
 411     coeffReader->SetFileName(s);
 412     if(m_Verbose) std::cout<<"Reading initial coefficients from file: "<<s<<"..."<<std::endl;
 413     coeffReader->Update();
 414     typename CoefficientImageType::Pointer coefficientImage=coeffReader->GetOutput();
 415     this->SetInitialParameters(coefficientImage, params);
 416   }
 417
 418   template < class TTransform, class TImage >
 419   void
 420   ShapedBLUTSpatioTemporalDeformableTransformInitializer<TTransform, TImage >
 421   ::SetInitialParameters(const typename CoefficientImageType::Pointer coefficientImage, ParametersType& params)
 422   {
 423     // Keep a reference
 424     m_Parameters=&params;
 425
 426     // Resample
 427     typedef ResampleBSplineDeformableTransformImageFilter<CoefficientImageType,CoefficientImageType> ResamplerType;
 428     typename ResamplerType::Pointer resampler=ResamplerType::New();
 429     resampler->SetVerbose(m_Verbose);
 430     resampler->SetInput(coefficientImage);
 431     resampler->SetOutputParametersFromImage(m_Transform->GetCoefficientImage());
 432     if(m_Verbose) std::cout<<"Resampling initial coefficients..."<<std::endl;
 433     resampler->Update();
 434
 435     // Copy parameters into the existing array, I know its crappy
 436     typedef itk::ImageRegionConstIterator<CoefficientImageType> Iterator;
 437     Iterator it (resampler->GetOutput(), resampler->GetOutput()->GetLargestPossibleRegion() );
 438     it.GoToBegin();
 439     unsigned int i=0;
 440     while(! it.IsAtEnd())
 441       {
 442         for (unsigned int j=0; j<InputDimension-1;j++)
 443           params[i+j]=it.Get()[j];
 444
 445         ++it;
 446         i+=InputDimension-1;
 447       }
 448
 449     // JV pass the reference to the array !!
 450     // JV Update the padded version !!
 451     this->m_Transform->SetParameters(params);
 452   }
 453
 454
 455   //   template < class TTransform, class TImage >
 456   //   void
 457   //   ShapedBLUTSpatioTemporalDeformableTransformInitializer<TTransform, TImage >
 458   //   ::SetInitialPaddedParameters( const std::string& s, ParametersType& params)
 459   //   {
 460   //     //---------------------------------------
 461   //     // Read Initial parameters
 462   //     //---------------------------------------
 463   //     typedef itk::ImageFileReader<CoefficientImageType> CoefficientReaderType;
 464   //     typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New();
 465   //     coeffReader->SetFileName(s);
 466   //     if(m_Verbose) std::cout<<"Reading initial padded coefficients from file: "<<s<<"..."<<std::endl;
 467   //     coeffReader->Update();
 468   //     typename CoefficientImageType::Pointer coefficientImage=coeffReader->GetOutput();
 469   //     this->SetInitialPaddedParameters(coefficientImage, params);
 470   //   }
 471
 472   //   template < class TTransform, class TImage >
 473   //   void
 474   //   ShapedBLUTSpatioTemporalDeformableTransformInitializer<TTransform, TImage >
 475   //   ::SetInitialPaddedParameters(const typename CoefficientImageType::Pointer coefficientImage, ParametersType& params)
 476   //   {
 477   //     // Resample
 478   //     typedef ResampleBSplineSpatioTemporalDeformableTransformImageFilter<CoefficientImageType,CoefficientImageType> ResamplerType;
 479   //     typename ResamplerType::Pointer resampler=ResamplerType::New();
 480   //     resampler->SetInput(coefficientImage);
 481   //     resampler->SetOutputParametersFromImage(m_Transform->GetCoefficientImage());
 482   //     if(m_Verbose) std::cout<<"Resampling initial padded coefficients..."<<std::endl;
 483   //     resampler->Update();
 484
 485   //     // Copy parameters into the existing array, I know its crappy
 486   //     typedef itk::ImageRegionConstIterator<CoefficientImageType> Iterator;
 487   //     Iterator it (resampler->GetOutput(), resampler->GetOutput()->GetLargestPossibleRegion() );
 488   //     it.GoToBegin();
 489   //     unsigned int i=0;
 490   //     while(! it.IsAtEnd())
 491   //       {
 492   //    for (unsigned int j=0; j<InputDimension;j++)
 493   //      params[i+j]=it.Get()[j];
 494
 495   //    ++it;
 496   //    i+=InputDimension;
 497   //       }
 498   //   }
 499
 500
 501 }  // namespace itk
 502
 503 #endif