1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ======================================================================-====*/
18 /* =========================================================================
20 @file itkBSplineWeightsCalculator.h
21 @author David Sarrut <david.sarrut@creatis.insa-lyon.fr>
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
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.
34 ========================================================================= */
36 #ifndef ITKBSPLINEWEIGHTSCALCULATOR_TXX
37 #define ITKBSPLINEWEIGHTSCALCULATOR_TXX
39 //====================================================================
40 template<int VDimension>
41 typename Index<VDimension>::IndexValueType Index2Offset(const Index<VDimension> & index, const Index<VDimension> & offsetTable) {
43 for(int l=1; l<VDimension; l++) {
44 v += index[l] * offsetTable[l];
48 //====================================================================
50 //====================================================================
51 template<class TCoefficientType, int VDimension>
52 BSplineWeightsCalculator<TCoefficientType, VDimension>::
53 BSplineWeightsCalculator() {
54 mWeightsAreUpToDate = false;
56 //====================================================================
58 //====================================================================
59 template<class TCoefficientType, int VDimension>
60 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
61 SetSplineOrder(int splineOrder) {
63 for(int l=0; l<VDimension; l++) temp[l] = splineOrder;
64 SetSplineOrders(temp);
66 //====================================================================
68 //====================================================================
69 template<class TCoefficientType, int VDimension>
70 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
71 SetSplineOrders(const SizeType & splineOrder) {
72 // Compute support size
74 for(int l=0; l<VDimension; l++) {
75 mSplineOrders[l] = splineOrder[l];
76 mSplineSupport[l] = mSplineOrders[l]+1;
77 mSupportSize *= mSplineSupport[l];
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];
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;
98 mWeightsAreUpToDate = false;
100 //====================================================================
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;
109 //====================================================================
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;
118 //====================================================================
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];
128 //====================================================================
130 //====================================================================
131 template<class T> inline T factorial(T rhs) {
133 for(T x=(T)1; x<=rhs; ++x) {
138 //====================================================================
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));
147 //====================================================================
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.
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;
166 mBasisFunctionCoefficientsMatrix[order][i][j] *= BinomialCoefficient(s-1, i);
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;
175 mBasisFunctionCoefficientsMatrixAreUpToDate[order] = true;
177 //====================================================================
179 //====================================================================
180 template<class TCoefficientType, int VDimension>
181 TCoefficientType BSplineWeightsCalculator<TCoefficientType, VDimension>::
182 BSplineEvaluate(int order, int k, double e) {
183 // Evaluate a BSpline
185 TCoefficientType v = 0.0;
186 for(int p=0; p<s; p++) {
187 v += pow(e,order-p)*GetInitialWeights(order)[p][k];
191 //====================================================================
193 //====================================================================
194 template<class TCoefficientType, int VDimension>
195 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
196 ComputeSampledWeights1D(std::vector<std::vector<TCoefficientType> > & w, int order, int sampling) {
199 for(int k=0; k<s; k++) w[k].resize(sampling);
200 double offset = 1.0/sampling;
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
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);
212 if (fabs(1.0-e)<=1e-6) e = 1.0; // (for even order)
213 if (e>=1.0) e = e-1.0;
216 //====================================================================
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]);
227 mWeightsAreUpToDate = true;
229 //====================================================================
231 //====================================================================
232 template<class TCoefficientType, int VDimension>
233 void BSplineWeightsCalculator<TCoefficientType, VDimension>::
234 ComputeTensorProducts() {
235 // Initial BSpline samples weights
236 ComputeSampledWeights();
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];
244 // Tensor product initialisation
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;
253 for(int a=0; a<nbPositions; a++) { // Loop over sampling positions
254 mTensorProducts[a].resize(mSupportSize);
256 for(int k=0; k<mSupportSize; k++) { // Loop over support positions
257 TCoefficientType B = 1.0;
259 for(int l=0; l<VDimension; l++) { // loop for tensor product
260 B *= mWeights[l][mSupportIndex[k][l]][mPositionIndex[l]];
262 mTensorProducts[a][k] = B;
265 // Next sample Position index
270 if (mPositionIndex[l] == (int)mSamplingFactors[l]) {
271 mPositionIndex[l] = 0;
278 //====================================================================
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]);
287 //====================================================================
289 #endif /* end #define ITKBSPLINEWEIGHTSCALCULATOR_TXX */