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