Edit comment text.
[synfig.git] / ETL / trunk / ETL / _bezier.h
1 /*! ========================================================================
2 ** Extended Template Library
3 ** Bezier Template Class Implementation
4 ** $Id$
5 **
6 ** Copyright (c) 2002 Robert B. Quattlebaum Jr.
7 ** Copyright (c) 2007 Chris Moore
8 **
9 ** This package is free software; you can redistribute it and/or
10 ** modify it under the terms of the GNU General Public License as
11 ** published by the Free Software Foundation; either version 2 of
12 ** the License, or (at your option) any later version.
13 **
14 ** This package is distributed in the hope that it will be useful,
15 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
16 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17 ** General Public License for more details.
18 **
19 ** === N O T E S ===========================================================
20 **
21 ** This is an internal header file, included by other ETL headers.
22 ** You should not attempt to use it directly.
23 **
24 ** ========================================================================= */
25
26 /* === S T A R T =========================================================== */
27
28 #ifndef __ETL_BEZIER_H
29 #define __ETL_BEZIER_H
30
31 /* === H E A D E R S ======================================================= */
32
33 #include "_curve_func.h"
34 #include <cmath>                                // for ldexp
35 // #include <ETL/fixed>                 // not used
36
37 /* === M A C R O S ========================================================= */
38
39 #define MAXDEPTH 64     /*  Maximum depth for recursion */
40
41 /* take binary sign of a, either -1, or 1 if >= 0 */
42 #define SGN(a)          (((a)<0) ? -1 : 1)
43
44 /* find minimum of a and b */
45 #ifndef MIN
46 #define MIN(a,b)        (((a)<(b))?(a):(b))
47 #endif
48
49 /* find maximum of a and b */
50 #ifndef MAX
51 #define MAX(a,b)        (((a)>(b))?(a):(b))
52 #endif
53
54 #define BEZIER_EPSILON  (ldexp(1.0,-MAXDEPTH-1)) /*Flatness control value */
55 //#define       BEZIER_EPSILON  0.00005 /*Flatness control value */
56 #define DEGREE  3                       /*  Cubic Bezier curve          */
57 #define W_DEGREE 5                      /*  Degree of eqn to find roots of */
58
59 /* === T Y P E D E F S ===================================================== */
60
61 /* === C L A S S E S & S T R U C T S ======================================= */
62
63 _ETL_BEGIN_NAMESPACE
64
65 template<typename V,typename T> class bezier;
66
67 //! Cubic Bezier Curve Base Class
68 // This generic implementation uses the DeCasteljau algorithm.
69 // Works for just about anything that has an affine combination function
70 template <typename V,typename T=float>
71 class bezier_base : public std::unary_function<T,V>
72 {
73 public:
74         typedef V value_type;
75         typedef T time_type;
76
77 private:
78         value_type a,b,c,d;
79         time_type r,s;
80
81 protected:
82         affine_combo<value_type,time_type> affine_func;
83
84 public:
85         bezier_base():r(0.0),s(1.0) { }
86         bezier_base(
87                 const value_type &a, const value_type &b, const value_type &c, const value_type &d,
88                 const time_type &r=0.0, const time_type &s=1.0):
89                 a(a),b(b),c(c),d(d),r(r),s(s) { sync(); }
90
91         void sync()
92         {
93         }
94
95         value_type
96         operator()(time_type t)const
97         {
98                 t=(t-r)/(s-r);
99                 return
100                 affine_func(
101                         affine_func(
102                                 affine_func(a,b,t),
103                                 affine_func(b,c,t)
104                         ,t),
105                         affine_func(
106                                 affine_func(b,c,t),
107                                 affine_func(c,d,t)
108                         ,t)
109                 ,t);
110         }
111
112         /*
113         void evaluate(time_type t, value_type &f, value_type &df) const
114         {
115                 t=(t-r)/(s-r);
116
117                 value_type p1 = affine_func(
118                                                         affine_func(a,b,t),
119                                                         affine_func(b,c,t)
120                                                         ,t);
121                 value_type p2 = affine_func(
122                                                         affine_func(b,c,t),
123                                                         affine_func(c,d,t)
124                                                 ,t);
125
126                 f = affine_func(p1,p2,t);
127                 df = (p2-p1)*3;
128         }
129         */
130
131         void set_rs(time_type new_r, time_type new_s) { r=new_r; s=new_s; }
132         void set_r(time_type new_r) { r=new_r; }
133         void set_s(time_type new_s) { s=new_s; }
134         const time_type &get_r()const { return r; }
135         const time_type &get_s()const { return s; }
136         time_type get_dt()const { return s-r; }
137
138         bool intersect_hull(const bezier_base<value_type,time_type> &x)const
139         {
140                 return 0;
141         }
142
143         //! Bezier curve intersection function
144         /*! Calculates the time of intersection
145         **      for the calling curve.
146         **
147         **      I still have not figured out a good generic
148         **      method of doing this for a bi-infinite
149         **      cubic bezier curve calculated with the DeCasteljau
150         **      algorithm.
151         **
152         **      One method, although it does not work for the
153         **      entire bi-infinite curve, is to iteratively
154         **      intersect the hulls. However, we would only detect
155         **      intersections that occur between R and S.
156         **
157         **      It is entirely possible that a new construct similar
158         **      to the affine combination function will be necessary
159         **      for this to work properly.
160         **
161         **      For now, this function is BROKEN. (although it works
162         **      for the floating-point specializations, using newton's method)
163         */
164         time_type intersect(const bezier_base<value_type,time_type> &x, time_type near=0.0)const
165         {
166                 return 0;
167         }
168
169         /* subdivide at some time t into 2 separate curves left and right
170
171                 b0 l1
172                 *               0+1 l2
173                 b1              *               1+2*1+2 l3
174                 *               1+2             *                       0+3*1+3*2+3 l4,r1
175                 b2              *               1+2*2+2 r2      *
176                 *               2+3     r3      *
177                 b3 r4   *
178                 *
179
180                 0.1 2.3 ->      0.1 2 3 4 5.6
181         */
182 /*      void subdivide(bezier_base *left, bezier_base *right, const time_type &time = (time_type)0.5) const
183         {
184                 time_type t = (time-r)/(s-r);
185                 bezier_base lt,rt;
186
187                 value_type temp;
188
189                 //1st stage points to keep
190                 lt.a = a;
191                 rt.d = d;
192
193                 //2nd stage calc
194                 lt.b = affine_func(a,b,t);
195                 temp = affine_func(b,c,t);
196                 rt.c = affine_func(c,d,t);
197
198                 //3rd stage calc
199                 lt.c = affine_func(lt.b,temp,t);
200                 rt.b = affine_func(temp,rt.c,t);
201
202                 //last stage calc
203                 lt.d = rt.a = affine_func(lt.c,rt.b,t);
204
205                 //set the time range for l,r (the inside values should be 1, 0 respectively)
206                 lt.r = r;
207                 rt.s = s;
208
209                 //give back the curves
210                 if(left) *left = lt;
211                 if(right) *right = rt;
212         }
213         */
214         value_type &
215         operator[](int i)
216         { return (&a)[i]; }
217
218         const value_type &
219         operator[](int i) const
220         { return (&a)[i]; }
221 };
222
223
224 #if 1
225 // Fast float implementation of a cubic bezier curve
226 template <>
227 class bezier_base<float,float> : public std::unary_function<float,float>
228 {
229 public:
230         typedef float value_type;
231         typedef float time_type;
232 private:
233         affine_combo<value_type,time_type> affine_func;
234         value_type a,b,c,d;
235         time_type r,s;
236
237         value_type _coeff[4];
238         time_type drs; // reciprocal of (s-r)
239 public:
240         bezier_base():r(0.0),s(1.0),drs(1.0) { }
241         bezier_base(
242                 const value_type &a, const value_type &b, const value_type &c, const value_type &d,
243                 const time_type &r=0.0, const time_type &s=1.0):
244                 a(a),b(b),c(c),d(d),r(r),s(s),drs(1.0/(s-r)) { sync(); }
245
246         void sync()
247         {
248 //              drs=1.0/(s-r);
249                 _coeff[0]=                 a;
250                 _coeff[1]=           b*3 - a*3;
251                 _coeff[2]=     c*3 - b*6 + a*3;
252                 _coeff[3]= d - c*3 + b*3 - a;
253         }
254
255         // Cost Summary: 4 products, 3 sums, and 1 difference.
256         inline value_type
257         operator()(time_type t)const
258         { t-=r; t*=drs; return _coeff[0]+(_coeff[1]+(_coeff[2]+(_coeff[3])*t)*t)*t; }
259
260         void set_rs(time_type new_r, time_type new_s) { r=new_r; s=new_s; drs=1.0/(s-r); }
261         void set_r(time_type new_r) { r=new_r; drs=1.0/(s-r); }
262         void set_s(time_type new_s) { s=new_s; drs=1.0/(s-r); }
263         const time_type &get_r()const { return r; }
264         const time_type &get_s()const { return s; }
265         time_type get_dt()const { return s-r; }
266
267         //! Bezier curve intersection function
268         /*! Calculates the time of intersection
269         **      for the calling curve.
270         */
271         time_type intersect(const bezier_base<value_type,time_type> &x, time_type t=0.0,int i=15)const
272         {
273                 //BROKEN - the time values of the 2 curves should be independent
274                 value_type system[4];
275                 system[0]=_coeff[0]-x._coeff[0];
276                 system[1]=_coeff[1]-x._coeff[1];
277                 system[2]=_coeff[2]-x._coeff[2];
278                 system[3]=_coeff[3]-x._coeff[3];
279
280                 t-=r;
281                 t*=drs;
282
283                 // Newton's method
284                 // Inner loop cost summary: 7 products, 5 sums, 1 difference
285                 for(;i;i--)
286                         t-= (system[0]+(system[1]+(system[2]+(system[3])*t)*t)*t)/
287                                 (system[1]+(system[2]*2+(system[3]*3)*t)*t);
288
289                 t*=(s-r);
290                 t+=r;
291
292                 return t;
293         }
294
295         value_type &
296         operator[](int i)
297         { return (&a)[i]; }
298
299         const value_type &
300         operator[](int i) const
301         { return (&a)[i]; }
302 };
303
304
305 // Fast double implementation of a cubic bezier curve
306 template <>
307 class bezier_base<double,float> : public std::unary_function<float,double>
308 {
309 public:
310         typedef double value_type;
311         typedef float time_type;
312 private:
313         affine_combo<value_type,time_type> affine_func;
314         value_type a,b,c,d;
315         time_type r,s;
316
317         value_type _coeff[4];
318         time_type drs; // reciprocal of (s-r)
319 public:
320         bezier_base():r(0.0),s(1.0),drs(1.0) { }
321         bezier_base(
322                 const value_type &a, const value_type &b, const value_type &c, const value_type &d,
323                 const time_type &r=0.0, const time_type &s=1.0):
324                 a(a),b(b),c(c),d(d),r(r),s(s),drs(1.0/(s-r)) { sync(); }
325
326         void sync()
327         {
328 //              drs=1.0/(s-r);
329                 _coeff[0]=                 a;
330                 _coeff[1]=           b*3 - a*3;
331                 _coeff[2]=     c*3 - b*6 + a*3;
332                 _coeff[3]= d - c*3 + b*3 - a;
333         }
334
335         // 4 products, 3 sums, and 1 difference.
336         inline value_type
337         operator()(time_type t)const
338         { t-=r; t*=drs; return _coeff[0]+(_coeff[1]+(_coeff[2]+(_coeff[3])*t)*t)*t; }
339
340         void set_rs(time_type new_r, time_type new_s) { r=new_r; s=new_s; drs=1.0/(s-r); }
341         void set_r(time_type new_r) { r=new_r; drs=1.0/(s-r); }
342         void set_s(time_type new_s) { s=new_s; drs=1.0/(s-r); }
343         const time_type &get_r()const { return r; }
344         const time_type &get_s()const { return s; }
345         time_type get_dt()const { return s-r; }
346
347         //! Bezier curve intersection function
348         /*! Calculates the time of intersection
349         **      for the calling curve.
350         */
351         time_type intersect(const bezier_base<value_type,time_type> &x, time_type t=0.0,int i=15)const
352         {
353                 //BROKEN - the time values of the 2 curves should be independent
354                 value_type system[4];
355                 system[0]=_coeff[0]-x._coeff[0];
356                 system[1]=_coeff[1]-x._coeff[1];
357                 system[2]=_coeff[2]-x._coeff[2];
358                 system[3]=_coeff[3]-x._coeff[3];
359
360                 t-=r;
361                 t*=drs;
362
363                 // Newton's method
364                 // Inner loop: 7 products, 5 sums, 1 difference
365                 for(;i;i--)
366                         t-= (system[0]+(system[1]+(system[2]+(system[3])*t)*t)*t)/
367                                 (system[1]+(system[2]*2+(system[3]*3)*t)*t);
368
369                 t*=(s-r);
370                 t+=r;
371
372                 return t;
373         }
374
375         value_type &
376         operator[](int i)
377         { return (&a)[i]; }
378
379         const value_type &
380         operator[](int i) const
381         { return (&a)[i]; }
382 };
383
384 //#ifdef __FIXED__
385
386 // Fast double implementation of a cubic bezier curve
387 /*
388 template <>
389 template <class T,unsigned int FIXED_BITS>
390 class bezier_base<fixed_base<T,FIXED_BITS> > : std::unary_function<fixed_base<T,FIXED_BITS>,fixed_base<T,FIXED_BITS> >
391 {
392 public:
393         typedef fixed_base<T,FIXED_BITS> value_type;
394         typedef fixed_base<T,FIXED_BITS> time_type;
395
396 private:
397         affine_combo<value_type,time_type> affine_func;
398         value_type a,b,c,d;
399         time_type r,s;
400
401         value_type _coeff[4];
402         time_type drs; // reciprocal of (s-r)
403 public:
404         bezier_base():r(0.0),s(1.0),drs(1.0) { }
405         bezier_base(
406                 const value_type &a, const value_type &b, const value_type &c, const value_type &d,
407                 const time_type &r=0, const time_type &s=1):
408                 a(a),b(b),c(c),d(d),r(r),s(s),drs(1.0/(s-r)) { sync(); }
409
410         void sync()
411         {
412                 drs=time_type(1)/(s-r);
413                 _coeff[0]=                 a;
414                 _coeff[1]=           b*3 - a*3;
415                 _coeff[2]=     c*3 - b*6 + a*3;
416                 _coeff[3]= d - c*3 + b*3 - a;
417         }
418
419         // 4 products, 3 sums, and 1 difference.
420         inline value_type
421         operator()(time_type t)const
422         { t-=r; t*=drs; return _coeff[0]+(_coeff[1]+(_coeff[2]+(_coeff[3])*t)*t)*t; }
423
424         void set_rs(time_type new_r, time_type new_s) { r=new_r; s=new_s; drs=time_type(1)/(s-r); }
425         void set_r(time_type new_r) { r=new_r; drs=time_type(1)/(s-r); }
426         void set_s(time_type new_s) { s=new_s; drs=time_type(1)/(s-r); }
427         const time_type &get_r()const { return r; }
428         const time_type &get_s()const { return s; }
429         time_type get_dt()const { return s-r; }
430
431         //! Bezier curve intersection function
432         //! Calculates the time of intersection
433         //      for the calling curve.
434         //
435         time_type intersect(const bezier_base<value_type,time_type> &x, time_type t=0,int i=15)const
436         {
437                 value_type system[4];
438                 system[0]=_coeff[0]-x._coeff[0];
439                 system[1]=_coeff[1]-x._coeff[1];
440                 system[2]=_coeff[2]-x._coeff[2];
441                 system[3]=_coeff[3]-x._coeff[3];
442
443                 t-=r;
444                 t*=drs;
445
446                 // Newton's method
447                 // Inner loop: 7 products, 5 sums, 1 difference
448                 for(;i;i--)
449                         t-=(time_type) ( (system[0]+(system[1]+(system[2]+(system[3])*t)*t)*t)/
450                                 (system[1]+(system[2]*2+(system[3]*3)*t)*t) );
451
452                 t*=(s-r);
453                 t+=r;
454
455                 return t;
456         }
457
458         value_type &
459         operator[](int i)
460         { return (&a)[i]; }
461
462         const value_type &
463         operator[](int i) const
464         { return (&a)[i]; }
465 };
466 */
467 //#endif
468
469 #endif
470
471
472
473 template <typename V, typename T>
474 class bezier_iterator
475 {
476 public:
477
478         struct iterator_category {};
479         typedef V value_type;
480         typedef T difference_type;
481         typedef V reference;
482
483 private:
484         difference_type t;
485         difference_type dt;
486         bezier_base<V,T>        curve;
487
488 public:
489
490 /*
491         reference
492         operator*(void)const { return curve(t); }
493         const surface_iterator&
494
495         operator++(void)
496         { t+=dt; return &this; }
497
498         const surface_iterator&
499         operator++(int)
500         { hermite_iterator _tmp=*this; t+=dt; return _tmp; }
501
502         const surface_iterator&
503         operator--(void)
504         { t-=dt; return &this; }
505
506         const surface_iterator&
507         operator--(int)
508         { hermite_iterator _tmp=*this; t-=dt; return _tmp; }
509
510
511         surface_iterator
512         operator+(difference_type __n) const
513         { return surface_iterator(data+__n[0]+__n[1]*pitch,pitch); }
514
515         surface_iterator
516         operator-(difference_type __n) const
517         { return surface_iterator(data-__n[0]-__n[1]*pitch,pitch); }
518 */
519
520 };
521
522 template <typename V,typename T=float>
523 class bezier : public bezier_base<V,T>
524 {
525 public:
526         typedef V value_type;
527         typedef T time_type;
528         typedef float distance_type;
529         typedef bezier_iterator<V,T> iterator;
530         typedef bezier_iterator<V,T> const_iterator;
531
532         distance_func<value_type> dist;
533
534         using bezier_base<V,T>::get_r;
535         using bezier_base<V,T>::get_s;
536         using bezier_base<V,T>::get_dt;
537
538 public:
539         bezier() { }
540         bezier(const value_type &a, const value_type &b, const value_type &c, const value_type &d):
541                 bezier_base<V,T>(a,b,c,d) { }
542
543
544         const_iterator begin()const;
545         const_iterator end()const;
546
547         time_type find_closest(bool fast, const value_type& x, int i=7)const
548         {
549             if (!fast)
550             {
551                         value_type array[4] = {
552                                 bezier<V,T>::operator[](0),
553                                 bezier<V,T>::operator[](1),
554                                 bezier<V,T>::operator[](2),
555                                 bezier<V,T>::operator[](3)};
556                         return NearestPointOnCurve(x, array);
557             }
558             else
559             {
560                         time_type r(0), s(1);
561                         float t((r+s)*0.5); /* half way between r and s */
562
563                         for(;i;i--)
564                         {
565                                 // compare 33% of the way between r and s with 67% of the way between r and s
566                                 if(dist(operator()((s-r)*(1.0/3.0)+r), x) <
567                                    dist(operator()((s-r)*(2.0/3.0)+r), x))
568                                         s=t;
569                                 else
570                                         r=t;
571                                 t=((r+s)*0.5);
572                         }
573                         return t;
574                 }
575         }
576
577         distance_type find_distance(time_type r, time_type s, int steps=7)const
578         {
579                 const time_type inc((s-r)/steps);
580                 if (!inc) return 0;
581                 distance_type ret(0);
582                 value_type last(operator()(r));
583
584                 for(r+=inc;r<s;r+=inc)
585                 {
586                         const value_type n(operator()(r));
587                         ret+=dist.uncook(dist(last,n));
588                         last=n;
589                 }
590                 ret+=dist.uncook(dist(last,operator()(r)))*(s-(r-inc))/inc;
591
592                 return ret;
593         }
594
595         distance_type length()const { return find_distance(get_r(),get_s()); }
596
597         /* subdivide at some time t into 2 separate curves left and right
598
599                 b0 l1
600                 *               0+1 l2
601                 b1              *               1+2*1+2 l3
602                 *               1+2             *                       0+3*1+3*2+3 l4,r1
603                 b2              *               1+2*2+2 r2      *
604                 *               2+3     r3      *
605                 b3 r4   *
606                 *
607
608                 0.1 2.3 ->      0.1 2 3 4 5.6
609         */
610         void subdivide(bezier *left, bezier *right, const time_type &time = (time_type)0.5) const
611         {
612                 time_type t=(time-get_r())/get_dt();
613                 bezier lt,rt;
614
615                 value_type temp;
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]);
620
621                 //1st stage points to keep
622                 lt[0] = a;
623                 rt[3] = d;
624
625                 //2nd stage calc
626                 lt[1] = affine_func(a,b,t);
627                 temp = affine_func(b,c,t);
628                 rt[2] = affine_func(c,d,t);
629
630                 //3rd stage calc
631                 lt[2] = affine_func(lt[1],temp,t);
632                 rt[1] = affine_func(temp,rt[2],t);
633
634                 //last stage calc
635                 lt[3] = rt[0] = affine_func(lt[2],rt[1],t);
636
637                 //set the time range for l,r (the inside values should be 1, 0 respectively)
638                 lt.set_r(get_r());
639                 rt.set_s(get_s());
640
641                 lt.sync();
642                 rt.sync();
643
644                 //give back the curves
645                 if(left) *left = lt;
646                 if(right) *right = rt;
647         }
648
649
650         void evaluate(time_type t, value_type &f, value_type &df) const
651         {
652                 t=(t-get_r())/get_dt();
653
654                 const value_type& a((*this)[0]);
655                 const value_type& b((*this)[1]);
656                 const value_type& c((*this)[2]);
657                 const value_type& d((*this)[3]);
658
659                 const value_type p1 = affine_func(
660                                                         affine_func(a,b,t),
661                                                         affine_func(b,c,t)
662                                                         ,t);
663                 const value_type p2 = affine_func(
664                                                         affine_func(b,c,t),
665                                                         affine_func(c,d,t)
666                                                 ,t);
667
668                 f = affine_func(p1,p2,t);
669                 df = (p2-p1)*3;
670         }
671
672 private:
673         /*
674          *  Bezier :
675          *      Evaluate a Bezier curve at a particular parameter value
676          *      Fill in control points for resulting sub-curves if "Left" and
677          *      "Right" are non-null.
678          *
679          *    int                       degree;         Degree of bezier curve
680          *    value_type        *VT;            Control pts
681          *    time_type         t;                      Parameter value
682          *    value_type        *Left;          RETURN left half ctl pts
683          *    value_type        *Right;         RETURN right half ctl pts
684          */
685         static value_type Bezier(value_type *VT, int degree, time_type t, value_type *Left, value_type *Right)
686         {
687                 int             i, j;           /* Index variables      */
688                 value_type      Vtemp[W_DEGREE+1][W_DEGREE+1];
689
690                 /* Copy control points  */
691                 for (j = 0; j <= degree; j++)
692                         Vtemp[0][j] = VT[j];
693
694                 /* Triangle computation */
695                 for (i = 1; i <= degree; i++)
696                         for (j =0 ; j <= degree - i; j++)
697                         {
698                                 Vtemp[i][j][0] = (1.0 - t) * Vtemp[i-1][j][0] + t * Vtemp[i-1][j+1][0];
699                                 Vtemp[i][j][1] = (1.0 - t) * Vtemp[i-1][j][1] + t * Vtemp[i-1][j+1][1];
700                         }
701
702                 if (Left != NULL)
703                         for (j = 0; j <= degree; j++)
704                                 Left[j]  = Vtemp[j][0];
705
706                 if (Right != NULL)
707                         for (j = 0; j <= degree; j++)
708                                 Right[j] = Vtemp[degree-j][j];
709
710                 return (Vtemp[degree][0]);
711         }
712
713         /*
714          * CrossingCount :
715          *      Count the number of times a Bezier control polygon
716          *      crosses the 0-axis. This number is >= the number of roots.
717          *
718          *    value_type        *VT;                    Control pts of Bezier curve
719          */
720         static int CrossingCount(value_type *VT)
721         {
722                 int     i;
723                 int     n_crossings = 0;        /*  Number of zero-crossings    */
724                 int             sign, old_sign;         /*  Sign of coefficients                */
725
726                 sign = old_sign = SGN(VT[0][1]);
727                 for (i = 1; i <= W_DEGREE; i++)
728                 {
729                         sign = SGN(VT[i][1]);
730                         if (sign != old_sign) n_crossings++;
731                         old_sign = sign;
732                 }
733
734                 return n_crossings;
735         }
736
737         /*
738          *  ControlPolygonFlatEnough :
739          *      Check if the control polygon of a Bezier curve is flat enough
740          *      for recursive subdivision to bottom out.
741          *
742          *    value_type        *VT;            Control points
743          */
744         static int ControlPolygonFlatEnough(value_type *VT)
745         {
746                 int                     i;                                      /* Index variable                                       */
747                 distance_type   distance[W_DEGREE];     /* Distances from pts to line           */
748                 distance_type   max_distance_above;     /* maximum of these                                     */
749                 distance_type   max_distance_below;
750                 time_type               intercept_1, intercept_2, left_intercept, right_intercept;
751                 distance_type   a, b, c;                        /* Coefficients of implicit                     */
752                 /* eqn for line from VT[0]-VT[deg]                      */
753                 /* Find the  perpendicular distance                     */
754                 /* from each interior control point to          */
755                 /* line connecting VT[0] and VT[W_DEGREE]       */
756                 {
757                         distance_type   abSquared;
758
759                         /* Derive the implicit equation for line connecting first *
760                          *  and last control points */
761                         a = VT[0][1] - VT[W_DEGREE][1];
762                         b = VT[W_DEGREE][0] - VT[0][0];
763                         c = VT[0][0] * VT[W_DEGREE][1] - VT[W_DEGREE][0] * VT[0][1];
764
765                         abSquared = (a * a) + (b * b);
766
767                         for (i = 1; i < W_DEGREE; i++)
768                         {
769                                 /* Compute distance from each of the points to that line        */
770                                 distance[i] = a * VT[i][0] + b * VT[i][1] + c;
771                                 if (distance[i] > 0.0) distance[i] =  (distance[i] * distance[i]) / abSquared;
772                                 if (distance[i] < 0.0) distance[i] = -(distance[i] * distance[i]) / abSquared;
773                         }
774                 }
775
776                 /* Find the largest distance */
777                 max_distance_above = max_distance_below = 0.0;
778
779                 for (i = 1; i < W_DEGREE; i++)
780                 {
781                         if (distance[i] < 0.0) max_distance_below = MIN(max_distance_below, distance[i]);
782                         if (distance[i] > 0.0) max_distance_above = MAX(max_distance_above, distance[i]);
783                 }
784
785                 /* Implicit equation for "above" line */
786                 intercept_1 = -(c + max_distance_above)/a;
787
788                 /*  Implicit equation for "below" line */
789                 intercept_2 = -(c + max_distance_below)/a;
790
791                 /* Compute intercepts of bounding box   */
792                 left_intercept = MIN(intercept_1, intercept_2);
793                 right_intercept = MAX(intercept_1, intercept_2);
794
795                 return 0.5 * (right_intercept-left_intercept) < BEZIER_EPSILON ? 1 : 0;
796         }
797
798         /*
799          *  ComputeXIntercept :
800          *      Compute intersection of chord from first control point to last
801          *  with 0-axis.
802          *
803          *    value_type        *VT;                    Control points
804          */
805         static time_type ComputeXIntercept(value_type *VT)
806         {
807                 distance_type YNM = VT[W_DEGREE][1] - VT[0][1];
808                 return (YNM*VT[0][0] - (VT[W_DEGREE][0] - VT[0][0])*VT[0][1]) / YNM;
809         }
810
811         /*
812          *  FindRoots :
813          *      Given a 5th-degree equation in Bernstein-Bezier form, find
814          *      all of the roots in the interval [0, 1].  Return the number
815          *      of roots found.
816          *
817          *    value_type        *w;                             The control points
818          *    time_type         *t;                             RETURN candidate t-values
819          *    int                       depth;                  The depth of the recursion
820          */
821         static int FindRoots(value_type *w, time_type *t, int depth)
822         {
823                 int             i;
824                 value_type      Left[W_DEGREE+1];       /* New left and right   */
825                 value_type      Right[W_DEGREE+1];      /* control polygons             */
826                 int             left_count;                     /* Solution count from  */
827                 int                     right_count;            /* children                             */
828                 time_type       left_t[W_DEGREE+1];     /* Solutions from kids  */
829                 time_type       right_t[W_DEGREE+1];
830
831                 switch (CrossingCount(w))
832                 {
833                         case 0 :
834                         {       /* No solutions here    */
835                                 return 0;
836                         }
837                         case 1 :
838                         {       /* Unique solution      */
839                                 /* Stop recursion when the tree is deep enough          */
840                                 /* if deep enough, return 1 solution at midpoint        */
841                                 if (depth >= MAXDEPTH)
842                                 {
843                                         t[0] = (w[0][0] + w[W_DEGREE][0]) / 2.0;
844                                         return 1;
845                                 }
846                                 if (ControlPolygonFlatEnough(w))
847                                 {
848                                         t[0] = ComputeXIntercept(w);
849                                         return 1;
850                                 }
851                                 break;
852                         }
853                 }
854
855                 /* Otherwise, solve recursively after   */
856                 /* subdividing control polygon                  */
857                 Bezier(w, W_DEGREE, 0.5, Left, Right);
858                 left_count  = FindRoots(Left,  left_t,  depth+1);
859                 right_count = FindRoots(Right, right_t, depth+1);
860
861                 /* Gather solutions together    */
862                 for (i = 0; i < left_count;  i++) t[i] = left_t[i];
863                 for (i = 0; i < right_count; i++) t[i+left_count] = right_t[i];
864
865                 /* Send back total number of solutions  */
866                 return (left_count+right_count);
867         }
868
869         /*
870          *  ConvertToBezierForm :
871          *              Given a point and a Bezier curve, generate a 5th-degree
872          *              Bezier-format equation whose solution finds the point on the
873          *      curve nearest the user-defined point.
874          *
875          *    value_type&       P;                              The point to find t for
876          *    value_type        *VT;                    The control points
877          */
878         static void ConvertToBezierForm(const value_type& P, value_type *VT, value_type w[W_DEGREE+1])
879         {
880                 int     i, j, k, m, n, ub, lb;
881                 int     row, column;                            /* Table indices                                */
882                 value_type      c[DEGREE+1];                    /* VT(i)'s - P                                  */
883                 value_type      d[DEGREE];                              /* VT(i+1) - VT(i)                              */
884                 distance_type   cdTable[3][4];          /* Dot product of c, d                  */
885                 static distance_type z[3][4] = {        /* Precomputed "z" for cubics   */
886                         {1.0, 0.6, 0.3, 0.1},
887                         {0.4, 0.6, 0.6, 0.4},
888                         {0.1, 0.3, 0.6, 1.0}};
889
890                 /* Determine the c's -- these are vectors created by subtracting */
891                 /* point P from each of the control points                                               */
892                 for (i = 0; i <= DEGREE; i++)
893                         c[i] = VT[i] - P;
894
895                 /* Determine the d's -- these are vectors created by subtracting */
896                 /* each control point from the next                                                              */
897                 for (i = 0; i <= DEGREE - 1; i++)
898                         d[i] = (VT[i+1] - VT[i]) * 3.0;
899
900                 /* Create the c,d table -- this is a table of dot products of the */
901                 /* c's and d's                                                                                                    */
902                 for (row = 0; row <= DEGREE - 1; row++)
903                         for (column = 0; column <= DEGREE; column++)
904                                 cdTable[row][column] = d[row] * c[column];
905
906                 /* Now, apply the z's to the dot products, on the skew diagonal */
907                 /* Also, set up the x-values, making these "points"                             */
908                 for (i = 0; i <= W_DEGREE; i++)
909                 {
910                         w[i][0] = (distance_type)(i) / W_DEGREE;
911                         w[i][1] = 0.0;
912                 }
913
914                 n = DEGREE;
915                 m = DEGREE-1;
916                 for (k = 0; k <= n + m; k++)
917                 {
918                         lb = MAX(0, k - m);
919                         ub = MIN(k, n);
920                         for (i = lb; i <= ub; i++)
921                         {
922                                 j = k - i;
923                                 w[i+j][1] += cdTable[j][i] * z[j][i];
924                         }
925                 }
926         }
927
928         /*
929          *  NearestPointOnCurve :
930          *      Compute the parameter value of the point on a Bezier
931          *              curve segment closest to some arbitrary, user-input point.
932          *              Return the point on the curve at that parameter value.
933          *
934          *    value_type&       P;                      The user-supplied point
935          *    value_type        *VT;            Control points of cubic Bezier
936          */
937         static time_type NearestPointOnCurve(const value_type& P, value_type VT[4])
938         {
939                 value_type      w[W_DEGREE+1];                  /* Ctl pts of 5th-degree curve  */
940                 time_type       t_candidate[W_DEGREE];  /* Possible roots                                */
941                 int             n_solutions;                    /* Number of roots found                 */
942                 time_type       t;                                              /* Parameter value of closest pt */
943
944                 /*  Convert problem to 5th-degree Bezier form */
945                 ConvertToBezierForm(P, VT, w);
946
947                 /* Find all possible roots of 5th-degree equation */
948                 n_solutions = FindRoots(w, t_candidate, 0);
949
950                 /* Compare distances of P to all candidates, and to t=0, and t=1 */
951                 {
952                         distance_type   dist, new_dist;
953                         value_type              p, v;
954                         int                             i;
955
956                         /* Check distance to beginning of curve, where t = 0    */
957                         dist = (P - VT[0]).mag_squared();
958                         t = 0.0;
959
960                         /* Find distances for candidate points  */
961                         for (i = 0; i < n_solutions; i++)
962                         {
963                                 p = Bezier(VT, DEGREE, t_candidate[i], (value_type *)NULL, (value_type *)NULL);
964                                 new_dist = (P - p).mag_squared();
965                                 if (new_dist < dist)
966                                 {
967                                         dist = new_dist;
968                                         t = t_candidate[i];
969                                 }
970                         }
971
972                         /* Finally, look at distance to end point, where t = 1.0 */
973                         new_dist = (P - VT[DEGREE]).mag_squared();
974                         if (new_dist < dist)
975                         {
976                                 dist = new_dist;
977                                 t = 1.0;
978                         }
979                 }
980
981                 /*  Return the point on the curve at parameter value t */
982                 return t;
983         }
984 };
985
986 _ETL_END_NAMESPACE
987
988 /* === E X T E R N S ======================================================= */
989
990 /* === E N D =============================================================== */
991
992 #endif