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  * $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
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
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 //=============================================================================
void bsplineBasisFunctions(std::vector< Scalar > &_N, const Vec2i &_span, Scalar _t, const std::vector< Scalar > &_knots)
Evaluate basis functions in a span.
VectorT< signed int, 2 > Vec2i
Definition: VectorT.hh:104
Vec2i bsplineSpan(Scalar _t, int _degree, const std::vector< Scalar > &_knots)
Find the span of a parameter value.
Definition: BSplineBasis.cc:77
Scalar bsplineBasisDerivative(int _i, int _degree, Scalar _t, int _der, const std::vector< Scalar > &_knots)
Compute derivative of a single basis function.
Namespace providing different geometric functions concerning angles.
Definition: DBSCANT.cc:51
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.
Scalar bsplineBasisFunction(int _i, int _degree, Scalar _t, const std::vector< Scalar > &_knots)
Evaluate a single basis function.