itk/itkBSplineWeightsCalculator.txx

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