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