]> Creatis software - clitk.git/blob - itk/itkBSplineWeightsCalculator.txx
Put console back for windows
[clitk.git] / 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 {
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