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