1 /*! ========================================================================
2 ** Extended Template Library
3 ** Bezier Template Class Implementation
6 ** Copyright (c) 2002 Robert B. Quattlebaum Jr.
8 ** This package is free software; you can redistribute it and/or
9 ** modify it under the terms of the GNU General Public License as
10 ** published by the Free Software Foundation; either version 2 of
11 ** the License, or (at your option) any later version.
13 ** This package is distributed in the hope that it will be useful,
14 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
15 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 ** General Public License for more details.
18 ** === N O T E S ===========================================================
20 ** This is an internal header file, included by other ETL headers.
21 ** You should not attempt to use it directly.
23 ** ========================================================================= */
25 /* === S T A R T =========================================================== */
27 #ifndef __ETL_BEZIER_H
28 #define __ETL_BEZIER_H
30 /* === H E A D E R S ======================================================= */
32 #include "_curve_func.h"
35 /* === M A C R O S ========================================================= */
37 /* === T Y P E D E F S ===================================================== */
39 /* === C L A S S E S & S T R U C T S ======================================= */
43 template<typename V,typename T> class bezier;
45 //! Cubic Bezier Curve Base Class
46 // This generic implementation uses the DeCasteljau algorithm.
47 // Works for just about anything that has an affine combination function
48 template <typename V,typename T=float>
49 class bezier_base : public std::unary_function<T,V>
60 affine_combo<value_type,time_type> affine_func;
63 bezier_base():r(0.0),s(1.0) { }
65 const value_type &a, const value_type &b, const value_type &c, const value_type &d,
66 const time_type &r=0.0, const time_type &s=1.0):
67 a(a),b(b),c(c),d(d),r(r),s(s) { sync(); }
74 operator()(time_type t)const
91 void evaluate(time_type t, value_type &f, value_type &df) const
95 value_type p1 = affine_func(
99 value_type p2 = affine_func(
104 f = affine_func(p1,p2,t);
109 void set_rs(time_type new_r, time_type new_s) { r=new_r; s=new_s; }
110 void set_r(time_type new_r) { r=new_r; }
111 void set_s(time_type new_s) { s=new_s; }
112 const time_type &get_r()const { return r; }
113 const time_type &get_s()const { return s; }
114 time_type get_dt()const { return s-r; }
116 bool intersect_hull(const bezier_base<value_type,time_type> &x)const
121 //! Bezier curve intersection function
122 /*! Calculates the time of intersection
123 ** for the calling curve.
125 ** I still have not figured out a good generic
126 ** method of doing this for a bi-infinite
127 ** cubic bezier curve calculated with the DeCasteljau
130 ** One method, although it does not work for the
131 ** entire bi-infinite curve, is to iteratively
132 ** intersect the hulls. However, we would only detect
133 ** intersections that occur between R and S.
135 ** It is entirely possible that a new construct similar
136 ** to the affine combination function will be necessary
137 ** for this to work properly.
139 ** For now, this function is BROKEN. (although it works
140 ** for the floating-point specializations, using newton's method)
142 time_type intersect(const bezier_base<value_type,time_type> &x, time_type near=0.0)const
147 /* subdivide at some time t into 2 separate curves left and right
152 * 1+2 * 0+3*1+3*2+3 l4,r1
158 0.1 2.3 -> 0.1 2 3 4 5.6
160 /* void subdivide(bezier_base *left, bezier_base *right, const time_type &time = (time_type)0.5) const
162 time_type t = (time-r)/(s-r);
167 //1st stage points to keep
172 lt.b = affine_func(a,b,t);
173 temp = affine_func(b,c,t);
174 rt.c = affine_func(c,d,t);
177 lt.c = affine_func(lt.b,temp,t);
178 rt.b = affine_func(temp,rt.c,t);
181 lt.d = rt.a = affine_func(lt.c,rt.b,t);
183 //set the time range for l,r (the inside values should be 1, 0 respectively)
187 //give back the curves
189 if(right) *right = rt;
197 operator[](int i) const
203 // Fast float implementation of a cubic bezier curve
205 class bezier_base<float,float> : public std::unary_function<float,float>
208 typedef float value_type;
209 typedef float time_type;
211 affine_combo<value_type,time_type> affine_func;
215 value_type _coeff[4];
216 time_type drs; // reciprocal of (s-r)
218 bezier_base():r(0.0),s(1.0),drs(1.0) { }
220 const value_type &a, const value_type &b, const value_type &c, const value_type &d,
221 const time_type &r=0.0, const time_type &s=1.0):
222 a(a),b(b),c(c),d(d),r(r),s(s),drs(1.0/(s-r)) { sync(); }
228 _coeff[1]= b*3 - a*3;
229 _coeff[2]= c*3 - b*6 + a*3;
230 _coeff[3]= d - c*3 + b*3 - a;
233 // Cost Summary: 4 products, 3 sums, and 1 difference.
235 operator()(time_type t)const
236 { t-=r; t*=drs; return _coeff[0]+(_coeff[1]+(_coeff[2]+(_coeff[3])*t)*t)*t; }
238 void set_rs(time_type new_r, time_type new_s) { r=new_r; s=new_s; drs=1.0/(s-r); }
239 void set_r(time_type new_r) { r=new_r; drs=1.0/(s-r); }
240 void set_s(time_type new_s) { s=new_s; drs=1.0/(s-r); }
241 const time_type &get_r()const { return r; }
242 const time_type &get_s()const { return s; }
243 time_type get_dt()const { return s-r; }
245 //! Bezier curve intersection function
246 /*! Calculates the time of intersection
247 ** for the calling curve.
249 time_type intersect(const bezier_base<value_type,time_type> &x, time_type t=0.0,int i=15)const
251 //BROKEN - the time values of the 2 curves should be independent
252 value_type system[4];
253 system[0]=_coeff[0]-x._coeff[0];
254 system[1]=_coeff[1]-x._coeff[1];
255 system[2]=_coeff[2]-x._coeff[2];
256 system[3]=_coeff[3]-x._coeff[3];
262 // Inner loop cost summary: 7 products, 5 sums, 1 difference
264 t-= (system[0]+(system[1]+(system[2]+(system[3])*t)*t)*t)/
265 (system[1]+(system[2]*2+(system[3]*3)*t)*t);
278 operator[](int i) const
283 // Fast double implementation of a cubic bezier curve
285 class bezier_base<double,float> : public std::unary_function<float,double>
288 typedef double value_type;
289 typedef float time_type;
291 affine_combo<value_type,time_type> affine_func;
295 value_type _coeff[4];
296 time_type drs; // reciprocal of (s-r)
298 bezier_base():r(0.0),s(1.0),drs(1.0) { }
300 const value_type &a, const value_type &b, const value_type &c, const value_type &d,
301 const time_type &r=0.0, const time_type &s=1.0):
302 a(a),b(b),c(c),d(d),r(r),s(s),drs(1.0/(s-r)) { sync(); }
308 _coeff[1]= b*3 - a*3;
309 _coeff[2]= c*3 - b*6 + a*3;
310 _coeff[3]= d - c*3 + b*3 - a;
313 // 4 products, 3 sums, and 1 difference.
315 operator()(time_type t)const
316 { t-=r; t*=drs; return _coeff[0]+(_coeff[1]+(_coeff[2]+(_coeff[3])*t)*t)*t; }
318 void set_rs(time_type new_r, time_type new_s) { r=new_r; s=new_s; drs=1.0/(s-r); }
319 void set_r(time_type new_r) { r=new_r; drs=1.0/(s-r); }
320 void set_s(time_type new_s) { s=new_s; drs=1.0/(s-r); }
321 const time_type &get_r()const { return r; }
322 const time_type &get_s()const { return s; }
323 time_type get_dt()const { return s-r; }
325 //! Bezier curve intersection function
326 /*! Calculates the time of intersection
327 ** for the calling curve.
329 time_type intersect(const bezier_base<value_type,time_type> &x, time_type t=0.0,int i=15)const
331 //BROKEN - the time values of the 2 curves should be independent
332 value_type system[4];
333 system[0]=_coeff[0]-x._coeff[0];
334 system[1]=_coeff[1]-x._coeff[1];
335 system[2]=_coeff[2]-x._coeff[2];
336 system[3]=_coeff[3]-x._coeff[3];
342 // Inner loop: 7 products, 5 sums, 1 difference
344 t-= (system[0]+(system[1]+(system[2]+(system[3])*t)*t)*t)/
345 (system[1]+(system[2]*2+(system[3]*3)*t)*t);
358 operator[](int i) const
364 // Fast double implementation of a cubic bezier curve
367 template <class T,unsigned int FIXED_BITS>
368 class bezier_base<fixed_base<T,FIXED_BITS> > : std::unary_function<fixed_base<T,FIXED_BITS>,fixed_base<T,FIXED_BITS> >
371 typedef fixed_base<T,FIXED_BITS> value_type;
372 typedef fixed_base<T,FIXED_BITS> time_type;
375 affine_combo<value_type,time_type> affine_func;
379 value_type _coeff[4];
380 time_type drs; // reciprocal of (s-r)
382 bezier_base():r(0.0),s(1.0),drs(1.0) { }
384 const value_type &a, const value_type &b, const value_type &c, const value_type &d,
385 const time_type &r=0, const time_type &s=1):
386 a(a),b(b),c(c),d(d),r(r),s(s),drs(1.0/(s-r)) { sync(); }
390 drs=time_type(1)/(s-r);
392 _coeff[1]= b*3 - a*3;
393 _coeff[2]= c*3 - b*6 + a*3;
394 _coeff[3]= d - c*3 + b*3 - a;
397 // 4 products, 3 sums, and 1 difference.
399 operator()(time_type t)const
400 { t-=r; t*=drs; return _coeff[0]+(_coeff[1]+(_coeff[2]+(_coeff[3])*t)*t)*t; }
402 void set_rs(time_type new_r, time_type new_s) { r=new_r; s=new_s; drs=time_type(1)/(s-r); }
403 void set_r(time_type new_r) { r=new_r; drs=time_type(1)/(s-r); }
404 void set_s(time_type new_s) { s=new_s; drs=time_type(1)/(s-r); }
405 const time_type &get_r()const { return r; }
406 const time_type &get_s()const { return s; }
407 time_type get_dt()const { return s-r; }
409 //! Bezier curve intersection function
410 //! Calculates the time of intersection
411 // for the calling curve.
413 time_type intersect(const bezier_base<value_type,time_type> &x, time_type t=0,int i=15)const
415 value_type system[4];
416 system[0]=_coeff[0]-x._coeff[0];
417 system[1]=_coeff[1]-x._coeff[1];
418 system[2]=_coeff[2]-x._coeff[2];
419 system[3]=_coeff[3]-x._coeff[3];
425 // Inner loop: 7 products, 5 sums, 1 difference
427 t-=(time_type) ( (system[0]+(system[1]+(system[2]+(system[3])*t)*t)*t)/
428 (system[1]+(system[2]*2+(system[3]*3)*t)*t) );
441 operator[](int i) const
451 template <typename V, typename T>
452 class bezier_iterator
456 struct iterator_category {};
457 typedef V value_type;
458 typedef T difference_type;
464 bezier_base<V,T> curve;
470 operator*(void)const { return curve(t); }
471 const surface_iterator&
474 { t+=dt; return &this; }
476 const surface_iterator&
478 { hermite_iterator _tmp=*this; t+=dt; return _tmp; }
480 const surface_iterator&
482 { t-=dt; return &this; }
484 const surface_iterator&
486 { hermite_iterator _tmp=*this; t-=dt; return _tmp; }
490 operator+(difference_type __n) const
491 { return surface_iterator(data+__n[0]+__n[1]*pitch,pitch); }
494 operator-(difference_type __n) const
495 { return surface_iterator(data-__n[0]-__n[1]*pitch,pitch); }
500 template <typename V,typename T=float>
501 class bezier : public bezier_base<V,T>
504 typedef V value_type;
506 typedef float distance_type;
507 typedef bezier_iterator<V,T> iterator;
508 typedef bezier_iterator<V,T> const_iterator;
510 distance_func<value_type> dist;
512 using bezier_base<V,T>::get_r;
513 using bezier_base<V,T>::get_s;
514 using bezier_base<V,T>::get_dt;
518 bezier(const value_type &a, const value_type &b, const value_type &c, const value_type &d):
519 bezier_base<V,T>(a,b,c,d) { }
522 const_iterator begin()const;
523 const_iterator end()const;
525 time_type find_closest(const value_type& x, int i=7, time_type r=(0), time_type s=(1))const
530 if(dist(operator()((s-r)*(1.0/3.0)+r),x) < dist(operator()((s-r)*(2.0/3.0)+r),x))
540 distance_type find_distance(time_type r, time_type s, int steps=7)const
542 const time_type inc((s-r)/steps);
543 distance_type ret(0);
544 value_type last(operator()(r));
546 for(r+=inc;r<s;r+=inc)
548 const value_type n(operator()(r));
549 ret+=dist.uncook(dist(last,n));
552 ret+=dist.uncook(dist(last,operator()(r)))*(s-(r-inc))/inc;
557 distance_type length()const { return find_distance(get_r(),get_s()); }
559 /* subdivide at some time t into 2 separate curves left and right
564 * 1+2 * 0+3*1+3*2+3 l4,r1
570 0.1 2.3 -> 0.1 2 3 4 5.6
572 void subdivide(bezier *left, bezier *right, const time_type &time = (time_type)0.5) const
574 time_type t=(t-get_r())/get_dt();
578 const value_type& a((*this)[0]);
579 const value_type& b((*this)[1]);
580 const value_type& c((*this)[2]);
581 const value_type& d((*this)[3]);
583 //1st stage points to keep
588 lt[1] = affine_func(a,b,t);
589 temp = affine_func(b,c,t);
590 rt[2] = affine_func(c,d,t);
593 lt[2] = affine_func(lt[1],temp,t);
594 rt[1] = affine_func(temp,rt[2],t);
597 lt[3] = rt[0] = affine_func(lt[2],rt[1],t);
599 //set the time range for l,r (the inside values should be 1, 0 respectively)
606 //give back the curves
608 if(right) *right = rt;
612 void evaluate(time_type t, value_type &f, value_type &df) const
614 t=(t-get_r())/get_dt();
616 const value_type& a((*this)[0]);
617 const value_type& b((*this)[1]);
618 const value_type& c((*this)[2]);
619 const value_type& d((*this)[3]);
621 const value_type p1 = affine_func(
625 const value_type p2 = affine_func(
630 f = affine_func(p1,p2,t);
637 /* === E X T E R N S ======================================================= */
639 /* === E N D =============================================================== */