Developer Documentation

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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 * $Revision$ *
45 * $Author$ *
46 * $Date$ *
47 * *
48 \*===========================================================================*/
49
50
51
52
53 //=============================================================================
54 //
55 // CLASS Geometry - IMPLEMENTATION
56 //
57 //=============================================================================
58
59 #define ACG_BSPLINEBASIS_C
60
61 //== INCLUDES =================================================================
62
63 #include "BSplineBasis.hh"
64
65 //----------------------------------------------------------------------------
66
67
68 namespace ACG {
69
70
71 //== IMPLEMENTATION ==========================================================
72
73
74
75 template<typename Scalar>
76 Vec2i
77 bsplineSpan(Scalar _t,
78 int _degree,
79 const std::vector<Scalar>& _knots)
80 {
81 // binary search
82
83 int lo = _degree;
84 int hi = _knots.size() - 1 - _degree;
85
86 Scalar upperBound = _knots[hi];
87
88 if (_t >= upperBound)
89 return Vec2i(hi-1-_degree, hi - 1);
90
91 int mid = (lo + hi) >> 1;
92
93 while (_t < _knots[mid] || _t >= _knots[mid + 1])
94 {
95 if (_t < _knots[mid])
96 hi = mid;
97 else
98 lo = mid;
99
100 mid = (lo + hi) >> 1;
101 }
102
103 return Vec2i(mid - _degree, mid);
104
105
106
107 // linear search:
108 //
109 // int i = 0;
110 //
111 // if (_t >= upperBound)
112 // i = _knots.size() - _degree - 2;
113 // else
114 // {
115 // while (_t >= _knots[i]) i++;
116 // while (_t < _knots[i]) i--;
117 // }
118 //
119 // return Vec2i(i - _degree, i);
120 }
121
122
123 template<typename Scalar>
124 void
125 bsplineBasisFunctions( std::vector<Scalar>& _N,
126 const Vec2i& _span,
127 Scalar _t,
128 const std::vector<Scalar>& _knots)
129 {
130 // inverted triangular scheme
131 // "The NURBS Book" : Algorithm A2.2 p70
132
133
134 // degree
135 int p = _span[1] - _span[0];
136
137 int i = _span[1];
138
139
140 _N[0] = 1.0;
141
142 // alloc temp buffer
143 static std::vector<Scalar> left(p+1);
144 static std::vector<Scalar> right(p+1);
145
146 if (left.size() < size_t(p+1))
147 {
148 left.resize(p+1);
149 right.resize(p+1);
150 }
151
152 // compute
153 for (int j = 1; j <= p; ++j)
154 {
155 left[j] = _t - _knots[i + 1 - j];
156 right[j] = _knots[i + j] - _t;
157
158 Scalar saved = 0.0;
159
160 for (int r = 0; r < j; ++r)
161 {
162 Scalar tmp = _N[r] / (right[r + 1] + left[j - r]);
163 _N[r] = saved + right[r + 1] * tmp;
164 saved = left[j - r] * tmp;
165 }
166 _N[j] = saved;
167 }
168 }
169
170
171 template<typename Scalar>
172 void
173 bsplineBasisDerivatives( std::vector<Scalar>& _ders,
174 const Vec2i& _span,
175 Scalar _t,
176 int _der,
177 const std::vector<Scalar>& _knots,
178 std::vector<Scalar>* _functionVals )
179 {
180 // The NURBS Book p72 algorithm A2.3
181
182 int p = _span[1] - _span[0];
183 int p1 = p+1;
184 int i = _span[1];
185
186
187 // alloc temp arrays
188 static std::vector< std::vector<Scalar> > ndu(p1);
189 static std::vector<Scalar> left(p1);
190 static std::vector<Scalar> right(p1);
191 static std::vector<Scalar> a(2*p1);
192
193 if (ndu[0].size() < size_t(p1))
194 {
195 ndu.resize(p1);
196 for (int i = 0; i < p1; ++i)
197 ndu[i].resize(p1);
198
199 left.resize(p1);
200 right.resize(p1);
201
202 a.resize(2*p1);
203 }
204
205
206
207
208 ndu[0][0] = 1.0;
209
210 for (int j = 1; j <= p; ++j)
211 {
212 left[j]= _t - _knots[i + 1 - j];
213 right[j]= _knots[i + j] - _t;
214 Scalar saved = 0.0;
215
216 for (int r = 0; r < j; ++r)
217 {
218 // lower triangle
219 ndu[j][r] = right[r+1] + left[j-r];
220 Scalar tmp = ndu[r][j-1] / ndu[j][r];
221
222 ndu[r][j] = saved + right[r+1] * tmp;
223 saved = left[j-r] * tmp;
224 }
225 ndu[j][j] = saved;
226 }
227
228 // load the basis functions
229 if (_functionVals)
230 {
231 for (int j = 0; j <= p; ++j)
232 (*_functionVals)[j] = ndu[j][p];
233 }
234
235 // compute derivatives
236
237
238
239 for (int r = 0; r <= p; ++r)
240 {
241 int s1 = 0, s2 = p1; // s1, s2: row offsets of linearized 2d array a[2][p+1]
242 a[0] = 1.0;
243
244 for (int k = 1; k <= _der; ++k)
245 {
246 Scalar d = 0.0;
247 int rk = r - k, pk = p - k;
248
249 if (r >= k)
250 {
251 a[s2 + 0] = a[s1 + 0] / ndu[pk+1][rk];
252 d = a[s2] * ndu[rk][pk];
253 }
254
255 int j1 = (rk >= -1) ? 1 : -rk;
256 int j2 = (r - 1 <= pk) ? k-1 : p-r;
257
258 for (int j = j1; j <= j2; ++j)
259 {
260 a[s2 + j] = (a[s1 + j] - a[s1 + j-1]) / ndu[pk+1][rk+j];
261 d += a[s2 + j] * ndu[rk+j][pk];
262 }
263
264 if (r <= pk)
265 {
266 a[s2 + k] = -a[s1 + k-1] / ndu[pk+1][r];
267 d += a[s2 + k] * ndu[r][pk];
268 }
269
270 if (k == _der)
271 _ders[r] = d;
272
273 std::swap(s1, s2); // switch rows
274 }
275 }
276
277 // correct factors
278
279 int r = p;
280
281 for (int k = 1; k <= _der; ++k)
282 {
283 Scalar rf = Scalar(r);
284 for (int j = 0; j <= p; ++j)
285 _ders[j] *= rf;
286
287 r *= (p - k);
288 }
289 }
290
291
292
293 template<typename Scalar>
294 Scalar
295 bsplineBasisFunction(int _i,
296 int _degree,
297 Scalar _t,
298 const std::vector<Scalar>& _knots)
299 {
300 int m = _knots.size() - 1;
301
302 // Mansfield Cox deBoor recursion
303 if ((_i == 0 && _t == _knots[0]) || (_i == m - _degree - 1 && _t == _knots[m]))
304 return 1.0;
305
306 if (_degree == 0) {
307 if (_t >= _knots[_i] && _t < _knots[_i + 1])
308 return 1.0;
309 else
310 return 0.0;
311 }
312
313 double Nin1 = basisFunction(_i, _degree - 1, _t, _knots);
314 double Nin2 = basisFunction(_i + 1, _degree - 1, _t, _knots);
315
316 double fac1 = 0;
317 // if ((_knotvector(_i+_n) - _knotvector(_i)) > 0.000001 )
318 if ((_knots[_i + _degree] - _knots[_i]) != 0)
319 fac1 = (_t - _knots[_i]) / (_knots[_i + _degree] - _knots[_i]);
320
321 double fac2 = 0;
322 // if ( (_knotvector(_i+1+_n) - _knotvector(_i+1)) > 0.000001 )
323 if ((_knots[_i + 1 + _degree] - _knots[_i + 1]) != 0)
324 fac2 = (_knots[_i + 1 + _degree] - _t) / (_knots[_i + 1 + _degree] - _knots[_i + 1]);
325
326 // std::cout << "Nin1 = " << Nin1 << ", Nin2 = " << Nin2 << ", fac1 = " << fac1 << ", fac2 = " << fac2 << std::endl;
327
328 return (fac1*Nin1 + fac2*Nin2);
329 }
330
331
332 template<typename Scalar>
333 Scalar
334 bsplineBasisDerivative(int _i,
335 int _degree,
336 Scalar _t,
337 int _der,
338 const std::vector<Scalar>& _knots)
339 {
340 assert(_degree >= 0);
341 assert(_i >= 0);
342
343
344 if (_der == 0)
345 return basisFunction(_i, _degree, _t, _knots);
346
347 Scalar Nin1 = derivativeBasisFunction(_i, _degree - 1, _t, _der - 1, _knots);
348 Scalar Nin2 = derivativeBasisFunction(_i + 1, _degree - 1, _t, _der - 1, _knots);
349
350 Scalar fac1 = 0;
351 if (std::abs(_knots[_i + _degree] - _knots[_i]) > 1e-6)
352 fac1 = Scalar(_degree) / (_knots[_i + _degree] - _knots[_i]);
353
354 Scalar fac2 = 0;
355 if (std::abs(_knots[_i + _degree + 1] - _knots[_i + 1]) > 1e-6)
356 fac2 = Scalar(_degree) / (_knots[_i + _degree + 1] - _knots[_i + 1]);
357
358 return (fac1*Nin1 - fac2*Nin2);
359 }
360
361 //=============================================================================
362 } // namespace ACG
363 //=============================================================================
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:125
ACG
Namespace providing different geometric functions concerning angles.
Definition: DBSCANT.cc:51
ACG::Vec2i
VectorT< signed int, 2 > Vec2i
Definition: VectorT.hh:104
ACG::bsplineSpan
Vec2i bsplineSpan(Scalar _t, int _degree, const std::vector< Scalar > &_knots)
Find the span of a parameter value.
Definition: BSplineBasis.cc:77
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:173
OpenMesh::VectorT
Definition: Vector11T.hh:83
ACG::bsplineBasisFunction
Scalar bsplineBasisFunction(int _i, int _degree, Scalar _t, const std::vector< Scalar > &_knots)
Evaluate a single basis function.
Definition: BSplineBasis.cc:295
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:334