itk/clitkVectorBSplineInterpolateImageFunctionWithLUT.txx

   1 #ifndef _CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX
   2 #define _CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX
   3
   4 /* =========================================================================
   5
   6 @file   itkBSplineInterpolateImageFunctionWithLUT.txx
   7 @author jefvdmb@gmail.com
   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 clitk
  22 {
  23   //====================================================================
  24   template <class TImageType, class TCoordRep, class TCoefficientType>
  25   VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
  26   VectorBSplineInterpolateImageFunctionWithLUT():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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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
 144     if (!inputData) return;
 145     UpdateWeightsProperties();
 146
 147   }
 148
 149   //====================================================================
 150   template <class TImageType, class TCoordRep, class TCoefficientType>
 151   void VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
 152   UpdateWeightsProperties() {
 153
 154     // Compute Memory offset inside coefficients images (for looping over coefficients)
 155     static const unsigned int d = TImageType::ImageDimension;
 156     mInputMemoryOffset[0] = 1;
 157     for(unsigned int l=1; l<d; l++) {
 158       mInputMemoryOffset[l] =
 159         mInputMemoryOffset[l-1]*this->m_Coefficients->GetLargestPossibleRegion().GetSize(l-1);
 160     }
 161
 162     //JV Put here?
 163     if (!mWeightsAreUpToDate)
 164       {
 165         // Compute mSupportOffset according to input size
 166         mSupportOffset.resize(mSupportSize);
 167         mSupportIndex.resize(mSupportSize);
 168         for(unsigned int l=0; l<d; l++) mSupportIndex[0][l] = 0;
 169         for(unsigned int k=0; k<mSupportSize; k++) {
 170           // Get memory offset
 171           mSupportOffset[k] = itk::Index2Offset<TImageType::ImageDimension>(mSupportIndex[k], mInputMemoryOffset);
 172           // next coefficient index
 173           if (k != mSupportSize-1) {
 174             for(unsigned int l=0; l<d; l++) mSupportIndex[k+1][l] = mSupportIndex[k][l];
 175             int l=0;
 176             bool stop = false;
 177             while (!stop) {
 178               mSupportIndex[k+1][l]++;
 179               if (static_cast<unsigned int>(mSupportIndex[k+1][l]) == mSupport[l]) {
 180                 mSupportIndex[k+1][l] = 0;
 181                 l++;
 182               }
 183               else stop = true;
 184             }
 185           }
 186         }
 187
 188         //  // Check
 189         //   for(unsigned int l=0; l<d; l++) {
 190         //     if (mOutputSpacing[l] == -1) {
 191         //       std::cerr << "Please use SetOutputSpacing before using BLUT interpolator" << std::endl;
 192         //       exit(0);
 193         //     }
 194         //   }
 195
 196         // Compute BSpline weights if not up to date
 197         //if (!mWeightsAreUpToDate)
 198         UpdatePrecomputedWeights();
 199       }
 200     // Initialisation
 201     //   *mIntrinsecErrorMax = *mIntrinsecError = 0.0;
 202     //   *mNumberOfError = 0;
 203   }
 204   //====================================================================
 205
 206   //====================================================================
 207   template <class TImageType, class TCoordRep, class TCoefficientType>
 208   void VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
 209   UpdatePrecomputedWeights() {
 210     //   mLUTTimer.Reset();
 211     //   mLUTTimer.Start();
 212     mWeightsCalculator.SetSplineOrders(mSplineOrders);
 213     mWeightsCalculator.SetSamplingFactors(mSamplingFactors);
 214     mWeightsCalculator.ComputeTensorProducts();
 215     mWeightsAreUpToDate = true;
 216     //DS
 217     //   coef = new TCoefficientType[mSupportSize];
 218     //     mCorrectedSupportIndex.resize(mSupportSize);
 219     //     mCorrectedSupportOffset.resize(mSupportSize);
 220     //  mLUTTimer.Stop();
 221     //   mLUTTimer.Print("LUT      \t");
 222   }
 223   //====================================================================
 224
 225   //====================================================================
 226   template <class TImageType, class TCoordRep, class TCoefficientType>
 227   typename VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::IndexType
 228   VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
 229   GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const {
 230
 231     /*
 232       WARNING : sometimes, a floating number x could not really be
 233       represented in memory. In this case, the difference between x and
 234       floor(x) can be almost 1 (instead of 0).  So we take into account
 235       such special case, otherwise it could lead to annoying results.
 236     */
 237     //  static const TCoefficientType tiny = 1.0e-7;
 238     IndexType index;
 239
 240     for(int l=0; l<TImageType::ImageDimension; l++)
 241       {
 242         //  bool mChange = false;
 243
 244         // Compute t1 = distance to floor
 245         TCoefficientType t1 = x[l]- vcl_floor(x[l]);
 246
 247         // Compute index in precomputed weights table
 248         TCoefficientType t2 = mSamplingFactors[l]*t1;
 249         index[l] = (IndexValueType)lrint(t2);
 250
 251         // For even order : test if too close to 0.5 (but lower). In this
 252         // case : take the next coefficient
 253         if (!(mSplineOrders[l] & 1)) {
 254           if (t1<0.5) {
 255             if (mSamplingFactors[l] & 1) {
 256               if (index[l] ==  (int) mSamplingFactors[l]/2+1) EvaluateIndex[l] = EvaluateIndex[l]+1;
 257             }
 258
 259             else if (index[l] == (int) mSamplingFactors[l]/2) EvaluateIndex[l] = EvaluateIndex[l]+1;
 260           }
 261         }
 262
 263         // Statistics (to be removed)
 264         /*
 265          *mIntrinsecError += fabs(index[l]-t2);
 266          (*mNumberOfError)++;
 267          if (fabs(index[l]-t2)> *mIntrinsecErrorMax) *mIntrinsecErrorMax = fabs(index[l]-t2);
 268         */
 269
 270         // When to close to 1, take the next coefficient for odd order, but
 271         // only change index for odd
 272         if (index[l] == (int)mSamplingFactors[l]) {
 273           index[l] = 0;
 274           if (mSplineOrders[l] & 1) EvaluateIndex[l] = EvaluateIndex[l]+1;
 275         }
 276       }
 277
 278     // The end
 279     return index;
 280   }
 281
 282
 283   //====================================================================
 284
 285   //====================================================================
 286   template <class TImageType, class TCoordRep, class TCoefficientType>
 287   typename VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::OutputType
 288   VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
 289   EvaluateAtContinuousIndex(const ContinuousIndexType & x) const {
 290
 291     // For shorter coding
 292     static const unsigned int d = TImageType::ImageDimension;
 293
 294     // Compute the index of the first interpolation coefficient in the coefficient image
 295     IndexType evaluateIndex;
 296     long indx;
 297     for (unsigned int l=0; l<d; l++)  {
 298       if (mSplineOrders[l] & 1) {  // Use this index calculation for odd splineOrder (like cubic)
 299         indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
 300         evaluateIndex[l] = indx;
 301       }
 302       else { // Use this index calculation for even splineOrder
 303         if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]);
 304         else {
 305           indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
 306           evaluateIndex[l] = indx;
 307         }
 308       }
 309     }
 310
 311     // Compute index of precomputed weights and get pointer to first weights
 312     const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex);
 313     const TCoefficientType * pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
 314
 315     // Check boundaries
 316     bool boundaryCase = false;
 317     for (unsigned int l=0; l<d; l++) {
 318       if ((evaluateIndex[l] < 0) ||
 319           (evaluateIndex[l]+mSupport[l]) >= this->m_Coefficients->GetLargestPossibleRegion().GetSize(l)) {
 320         boundaryCase = true;
 321       }
 322     }
 323
 324     // Pointer to support offset
 325     const int * psupport;
 326
 327     // Special case for boundary (to be changed ....)
 328     std::vector<int> correctedSupportOffset;
 329     if (boundaryCase) {
 330       //    return -1000;
 331       //std::vector<TCoefficientType> coef(mSupportSize);
 332       // DD(EvaluateIndex);
 333       //std::vector<int> CorrectedSupportOffset;//(mSupportSize);
 334       std::vector<IndexType> correctedSupportIndex;//(mSupportSize);
 335       correctedSupportIndex.resize(mSupportSize);
 336       correctedSupportOffset.resize(mSupportSize);
 337       for(unsigned int i=0; i<mSupportSize; i++) {
 338         // DD(mSupportIndex[i]);
 339         for (unsigned int l=0; l<d; l++) {
 340           long a = mSupportIndex[i][l] + evaluateIndex[l];
 341           long b = this->m_Coefficients->GetLargestPossibleRegion().GetSize(l);
 342           // DD(a);
 343           // DD(b);
 344           long d2 = 2 * b - 2;
 345           if (a < 0) {
 346             correctedSupportIndex[i][l] = -a - d2*(-a/d2) - evaluateIndex[l];//mSupportIndex[i][l]-a;
 347           }
 348           else {
 349             if (a>=b) {
 350               correctedSupportIndex[i][l] = d2 - a - evaluateIndex[l];
 351             }
 352             else {
 353               correctedSupportIndex[i][l] = mSupportIndex[i][l]; //a - d2*(a/d2) - EvaluateIndex[l];
 354             }
 355             /*
 356               if (a>=b) {
 357               correctedSupportIndex[i][l] = d2 - a - EvaluateIndex[l];//mSupportIndex[i][l] - (a-(b-1));
 358               }
 359               else {
 360               correctedSupportIndex[i][l] = mSupportIndex[i][l];
 361               }
 362             */
 363           }
 364         }
 365         // DD(correctedSupportIndex[i]);
 366         correctedSupportOffset[i] = itk::Index2Offset<TImageType::ImageDimension>(correctedSupportIndex[i], mInputMemoryOffset);
 367       }
 368       // for (unsigned int l=0; l<d; l++) {
 369       //       EvaluateIndex[l] = EvaluateIndex[l] + correctedSupportIndex[0][l];
 370       //     }
 371       // DD(EvaluateIndex);
 372       psupport = &correctedSupportOffset[0];
 373     }
 374     else {
 375       psupport = &mSupportOffset[0];
 376     }
 377
 378     // Get pointer to first coefficient. Even if in some boundary cases,
 379     // EvaluateIndex is out of the coefficient image,
 380     //JV
 381     const CoefficientImagePixelType * pcoef = &(this->m_Coefficients->GetPixel(evaluateIndex));
 382
 383     // Main loop over BSpline support
 384     CoefficientImagePixelType interpolated = 0.0;
 385     for (unsigned int p=0; p<mSupportSize; p++) {
 386       interpolated += pcoef[*psupport] * (*pweights);
 387       ++psupport;
 388       ++pweights;
 389     }
 390
 391     // Return interpolated value
 392     return(interpolated);
 393   }
 394   //====================================================================
 395
 396
 397   //====================================================================
 398   template <class TImageType, class TCoordRep, class TCoefficientType>
 399   void
 400   VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
 401   EvaluateWeightsAtContinuousIndex(const ContinuousIndexType&  x,  const TCoefficientType ** pweights, IndexType & evaluateIndex)const {
 402
 403     // JV Compute BSpline weights if not up to date! Problem const: pass image as last
 404     //  if (!mWeightsAreUpToDate) UpdatePrecomputedWeights();
 405
 406     // For shorter coding
 407     static const unsigned int d = TImageType::ImageDimension;
 408
 409     // Compute the index of the first interpolation coefficient in the coefficient image
 410     //IndexType evaluateIndex;
 411     long indx;
 412     for (unsigned int l=0; l<d; l++)  {
 413       if (mSplineOrders[l] & 1) {  // Use this index calculation for odd splineOrder (like cubic)
 414         indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
 415         evaluateIndex[l] = indx;
 416       }
 417       else { // Use this index calculation for even splineOrder
 418         if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]);
 419         else {
 420           indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
 421           evaluateIndex[l] = indx;
 422         }
 423       }
 424     }
 425
 426     // Compute index of precomputed weights and get pointer to first weights
 427     const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex);
 428     *pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
 429
 430   }
 431   //====================================================================
 432
 433
 434
 435
 436 } //end namespace
 437
 438 #endif //_CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX