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