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