itk/itkBSplineWeightsCalculator.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 /* =========================================================================
  19
  20   @file   itkBSplineWeightsCalculator.h
  21   @author David Sarrut <david.sarrut@creatis.insa-lyon.fr>
  22
  23   Copyright (c)
  24   * CREATIS (Centre de Recherche et d'Applications en Traitement de l'Image).
  25     All rights reserved. See Doc/License.txt or
  26     http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
  27   * Léon Bérard cancer center, 28 rue Laënnec, 69373 Lyon cedex 08, France
  28   * http://www.creatis.insa-lyon.fr/rio
  29
  30   This software is distributed WITHOUT ANY WARRANTY; without even the
  31   implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  32   PURPOSE.  See the above copyright notices for more information.
  33
  34 ========================================================================= */
  35
  36 #ifndef ITKBSPLINEWEIGHTSCALCULATOR_TXX
  37 #define ITKBSPLINEWEIGHTSCALCULATOR_TXX
  38
  39 //====================================================================
  40 template<int VDimension>
  41 typename Index<VDimension>::IndexValueType Index2Offset(const Index<VDimension> & index, const Index<VDimension> & offsetTable)
  42 {
  43   long v = index[0];
  44   for(int l=1; l<VDimension; l++) {
  45     v += index[l] * offsetTable[l];
  46   }
  47   return v;
  48 }
  49 //====================================================================
  50
  51 //====================================================================
  52 template<class TCoefficientType, int VDimension>
  53 BSplineWeightsCalculator<TCoefficientType, VDimension>::
  54 BSplineWeightsCalculator()
  55 {
  56   mWeightsAreUpToDate = false;
  57 }
  58 //====================================================================
  59
  60 //====================================================================
  61 template<class TCoefficientType, int VDimension>
  62 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
  63 SetSplineOrder(int splineOrder)
  64 {
  65   SizeType temp;
  66   for(int l=0; l<VDimension; l++) temp[l] = splineOrder;
  67   SetSplineOrders(temp);
  68 }
  69 //====================================================================
  70
  71 //====================================================================
  72 template<class TCoefficientType, int VDimension>
  73 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
  74 SetSplineOrders(const SizeType & splineOrder)
  75 {
  76   // Compute support size
  77   mSupportSize = 1;
  78   for(int l=0; l<VDimension; l++) {
  79     mSplineOrders[l] = splineOrder[l];
  80     mSplineSupport[l] = mSplineOrders[l]+1;
  81     mSupportSize *= mSplineSupport[l];
  82   }
  83
  84   // Compute indexes of points in support
  85   mSupportIndex.resize(mSupportSize);
  86   for(int l=0; l<VDimension; l++) mSupportIndex[0][l] = 0;
  87   for(int k=0; k<mSupportSize; k++) {
  88     if (k != mSupportSize-1) {
  89       for(int l=0; l<VDimension; l++) mSupportIndex[k+1][l] = mSupportIndex[k][l];
  90       int l=0;
  91       bool stop = false;
  92       while (!stop) {
  93         mSupportIndex[k+1][l]++;
  94         if (mSupportIndex[k+1][l] == (int)mSplineSupport[l]) { //Ds supportindex=int support=uint
  95           mSupportIndex[k+1][l] = 0;
  96           l++;
  97         } else stop = true;
  98       }
  99     }
 100   }
 101   mWeightsAreUpToDate = false;
 102 }
 103 //====================================================================
 104
 105 //====================================================================
 106 template<class TCoefficientType, int VDimension>
 107 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
 108 SetSamplingFactor(int sampling)
 109 {
 110   for(int l=0; l<VDimension; l++) mSamplingFactors[l] = sampling;
 111   mWeightsAreUpToDate = false;
 112 }
 113 //====================================================================
 114
 115 ///====================================================================
 116 template<class TCoefficientType, int VDimension>
 117 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
 118 SetSamplingFactors(const SizeType & sampling)
 119 {
 120   for(int l=0; l<VDimension; l++) mSamplingFactors[l] = sampling[l];
 121   mWeightsAreUpToDate = false;
 122 }
 123 //====================================================================
 124
 125 //====================================================================
 126 template<class TCoefficientType, int VDimension>
 127 typename BSplineWeightsCalculator<TCoefficientType, VDimension>::InitialWeightsType &
 128 BSplineWeightsCalculator<TCoefficientType, VDimension>::
 129 GetInitialWeights(int order)
 130 {
 131   if (!mBasisFunctionCoefficientsMatrixAreUpToDate[order]) ComputeBasisFunctionCoefficientsMatrix(order);
 132   return mBasisFunctionCoefficientsMatrix[order];
 133 }
 134 //====================================================================
 135
 136 //====================================================================
 137 template<class T> inline T factorial(T rhs)
 138 {
 139   T lhs = (T)1;
 140   for(T x=(T)1; x<=rhs; ++x) {
 141     lhs *= x;
 142   }
 143   return lhs;
 144 }
 145 //====================================================================
 146
 147 //====================================================================
 148 template<class TCoefficientType, int VDimension>
 149 double BSplineWeightsCalculator<TCoefficientType, VDimension>::
 150 BinomialCoefficient(int i, int j)
 151 {
 152   double f = (factorial(i))/(factorial(j) * factorial(i-j));
 153   return f;
 154 }
 155 //====================================================================
 156
 157 //====================================================================
 158 template<class TCoefficientType, int VDimension>
 159 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
 160 ComputeBasisFunctionCoefficientsMatrix(int order)
 161 {
 162   // Compute the sxs matrix of coefficients used to build the different
 163   // polynomials. With s the support is order+1.
 164   int s = order+1;
 165   mBasisFunctionCoefficientsMatrix[order] = InitialWeightsType();
 166   mBasisFunctionCoefficientsMatrix[order].resize(s);
 167   for(int i=0; i<s; i++) {
 168     mBasisFunctionCoefficientsMatrix[order][i].resize(s);
 169     for(int j=0; j<s; j++) {
 170       mBasisFunctionCoefficientsMatrix[order][i][j] = 0.0;
 171       for(int m=j; m<s; m++) {
 172         double a = pow((double)(s-(m+1)),i) * pow((double)-1,m-j) * BinomialCoefficient(s, m-j);
 173         mBasisFunctionCoefficientsMatrix[order][i][j] += a;
 174       }
 175       mBasisFunctionCoefficientsMatrix[order][i][j] *= BinomialCoefficient(s-1, i);
 176     }
 177   }
 178   int f = factorial(order);
 179   for(int i=0; i<s; i++) {
 180     for(int j=0; j<s; j++) {
 181       mBasisFunctionCoefficientsMatrix[order][i][j] /= f;
 182     }
 183   }
 184   mBasisFunctionCoefficientsMatrixAreUpToDate[order] = true;
 185 }
 186 //====================================================================
 187
 188 //====================================================================
 189 template<class TCoefficientType, int VDimension>
 190 TCoefficientType BSplineWeightsCalculator<TCoefficientType, VDimension>::
 191 BSplineEvaluate(int order, int k, double e)
 192 {
 193   // Evaluate a BSpline
 194   int s=order+1;
 195   TCoefficientType v = 0.0;
 196   for(int p=0; p<s; p++) {
 197     v += pow(e,order-p)*GetInitialWeights(order)[p][k];
 198   }
 199   return v;
 200 }
 201 //====================================================================
 202
 203 //====================================================================
 204 template<class TCoefficientType, int VDimension>
 205 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
 206 ComputeSampledWeights1D(std::vector<std::vector<TCoefficientType> > & w, int order, int sampling)
 207 {
 208   int s = order+1;
 209   w.resize(s);
 210   for(int k=0; k<s; k++) w[k].resize(sampling);
 211   double offset = 1.0/sampling;
 212   double e=0.0;
 213   // For odd order (like cubic) start with e=0
 214   if (order & 1) e = 0;
 215   // or even order (like quadratic), shift by 0.5
 216   else  e = +0.5;
 217   for(int a=0; a<sampling; a++) {   // loop over positions
 218     // std::cout << a << " = " << e << std::endl;
 219     for(int k=0; k<s; k++) {        // loop over spline support
 220       w[k][a] = BSplineEvaluate(order, k, e);
 221     }
 222     e += offset;
 223     if (fabs(1.0-e)<=1e-6) e = 1.0; // (for even order)
 224     if (e>=1.0) e = e-1.0;
 225   }
 226 }
 227 //====================================================================
 228
 229 //====================================================================
 230 template<class TCoefficientType, int VDimension>
 231 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
 232 ComputeSampledWeights()
 233 {
 234   mWeights.resize(VDimension);
 235   // Loop over dimension to compute weights
 236   for(int l=0; l<VDimension; l++) {
 237     ComputeSampledWeights1D(mWeights[l], mSplineOrders[l], mSamplingFactors[l]);
 238   }
 239   mWeightsAreUpToDate = true;
 240 }
 241 //====================================================================
 242
 243 //====================================================================
 244 template<class TCoefficientType, int VDimension>
 245 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
 246 ComputeTensorProducts()
 247 {
 248   // Initial BSpline samples weights
 249   ComputeSampledWeights();
 250
 251   // tensor product memory offsets
 252   mTensorProductMemoryOffset[0] = 1;
 253   for(int l=1; l<VDimension; l++) {
 254     mTensorProductMemoryOffset[l] = mTensorProductMemoryOffset[l-1]*mSamplingFactors[l-1];
 255   }
 256
 257   // Tensor product initialisation
 258   int nbPositions = 1;
 259   for(int l=0; l<VDimension; l++) nbPositions *= mSamplingFactors[l];
 260   // int nb = mSupportSize * nbPositions;
 261   mTensorProducts.resize(nbPositions);
 262   IndexType mPositionIndex;
 263   for(int l=0; l<VDimension; l++) mPositionIndex[l] = 0;
 264
 265   // Tensor product
 266   for(int a=0; a<nbPositions; a++) { // Loop over sampling positions
 267     mTensorProducts[a].resize(mSupportSize);
 268
 269     for(int k=0; k<mSupportSize; k++) { // Loop over support positions
 270       TCoefficientType B = 1.0;
 271
 272       for(int l=0; l<VDimension; l++) { // loop for tensor product
 273         B *= mWeights[l][mSupportIndex[k][l]][mPositionIndex[l]];
 274       }
 275       mTensorProducts[a][k] = B;
 276     }
 277
 278     // Next sample Position index
 279     int l=0;
 280     bool stop = false;
 281     while (!stop) {
 282       mPositionIndex[l]++;
 283       if (mPositionIndex[l] == (int)mSamplingFactors[l]) {
 284         mPositionIndex[l] = 0;
 285         l++;
 286       } else stop = true;
 287     }
 288   }
 289 }
 290 //====================================================================
 291
 292 //====================================================================
 293 template<class TCoefficientType, int VDimension>
 294 const TCoefficientType * BSplineWeightsCalculator<TCoefficientType, VDimension>::
 295 GetFirstTensorProduct(const IndexType & index) const
 296 {
 297   int i = Index2Offset<VDimension>(index, mTensorProductMemoryOffset);
 298   return &(mTensorProducts[i][0]);
 299 }
 300 //====================================================================
 301
 302 #endif /* end #define ITKBSPLINEWEIGHTSCALCULATOR_TXX */
 303