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