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://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 /* =========================================================================
  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   long v = index[0];
  43   for(int l=1; l<VDimension; l++) {
  44     v += index[l] * offsetTable[l];
  45   }
  46   return v;
  47 }
  48 //====================================================================
  49
  50 //====================================================================
  51 template<class TCoefficientType, int VDimension>
  52 BSplineWeightsCalculator<TCoefficientType, VDimension>::
  53 BSplineWeightsCalculator() {
  54   mWeightsAreUpToDate = false;
  55 }
  56 //====================================================================
  57
  58 //====================================================================
  59 template<class TCoefficientType, int VDimension>
  60 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
  61 SetSplineOrder(int splineOrder) {
  62   SizeType temp;
  63   for(int l=0; l<VDimension; l++) temp[l] = splineOrder;
  64   SetSplineOrders(temp);
  65 }
  66 //====================================================================
  67
  68 //====================================================================
  69 template<class TCoefficientType, int VDimension>
  70 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
  71 SetSplineOrders(const SizeType & splineOrder) {
  72   // Compute support size
  73   mSupportSize = 1;
  74   for(int l=0; l<VDimension; l++) {
  75     mSplineOrders[l] = splineOrder[l];
  76     mSplineSupport[l] = mSplineOrders[l]+1;
  77     mSupportSize *= mSplineSupport[l];
  78   }
  79
  80   // Compute indexes of points in support
  81   mSupportIndex.resize(mSupportSize);
  82   for(int l=0; l<VDimension; l++) mSupportIndex[0][l] = 0;
  83   for(int k=0; k<mSupportSize; k++) {
  84     if (k != mSupportSize-1) {
  85       for(int l=0; l<VDimension; l++) mSupportIndex[k+1][l] = mSupportIndex[k][l];
  86       int l=0;
  87       bool stop = false;
  88       while (!stop) {
  89         mSupportIndex[k+1][l]++;
  90         if (mSupportIndex[k+1][l] == (int)mSplineSupport[l]) { //Ds supportindex=int support=uint
  91           mSupportIndex[k+1][l] = 0;
  92           l++;
  93         }
  94         else stop = true;
  95       }
  96     }
  97   }
  98   mWeightsAreUpToDate = false;
  99 }
 100 //====================================================================
 101
 102 //====================================================================
 103 template<class TCoefficientType, int VDimension>
 104 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
 105 SetSamplingFactor(int sampling) {
 106   for(int l=0; l<VDimension; l++) mSamplingFactors[l] = sampling;
 107   mWeightsAreUpToDate = false;
 108 }
 109 //====================================================================
 110
 111 ///====================================================================
 112 template<class TCoefficientType, int VDimension>
 113 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
 114 SetSamplingFactors(const SizeType & sampling) {
 115   for(int l=0; l<VDimension; l++) mSamplingFactors[l] = sampling[l];
 116   mWeightsAreUpToDate = false;
 117 }
 118 //====================================================================
 119
 120 //====================================================================
 121 template<class TCoefficientType, int VDimension>
 122 typename BSplineWeightsCalculator<TCoefficientType, VDimension>::InitialWeightsType &
 123 BSplineWeightsCalculator<TCoefficientType, VDimension>::
 124 GetInitialWeights(int order) {
 125   if (!mBasisFunctionCoefficientsMatrixAreUpToDate[order]) ComputeBasisFunctionCoefficientsMatrix(order);
 126   return mBasisFunctionCoefficientsMatrix[order];
 127 }
 128 //====================================================================
 129
 130 //====================================================================
 131 template<class T> inline T factorial(T rhs) {
 132   T lhs = (T)1;
 133   for(T x=(T)1; x<=rhs; ++x) {
 134     lhs *= x;
 135   }
 136   return lhs;
 137 }
 138 //====================================================================
 139
 140 //====================================================================
 141 template<class TCoefficientType, int VDimension>
 142 double BSplineWeightsCalculator<TCoefficientType, VDimension>::
 143 BinomialCoefficient(int i, int j) {
 144   double f = (factorial(i))/(factorial(j) * factorial(i-j));
 145   return f;
 146 }
 147 //====================================================================
 148
 149 //====================================================================
 150 template<class TCoefficientType, int VDimension>
 151 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
 152 ComputeBasisFunctionCoefficientsMatrix(int order) {
 153   // Compute the sxs matrix of coefficients used to build the different
 154   // polynomials. With s the support is order+1.
 155   int s = order+1;
 156   mBasisFunctionCoefficientsMatrix[order] = InitialWeightsType();
 157   mBasisFunctionCoefficientsMatrix[order].resize(s);
 158   for(int i=0; i<s; i++) {
 159     mBasisFunctionCoefficientsMatrix[order][i].resize(s);
 160     for(int j=0; j<s; j++) {
 161       mBasisFunctionCoefficientsMatrix[order][i][j] = 0.0;
 162       for(int m=j; m<s; m++) {
 163         double a = pow((double)(s-(m+1)),i) * pow((double)-1,m-j) * BinomialCoefficient(s, m-j);
 164         mBasisFunctionCoefficientsMatrix[order][i][j] += a;
 165       }
 166       mBasisFunctionCoefficientsMatrix[order][i][j] *= BinomialCoefficient(s-1, i);
 167     }
 168   }
 169   int f = factorial(order);
 170   for(int i=0; i<s; i++) {
 171     for(int j=0; j<s; j++) {
 172       mBasisFunctionCoefficientsMatrix[order][i][j] /= f;
 173     }
 174   }
 175   mBasisFunctionCoefficientsMatrixAreUpToDate[order] = true;
 176 }
 177 //====================================================================
 178
 179 //====================================================================
 180 template<class TCoefficientType, int VDimension>
 181 TCoefficientType BSplineWeightsCalculator<TCoefficientType, VDimension>::
 182 BSplineEvaluate(int order, int k, double e) {
 183   // Evaluate a BSpline
 184   int s=order+1;
 185   TCoefficientType v = 0.0;
 186   for(int p=0; p<s; p++) {
 187     v += pow(e,order-p)*GetInitialWeights(order)[p][k];
 188   }
 189   return v;
 190 }
 191 //====================================================================
 192
 193 //====================================================================
 194 template<class TCoefficientType, int VDimension>
 195 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
 196 ComputeSampledWeights1D(std::vector<std::vector<TCoefficientType> > & w, int order, int sampling) {
 197   int s = order+1;
 198   w.resize(s);
 199   for(int k=0; k<s; k++) w[k].resize(sampling);
 200   double offset = 1.0/sampling;
 201   double e=0.0;
 202   // For odd order (like cubic) start with e=0
 203   if (order & 1) e = 0;
 204   // or even order (like quadratic), shift by 0.5
 205   else  e = +0.5;
 206   for(int a=0; a<sampling; a++) {   // loop over positions
 207     // std::cout << a << " = " << e << std::endl;
 208     for(int k=0; k<s; k++) {        // loop over spline support
 209       w[k][a] = BSplineEvaluate(order, k, e);
 210     }
 211     e += offset;
 212     if (fabs(1.0-e)<=1e-6) e = 1.0; // (for even order)
 213     if (e>=1.0) e = e-1.0;
 214   }
 215 }
 216 //====================================================================
 217
 218 //====================================================================
 219 template<class TCoefficientType, int VDimension>
 220 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
 221 ComputeSampledWeights() {
 222   mWeights.resize(VDimension);
 223   // Loop over dimension to compute weights
 224   for(int l=0; l<VDimension; l++) {
 225     ComputeSampledWeights1D(mWeights[l], mSplineOrders[l], mSamplingFactors[l]);
 226   }
 227   mWeightsAreUpToDate = true;
 228 }
 229 //====================================================================
 230
 231 //====================================================================
 232 template<class TCoefficientType, int VDimension>
 233 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
 234 ComputeTensorProducts() {
 235   // Initial BSpline samples weights
 236   ComputeSampledWeights();
 237
 238   // tensor product memory offsets
 239   mTensorProductMemoryOffset[0] = 1;
 240   for(int l=1; l<VDimension; l++) {
 241     mTensorProductMemoryOffset[l] = mTensorProductMemoryOffset[l-1]*mSamplingFactors[l-1];
 242   }
 243
 244   // Tensor product initialisation
 245   int nbPositions = 1;
 246   for(int l=0; l<VDimension; l++) nbPositions *= mSamplingFactors[l];
 247   // int nb = mSupportSize * nbPositions;
 248   mTensorProducts.resize(nbPositions);
 249   IndexType mPositionIndex;
 250   for(int l=0; l<VDimension; l++) mPositionIndex[l] = 0;
 251
 252   // Tensor product
 253   for(int a=0; a<nbPositions; a++) { // Loop over sampling positions
 254     mTensorProducts[a].resize(mSupportSize);
 255
 256     for(int k=0; k<mSupportSize; k++) { // Loop over support positions
 257       TCoefficientType B = 1.0;
 258
 259       for(int l=0; l<VDimension; l++) { // loop for tensor product
 260         B *= mWeights[l][mSupportIndex[k][l]][mPositionIndex[l]];
 261       }
 262       mTensorProducts[a][k] = B;
 263     }
 264
 265     // Next sample Position index
 266     int l=0;
 267     bool stop = false;
 268     while (!stop) {
 269       mPositionIndex[l]++;
 270       if (mPositionIndex[l] == (int)mSamplingFactors[l]) {
 271         mPositionIndex[l] = 0;
 272         l++;
 273       }
 274       else stop = true;
 275     }
 276   }
 277 }
 278 //====================================================================
 279
 280 //====================================================================
 281 template<class TCoefficientType, int VDimension>
 282 const TCoefficientType * BSplineWeightsCalculator<TCoefficientType, VDimension>::
 283 GetFirstTensorProduct(const IndexType & index) const {
 284   int i = Index2Offset<VDimension>(index, mTensorProductMemoryOffset);
 285   return &(mTensorProducts[i][0]);
 286 }
 287 //====================================================================
 288
 289 #endif /* end #define ITKBSPLINEWEIGHTSCALCULATOR_TXX */
 290