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