1 /*! ========================================================================
2 ** Extended Template Library
3 ** Bezier Template Class Implementation
6 ** Copyright (c) 2002 Robert B. Quattlebaum Jr.
7 ** Copyright (c) 2007 Chris Moore
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.
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.
19 ** === N O T E S ===========================================================
21 ** This is an internal header file, included by other ETL headers.
22 ** You should not attempt to use it directly.
24 ** ========================================================================= */
26 /* === S T A R T =========================================================== */
28 #ifndef __ETL_BEZIER_H
29 #define __ETL_BEZIER_H
31 /* === H E A D E R S ======================================================= */
33 #include "_curve_func.h"
36 /* === M A C R O S ========================================================= */
38 #define MAXDEPTH 64 /* Maximum depth for recursion */
40 /* take binary sign of a, either -1, or 1 if >= 0 */
41 #define SGN(a) (((a)<0) ? -1 : 1)
43 /* find minimum of a and b */
45 #define MIN(a,b) (((a)<(b))?(a):(b))
48 /* find maximum of a and b */
50 #define MAX(a,b) (((a)>(b))?(a):(b))
53 #define BEZIER_EPSILON (ldexp(1.0,-MAXDEPTH-1)) /*Flatness control value */
54 //#define BEZIER_EPSILON 0.00005 /*Flatness control value */
55 #define DEGREE 3 /* Cubic Bezier curve */
56 #define W_DEGREE 5 /* Degree of eqn to find roots of */
58 /* === T Y P E D E F S ===================================================== */
60 /* === C L A S S E S & S T R U C T S ======================================= */
64 template<typename V,typename T> class bezier;
66 //! Cubic Bezier Curve Base Class
67 // This generic implementation uses the DeCasteljau algorithm.
68 // Works for just about anything that has an affine combination function
69 template <typename V,typename T=float>
70 class bezier_base : public std::unary_function<T,V>
81 affine_combo<value_type,time_type> affine_func;
84 bezier_base():r(0.0),s(1.0) { }
86 const value_type &a, const value_type &b, const value_type &c, const value_type &d,
87 const time_type &r=0.0, const time_type &s=1.0):
88 a(a),b(b),c(c),d(d),r(r),s(s) { sync(); }
95 operator()(time_type t)const
112 void evaluate(time_type t, value_type &f, value_type &df) const
116 value_type p1 = affine_func(
120 value_type p2 = affine_func(
125 f = affine_func(p1,p2,t);
130 void set_rs(time_type new_r, time_type new_s) { r=new_r; s=new_s; }
131 void set_r(time_type new_r) { r=new_r; }
132 void set_s(time_type new_s) { s=new_s; }
133 const time_type &get_r()const { return r; }
134 const time_type &get_s()const { return s; }
135 time_type get_dt()const { return s-r; }
137 bool intersect_hull(const bezier_base<value_type,time_type> &x)const
142 //! Bezier curve intersection function
143 /*! Calculates the time of intersection
144 ** for the calling curve.
146 ** I still have not figured out a good generic
147 ** method of doing this for a bi-infinite
148 ** cubic bezier curve calculated with the DeCasteljau
151 ** One method, although it does not work for the
152 ** entire bi-infinite curve, is to iteratively
153 ** intersect the hulls. However, we would only detect
154 ** intersections that occur between R and S.
156 ** It is entirely possible that a new construct similar
157 ** to the affine combination function will be necessary
158 ** for this to work properly.
160 ** For now, this function is BROKEN. (although it works
161 ** for the floating-point specializations, using newton's method)
163 time_type intersect(const bezier_base<value_type,time_type> &x, time_type near=0.0)const
168 /* subdivide at some time t into 2 separate curves left and right
173 * 1+2 * 0+3*1+3*2+3 l4,r1
179 0.1 2.3 -> 0.1 2 3 4 5.6
181 /* void subdivide(bezier_base *left, bezier_base *right, const time_type &time = (time_type)0.5) const
183 time_type t = (time-r)/(s-r);
188 //1st stage points to keep
193 lt.b = affine_func(a,b,t);
194 temp = affine_func(b,c,t);
195 rt.c = affine_func(c,d,t);
198 lt.c = affine_func(lt.b,temp,t);
199 rt.b = affine_func(temp,rt.c,t);
202 lt.d = rt.a = affine_func(lt.c,rt.b,t);
204 //set the time range for l,r (the inside values should be 1, 0 respectively)
208 //give back the curves
210 if(right) *right = rt;
218 operator[](int i) const
224 // Fast float implementation of a cubic bezier curve
226 class bezier_base<float,float> : public std::unary_function<float,float>
229 typedef float value_type;
230 typedef float time_type;
232 affine_combo<value_type,time_type> affine_func;
236 value_type _coeff[4];
237 time_type drs; // reciprocal of (s-r)
239 bezier_base():r(0.0),s(1.0),drs(1.0) { }
241 const value_type &a, const value_type &b, const value_type &c, const value_type &d,
242 const time_type &r=0.0, const time_type &s=1.0):
243 a(a),b(b),c(c),d(d),r(r),s(s),drs(1.0/(s-r)) { sync(); }
249 _coeff[1]= b*3 - a*3;
250 _coeff[2]= c*3 - b*6 + a*3;
251 _coeff[3]= d - c*3 + b*3 - a;
254 // Cost Summary: 4 products, 3 sums, and 1 difference.
256 operator()(time_type t)const
257 { t-=r; t*=drs; return _coeff[0]+(_coeff[1]+(_coeff[2]+(_coeff[3])*t)*t)*t; }
259 void set_rs(time_type new_r, time_type new_s) { r=new_r; s=new_s; drs=1.0/(s-r); }
260 void set_r(time_type new_r) { r=new_r; drs=1.0/(s-r); }
261 void set_s(time_type new_s) { s=new_s; drs=1.0/(s-r); }
262 const time_type &get_r()const { return r; }
263 const time_type &get_s()const { return s; }
264 time_type get_dt()const { return s-r; }
266 //! Bezier curve intersection function
267 /*! Calculates the time of intersection
268 ** for the calling curve.
270 time_type intersect(const bezier_base<value_type,time_type> &x, time_type t=0.0,int i=15)const
272 //BROKEN - the time values of the 2 curves should be independent
273 value_type system[4];
274 system[0]=_coeff[0]-x._coeff[0];
275 system[1]=_coeff[1]-x._coeff[1];
276 system[2]=_coeff[2]-x._coeff[2];
277 system[3]=_coeff[3]-x._coeff[3];
283 // Inner loop cost summary: 7 products, 5 sums, 1 difference
285 t-= (system[0]+(system[1]+(system[2]+(system[3])*t)*t)*t)/
286 (system[1]+(system[2]*2+(system[3]*3)*t)*t);
299 operator[](int i) const
304 // Fast double implementation of a cubic bezier curve
306 class bezier_base<double,float> : public std::unary_function<float,double>
309 typedef double value_type;
310 typedef float time_type;
312 affine_combo<value_type,time_type> affine_func;
316 value_type _coeff[4];
317 time_type drs; // reciprocal of (s-r)
319 bezier_base():r(0.0),s(1.0),drs(1.0) { }
321 const value_type &a, const value_type &b, const value_type &c, const value_type &d,
322 const time_type &r=0.0, const time_type &s=1.0):
323 a(a),b(b),c(c),d(d),r(r),s(s),drs(1.0/(s-r)) { sync(); }
329 _coeff[1]= b*3 - a*3;
330 _coeff[2]= c*3 - b*6 + a*3;
331 _coeff[3]= d - c*3 + b*3 - a;
334 // 4 products, 3 sums, and 1 difference.
336 operator()(time_type t)const
337 { t-=r; t*=drs; return _coeff[0]+(_coeff[1]+(_coeff[2]+(_coeff[3])*t)*t)*t; }
339 void set_rs(time_type new_r, time_type new_s) { r=new_r; s=new_s; drs=1.0/(s-r); }
340 void set_r(time_type new_r) { r=new_r; drs=1.0/(s-r); }
341 void set_s(time_type new_s) { s=new_s; drs=1.0/(s-r); }
342 const time_type &get_r()const { return r; }
343 const time_type &get_s()const { return s; }
344 time_type get_dt()const { return s-r; }
346 //! Bezier curve intersection function
347 /*! Calculates the time of intersection
348 ** for the calling curve.
350 time_type intersect(const bezier_base<value_type,time_type> &x, time_type t=0.0,int i=15)const
352 //BROKEN - the time values of the 2 curves should be independent
353 value_type system[4];
354 system[0]=_coeff[0]-x._coeff[0];
355 system[1]=_coeff[1]-x._coeff[1];
356 system[2]=_coeff[2]-x._coeff[2];
357 system[3]=_coeff[3]-x._coeff[3];
363 // Inner loop: 7 products, 5 sums, 1 difference
365 t-= (system[0]+(system[1]+(system[2]+(system[3])*t)*t)*t)/
366 (system[1]+(system[2]*2+(system[3]*3)*t)*t);
379 operator[](int i) const
385 // Fast double implementation of a cubic bezier curve
388 template <class T,unsigned int FIXED_BITS>
389 class bezier_base<fixed_base<T,FIXED_BITS> > : std::unary_function<fixed_base<T,FIXED_BITS>,fixed_base<T,FIXED_BITS> >
392 typedef fixed_base<T,FIXED_BITS> value_type;
393 typedef fixed_base<T,FIXED_BITS> time_type;
396 affine_combo<value_type,time_type> affine_func;
400 value_type _coeff[4];
401 time_type drs; // reciprocal of (s-r)
403 bezier_base():r(0.0),s(1.0),drs(1.0) { }
405 const value_type &a, const value_type &b, const value_type &c, const value_type &d,
406 const time_type &r=0, const time_type &s=1):
407 a(a),b(b),c(c),d(d),r(r),s(s),drs(1.0/(s-r)) { sync(); }
411 drs=time_type(1)/(s-r);
413 _coeff[1]= b*3 - a*3;
414 _coeff[2]= c*3 - b*6 + a*3;
415 _coeff[3]= d - c*3 + b*3 - a;
418 // 4 products, 3 sums, and 1 difference.
420 operator()(time_type t)const
421 { t-=r; t*=drs; return _coeff[0]+(_coeff[1]+(_coeff[2]+(_coeff[3])*t)*t)*t; }
423 void set_rs(time_type new_r, time_type new_s) { r=new_r; s=new_s; drs=time_type(1)/(s-r); }
424 void set_r(time_type new_r) { r=new_r; drs=time_type(1)/(s-r); }
425 void set_s(time_type new_s) { s=new_s; drs=time_type(1)/(s-r); }
426 const time_type &get_r()const { return r; }
427 const time_type &get_s()const { return s; }
428 time_type get_dt()const { return s-r; }
430 //! Bezier curve intersection function
431 //! Calculates the time of intersection
432 // for the calling curve.
434 time_type intersect(const bezier_base<value_type,time_type> &x, time_type t=0,int i=15)const
436 value_type system[4];
437 system[0]=_coeff[0]-x._coeff[0];
438 system[1]=_coeff[1]-x._coeff[1];
439 system[2]=_coeff[2]-x._coeff[2];
440 system[3]=_coeff[3]-x._coeff[3];
446 // Inner loop: 7 products, 5 sums, 1 difference
448 t-=(time_type) ( (system[0]+(system[1]+(system[2]+(system[3])*t)*t)*t)/
449 (system[1]+(system[2]*2+(system[3]*3)*t)*t) );
462 operator[](int i) const
472 template <typename V, typename T>
473 class bezier_iterator
477 struct iterator_category {};
478 typedef V value_type;
479 typedef T difference_type;
485 bezier_base<V,T> curve;
491 operator*(void)const { return curve(t); }
492 const surface_iterator&
495 { t+=dt; return &this; }
497 const surface_iterator&
499 { hermite_iterator _tmp=*this; t+=dt; return _tmp; }
501 const surface_iterator&
503 { t-=dt; return &this; }
505 const surface_iterator&
507 { hermite_iterator _tmp=*this; t-=dt; return _tmp; }
511 operator+(difference_type __n) const
512 { return surface_iterator(data+__n[0]+__n[1]*pitch,pitch); }
515 operator-(difference_type __n) const
516 { return surface_iterator(data-__n[0]-__n[1]*pitch,pitch); }
521 template <typename V,typename T=float>
522 class bezier : public bezier_base<V,T>
525 typedef V value_type;
527 typedef float distance_type;
528 typedef bezier_iterator<V,T> iterator;
529 typedef bezier_iterator<V,T> const_iterator;
531 distance_func<value_type> dist;
533 using bezier_base<V,T>::get_r;
534 using bezier_base<V,T>::get_s;
535 using bezier_base<V,T>::get_dt;
539 bezier(const value_type &a, const value_type &b, const value_type &c, const value_type &d):
540 bezier_base<V,T>(a,b,c,d) { }
543 const_iterator begin()const;
544 const_iterator end()const;
546 time_type find_closest(bool fast, const value_type& x, int i=7)const
550 value_type array[4] = {
551 bezier<V,T>::operator[](0),
552 bezier<V,T>::operator[](1),
553 bezier<V,T>::operator[](2),
554 bezier<V,T>::operator[](3)};
555 float t = NearestPointOnCurve(x, array);
556 return t > 0.999999 ? 0.999999 : t < 0.000001 ? 0.000001 : t;
560 time_type r(0), s(1);
561 float t((r+s)*0.5); /* half way between r and s */
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))
577 distance_type find_distance(time_type r, time_type s, int steps=7)const
579 const time_type inc((s-r)/steps);
580 distance_type ret(0);
581 value_type last(operator()(r));
583 for(r+=inc;r<s;r+=inc)
585 const value_type n(operator()(r));
586 ret+=dist.uncook(dist(last,n));
589 ret+=dist.uncook(dist(last,operator()(r)))*(s-(r-inc))/inc;
594 distance_type length()const { return find_distance(get_r(),get_s()); }
596 /* subdivide at some time t into 2 separate curves left and right
601 * 1+2 * 0+3*1+3*2+3 l4,r1
607 0.1 2.3 -> 0.1 2 3 4 5.6
609 void subdivide(bezier *left, bezier *right, const time_type &time = (time_type)0.5) const
611 time_type t=(time-get_r())/get_dt();
615 const value_type& a((*this)[0]);
616 const value_type& b((*this)[1]);
617 const value_type& c((*this)[2]);
618 const value_type& d((*this)[3]);
620 //1st stage points to keep
625 lt[1] = affine_func(a,b,t);
626 temp = affine_func(b,c,t);
627 rt[2] = affine_func(c,d,t);
630 lt[2] = affine_func(lt[1],temp,t);
631 rt[1] = affine_func(temp,rt[2],t);
634 lt[3] = rt[0] = affine_func(lt[2],rt[1],t);
636 //set the time range for l,r (the inside values should be 1, 0 respectively)
643 //give back the curves
645 if(right) *right = rt;
649 void evaluate(time_type t, value_type &f, value_type &df) const
651 t=(t-get_r())/get_dt();
653 const value_type& a((*this)[0]);
654 const value_type& b((*this)[1]);
655 const value_type& c((*this)[2]);
656 const value_type& d((*this)[3]);
658 const value_type p1 = affine_func(
662 const value_type p2 = affine_func(
667 f = affine_func(p1,p2,t);
674 * Evaluate a Bezier curve at a particular parameter value
675 * Fill in control points for resulting sub-curves if "Left" and
676 * "Right" are non-null.
678 * int degree; Degree of bezier curve
679 * value_type *VT; Control pts
680 * time_type t; Parameter value
681 * value_type *Left; RETURN left half ctl pts
682 * value_type *Right; RETURN right half ctl pts
684 static value_type Bezier(value_type *VT, int degree, time_type t, value_type *Left, value_type *Right)
686 int i, j; /* Index variables */
687 value_type Vtemp[W_DEGREE+1][W_DEGREE+1];
689 /* Copy control points */
690 for (j = 0; j <= degree; j++)
693 /* Triangle computation */
694 for (i = 1; i <= degree; i++)
695 for (j =0 ; j <= degree - i; j++)
697 Vtemp[i][j][0] = (1.0 - t) * Vtemp[i-1][j][0] + t * Vtemp[i-1][j+1][0];
698 Vtemp[i][j][1] = (1.0 - t) * Vtemp[i-1][j][1] + t * Vtemp[i-1][j+1][1];
702 for (j = 0; j <= degree; j++)
703 Left[j] = Vtemp[j][0];
706 for (j = 0; j <= degree; j++)
707 Right[j] = Vtemp[degree-j][j];
709 return (Vtemp[degree][0]);
714 * Count the number of times a Bezier control polygon
715 * crosses the 0-axis. This number is >= the number of roots.
717 * value_type *VT; Control pts of Bezier curve
719 static int CrossingCount(value_type *VT)
722 int n_crossings = 0; /* Number of zero-crossings */
723 int sign, old_sign; /* Sign of coefficients */
725 sign = old_sign = SGN(VT[0][1]);
726 for (i = 1; i <= W_DEGREE; i++)
728 sign = SGN(VT[i][1]);
729 if (sign != old_sign) n_crossings++;
737 * ControlPolygonFlatEnough :
738 * Check if the control polygon of a Bezier curve is flat enough
739 * for recursive subdivision to bottom out.
741 * value_type *VT; Control points
743 static int ControlPolygonFlatEnough(value_type *VT)
745 int i; /* Index variable */
746 distance_type distance[W_DEGREE]; /* Distances from pts to line */
747 distance_type max_distance_above; /* maximum of these */
748 distance_type max_distance_below;
749 time_type intercept_1, intercept_2, left_intercept, right_intercept;
750 distance_type a, b, c; /* Coefficients of implicit */
751 /* eqn for line from VT[0]-VT[deg] */
752 /* Find the perpendicular distance */
753 /* from each interior control point to */
754 /* line connecting VT[0] and VT[W_DEGREE] */
756 distance_type abSquared;
758 /* Derive the implicit equation for line connecting first *
759 * and last control points */
760 a = VT[0][1] - VT[W_DEGREE][1];
761 b = VT[W_DEGREE][0] - VT[0][0];
762 c = VT[0][0] * VT[W_DEGREE][1] - VT[W_DEGREE][0] * VT[0][1];
764 abSquared = (a * a) + (b * b);
766 for (i = 1; i < W_DEGREE; i++)
768 /* Compute distance from each of the points to that line */
769 distance[i] = a * VT[i][0] + b * VT[i][1] + c;
770 if (distance[i] > 0.0) distance[i] = (distance[i] * distance[i]) / abSquared;
771 if (distance[i] < 0.0) distance[i] = -(distance[i] * distance[i]) / abSquared;
775 /* Find the largest distance */
776 max_distance_above = max_distance_below = 0.0;
778 for (i = 1; i < W_DEGREE; i++)
780 if (distance[i] < 0.0) max_distance_below = MIN(max_distance_below, distance[i]);
781 if (distance[i] > 0.0) max_distance_above = MAX(max_distance_above, distance[i]);
784 /* Implicit equation for "above" line */
785 intercept_1 = -(c + max_distance_above)/a;
787 /* Implicit equation for "below" line */
788 intercept_2 = -(c + max_distance_below)/a;
790 /* Compute intercepts of bounding box */
791 left_intercept = MIN(intercept_1, intercept_2);
792 right_intercept = MAX(intercept_1, intercept_2);
794 return 0.5 * (right_intercept-left_intercept) < BEZIER_EPSILON ? 1 : 0;
798 * ComputeXIntercept :
799 * Compute intersection of chord from first control point to last
802 * value_type *VT; Control points
804 static time_type ComputeXIntercept(value_type *VT)
806 distance_type YNM = VT[W_DEGREE][1] - VT[0][1];
807 return (YNM*VT[0][0] - (VT[W_DEGREE][0] - VT[0][0])*VT[0][1]) / YNM;
812 * Given a 5th-degree equation in Bernstein-Bezier form, find
813 * all of the roots in the interval [0, 1]. Return the number
816 * value_type *w; The control points
817 * time_type *t; RETURN candidate t-values
818 * int depth; The depth of the recursion
820 static int FindRoots(value_type *w, time_type *t, int depth)
823 value_type Left[W_DEGREE+1]; /* New left and right */
824 value_type Right[W_DEGREE+1]; /* control polygons */
825 int left_count; /* Solution count from */
826 int right_count; /* children */
827 time_type left_t[W_DEGREE+1]; /* Solutions from kids */
828 time_type right_t[W_DEGREE+1];
830 switch (CrossingCount(w))
833 { /* No solutions here */
837 { /* Unique solution */
838 /* Stop recursion when the tree is deep enough */
839 /* if deep enough, return 1 solution at midpoint */
840 if (depth >= MAXDEPTH)
842 t[0] = (w[0][0] + w[W_DEGREE][0]) / 2.0;
845 if (ControlPolygonFlatEnough(w))
847 t[0] = ComputeXIntercept(w);
854 /* Otherwise, solve recursively after */
855 /* subdividing control polygon */
856 Bezier(w, W_DEGREE, 0.5, Left, Right);
857 left_count = FindRoots(Left, left_t, depth+1);
858 right_count = FindRoots(Right, right_t, depth+1);
860 /* Gather solutions together */
861 for (i = 0; i < left_count; i++) t[i] = left_t[i];
862 for (i = 0; i < right_count; i++) t[i+left_count] = right_t[i];
864 /* Send back total number of solutions */
865 return (left_count+right_count);
869 * ConvertToBezierForm :
870 * Given a point and a Bezier curve, generate a 5th-degree
871 * Bezier-format equation whose solution finds the point on the
872 * curve nearest the user-defined point.
874 * value_type& P; The point to find t for
875 * value_type *VT; The control points
877 static void ConvertToBezierForm(const value_type& P, value_type *VT, value_type w[W_DEGREE+1])
879 int i, j, k, m, n, ub, lb;
880 int row, column; /* Table indices */
881 value_type c[DEGREE+1]; /* VT(i)'s - P */
882 value_type d[DEGREE]; /* VT(i+1) - VT(i) */
883 distance_type cdTable[3][4]; /* Dot product of c, d */
884 static distance_type z[3][4] = { /* Precomputed "z" for cubics */
885 {1.0, 0.6, 0.3, 0.1},
886 {0.4, 0.6, 0.6, 0.4},
887 {0.1, 0.3, 0.6, 1.0}};
889 /* Determine the c's -- these are vectors created by subtracting */
890 /* point P from each of the control points */
891 for (i = 0; i <= DEGREE; i++)
894 /* Determine the d's -- these are vectors created by subtracting */
895 /* each control point from the next */
896 for (i = 0; i <= DEGREE - 1; i++)
897 d[i] = (VT[i+1] - VT[i]) * 3.0;
899 /* Create the c,d table -- this is a table of dot products of the */
901 for (row = 0; row <= DEGREE - 1; row++)
902 for (column = 0; column <= DEGREE; column++)
903 cdTable[row][column] = d[row] * c[column];
905 /* Now, apply the z's to the dot products, on the skew diagonal */
906 /* Also, set up the x-values, making these "points" */
907 for (i = 0; i <= W_DEGREE; i++)
909 w[i][0] = (distance_type)(i) / W_DEGREE;
915 for (k = 0; k <= n + m; k++)
919 for (i = lb; i <= ub; i++)
922 w[i+j][1] += cdTable[j][i] * z[j][i];
928 * NearestPointOnCurve :
929 * Compute the parameter value of the point on a Bezier
930 * curve segment closest to some arbitrary, user-input point.
931 * Return the point on the curve at that parameter value.
933 * value_type& P; The user-supplied point
934 * value_type *VT; Control points of cubic Bezier
936 static time_type NearestPointOnCurve(const value_type& P, value_type VT[4])
938 value_type w[W_DEGREE+1]; /* Ctl pts of 5th-degree curve */
939 time_type t_candidate[W_DEGREE]; /* Possible roots */
940 int n_solutions; /* Number of roots found */
941 time_type t; /* Parameter value of closest pt */
943 /* Convert problem to 5th-degree Bezier form */
944 ConvertToBezierForm(P, VT, w);
946 /* Find all possible roots of 5th-degree equation */
947 n_solutions = FindRoots(w, t_candidate, 0);
949 /* Compare distances of P to all candidates, and to t=0, and t=1 */
951 distance_type dist, new_dist;
955 /* Check distance to beginning of curve, where t = 0 */
956 dist = (P - VT[0]).mag_squared();
959 /* Find distances for candidate points */
960 for (i = 0; i < n_solutions; i++)
962 p = Bezier(VT, DEGREE, t_candidate[i], (value_type *)NULL, (value_type *)NULL);
963 new_dist = (P - p).mag_squared();
971 /* Finally, look at distance to end point, where t = 1.0 */
972 new_dist = (P - VT[DEGREE]).mag_squared();
980 /* Return the point on the curve at parameter value t */
987 /* === E X T E R N S ======================================================= */
989 /* === E N D =============================================================== */