itk/clitkVectorBSplineInterpolateImageFunctionWithLUT.txx

   1 #ifndef _CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX
   2 #define _CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX
   3 /* =========================================================================
   4
   5 @file   itkBSplineInterpolateImageFunctionWithLUT.txx
   6 @author jefvdmb@gmail.com
   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 clitk
  21 {
  22   //====================================================================
  23   template <class TImageType, class TCoordRep, class TCoefficientType>
  24   VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
  25   VectorBSplineInterpolateImageFunctionWithLUT():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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<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
 143     if (!inputData) return;
 144     UpdateWeightsProperties();
 145
 146   }
 147
 148   //====================================================================
 149   template <class TImageType, class TCoordRep, class TCoefficientType>
 150   void VectorBSplineInterpolateImageFunctionWithLUT<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] = itk::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 VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::IndexType
 227   VectorBSplineInterpolateImageFunctionWithLUT<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 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::OutputType
 287   VectorBSplineInterpolateImageFunctionWithLUT<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     //JV
 380     const CoefficientImagePixelType * pcoef = &(this->m_Coefficients->GetPixel(evaluateIndex));
 381
 382     // Main loop over BSpline support
 383     CoefficientImagePixelType interpolated = 0.0;
 384     for (unsigned int p=0; p<mSupportSize; p++) {
 385       interpolated += pcoef[*psupport] * (*pweights);
 386       ++psupport;
 387       ++pweights;
 388     }
 389
 390     // Return interpolated value
 391     return(interpolated);
 392   }
 393   //====================================================================
 394
 395
 396   //====================================================================
 397   template <class TImageType, class TCoordRep, class TCoefficientType>
 398   void
 399   VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
 400   EvaluateWeightsAtContinuousIndex(const ContinuousIndexType&  x,  const TCoefficientType ** pweights, IndexType & evaluateIndex)const {
 401
 402     // JV Compute BSpline weights if not up to date! Problem const: pass image as last
 403     //  if (!mWeightsAreUpToDate) UpdatePrecomputedWeights();
 404
 405     // For shorter coding
 406     static const unsigned int d = TImageType::ImageDimension;
 407
 408     // Compute the index of the first interpolation coefficient in the coefficient image
 409     //IndexType evaluateIndex;
 410     long indx;
 411     for (unsigned int l=0; l<d; l++)  {
 412       if (mSplineOrders[l] & 1) {  // Use this index calculation for odd splineOrder (like cubic)
 413         indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
 414         evaluateIndex[l] = indx;
 415       }
 416       else { // Use this index calculation for even splineOrder
 417         if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]);
 418         else {
 419           indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
 420           evaluateIndex[l] = indx;
 421         }
 422       }
 423     }
 424
 425     // Compute index of precomputed weights and get pointer to first weights
 426     const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex);
 427     *pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
 428
 429   }
 430   //====================================================================
 431
 432
 433
 434
 435 } //end namespace
 436
 437 #endif //_CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX