1 /* =========================================================================
3 @file itkBSplineWeightsCalculator.h
4 @author David Sarrut <david.sarrut@creatis.insa-lyon.fr>
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
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.
17 ========================================================================= */
19 #ifndef ITKBSPLINEWEIGHTSCALCULATOR_TXX
20 #define ITKBSPLINEWEIGHTSCALCULATOR_TXX
22 //====================================================================
23 template<int VDimension>
24 typename Index<VDimension>::IndexValueType Index2Offset(const Index<VDimension> & index, const Index<VDimension> & offsetTable) {
26 for(int l=1; l<VDimension; l++) {
27 v += index[l] * offsetTable[l];
31 //====================================================================
33 //====================================================================
34 template<class TCoefficientType, int VDimension>
35 BSplineWeightsCalculator<TCoefficientType, VDimension>::
36 BSplineWeightsCalculator() {
37 mWeightsAreUpToDate = false;
39 //====================================================================
41 //====================================================================
42 template<class TCoefficientType, int VDimension>
43 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
44 SetSplineOrder(int splineOrder) {
46 for(int l=0; l<VDimension; l++) temp[l] = splineOrder;
47 SetSplineOrders(temp);
49 //====================================================================
51 //====================================================================
52 template<class TCoefficientType, int VDimension>
53 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
54 SetSplineOrders(const SizeType & splineOrder) {
55 // Compute support size
57 for(int l=0; l<VDimension; l++) {
58 mSplineOrders[l] = splineOrder[l];
59 mSplineSupport[l] = mSplineOrders[l]+1;
60 mSupportSize *= mSplineSupport[l];
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];
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;
81 mWeightsAreUpToDate = false;
83 //====================================================================
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;
92 //====================================================================
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;
101 //====================================================================
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];
111 //====================================================================
113 //====================================================================
114 template<class T> inline T factorial(T rhs) {
116 for(T x=(T)1; x<=rhs; ++x) {
121 //====================================================================
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));
130 //====================================================================
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.
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;
149 mBasisFunctionCoefficientsMatrix[order][i][j] *= BinomialCoefficient(s-1, i);
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;
158 mBasisFunctionCoefficientsMatrixAreUpToDate[order] = true;
160 //====================================================================
162 //====================================================================
163 template<class TCoefficientType, int VDimension>
164 TCoefficientType BSplineWeightsCalculator<TCoefficientType, VDimension>::
165 BSplineEvaluate(int order, int k, double e) {
166 // Evaluate a BSpline
168 TCoefficientType v = 0.0;
169 for(int p=0; p<s; p++) {
170 v += pow(e,order-p)*GetInitialWeights(order)[p][k];
174 //====================================================================
176 //====================================================================
177 template<class TCoefficientType, int VDimension>
178 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
179 ComputeSampledWeights1D(std::vector<std::vector<TCoefficientType> > & w, int order, int sampling) {
182 for(int k=0; k<s; k++) w[k].resize(sampling);
183 double offset = 1.0/sampling;
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
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);
195 if (fabs(1.0-e)<=1e-6) e = 1.0; // (for even order)
196 if (e>=1.0) e = e-1.0;
199 //====================================================================
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]);
210 mWeightsAreUpToDate = true;
212 //====================================================================
214 //====================================================================
215 template<class TCoefficientType, int VDimension>
216 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
217 ComputeTensorProducts() {
218 // Initial BSpline samples weights
219 ComputeSampledWeights();
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];
227 // Tensor product initialisation
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;
236 for(int a=0; a<nbPositions; a++) { // Loop over sampling positions
237 mTensorProducts[a].resize(mSupportSize);
239 for(int k=0; k<mSupportSize; k++) { // Loop over support positions
240 TCoefficientType B = 1.0;
242 for(int l=0; l<VDimension; l++) { // loop for tensor product
243 B *= mWeights[l][mSupportIndex[k][l]][mPositionIndex[l]];
245 mTensorProducts[a][k] = B;
248 // Next sample Position index
253 if (mPositionIndex[l] == (int)mSamplingFactors[l]) {
254 mPositionIndex[l] = 0;
261 //====================================================================
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]);
270 //====================================================================
272 #endif /* end #define ITKBSPLINEWEIGHTSCALCULATOR_TXX */