itk/itkBSplineInterpolateImageFunctionWithLUT.txx

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