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