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