Developer Documentation
BSplineCurveT_impl.hh
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 // CLASS BSplineCurveT - IMPLEMENTATION
48 // Author: Ellen Dekkers <dekkers@cs.rwth-aachen.de>
49 //
50 //=============================================================================
51 
52 #define BSPLINECURVE_BSPLINECURVET_C
53 
54 //== INCLUDES =================================================================
55 
56 #include <OpenMesh/Core/Geometry/VectorT.hh>
57 
58 #include <iostream>
59 #include <fstream>
60 
61 #include "BSplineCurve.hh"
62 
63 #include <cfloat>
64 #include <ACG/Geometry/Algorithms.hh>
65 #include <ACG/Math/BSplineBasis.hh>
66 
67 //== NAMESPACES ===============================================================
68 
69 namespace ACG {
70 
71 //== IMPLEMENTATION ==========================================================
72 
73 
74 template <class PointT>
76 BSplineCurveT( unsigned int _degree )
77 : degree_(_degree),
78  autocompute_knotvector_(true),
79  fix_number_control_points_(false),
80  ref_count_cpselections_(0),
81  ref_count_eselections_(0)
82 {
83 }
84 
85 //-----------------------------------------------------------------------------
86 
87 template <class PointT>
89 BSplineCurveT( const BSplineCurveT& _curve )
90 {
91  //copy control points
92  control_polygon_ = _curve.control_polygon_;
93 
94  //copy knotvector
95  knotvector_ = _curve.knotvector_;
96 
97  degree_ = _curve.degree_;
98  autocompute_knotvector_ = _curve.autocompute_knotvector_;
99  fix_number_control_points_ = _curve.fix_number_control_points_;
100 
101  // copy properties
102  cpselections_ = _curve.cpselections_;
103  eselections_ = _curve.eselections_;
104 
105  // copy property reference counter
106  ref_count_cpselections_ = _curve.ref_count_cpselections_;
107  ref_count_eselections_ = _curve.ref_count_eselections_;
108 
109 }
110 
111 //-----------------------------------------------------------------------------
112 
113 template <class PointT>
114 template <class PropT>
115 void
117 request_prop( unsigned int& _ref_count, PropT& _prop)
118 {
119  if(_ref_count == 0)
120  {
121  _ref_count = 1;
122  // always use vertex size!!!
123  _prop.resize(n_control_points());
124  }
125  else ++_ref_count;
126 }
127 
128 //-----------------------------------------------------------------------------
129 
130 template <class PointT>
131 template <class PropT>
132 void
134 release_prop( unsigned int& _ref_count, PropT& _prop)
135 {
136  if( _ref_count <= 1)
137  {
138  _ref_count = 0;
139  _prop.clear();
140  }
141  else --_ref_count;
142 }
143 
144 //-----------------------------------------------------------------------------
145 
146 template <class PointT>
147 void
149 add_control_point(const Point& _cp)
150 {
151  if (fix_number_control_points_)
152  return;
153 
154  // add new point
155  control_polygon_.push_back(_cp);
156 
157  // update knot vector
158  if (autocompute_knotvector_)
159  knotvector_.createKnots(degree_, control_polygon_.size());
160 
161  // add available properties
162  if( controlpoint_selections_available())
163  cpselections_.push_back( false);
164 
165  if( edge_selections_available())
166  eselections_.push_back(false);
167 }
168 
169 //-----------------------------------------------------------------------------
170 
171 template <class PointT>
172 void
174 insert_control_point(int _idx, const Point& _cp)
175 {
176  if (fix_number_control_points_)
177  return;
178 
179  assert(_idx < (int)control_polygon_.size());
180  control_polygon_.insert(control_polygon_.begin()+_idx, _cp);
181 
182  // update knot vector
183  if (autocompute_knotvector_)
184  knotvector_.createKnots(degree_, control_polygon_.size());
185  else{
186  // compute knot in between its wo neighbors
187  double knot = ( knotvector_(_idx-1) + knotvector_(_idx+1) ) / 2.0;
188  knotvector_.insertKnot(_idx, knot);
189  }
190 
191  // add available properties
192  if( controlpoint_selections_available())
193  cpselections_.insert(cpselections_.begin()+_idx, false);
194 
195  if( edge_selections_available())
196  eselections_.insert(eselections_.begin()+_idx, false);
197 }
198 
199 
200 //-----------------------------------------------------------------------------
201 
202 template <class PointT>
203 void
206 {
207  if (fix_number_control_points_)
208  return;
209 
210  assert(_idx < (int)control_polygon_.size());
211  control_polygon_.erase(control_polygon_.begin()+_idx);
212 
213  // update knot vector
214  if (autocompute_knotvector_)
215  knotvector_.createKnots(degree_, control_polygon_.size());
216  else{
217  // delete knot at given index
218  knotvector_.deleteKnot(_idx);
219  }
220 
221  // add available properties
222  if( controlpoint_selections_available())
223  cpselections_.erase(cpselections_.begin()+_idx);
224 
225  if( edge_selections_available())
226  eselections_.erase(eselections_.begin()+_idx);
227 }
228 
229 //-----------------------------------------------------------------------------
230 
231 template <class PointT>
232 void
234 set_control_point(int _idx, const Point& _cp)
235 {
236  assert(_idx < (int)control_polygon_.size());
237  control_polygon_[_idx] = _cp;
238 }
239 
240 //-----------------------------------------------------------------------------
241 
242 template <class PointT>
243 void
245 set_control_polygon(std::vector< Point> & _control_polygon)
246 {
247  control_polygon_.clear();
248 
249  for (unsigned int i = 0; i < _control_polygon.size(); ++i)
250  control_polygon_.push_back(_control_polygon[i]);
251 
252  cpselections_.clear();
253  if( controlpoint_selections_available())
254  cpselections_.resize(control_polygon_.size(), false);
255 
256  eselections_.clear();
257  if( edge_selections_available())
258  eselections_.resize(control_polygon_.size(), false);
259 }
260 
261 //-----------------------------------------------------------------------------
262 
263 template <class PointT>
264 void
267 {
268  control_polygon_.clear();
269 
270  // reset properties
271  cpselections_.clear();
272  eselections_.clear();
273 }
274 
275 //-----------------------------------------------------------------------------
276 
277 template <class PointT>
278 PointT
280 curvePoint(Scalar _u)
281 {
282  Scalar epsilon = 0.0000001;
283 
284  if (_u > upper() && _u < upper()+epsilon)
285  _u = upper();
286 
287  assert(_u >= lower() && _u <= upper());
288 
289  int p = degree(); // spline degree
290 
291  Vec2i span = ACG::bsplineSpan(_u, p, knotvector_.getKnotvector());
292 
293  // eval non-zero basis functions
294  std::vector<Scalar> N(p+1);
295  ACG::bsplineBasisFunctions(N, span, _u, knotvector_.getKnotvector());
296 
297  Point point = Point(0.0, 0.0, 0.0);
298 
299  for (int i = 0; i <= p; ++i)
300  point += get_control_point(span[0] + i) * N[i];
301 
302  return point;
303 }
304 
305 //-----------------------------------------------------------------------------
306 
307 template <class PointT>
308 PointT
310 derivativeCurvePoint(Scalar _u, unsigned int _der)
311 {
312  assert(_u >= lower() && _u <= upper());
313 
314  int p = degree(); // spline degree
315 
316  Vec2i span = ACG::bsplineSpan(_u, p, knotvector_.getKnotvector());
317 
318  // eval non-zero basis functions
319  std::vector<Scalar> dNdu(p+1);
320  ACG::bsplineBasisDerivatives(dNdu, span, _u, 1, knotvector_.getKnotvector(), 0);
321 
322  Point point = Point(0.0, 0.0, 0.0);
323 
324  for (int i = 0; i <= p; ++i)
325  point += get_control_point(i + span[0]) * dNdu[i];
326 
327  return point;
328 }
329 
330 //-----------------------------------------------------------------------------
331 
332 template <class PointT>
333 typename BSplineCurveT<PointT>::Scalar
335 basisFunction(int _i, int _n, Scalar _t)
336 {
337  int m = knotvector_.size() - 1;
338 
339  // Mansfield Cox deBoor recursion
340  if ((_i==0 && _t== knotvector_(0)) || (_i==m-_n-1 && _t==knotvector_(m)))
341  return 1.0;
342 
343  if (_n == 0){
344  if (_t >= knotvector_(_i) && _t < knotvector_(_i+1))
345  return 1.0;
346  else
347  return 0.0;
348  }
349 
350  typename BSplineCurveT<PointT>::Scalar Nin1 = basisFunction(_i, _n-1, _t);
351  typename BSplineCurveT<PointT>::Scalar Nin2 = basisFunction(_i+1, _n-1, _t);
352 
353  Scalar fac1 = 0;
354  if ((knotvector_(_i+_n)-knotvector_(_i)) !=0 )
355  fac1 = (_t - knotvector_(_i)) / (knotvector_(_i+_n) - knotvector_(_i)) ;
356 
357  Scalar fac2 = 0;
358  if ( (knotvector_(_i+1+_n)-knotvector_(_i+1)) !=0 )
359  fac2 = (knotvector_(_i+1+_n) - _t)/ (knotvector_(_i+1+_n) - knotvector_(_i+1));
360 
361 // std::cout << "N " << _i << " " << _n-1 << " = " << Nip1 << ", fac1 = " << fac1
362 // << ", N " << _i+1 << " " << _n-1 << " = " << Nip2 << ", fac2 = " << fac2
363 // << ", result = " << (fac1*Nip1 - fac2*Nip2)
364 // << std::endl;
365 
366  return (fac1*Nin1 + fac2*Nin2);
367 }
368 
369 //-----------------------------------------------------------------------------
370 
371 template <class PointT>
372 typename BSplineCurveT<PointT>::Scalar
374 derivativeBasisFunction(int _i, int _n, Scalar _t, int _der)
375 {
376  if (_der == 0)
377  return basisFunction(_i, _n, _t);
378 
379  typename BSplineCurveT<PointT>::Scalar Nin1 = derivativeBasisFunction(_i, _n-1, _t, _der-1);
380  typename BSplineCurveT<PointT>::Scalar Nin2 = derivativeBasisFunction(_i+1, _n-1, _t, _der-1);
381 // std::cout << "der = " << _der << ", i = " << _i << ", n = " << _n << ", t = " << _t << " ===> " << Nin1 << " , " << Nin2 << std::endl;
382 // std::cout << "Nin1 = " << Nin1 << ", Nin2 = " << Nin2 << std::endl;
383 
384  double fac1 = 0;
385  if ( (knotvector_(_i+_n)-knotvector_(_i)) !=0 )
386  fac1 = double(_n) / (knotvector_(_i+_n)-knotvector_(_i));
387 
388  double fac2 = 0;
389  if ( (knotvector_(_i+_n+1)-knotvector_(_i+1)) !=0 )
390  fac2 = double(_n) / (knotvector_(_i+_n+1)-knotvector_(_i+1));
391 
392 // std::cout << "fac1 = " << fac1 << ", fac2 = " << fac2 << std::endl;
393  return (fac1*Nin1 - fac2*Nin2);
394 }
395 
396 //-----------------------------------------------------------------------------
397 
398 template <class PointT>
399 std::vector<PointT>
401 deBoorAlgorithm( double _u)
402 {
403 // std::cout << "\n Evaluating via deBoor algorithm at " << _u << std::endl;
404  assert(_u >= lower() && _u <= upper());
405 
406  int n = degree();
407 
408  Point point = Point(0.0, 0.0, 0.0);
409 
410  ACG::Vec2i span_u = span(_u);
411 
412  std::vector<Point> allPoints;
413 
414  // control points in current iteration
415  std::vector<Point> controlPoints_r;
416 
417  // control points in last iteration
418  std::vector<Point> controlPoints_r1;
419  for (int i = span_u[0]; i <= span_u[1]; ++i)
420  controlPoints_r1.push_back(control_polygon_[i]);
421 
422  for (int r = 1; r <= n; ++r)
423  {
424  controlPoints_r.clear();
425 
426  for (int i = r; i <= n; ++i)
427  {
428  // compute alpha _i ^ n-r
429  double alpha = (_u - knotvector_(span_u[0]+i)) / (knotvector_(span_u[0]+i + n + 1 - r) - knotvector_(span_u[0]+i));
430  Point c = controlPoints_r1[i-r] * (1.0 - alpha) + controlPoints_r1[i-r+1] * alpha;
431 
432  // save
433  controlPoints_r.push_back(c);
434  allPoints.push_back(c);
435  }
436 
437  controlPoints_r1 = controlPoints_r;
438  }
439 
440  return allPoints;
441 }
442 
443 //-----------------------------------------------------------------------------
444 
445 template <class PointT>
446 typename BSplineCurveT<PointT>::Scalar
448 lower() const
449 {
450  return knotvector_(degree());
451 }
452 
453 //-----------------------------------------------------------------------------
454 
455 template <class PointT>
456 typename BSplineCurveT<PointT>::Scalar
458 upper() const
459 {
460  return knotvector_(knotvector_.size() - 1 - degree());
461 }
462 
463 //-----------------------------------------------------------------------------
464 
465 template <class PointT>
468 span(double _t)
469 {
470  unsigned int i(0);
471 
472  if ( _t >= upper())
473  i = n_control_points() - 1;
474  else {
475  while (_t >= knotvector_(i)) i++;
476  while (_t < knotvector_(i)) i--;
477  }
478 
479  return Vec2i(i-degree(), i);
480 }
481 
482 //-----------------------------------------------------------------------------
483 
484 template <class PointT>
487 interval(double _t)
488 {
489  unsigned int i(0);
490  while (_t >= knotvector_(i)) i++;
491  while (_t < knotvector_(i)) i--;
492 
493  return Vec2i(i, i+1);
494 }
495 
496 //-----------------------------------------------------------------------------
497 
498 template <class PointT>
499 void
501 print() const
502 {
503  std::cerr << "****** BSplineInfo ******\n";
504 
505  std::cerr << "#control_points: " << control_polygon_.size() << std::endl;
506  for (unsigned int i = 0; i < control_polygon_.size(); ++i)
507  std::cerr << get_control_point(i) << std::endl;
508 }
509 
510 //-----------------------------------------------------------------------------
511 
512 template <class PointT>
513 void
515 insertKnot(double _u)
516 {
517  // insert a knot at parameter u and recompute controlpoint s.t. curve stays the same
518 
519 // std::cout << "Lower: " << lower() << ", upper: " << upper() << std::endl;
520  assert(_u >= lower() && _u <= upper());
521 
522  ACG::Vec2i span_u = span(_u);
523 
524  std::vector<Point> updateControlPoints;
525 
526  // keep control points that are not affected by knot insertion
527  for (int i = 0; i <= span_u[0]; ++i)
528  updateControlPoints.push_back(get_control_point(i));
529 
530  // add control points in second column of triangular array
531  std::vector<Vec3d> cp = deBoorAlgorithm(_u);
532  for (unsigned int i = 0; i < degree_; ++i)
533  updateControlPoints.push_back(cp[i]);
534 
535  // add remaining control points not affected by knot insertion
536  for (unsigned int i = span_u[1]; i < n_control_points(); ++i)
537  updateControlPoints.push_back(get_control_point(i));
538 
539 
540  // find out where to insert the new knot
541  int index = 0;
542  for (unsigned int i = 0; i < knotvector_.size(); ++i)
543  if (_u > knotvector_(i))
544  ++index;
545 
546  // update knot vector
547  knotvector_.insertKnot(index, _u);
548 
549  // update control points
550  set_control_polygon(updateControlPoints);
551 }
552 
553 //-----------------------------------------------------------------------------
554 
555 template <class PointT>
556 void
558 set_knots(std::vector< Scalar > _knots)
559 {
560  autocompute_knotvector(false);
561 
562  // set the knotvector
563  knotvector_.setKnotvector(_knots);
564 }
565 
566 //-----------------------------------------------------------------------------
567 
568 template <class PointT>
569 void
572 {
573  // reverse vector of controlpoints
574  std::vector<Point> reversed_control_polygon;
575  for (int i = n_control_points()-1; i > -1; --i)
576  reversed_control_polygon.push_back( get_control_point(i) );
577 
578  // compute new knot vector
579  int num_knots = n_knots();
580  std::vector< Scalar > reversed_knotvector(num_knots) ;
581 
582  int m = num_knots-1;
583  int p = degree();
584 
585  double a = get_knot(0);
586  double b = get_knot(num_knots-1);
587 
588  for (int i = 0; i <= p; ++i)
589  reversed_knotvector[i] = get_knot(i);
590 
591  for (int i = m-p-1; i <= m; ++i)
592  reversed_knotvector[i] = get_knot(i);
593 
594  for (int i = 1; i < m - 2*p; ++i)
595  reversed_knotvector[m - p - i] = - get_knot(p+i) + a + b;
596 
597  // reset control polygon
598  set_control_polygon(reversed_control_polygon);
599 
600  // reset knot vector
601  set_knots(reversed_knotvector);
602 }
603 
604 //-----------------------------------------------------------------------------
605 
606 //=============================================================================
607 } // namespace ACG
608 //=============================================================================
VectorT< signed int, 2 > Vec2i
Definition: VectorT.hh:98
Point curvePoint(Scalar _u)
unsigned int n_knots() const
Returns the number of knots.
Point derivativeCurvePoint(Scalar _u, unsigned int _der)
void set_control_point(int _idx, const Point &_cp)
reset a control point
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.
Namespace providing different geometric functions concerning angles.
void add_control_point(const Point &_cp)
add a control point
void insert_control_point(int _idx, const Point &_cp)
insert a control point at given index
void reset_control_polygon()
Clears the control polygon.
unsigned int degree() const
Returns the spline degree.
Scalar upper() const
Returns the upper parameter.
void delete_control_point(int _idx)
delete control point at given index
Vec2i bsplineSpan(Scalar _t, int _degree, const std::vector< Scalar > &_knots)
Find the span of a parameter value.
Definition: BSplineBasis.cc:71
Scalar lower() const
Returns the lower parameter.
ACG::Vec2i span(double _t)
unsigned int n_control_points() const
Returns the number of control points.
Point & get_control_point(int _i)
get control point i
void bsplineBasisFunctions(std::vector< Scalar > &_N, const Vec2i &_span, Scalar _t, const std::vector< Scalar > &_knots)
Evaluate basis functions in a span.
Scalar derivativeBasisFunction(int _i, int _n, Scalar _t, int _der)
void print() const
print information string
void set_control_polygon(std::vector< Point > &_control_polygon)
set whole control polygon
Scalar basisFunction(int _i, int _n, Scalar _t)
ACG::Vec2i interval(double _t)
void reverse()
Reverses the curve.
void set_knots(std::vector< Scalar > _knots)
set the knotvector of the bspline curve
BSplineCurveT(unsigned int _degree=3)
Constructor.
Scalar get_knot(int _i)
get knot i
std::vector< Point > deBoorAlgorithm(double _u)
void insertKnot(double _u)
Inserts a new knot at parameter u.