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://www.centreleonberard.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 #include <itkMath.h>
  21
  22 namespace itk
  23 {
  24   //====================================================================
  25   template <class TImageType, class TCoordRep, class TCoefficientType>
  26   BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
  27   BSplineInterpolateImageFunctionWithLUT():Superclass() {
  28
  29     // Set default values
  30     SetLUTSamplingFactor(20);
  31     SetSplineOrder(3);
  32     mWeightsAreUpToDate = false;
  33   }
  34
  35   //====================================================================
  36
  37   //====================================================================
  38   template <class TImageType, class TCoordRep, class TCoefficientType>
  39   void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
  40   SetLUTSamplingFactor(const int& s) {
  41     for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s;
  42     mWeightsAreUpToDate = false;
  43   }
  44
  45   //====================================================================
  46
  47   //====================================================================
  48   template <class TImageType, class TCoordRep, class TCoefficientType>
  49   void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
  50   SetLUTSamplingFactors(const SizeType& s) {
  51     for(int i=0; i<TImageType::ImageDimension; i++) mSamplingFactors[i] = s[i];
  52     mWeightsAreUpToDate = false;
  53   }
  54
  55   //====================================================================
  56
  57   //====================================================================
  58   template <class TImageType, class TCoordRep, class TCoefficientType>
  59   void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
  60   SetSplineOrder(const unsigned int& SplineOrder) {
  61     Superclass::SetSplineOrder(SplineOrder);
  62     // Compute support and half size
  63     static const int d = TImageType::ImageDimension;
  64     for(int l=0; l<d; l++) {
  65       mSplineOrders[l]= SplineOrder;
  66       mSupport[l] = SplineOrder+1;
  67       if (mSupport[l] % 2 == 0) { // support is even
  68         mHalfSupport[l] = mSupport[l]/2-1;
  69       }
  70       else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
  71     }
  72     mSupportSize = 1;
  73     for(int l=0; l<d; l++) {
  74       mSupportSize *= mSupport[l];
  75     }
  76     mWeightsAreUpToDate = false;
  77   }
  78   //====================================================================
  79
  80   //JV
  81   //====================================================================
  82   template <class TImageType, class TCoordRep, class TCoefficientType>
  83   void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
  84   SetSplineOrders(const SizeType& SplineOrders) {
  85     mSplineOrders=SplineOrders;
  86
  87     // Compute support and half size
  88     static const int d = TImageType::ImageDimension;
  89     for(int l=0; l<d; l++) {
  90       mSupport[l] = mSplineOrders[l]+1;
  91       if (mSupport[l] % 2 == 0) { // support is even
  92         mHalfSupport[l] = mSupport[l]/2-1;
  93       }
  94       else mHalfSupport[l] = mSupport[l]/2; // support is odd (like cubic spline)
  95     }
  96     mSupportSize = 1;
  97     for(int l=0; l<d; l++) {
  98       mSupportSize *= mSupport[l];
  99     }
 100     mWeightsAreUpToDate = false;
 101   }
 102   //====================================================================
 103
 104   //====================================================================
 105   template <class TImageType, class TCoordRep, class TCoefficientType>
 106   void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
 107   SetInputImage(const TImageType * inputData) {
 108
 109     // JV  Call superclass (decomposition filter is executeed each time!)
 110     // JV Should call itkBSplineDecompositionFilterWithOBD to allow different order by dimension
 111     Superclass::SetInputImage(inputData);
 112
 113     // Update the weightproperties
 114     if (!inputData) return;
 115     UpdateWeightsProperties();
 116   }
 117
 118   //====================================================================
 119   template <class TImageType, class TCoordRep, class TCoefficientType>
 120   void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
 121   UpdateWeightsProperties() {
 122
 123     // Compute Memory offset inside coefficients images (for looping over coefficients)
 124     static const unsigned int d = TImageType::ImageDimension;
 125     mInputMemoryOffset[0] = 1;
 126     for(unsigned int l=1; l<d; l++) {
 127       mInputMemoryOffset[l] =
 128         mInputMemoryOffset[l-1]*this->m_Coefficients->GetLargestPossibleRegion().GetSize(l-1);
 129     }
 130
 131     // Compute mSupportOffset according to input size
 132     mSupportOffset.resize(mSupportSize);
 133     mSupportIndex.resize(mSupportSize);
 134     for(unsigned int l=0; l<d; l++) mSupportIndex[0][l] = 0;
 135     for(unsigned int k=0; k<mSupportSize; k++) {
 136       // Get memory offset
 137       mSupportOffset[k] = Index2Offset<TImageType::ImageDimension>(mSupportIndex[k], mInputMemoryOffset);
 138       // next coefficient index
 139       if (k != mSupportSize-1) {
 140         for(unsigned int l=0; l<d; l++) mSupportIndex[k+1][l] = mSupportIndex[k][l];
 141         int l=0;
 142         bool stop = false;
 143         while (!stop) {
 144           mSupportIndex[k+1][l]++;
 145           if (static_cast<unsigned int>(mSupportIndex[k+1][l]) == mSupport[l]) {
 146             mSupportIndex[k+1][l] = 0;
 147             l++;
 148           }
 149               else stop = true;
 150         }
 151       }
 152     }
 153
 154     // Compute BSpline weights if not up to date
 155     if (!mWeightsAreUpToDate)
 156       UpdatePrecomputedWeights();
 157   }
 158   //====================================================================
 159
 160   //====================================================================
 161   template <class TImageType, class TCoordRep, class TCoefficientType>
 162   void BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
 163   UpdatePrecomputedWeights() {
 164
 165     mWeightsCalculator.SetSplineOrders(mSplineOrders);
 166     mWeightsCalculator.SetSamplingFactors(mSamplingFactors);
 167     mWeightsCalculator.ComputeTensorProducts();
 168     mWeightsAreUpToDate = true;
 169   }
 170   //====================================================================
 171
 172   //====================================================================
 173   template <class TImageType, class TCoordRep, class TCoefficientType>
 174   typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::IndexType
 175   BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
 176   GetSampleIndexOfPixelPosition(const ContinuousIndexType & x, IndexType & EvaluateIndex) const {
 177
 178     /*
 179       WARNING : sometimes, a floating number x could not really be
 180       represented in memory. In this case, the difference between x and
 181       floor(x) can be almost 1 (instead of 0).  So we take into account
 182       such special case, otherwise it could lead to annoying results.
 183     */
 184     //  static const TCoefficientType tiny = 1.0e-7;
 185     IndexType index;
 186
 187     for(int l=0; l<TImageType::ImageDimension; l++)
 188       {
 189         // Compute t1 = distance to floor
 190         TCoefficientType t1 = x[l]- vcl_floor(x[l]);
 191
 192         // Compute index in precomputed weights table
 193         TCoefficientType t2 = mSamplingFactors[l]*t1;
 194         index[l] = (IndexValueType)Math::Round<TCoefficientType>(t2);
 195
 196         // For even order : test if too close to 0.5 (but lower). In this
 197         // case : take the next coefficient
 198         if (!(mSplineOrders[l] & 1)) {
 199           if (t1<0.5) {
 200             if (mSamplingFactors[l] & 1) {
 201               if (index[l] ==  (int) mSamplingFactors[l]/2+1) EvaluateIndex[l] = EvaluateIndex[l]+1;
 202             }
 203
 204             else if (index[l] == (int) mSamplingFactors[l]/2) EvaluateIndex[l] = EvaluateIndex[l]+1;
 205           }
 206         }
 207
 208         // When to close to 1, take the next coefficient for odd order, but
 209         // only change index for odd
 210         if (index[l] == (int)mSamplingFactors[l]) {
 211           index[l] = 0;
 212           if (mSplineOrders[l] & 1) EvaluateIndex[l] = EvaluateIndex[l]+1;
 213         }
 214       }
 215
 216     // The end
 217     return index;
 218   }
 219
 220
 221   //====================================================================
 222
 223   //====================================================================
 224   template <class TImageType, class TCoordRep, class TCoefficientType>
 225   typename BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::OutputType
 226   BSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
 227   EvaluateAtContinuousIndex(const ContinuousIndexType & x) const {
 228
 229     // JV Compute BSpline weights if not up to date! Problem const: pass image as last
 230     //  if (!mWeightsAreUpToDate) UpdatePrecomputedWeights();
 231
 232     // For shorter coding
 233     static const unsigned int d = TImageType::ImageDimension;
 234
 235     // Compute the index of the first interpolation coefficient in the coefficient image
 236     IndexType evaluateIndex;
 237     long indx;
 238     for (unsigned int l=0; l<d; l++)  {
 239       if (mSplineOrders[l] & 1) {  // Use this index calculation for odd splineOrder (like cubic)
 240         indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
 241         evaluateIndex[l] = indx;
 242       }
 243       else { // Use this index calculation for even splineOrder
 244         if (mSplineOrders[l] == 0) evaluateIndex[l] = Math::Round<typename ContinuousIndexType::ValueType>(x[l]);
 245         else {
 246           indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
 247           evaluateIndex[l] = indx;
 248         }
 249       }
 250     }
 251
 252     // Compute index of precomputed weights and get pointer to first weights
 253     const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex);
 254     const TCoefficientType * pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
 255
 256     // Check boundaries
 257     bool boundaryCase = false;
 258     for (unsigned int l=0; l<d; l++) {
 259       if ((evaluateIndex[l] < 0) ||
 260           (evaluateIndex[l]+mSupport[l]) >= this->m_Coefficients->GetLargestPossibleRegion().GetSize(l)) {
 261         boundaryCase = true;
 262       }
 263     }
 264
 265     // Pointer to support offset
 266     const int * psupport;
 267
 268     // Special case for boundary (to be changed ....)
 269     std::vector<int> correctedSupportOffset;
 270     if (boundaryCase) {
 271
 272       std::vector<IndexType> correctedSupportIndex;//(mSupportSize);
 273       correctedSupportIndex.resize(mSupportSize);
 274       correctedSupportOffset.resize(mSupportSize);
 275       for(unsigned int i=0; i<mSupportSize; i++) {
 276         for (unsigned int l=0; l<d; l++) {
 277           long a = mSupportIndex[i][l] + evaluateIndex[l];
 278           long b = this->m_Coefficients->GetLargestPossibleRegion().GetSize(l);
 279           long d2 = 2 * b - 2;
 280           if (a < 0) {
 281             correctedSupportIndex[i][l] = -a - d2*(-a/d2) - evaluateIndex[l];//mSupportIndex[i][l]-a;
 282           }
 283           else {
 284             if (a>=b) {
 285               correctedSupportIndex[i][l] = d2 - a - evaluateIndex[l];
 286             }
 287             else {
 288               correctedSupportIndex[i][l] = mSupportIndex[i][l]; //a - d2*(a/d2) - EvaluateIndex[l];
 289             }
 290             /*
 291               if (a>=b) {
 292               correctedSupportIndex[i][l] = d2 - a - EvaluateIndex[l];//mSupportIndex[i][l] - (a-(b-1));
 293               }
 294               else {
 295               correctedSupportIndex[i][l] = mSupportIndex[i][l];
 296               }
 297             */
 298           }
 299         }
 300         correctedSupportOffset[i] = itk::Index2Offset<TImageType::ImageDimension>(correctedSupportIndex[i], mInputMemoryOffset);
 301       }
 302       psupport = &correctedSupportOffset[0];
 303     }
 304     else {
 305       psupport = &mSupportOffset[0];
 306     }
 307
 308     // Get pointer to first coefficient. Even if in some boundary cases,
 309     // EvaluateIndex is out of the coefficient image,
 310     const TCoefficientType * pcoef = &(this->m_Coefficients->GetPixel(evaluateIndex));
 311
 312     // Main loop over BSpline support
 313     TCoefficientType interpolated = itk::NumericTraits<TCoefficientType>::Zero;
 314     for (unsigned int p=0; p<mSupportSize; p++) {
 315       interpolated += pcoef[*psupport] * (*pweights);
 316       ++psupport;
 317       ++pweights;
 318     }
 319
 320     // Return interpolated value
 321     return(interpolated);
 322   }
 323   //====================================================================
 324
 325 }
 326 #endif //_ITKINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX