itk/itkBSplineInterpolateImageFunctionWithLUT.txx

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