itk/itkBSplineInterpolateImageFunctionWithLUT.txx.original

   1 #ifndef _ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_
   2 #define _ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_
   3
   4 /* =========================================================================
   5
   6   @file   itkBSplineInterpolateImageFunctionWithLUT.txx
   7   @author David Sarrut <david.sarrut@creatis.insa-lyon.fr>
   8
   9   Copyright (c)
  10   * CREATIS (Centre de Recherche et d'Applications en Traitement de l'Image).
  11     All rights reserved. See Doc/License.txt or
  12     http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
  13   * L�on B�rard cancer center, 28 rue La�nnec, 69373 Lyon cedex 08, France
  14   * http://www.creatis.insa-lyon.fr/rio
  15
  16   This software is distributed WITHOUT ANY WARRANTY; without even the
  17   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  18   PURPOSE.  See the above copyright notices for more information.
  19
  20 ========================================================================= */
  21 #include "itkCastImageFilter.h"
  22
  23 //====================================================================
  24 template <class TImageType, class TCoordRep, class TCoefficientType>
  25 BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
  26 BSplineInterpolateImageFunctionWithLUT():Superclass() {
  27   // Set default values
  28   for(int i=0; i<TImageType::ImageDimension; i++) mOutputSpacing[i] = -1;
  29   SetLUTSamplingFactor(10);
  30   SetSplineOrder(3);
  31   mWeightsAreUpToDate = false;
  32   // Following need to be pointer beacause values are updated into
  33   // "const" function
  34   mIntrinsecError  = new double;
  35   mNumberOfError = new long;
  36   mIntrinsecErrorMax = new double;
  37   mInputIsCoef = false;
  38 }
  39 //====================================================================
  40
  41 //====================================================================
  42 template <class TImageType, class TCoordRep, class TCoefficientType>
  43 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
  44 SetLUTSamplingFactor(int s) {
  45   for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s;
  46   mWeightsAreUpToDate = false;
  47 }
  48
  49 //====================================================================
  50
  51 //====================================================================
  52 template <class TImageType, class TCoordRep, class TCoefficientType>
  53 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
  54 SetOutputSpacing(const SpacingType & s) {
  55   for(int i=0; i<TImageType::ImageDimension; i++)
  56     mOutputSpacing[i] = s[i];
  57   mWeightsAreUpToDate = false;
  58 }
  59 //====================================================================
  60
  61 //====================================================================
  62 template <class TImageType, class TCoordRep, class TCoefficientType>
  63 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
  64 SetSplineOrder(unsigned int SplineOrder) {
  65   Superclass::SetSplineOrder(SplineOrder);
  66   // Compute support and half size
  67   static const int d = TImageType::ImageDimension;
  68   for(int l=0; l<d; l++) {
  69     mSplineOrders[l]= SplineOrder;
  70     mSupport[l] = SplineOrder+1;
  71     if (mSupport[l] % 2 == 0) { // support is even
  72       mHalfSupport[l] = mSupport[l]/2-1;
  73     }
  74     else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
  75   }
  76   mSupportSize = 1;
  77   for(int l=0; l<d; l++) {
  78     mSupportSize *= mSupport[l];
  79   }
  80   mWeightsAreUpToDate = false;
  81 }
  82 //====================================================================
  83
  84 //JV
  85 //====================================================================
  86 template <class TImageType, class TCoordRep, class TCoefficientType>
  87 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
  88 SetSplineOrders(SizeType SplineOrders) {
  89   mSplineOrders=SplineOrders;
  90   DD(mSplineOrders);
  91 // Compute support and half size
  92   static const int d = TImageType::ImageDimension;
  93   for(int l=0; l<d; l++) {
  94     mSupport[l] = mSplineOrders[l]+1;
  95     if (mSupport[l] % 2 == 0) { // support is even
  96       mHalfSupport[l] = mSupport[l]/2-1;
  97     }
  98     else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
  99   }
 100   mSupportSize = 1;
 101   for(int l=0; l<d; l++) {
 102     mSupportSize *= mSupport[l];
 103   }
 104   mWeightsAreUpToDate = false;
 105 }
 106 //====================================================================
 107
 108 //====================================================================
 109 template <class TImageType, class TCoordRep, class TCoefficientType>
 110 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
 111 SetInputImage(const TImageType * inputData) {
 112
 113   // Call superclass SetInputImage and return if input is null
 114   // clitk::Timer t;
 115   //   t.Start();
 116   //ds
 117   // this->InterpolateImageFunction<TImageType,TCoordRep>::SetInputImage(inputData);
 118
 119   //==============================
 120   if (!mInputIsCoef) {
 121     Superclass::SetInputImage(inputData);
 122     //  t.Stop();
 123   }
 124   //==============================
 125   else {
 126     // If inputData is not float/double, we do not want to set
 127     // mInputIsCoef to true, but the compiler goes here ...
 128
 129     if (typeid(typename TImageType::PixelType) != typeid(TCoefficientType)) {
 130       std::cerr << "Input should be float or double to be set as 'coefficients' with SetInputImageIsCoefficient(true) and should be the same type than TCoefficientType" << std::endl;
 131       exit(0);
 132     }
 133     else {
 134       itkStaticConstMacro(ImageDimension, unsigned int, Superclass::ImageDimension);
 135       typedef Image<TCoefficientType, ImageDimension> CoefficientImageType;
 136       typedef itk::CastImageFilter< TImageType, CoefficientImageType> CastFilterType;
 137       typename CastFilterType::Pointer castFilter=CastFilterType::New();
 138       castFilter->SetInput(inputData);
 139       castFilter->Update();
 140       this->m_Coefficients=castFilter->GetOutput();
 141       if ( this->m_Coefficients.IsNotNull())
 142         this->m_DataLength = this->m_Coefficients->GetBufferedRegion().GetSize();
 143     }
 144
 145 //     //JV input should be double or float. To keep it generic: cast to the CoefficientImageType
 146 //     /** Dimension underlying input image. */
 147 //     itkStaticConstMacro(ImageDimension, unsigned int, Superclass::ImageDimension);
 148 //     /** Internal Coefficient typedef support */
 149 //     DD( ImageDimension);
 150 //     typedef Image<TCoefficientType,ImageDimension> CoefficientImageType;
 151 //     typedef itk::CastImageFilter< TImageType, CoefficientImageType> CastFilterType;
 152 //     typename CastFilterType::Pointer castFilter=CastFilterType::New();
 153 //     castFilter->SetInput(inputData);
 154 //     castFilter->Update();
 155 //      DS this->InterpolateImageFunction<TImageType,TCoordRep>::m_Coefficients = inputData;
 156
 157 //     //JV cast to the coeff type, don't tell my mother
 158 //     this->InterpolateImageFunction<TImageType,TCoordRep>::SetInputImage(inputData);
 159 //     this->m_Coefficients=castFilter->GetOutput();
 160 //     DS ??????    this->m_Coefficients = inputData;
 161
 162 //     //   //============================== old not generic solution
 163 //     //      this->InterpolateImageFunction<TImageType,TCoordRep>::SetInputImage(inputData);
 164 //          // DS this->InterpolateImageFunction<TImageType,TCoordRep>::m_Coefficients = inputData;
 165 //     //     this->m_Coefficients = inputData;
 166 //     //   //==============================
 167 //     if ( this->m_Coefficients.IsNotNull()) {
 168 //       this->m_DataLength = this->m_Coefficients->GetBufferedRegion().GetSize();
 169 //     }
 170   }
 171
 172   //   t.Stop();
 173   if (!inputData) return;
 174  //  mCoefficientTimer = t;
 175
 176   // Compute Memory offset inside coefficients images (for looping over coefficients)
 177   static const int d = TImageType::ImageDimension;
 178   mInputMemoryOffset[0] = 1;
 179   for(int l=1; l<d; l++) {
 180     mInputMemoryOffset[l] =
 181       mInputMemoryOffset[l-1]*this->m_Coefficients->GetLargestPossibleRegion().GetSize(l-1);
 182   }
 183
 184   // Compute mSupportOffset according to input size
 185   mSupportOffset.resize(mSupportSize);
 186   mSupportIndex.resize(mSupportSize);
 187   for(int l=0; l<d; l++) mSupportIndex[0][l] = 0;
 188   for(unsigned int k=0; k<mSupportSize; k++) {
 189     // Get memory offset
 190     mSupportOffset[k] = Index2Offset<TImageType::ImageDimension>(mSupportIndex[k], mInputMemoryOffset);
 191     // next coefficient index
 192     if (k != mSupportSize-1) {
 193       for(int l=0; l<d; l++) mSupportIndex[k+1][l] = mSupportIndex[k][l];
 194       int l=0;
 195       bool stop = false;
 196       while (!stop) {
 197         mSupportIndex[k+1][l]++;
 198         if (mSupportIndex[k+1][l] == (int)mSupport[l]) { //ds supportindex=uint et support=int
 199           mSupportIndex[k+1][l] = 0;
 200           l++;
 201         }
 202         else stop = true;
 203       }
 204     }
 205   }
 206
 207   // Check
 208   for(int l=0; l<d; l++) {
 209     if (mOutputSpacing[l] == -1) {
 210       std::cerr << "Please use SetOutputSpacing before using bslut interpolator" << std::endl;
 211       exit(0);
 212     }
 213   }
 214
 215   // Compute BSpline weights if not up to date
 216   if (!mWeightsAreUpToDate) UpdatePrecomputedWeights();
 217
 218   // Initialisation
 219   *mIntrinsecErrorMax = *mIntrinsecError = 0.0;
 220   *mNumberOfError = 0;
 221 }
 222 //====================================================================
 223
 224 //====================================================================
 225 template <class TImageType, class TCoordRep, class TCoefficientType>
 226 void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
 227 UpdatePrecomputedWeights() {
 228 //   mLUTTimer.Reset();
 229 //   mLUTTimer.Start();
 230   mWeightsCalculator.SetSplineOrders(mSplineOrders);
 231   mWeightsCalculator.SetSamplingFactors(mSamplingFactors);
 232   mWeightsCalculator.ComputeTensorProducts();
 233   mWeightsAreUpToDate = true;
 234   //DS
 235   //   coef = new TCoefficientType[mSupportSize];
 236   //     mCorrectedSupportIndex.resize(mSupportSize);
 237   //     mCorrectedSupportOffset.resize(mSupportSize);
 238  //  mLUTTimer.Stop();
 239 //   mLUTTimer.Print("LUT      \t");
 240 }
 241 //====================================================================
 242
 243 //====================================================================
 244 template <class TImageType, class TCoordRep, class TCoefficientType>
 245 typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::IndexType
 246 BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
 247 GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const {
 248
 249   /*
 250     WARNING : sometimes, a floating number x could not really be
 251     represented in memory. In this case, the difference between x and
 252     floor(x) can be almost 1 (instead of 0).  So we take into account
 253     such special case, otherwise it could lead to annoying results.
 254   */
 255   // static const TCoefficientType tiny = 1.0e-7;
 256   IndexType index;
 257
 258   for(int l=0; l<TImageType::ImageDimension; l++) {
 259   // bool mChange = false;
 260
 261   // Compute t1 = distance to floor
 262   TCoefficientType t1 = x[l]- vcl_floor(x[l]);
 263
 264   // Compute index in precomputed weights table
 265   TCoefficientType t2 = mSamplingFactors[l]*t1;
 266   index[l] = (IndexValueType)lrint(t2);
 267
 268   // For even order : test if too close to 0.5 (but lower). In this
 269   // case : take the next coefficient
 270   if (!(mSplineOrders[l] & 1)) {
 271     if (t1<0.5) {
 272       if (mSamplingFactors[l] & 1) {
 273         if (index[l] ==  (int) mSamplingFactors[l]/2+1) EvaluateIndex[l] = EvaluateIndex[l]+1;
 274       }
 275       else if (index[l] == (int) mSamplingFactors[l]/2) EvaluateIndex[l] = EvaluateIndex[l]+1;
 276     }
 277   }
 278
 279   // Statistics (to be removed)
 280   /*
 281   *mIntrinsecError += fabs(index[l]-t2);
 282   (*mNumberOfError)++;
 283   if (fabs(index[l]-t2)> *mIntrinsecErrorMax) *mIntrinsecErrorMax = fabs(index[l]-t2);
 284   */
 285
 286   // When to close to 1, take the next coefficient for ordd order, but
 287   // only change index for odd
 288     if (index[l] == (int)mSamplingFactors[l]) {
 289       index[l] = 0;
 290       if (mSplineOrders[l] & 1) EvaluateIndex[l] = EvaluateIndex[l]+1;
 291     }
 292   }
 293
 294   // The end
 295   return index;
 296 }
 297 //====================================================================
 298
 299 //====================================================================
 300 template <class TImageType, class TCoordRep, class TCoefficientType>
 301 typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::OutputType
 302 BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
 303 EvaluateAtContinuousIndex(const ContinuousIndexType & x) const {
 304
 305   // For shorter coding
 306   static const unsigned int d = TImageType::ImageDimension;
 307
 308   // Compute the indexe of the first interpolation coefficient in the
 309   // coefficient image
 310   IndexType EvaluateIndex;
 311   long indx;
 312   for (unsigned int l=0; l<d; l++)  {
 313     if (mSplineOrders[l] & 1) {  // Use this index calculation for odd splineOrder (like cubic)
 314       indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
 315       EvaluateIndex[l] = indx;
 316     }
 317     else { // Use this index calculation for even splineOrder
 318       if (mSplineOrders[l] == 0) EvaluateIndex[l] = (long)rint(x[l]);
 319       else {
 320         indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
 321         EvaluateIndex[l] = indx;
 322       }
 323     }
 324   }
 325
 326   // Compute index of precomputed weights and get pointer to first weights
 327   const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, EvaluateIndex);
 328   const TCoefficientType * pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
 329
 330   // Check boundaries
 331   bool mBoundaryCase = false;
 332   for (unsigned int l=0; l<d; l++) {
 333     if ((EvaluateIndex[l] < 0) ||
 334         (EvaluateIndex[l]+mSupport[l]) >= this->m_Coefficients->GetLargestPossibleRegion().GetSize(l)) {
 335       mBoundaryCase = true;
 336     }
 337   }
 338
 339   // Pointer to support offset
 340   const int * psupport;
 341
 342   // Special case for boundary (to be changed ....)
 343   std::vector<int> mCorrectedSupportOffset;//(mSupportSize);
 344   if (mBoundaryCase) {
 345     //    return -1000;
 346     //std::vector<TCoefficientType> coef(mSupportSize);
 347     // DD(EvaluateIndex);
 348      std::vector<IndexType> mCorrectedSupportIndex;//(mSupportSize);
 349     mCorrectedSupportIndex.resize(mSupportSize);
 350     mCorrectedSupportOffset.resize(mSupportSize);
 351     for(unsigned int i=0; i<mSupportSize; i++) {
 352       // DD(mSupportIndex[i]);
 353       for (unsigned int l=0; l<d; l++) {
 354         long a = mSupportIndex[i][l] + EvaluateIndex[l];
 355         long b = this->m_Coefficients->GetLargestPossibleRegion().GetSize(l);
 356         // DD(a);
 357         //      DD(b);
 358         long d2 = 2 * b - 2;
 359         if (a < 0) {
 360           mCorrectedSupportIndex[i][l] = -a - d2*(-a/d2) - EvaluateIndex[l];//mSupportIndex[i][l]-a;
 361         }
 362         else {
 363           if (a>=b) {
 364             mCorrectedSupportIndex[i][l] = d2 - a - EvaluateIndex[l];
 365           }
 366           else {
 367             mCorrectedSupportIndex[i][l] = mSupportIndex[i][l]; //a - d2*(a/d2) - EvaluateIndex[l];
 368           }
 369           /*
 370           if (a>=b) {
 371             mCorrectedSupportIndex[i][l] = d2 - a - EvaluateIndex[l];//mSupportIndex[i][l] - (a-(b-1));
 372           }
 373           else {
 374             mCorrectedSupportIndex[i][l] = mSupportIndex[i][l];
 375           }
 376           */
 377         }
 378       }
 379       // DD(mCorrectedSupportIndex[i]);
 380       mCorrectedSupportOffset[i] = Index2Offset<TImageType::ImageDimension>(mCorrectedSupportIndex[i], mInputMemoryOffset);
 381     }
 382     // for (unsigned int l=0; l<d; l++) {
 383     //       EvaluateIndex[l] = EvaluateIndex[l] + mCorrectedSupportIndex[0][l];
 384     //     }
 385     // DD(EvaluateIndex);
 386     psupport = &mCorrectedSupportOffset[0];
 387   }
 388   else {
 389     psupport = &mSupportOffset[0];
 390   }
 391
 392   // Get pointer to first coefficient. Even if in some boundary cases,
 393   // EvaluateIndex is out of the coefficient image,
 394   const TCoefficientType * pcoef = &(this->m_Coefficients->GetPixel(EvaluateIndex));
 395
 396   // Main loop over BSpline support
 397   TCoefficientType interpolated = 0.0;
 398   for (unsigned int p=0; p<mSupportSize; p++) {
 399     interpolated += pcoef[*psupport] * (*pweights);
 400     ++psupport;
 401     ++pweights;
 402   }
 403
 404   // Return interpolated value
 405   return(interpolated);
 406 }
 407 //====================================================================
 408
 409 #endif