registration/clitkBSplineDeformableTransformInitializer.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 __clitkBSplineDeformableTransformInitializer_txx
  19 #define __clitkBSplineDeformableTransformInitializer_txx
  20 #include "clitkBSplineDeformableTransformInitializer.h"
  21
  22 namespace clitk
  23 {
  24
  25
  26   template < class TTransform, class TImage >
  27   BSplineDeformableTransformInitializer<TTransform, TImage >
  28   ::BSplineDeformableTransformInitializer()
  29   {
  30     this->m_NumberOfControlPointsInsideTheImage.Fill( 5 );
  31     this->m_ControlPointSpacing.Fill(64.);
  32     this->m_ChosenSpacing.Fill(64.);
  33     m_NumberOfControlPointsIsGiven=false;
  34     m_ControlPointSpacingIsGiven=false;
  35     m_TransformRegionIsGiven=false;
  36     m_SamplingFactorIsGiven=false;
  37     m_SplineOrders.Fill(3);
  38     m_Parameters=NULL;
  39   }
  40
  41   template < class TTransform, class TImage >
  42   void
  43   BSplineDeformableTransformInitializer<TTransform, TImage >
  44   ::InitializeTransform()
  45   {
  46     // Sanity check:
  47     // The image is required for the region and spacing
  48     if( ! this->m_Image )
  49       {
  50         itkExceptionMacro( "Reference Image has not been set" );
  51         return;
  52       }
  53
  54     // A pointer to the transform is required
  55     if( ! this->m_Transform )
  56       {
  57         itkExceptionMacro( "Transform has not been set" );
  58         return;
  59       }
  60
  61     // If the image come from a filter, then update that filter.
  62     if( this->m_Image->GetSource() )
  63       {
  64         this->m_Image->GetSource()->Update();
  65       }
  66
  67
  68     //--------------------------------------
  69     // Spacing & Size on image
  70     //--------------------------------------
  71     SpacingType fixedImageSpacing = m_Image->GetSpacing();
  72     SizeType fixedImageSize=m_Image->GetLargestPossibleRegion().GetSize();
  73     typename RegionType::SizeType   gridBorderSize;
  74     typename RegionType::SizeType   totalGridSize;
  75
  76     // Only spacing given: adjust if necessary
  77     if (m_ControlPointSpacingIsGiven  && !m_NumberOfControlPointsIsGiven)
  78       {
  79         for(unsigned int r=0; r<InputDimension; r++)
  80           {
  81             // JV
  82             m_ChosenSpacing[r]= m_ControlPointSpacing[r];
  83             m_ControlPointSpacing[r]= ( round(m_ChosenSpacing[r]/fixedImageSpacing[r]) *fixedImageSpacing[r] ) ;
  84             m_NumberOfControlPointsInsideTheImage[r] = ceil( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] );
  85             if (  ( ceil( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] ) )
  86                  == ( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] ) )
  87               {
  88                 m_NumberOfControlPointsInsideTheImage[r] +=1;
  89                 DD("Adding control point");
  90               }
  91           }
  92       }
  93
  94     // Only number of CP given: adjust if necessary
  95     // JV -1 ?
  96     else if (!m_ControlPointSpacingIsGiven  && m_NumberOfControlPointsIsGiven)
  97       {
  98         for(unsigned int r=0; r<InputDimension; r++)
  99           {
 100             m_ChosenSpacing[r]=fixedImageSpacing[r]*( (double)(fixedImageSize[r])  /
 101                                                       (double)(m_NumberOfControlPointsInsideTheImage[r]) );
 102             m_ControlPointSpacing[r]= fixedImageSpacing[r]* ceil( (double)(fixedImageSize[r] - 1)  /
 103                                                                   (double)(m_NumberOfControlPointsInsideTheImage[r] - 1) );
 104           }
 105       }
 106
 107     // Both or none of Spacing and number of CP given: no adjustment adjust, just warnings
 108     else
 109       {
 110         for(unsigned int r=0; r<InputDimension; r++)
 111           {
 112             m_ChosenSpacing[r]= m_ControlPointSpacing[r];
 113             if (m_NumberOfControlPointsInsideTheImage[r]*m_ControlPointSpacing[r]< fixedImageSize[r]*fixedImageSpacing[r])
 114               {
 115                 std::cout<<"WARNING: Specified control point region ("<<m_NumberOfControlPointsInsideTheImage[r]*m_ControlPointSpacing[r]
 116                          <<"mm) does not cover the fixed image region ("<< fixedImageSize[r]*fixedImageSpacing[r]
 117                          <<"mm) for dimension "<<r<<"!" <<std::endl
 118                          <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
 119               }
 120             if (  fmod(m_ControlPointSpacing[r], fixedImageSpacing[r]) )
 121               {
 122                 std::cout<<"WARNING: Specified control point spacing for dimension "<<r
 123                          <<" does not allow exact representation of BLUT FFD!"<<std::endl
 124                          <<"Spacing ratio is non-integer: "<<m_ControlPointSpacing[r]/ fixedImageSpacing[r]<<std::endl
 125                          <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
 126               }
 127           }
 128       }
 129
 130     if (m_Verbose) std::cout<<"The chosen control point spacing "<<m_ChosenSpacing<<"..."<<std::endl;
 131     if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<m_ControlPointSpacing<<"..."<<std::endl;
 132     if (m_Verbose) std::cout<<"The number of (internal) control points is "<<m_NumberOfControlPointsInsideTheImage<<"..."<<std::endl;
 133
 134
 135     //--------------------------------------
 136     // Border size
 137     //--------------------------------------
 138     for(unsigned int r=0; r<InputDimension; r++) gridBorderSize[r]=m_SplineOrders[r]; // Border for spline order = 3 ( 1 lower, 2 upper )
 139     totalGridSize = m_NumberOfControlPointsInsideTheImage + gridBorderSize;
 140     if (m_Verbose) std::cout<<"The total control point grid size was set to "<<totalGridSize<<"..."<<std::endl;
 141     RegionType gridRegion;
 142     gridRegion.SetSize(totalGridSize);
 143
 144
 145     //--------------------------------------
 146     // Origin
 147     //--------------------------------------
 148     typedef typename TransformType::OriginType OriginType;
 149     OriginType fixedImageOrigin, gridOrigin;
 150
 151     // In case of non-zero index
 152     m_Image->TransformIndexToPhysicalPoint(m_Image->GetLargestPossibleRegion().GetIndex(),fixedImageOrigin);
 153     typename ImageType::DirectionType gridDirection = m_Image->GetDirection();
 154     SizeType orderShift;
 155
 156     // Shift depends on order
 157     for(unsigned int r=0; r<InputDimension; r++)
 158       {
 159         orderShift[r] = m_SplineOrders[r] / 2;
 160         gridOrigin[r] = fixedImageOrigin[r]- (double) orderShift[r]* m_ControlPointSpacing[r];
 161       }
 162     if (m_Verbose) std::cout<<"The control point grid origin was set to "<<gridOrigin<<"..."<<std::endl;
 163
 164
 165     //--------------------------------------
 166     // LUT sampling factors
 167     //--------------------------------------
 168     for (unsigned int i=0; i< InputDimension; i++)
 169       {
 170         if (m_Verbose) std::cout<<"For dimension "<<i<<", the ideal sampling factor (if integer) is a multitude of "
 171                                 << m_ControlPointSpacing[i]/ fixedImageSpacing[i]<<"..."<<std::endl;
 172         if ( !m_SamplingFactorIsGiven) m_SamplingFactors[i]= (int) ( m_ControlPointSpacing[i]/ fixedImageSpacing[i]);
 173         if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<m_SamplingFactors[i]<<"..."<<std::endl;
 174       }
 175
 176
 177     //--------------------------------------
 178     // Set
 179     //--------------------------------------
 180     this->m_Transform->SetSplineOrders(m_SplineOrders);
 181     this->m_Transform->SetGridRegion( gridRegion );
 182     this->m_Transform->SetGridOrigin( gridOrigin );
 183     this->m_Transform->SetGridSpacing( m_ControlPointSpacing );
 184     this->m_Transform->SetGridDirection( gridDirection );
 185     this->m_Transform->SetLUTSamplingFactors(m_SamplingFactors);
 186
 187   }
 188
 189   template < class TTransform, class TImage >
 190   void
 191   BSplineDeformableTransformInitializer<TTransform, TImage >
 192   ::SetInitialParameters( const std::string& s, ParametersType& params)
 193   {
 194     //---------------------------------------
 195     // Read Initial parameters
 196     //---------------------------------------
 197     typedef itk::ImageFileReader<CoefficientImageType> CoefficientReaderType;
 198     typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New();
 199     coeffReader->SetFileName(s);
 200     if(m_Verbose) std::cout<<"Reading initial coefficients from file: "<<s<<"..."<<std::endl;
 201     coeffReader->Update();
 202     typename CoefficientImageType::Pointer coefficientImage=coeffReader->GetOutput();
 203     this->SetInitialParameters(coefficientImage, params);
 204   }
 205
 206   template < class TTransform, class TImage >
 207   void
 208   BSplineDeformableTransformInitializer<TTransform, TImage >
 209   ::SetInitialParameters(const typename CoefficientImageType::Pointer coefficientImage, ParametersType& params)
 210   {
 211     // Keep a reference
 212     m_Parameters=&params;
 213
 214     // Resample
 215     typedef ResampleBSplineDeformableTransformImageFilter<CoefficientImageType,CoefficientImageType> ResamplerType;
 216     typename ResamplerType::Pointer resampler=ResamplerType::New();
 217     resampler->SetInput(coefficientImage);
 218     resampler->SetOutputParametersFromImage(m_Transform->GetCoefficientImage());
 219     if(m_Verbose) std::cout<<"Resampling initial coefficients..."<<std::endl;
 220     resampler->Update();
 221
 222     // Copy parameters into the existing array, I know its crappy
 223     typedef itk::ImageRegionConstIterator<CoefficientImageType> Iterator;
 224     Iterator it (resampler->GetOutput(), resampler->GetOutput()->GetLargestPossibleRegion() );
 225     it.GoToBegin();
 226     unsigned int i=0;
 227     while(! it.IsAtEnd())
 228       {
 229         for (unsigned int j=0; j<InputDimension;j++)
 230           params[i+j]=it.Get()[j];
 231
 232         ++it;
 233         i+=InputDimension;
 234       }
 235
 236     // JV pass the reference to the array !!
 237     this->m_Transform->SetParameters(params);
 238   }
 239
 240
 241
 242 }  // namespace itk
 243
 244 #endif