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)lrint(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
 249     // JV Compute BSpline weights if not up to date! Problem const: pass image as last
 250     //  if (!mWeightsAreUpToDate) UpdatePrecomputedWeights();
 251   // For shorter coding
 252   static const unsigned int d = TImageType::ImageDimension;
 253
 254   // Compute the index of the first interpolation coefficient in the coefficient image
 255   IndexType evaluateIndex;
 256   long indx;
 257   for (unsigned int l=0; l<d; l++)  {
 258     if (mSplineOrders[l] & 1) {  // Use this index calculation for odd splineOrder (like cubic)
 259       indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
 260       evaluateIndex[l] = indx;
 261     } else { // Use this index calculation for even splineOrder
 262       if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]);
 263       else {
 264         indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
 265         evaluateIndex[l] = indx;
 266       }
 267     }
 268   }
 269
 270   // Compute index of precomputed weights and get pointer to first weights
 271   const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex);
 272   const TCoefficientType * pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
 273
 274   // Check boundaries
 275   bool boundaryCase = false;
 276   for (unsigned int l=0; l<d; l++) {
 277     if ((evaluateIndex[l] < 0) ||
 278         (evaluateIndex[l]+mSupport[l]) >= this->m_Coefficients->GetLargestPossibleRegion().GetSize(l)) {
 279       boundaryCase = true;
 280     }
 281   }
 282
 283   // Pointer to support offset
 284   const int * psupport;
 285
 286   // Special case for boundary (to be changed ....)
 287   std::vector<int> correctedSupportOffset;
 288   if (boundaryCase) {
 289     std::vector<IndexType> correctedSupportIndex;//(mSupportSize);
 290     correctedSupportIndex.resize(mSupportSize);
 291     correctedSupportOffset.resize(mSupportSize);
 292     for(unsigned int i=0; i<mSupportSize; i++) {
 293       for (unsigned int l=0; l<d; l++) {
 294         long a = mSupportIndex[i][l] + evaluateIndex[l];
 295         long b = this->m_Coefficients->GetLargestPossibleRegion().GetSize(l);
 296         long d2 = 2 * b - 2;
 297         if (a < 0) {
 298           correctedSupportIndex[i][l] = -a - d2*(-a/d2) - evaluateIndex[l];//mSupportIndex[i][l]-a;
 299         } else {
 300           if (a>=b) {
 301             correctedSupportIndex[i][l] = d2 - a - evaluateIndex[l];
 302           } else {
 303             correctedSupportIndex[i][l] = mSupportIndex[i][l]; //a - d2*(a/d2) - EvaluateIndex[l];
 304           }
 305           /*
 306             if (a>=b) {
 307             correctedSupportIndex[i][l] = d2 - a - EvaluateIndex[l];//mSupportIndex[i][l] - (a-(b-1));
 308             }
 309             else {
 310             correctedSupportIndex[i][l] = mSupportIndex[i][l];
 311             }
 312           */
 313         }
 314       }
 315       // DD(correctedSupportIndex[i]);
 316       correctedSupportOffset[i] = itk::Index2Offset<TImageType::ImageDimension>(correctedSupportIndex[i], mInputMemoryOffset);
 317     }
 318     psupport = &correctedSupportOffset[0];
 319   } else {
 320     psupport = &mSupportOffset[0];
 321   }
 322
 323   // Get pointer to first coefficient. Even if in some boundary cases,
 324   // EvaluateIndex is out of the coefficient image,
 325   //JV
 326   const CoefficientImagePixelType * pcoef = &(this->m_Coefficients->GetPixel(evaluateIndex));
 327
 328   // Main loop over BSpline support
 329   CoefficientImagePixelType interpolated = itk::NumericTraits<CoefficientImagePixelType>::Zero;;
 330   for (unsigned int p=0; p<mSupportSize; p++) {
 331     interpolated += pcoef[*psupport] * (*pweights);
 332     ++psupport;
 333     ++pweights;
 334   }
 335
 336   // Return interpolated value
 337   return(interpolated);
 338 }
 339 //====================================================================
 340
 341
 342 //====================================================================
 343 template <class TImageType, class TCoordRep, class TCoefficientType>
 344 void
 345 VectorBSplineInterpolateImageFunctionWithLUT<TImageType,TCoordRep,TCoefficientType>::
 346 EvaluateWeightsAtContinuousIndex(const ContinuousIndexType&  x,  const TCoefficientType ** pweights, IndexType & evaluateIndex)const
 347 {
 348
 349   // JV Compute BSpline weights if not up to date! Problem const: pass image as last
 350   //  if (!mWeightsAreUpToDate) UpdatePrecomputedWeights();
 351
 352   // For shorter coding
 353   static const unsigned int d = TImageType::ImageDimension;
 354
 355   // Compute the index of the first interpolation coefficient in the coefficient image
 356   long indx;
 357   for (unsigned int l=0; l<d; l++)  {
 358     if (mSplineOrders[l] & 1) {  // Use this index calculation for odd splineOrder (like cubic)
 359       indx = (long)vcl_floor(x[l]) - mSplineOrders[l] / 2 ; //this->m_SplineOrder / 2;
 360       evaluateIndex[l] = indx;
 361     } else { // Use this index calculation for even splineOrder
 362       if (mSplineOrders[l] == 0) evaluateIndex[l] = (long)rint(x[l]);
 363       else {
 364         indx = (long)vcl_floor((x[l]+ 0.5)) - mSplineOrders[l] / 2; //this->m_SplineOrder / 2;
 365         evaluateIndex[l] = indx;
 366       }
 367     }
 368   }
 369
 370   // Compute index of precomputed weights and get pointer to first weights
 371   const IndexType weightIndex = GetSampleIndexOfPixelPosition(x, evaluateIndex);
 372   *pweights = mWeightsCalculator.GetFirstTensorProduct(weightIndex);
 373
 374 }
 375 //====================================================================
 376
 377
 378
 379
 380 } //end namespace
 381
 382 #endif //_CLITKVECTORBSPLINEINTERPOLATEIMAGEFUNCTIONWITHLUT_TXX