Developer Documentation

BSplineBasis.cc
1 /*===========================================================================*\
2 * *
3 * OpenFlipper *
4 * Copyright (c) 2001-2015, RWTH-Aachen University *
5 * Department of Computer Graphics and Multimedia *
6 * All rights reserved. *
7 * www.openflipper.org *
8 * *
9 *---------------------------------------------------------------------------*
10 * This file is part of OpenFlipper. *
11 *---------------------------------------------------------------------------*
12 * *
13 * Redistribution and use in source and binary forms, with or without *
14 * modification, are permitted provided that the following conditions *
15 * are met: *
16 * *
17 * 1. Redistributions of source code must retain the above copyright notice, *
18 * this list of conditions and the following disclaimer. *
19 * *
20 * 2. Redistributions in binary form must reproduce the above copyright *
21 * notice, this list of conditions and the following disclaimer in the *
22 * documentation and/or other materials provided with the distribution. *
23 * *
24 * 3. Neither the name of the copyright holder nor the names of its *
25 * contributors may be used to endorse or promote products derived from *
26 * this software without specific prior written permission. *
27 * *
28 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS *
29 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED *
30 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A *
31 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER *
32 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, *
33 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, *
34 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR *
35 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF *
36 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING *
37 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
38 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
39 * *
40 \*===========================================================================*/
41
42
43
44
45
46
47 //=============================================================================
48 //
49 // CLASS Geometry - IMPLEMENTATION
50 //
51 //=============================================================================
52
53 #define ACG_BSPLINEBASIS_C
54
55 //== INCLUDES =================================================================
56
57 #include "BSplineBasis.hh"
58
59 //----------------------------------------------------------------------------
60
61
62 namespace ACG {
63
64
65 //== IMPLEMENTATION ==========================================================
66
67
68
69 template<typename Scalar>
70 Vec2i
71 bsplineSpan(Scalar _t,
72 int _degree,
73 const std::vector<Scalar>& _knots)
74 {
75 // binary search
76
77 int lo = _degree;
78 int hi = _knots.size() - 1 - _degree;
79
80 Scalar upperBound = _knots[hi];
81
82 if (_t >= upperBound)
83 return Vec2i(hi-1-_degree, hi - 1);
84
85 int mid = (lo + hi) >> 1;
86
87 while (_t < _knots[mid] || _t >= _knots[mid + 1])
88 {
89 if (_t < _knots[mid])
90 hi = mid;
91 else
92 lo = mid;
93
94 mid = (lo + hi) >> 1;
95 }
96
97 return Vec2i(mid - _degree, mid);
98
99
100
101 // linear search:
102 //
103 // int i = 0;
104 //
105 // if (_t >= upperBound)
106 // i = _knots.size() - _degree - 2;
107 // else
108 // {
109 // while (_t >= _knots[i]) i++;
110 // while (_t < _knots[i]) i--;
111 // }
112 //
113 // return Vec2i(i - _degree, i);
114 }
115
116
117 template<typename Scalar>
118 void
119 bsplineBasisFunctions( std::vector<Scalar>& _N,
120 const Vec2i& _span,
121 Scalar _t,
122 const std::vector<Scalar>& _knots)
123 {
124 // inverted triangular scheme
125 // "The NURBS Book" : Algorithm A2.2 p70
126
127
128 // degree
129 int p = _span[1] - _span[0];
130
131 int i = _span[1];
132
133
134 _N[0] = 1.0;
135
136 // alloc temp buffer
137 static std::vector<Scalar> left(p+1);
138 static std::vector<Scalar> right(p+1);
139
140 if (left.size() < size_t(p+1))
141 {
142 left.resize(p+1);
143 right.resize(p+1);
144 }
145
146 // compute
147 for (int j = 1; j <= p; ++j)
148 {
149 left[j] = _t - _knots[i + 1 - j];
150 right[j] = _knots[i + j] - _t;
151
152 Scalar saved = 0.0;
153
154 for (int r = 0; r < j; ++r)
155 {
156 Scalar tmp = _N[r] / (right[r + 1] + left[j - r]);
157 _N[r] = saved + right[r + 1] * tmp;
158 saved = left[j - r] * tmp;
159 }
160 _N[j] = saved;
161 }
162 }
163
164
165 template<typename Scalar>
166 void
167 bsplineBasisDerivatives( std::vector<Scalar>& _ders,
168 const Vec2i& _span,
169 Scalar _t,
170 int _der,
171 const std::vector<Scalar>& _knots,
172 std::vector<Scalar>* _functionVals )
173 {
174 // The NURBS Book p72 algorithm A2.3
175
176 int p = _span[1] - _span[0];
177 int p1 = p+1;
178 int i = _span[1];
179
180
181 // alloc temp arrays
182 static std::vector< std::vector<Scalar> > ndu(p1);
183 static std::vector<Scalar> left(p1);
184 static std::vector<Scalar> right(p1);
185 static std::vector<Scalar> a(2*p1);
186
187 if (ndu[0].size() < size_t(p1))
188 {
189 ndu.resize(p1);
190 for (int i = 0; i < p1; ++i)
191 ndu[i].resize(p1);
192
193 left.resize(p1);
194 right.resize(p1);
195
196 a.resize(2*p1);
197 }
198
199
200
201
202 ndu[0][0] = 1.0;
203
204 for (int j = 1; j <= p; ++j)
205 {
206 left[j]= _t - _knots[i + 1 - j];
207 right[j]= _knots[i + j] - _t;
208 Scalar saved = 0.0;
209
210 for (int r = 0; r < j; ++r)
211 {
212 // lower triangle
213 ndu[j][r] = right[r+1] + left[j-r];
214 Scalar tmp = ndu[r][j-1] / ndu[j][r];
215
216 ndu[r][j] = saved + right[r+1] * tmp;
217 saved = left[j-r] * tmp;
218 }
219 ndu[j][j] = saved;
220 }
221
222 // load the basis functions
223 if (_functionVals)
224 {
225 for (int j = 0; j <= p; ++j)
226 (*_functionVals)[j] = ndu[j][p];
227 }
228
229 // compute derivatives
230
231
232
233 for (int r = 0; r <= p; ++r)
234 {
235 int s1 = 0, s2 = p1; // s1, s2: row offsets of linearized 2d array a[2][p+1]
236 a[0] = 1.0;
237
238 for (int k = 1; k <= _der; ++k)
239 {
240 Scalar d = 0.0;
241 int rk = r - k, pk = p - k;
242
243 if (r >= k)
244 {
245 a[s2 + 0] = a[s1 + 0] / ndu[pk+1][rk];
246 d = a[s2] * ndu[rk][pk];
247 }
248
249 int j1 = (rk >= -1) ? 1 : -rk;
250 int j2 = (r - 1 <= pk) ? k-1 : p-r;
251
252 for (int j = j1; j <= j2; ++j)
253 {
254 a[s2 + j] = (a[s1 + j] - a[s1 + j-1]) / ndu[pk+1][rk+j];
255 d += a[s2 + j] * ndu[rk+j][pk];
256 }
257
258 if (r <= pk)
259 {
260 a[s2 + k] = -a[s1 + k-1] / ndu[pk+1][r];
261 d += a[s2 + k] * ndu[r][pk];
262 }
263
264 if (k == _der)
265 _ders[r] = d;
266
267 std::swap(s1, s2); // switch rows
268 }
269 }
270
271 // correct factors
272
273 int r = p;
274
275 for (int k = 1; k <= _der; ++k)
276 {
277 Scalar rf = Scalar(r);
278 for (int j = 0; j <= p; ++j)
279 _ders[j] *= rf;
280
281 r *= (p - k);
282 }
283 }
284
285
286
287 template<typename Scalar>
288 Scalar
289 bsplineBasisFunction(int _i,
290 int _degree,
291 Scalar _t,
292 const std::vector<Scalar>& _knots)
293 {
294 int m = _knots.size() - 1;
295
296 // Mansfield Cox deBoor recursion
297 if ((_i == 0 && _t == _knots[0]) || (_i == m - _degree - 1 && _t == _knots[m]))
298 return 1.0;
299
300 if (_degree == 0) {
301 if (_t >= _knots[_i] && _t < _knots[_i + 1])
302 return 1.0;
303 else
304 return 0.0;
305 }
306
307 double Nin1 = basisFunction(_i, _degree - 1, _t, _knots);
308 double Nin2 = basisFunction(_i + 1, _degree - 1, _t, _knots);
309
310 double fac1 = 0;
311 // if ((_knotvector(_i+_n) - _knotvector(_i)) > 0.000001 )
312 if ((_knots[_i + _degree] - _knots[_i]) != 0)
313 fac1 = (_t - _knots[_i]) / (_knots[_i + _degree] - _knots[_i]);
314
315 double fac2 = 0;
316 // if ( (_knotvector(_i+1+_n) - _knotvector(_i+1)) > 0.000001 )
317 if ((_knots[_i + 1 + _degree] - _knots[_i + 1]) != 0)
318 fac2 = (_knots[_i + 1 + _degree] - _t) / (_knots[_i + 1 + _degree] - _knots[_i + 1]);
319
320 // std::cout << "Nin1 = " << Nin1 << ", Nin2 = " << Nin2 << ", fac1 = " << fac1 << ", fac2 = " << fac2 << std::endl;
321
322 return (fac1*Nin1 + fac2*Nin2);
323 }
324
325
326 template<typename Scalar>
327 Scalar
328 bsplineBasisDerivative(int _i,
329 int _degree,
330 Scalar _t,
331 int _der,
332 const std::vector<Scalar>& _knots)
333 {
334 assert(_degree >= 0);
335 assert(_i >= 0);
336
337
338 if (_der == 0)
339 return basisFunction(_i, _degree, _t, _knots);
340
341 Scalar Nin1 = derivativeBasisFunction(_i, _degree - 1, _t, _der - 1, _knots);
342 Scalar Nin2 = derivativeBasisFunction(_i + 1, _degree - 1, _t, _der - 1, _knots);
343
344 Scalar fac1 = 0;
345 if (std::abs(_knots[_i + _degree] - _knots[_i]) > 1e-6)
346 fac1 = Scalar(_degree) / (_knots[_i + _degree] - _knots[_i]);
347
348 Scalar fac2 = 0;
349 if (std::abs(_knots[_i + _degree + 1] - _knots[_i + 1]) > 1e-6)
350 fac2 = Scalar(_degree) / (_knots[_i + _degree + 1] - _knots[_i + 1]);
351
352 return (fac1*Nin1 - fac2*Nin2);
353 }
354
355 //=============================================================================
356 } // namespace ACG
357 //=============================================================================
ACG
Namespace providing different geometric functions concerning angles.
Definition: BaseObjectData.hh:68
ACG::bsplineSpan
Vec2i bsplineSpan(Scalar _t, int _degree, const std::vector< Scalar > &_knots)
Find the span of a parameter value.
Definition: BSplineBasis.cc:71
ACG::Vec2i
VectorT< signed int, 2 > Vec2i
Definition: VectorT.hh:98
ACG::bsplineBasisDerivatives
void bsplineBasisDerivatives(std::vector< Scalar > &_ders, const Vec2i &_span, Scalar _t, int _der, const std::vector< Scalar > &_knots, std::vector< Scalar > *_functionVals)
Compute derivatives of basis functions in a span.
Definition: BSplineBasis.cc:167
ACG::bsplineBasisFunction
Scalar bsplineBasisFunction(int _i, int _degree, Scalar _t, const std::vector< Scalar > &_knots)
Evaluate a single basis function.
Definition: BSplineBasis.cc:289
ACG::bsplineBasisFunctions
void bsplineBasisFunctions(std::vector< Scalar > &_N, const Vec2i &_span, Scalar _t, const std::vector< Scalar > &_knots)
Evaluate basis functions in a span.
Definition: BSplineBasis.cc:119
ACG::bsplineBasisDerivative
Scalar bsplineBasisDerivative(int _i, int _degree, Scalar _t, int _der, const std::vector< Scalar > &_knots)
Compute derivative of a single basis function.
Definition: BSplineBasis.cc:328
OpenMesh::VectorT
Definition: Vector11T.hh:83