X-Git-Url: https://git.pterodactylus.net/?a=blobdiff_plain;f=ETL%2Ftrunk%2FETL%2F_bezier.h;fp=ETL%2Ftrunk%2FETL%2F_bezier.h;h=0000000000000000000000000000000000000000;hb=a095981e18cc37a8ecc7cd237cc22b9c10329264;hp=fe6b17017f1e89b13a62f40efbc9f3ed67322d58;hpb=9459638ad6797b8139f1e9f0715c96076dbf0890;p=synfig.git diff --git a/ETL/trunk/ETL/_bezier.h b/ETL/trunk/ETL/_bezier.h deleted file mode 100644 index fe6b170..0000000 --- a/ETL/trunk/ETL/_bezier.h +++ /dev/null @@ -1,992 +0,0 @@ -/*! ======================================================================== -** Extended Template Library -** Bezier Template Class Implementation -** $Id$ -** -** Copyright (c) 2002 Robert B. Quattlebaum Jr. -** Copyright (c) 2007 Chris Moore -** -** This package is free software; you can redistribute it and/or -** modify it under the terms of the GNU General Public License as -** published by the Free Software Foundation; either version 2 of -** the License, or (at your option) any later version. -** -** This package is distributed in the hope that it will be useful, -** but WITHOUT ANY WARRANTY; without even the implied warranty of -** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -** General Public License for more details. -** -** === N O T E S =========================================================== -** -** This is an internal header file, included by other ETL headers. -** You should not attempt to use it directly. -** -** ========================================================================= */ - -/* === S T A R T =========================================================== */ - -#ifndef __ETL__BEZIER_H -#define __ETL__BEZIER_H - -/* === H E A D E R S ======================================================= */ - -#include "_curve_func.h" -#include // for ldexp -// #include // not used - -/* === M A C R O S ========================================================= */ - -#define MAXDEPTH 64 /* Maximum depth for recursion */ - -/* take binary sign of a, either -1, or 1 if >= 0 */ -#define SGN(a) (((a)<0) ? -1 : 1) - -/* find minimum of a and b */ -#ifndef MIN -#define MIN(a,b) (((a)<(b))?(a):(b)) -#endif - -/* find maximum of a and b */ -#ifndef MAX -#define MAX(a,b) (((a)>(b))?(a):(b)) -#endif - -#define BEZIER_EPSILON (ldexp(1.0,-MAXDEPTH-1)) /*Flatness control value */ -//#define BEZIER_EPSILON 0.00005 /*Flatness control value */ -#define DEGREE 3 /* Cubic Bezier curve */ -#define W_DEGREE 5 /* Degree of eqn to find roots of */ - -/* === T Y P E D E F S ===================================================== */ - -/* === C L A S S E S & S T R U C T S ======================================= */ - -_ETL_BEGIN_NAMESPACE - -template class bezier; - -//! Cubic Bezier Curve Base Class -// This generic implementation uses the DeCasteljau algorithm. -// Works for just about anything that has an affine combination function -template -class bezier_base : public std::unary_function -{ -public: - typedef V value_type; - typedef T time_type; - -private: - value_type a,b,c,d; - time_type r,s; - -protected: - affine_combo affine_func; - -public: - bezier_base():r(0.0),s(1.0) { } - bezier_base( - const value_type &a, const value_type &b, const value_type &c, const value_type &d, - const time_type &r=0.0, const time_type &s=1.0): - a(a),b(b),c(c),d(d),r(r),s(s) { sync(); } - - void sync() - { - } - - value_type - operator()(time_type t)const - { - t=(t-r)/(s-r); - return - affine_func( - affine_func( - affine_func(a,b,t), - affine_func(b,c,t) - ,t), - affine_func( - affine_func(b,c,t), - affine_func(c,d,t) - ,t) - ,t); - } - - /* - void evaluate(time_type t, value_type &f, value_type &df) const - { - t=(t-r)/(s-r); - - value_type p1 = affine_func( - affine_func(a,b,t), - affine_func(b,c,t) - ,t); - value_type p2 = affine_func( - affine_func(b,c,t), - affine_func(c,d,t) - ,t); - - f = affine_func(p1,p2,t); - df = (p2-p1)*3; - } - */ - - void set_rs(time_type new_r, time_type new_s) { r=new_r; s=new_s; } - void set_r(time_type new_r) { r=new_r; } - void set_s(time_type new_s) { s=new_s; } - const time_type &get_r()const { return r; } - const time_type &get_s()const { return s; } - time_type get_dt()const { return s-r; } - - bool intersect_hull(const bezier_base &x)const - { - return 0; - } - - //! Bezier curve intersection function - /*! Calculates the time of intersection - ** for the calling curve. - ** - ** I still have not figured out a good generic - ** method of doing this for a bi-infinite - ** cubic bezier curve calculated with the DeCasteljau - ** algorithm. - ** - ** One method, although it does not work for the - ** entire bi-infinite curve, is to iteratively - ** intersect the hulls. However, we would only detect - ** intersections that occur between R and S. - ** - ** It is entirely possible that a new construct similar - ** to the affine combination function will be necessary - ** for this to work properly. - ** - ** For now, this function is BROKEN. (although it works - ** for the floating-point specializations, using newton's method) - */ - time_type intersect(const bezier_base &x, time_type near=0.0)const - { - return 0; - } - - /* subdivide at some time t into 2 separate curves left and right - - b0 l1 - * 0+1 l2 - b1 * 1+2*1+2 l3 - * 1+2 * 0+3*1+3*2+3 l4,r1 - b2 * 1+2*2+2 r2 * - * 2+3 r3 * - b3 r4 * - * - - 0.1 2.3 -> 0.1 2 3 4 5.6 - */ -/* void subdivide(bezier_base *left, bezier_base *right, const time_type &time = (time_type)0.5) const - { - time_type t = (time-r)/(s-r); - bezier_base lt,rt; - - value_type temp; - - //1st stage points to keep - lt.a = a; - rt.d = d; - - //2nd stage calc - lt.b = affine_func(a,b,t); - temp = affine_func(b,c,t); - rt.c = affine_func(c,d,t); - - //3rd stage calc - lt.c = affine_func(lt.b,temp,t); - rt.b = affine_func(temp,rt.c,t); - - //last stage calc - lt.d = rt.a = affine_func(lt.c,rt.b,t); - - //set the time range for l,r (the inside values should be 1, 0 respectively) - lt.r = r; - rt.s = s; - - //give back the curves - if(left) *left = lt; - if(right) *right = rt; - } - */ - value_type & - operator[](int i) - { return (&a)[i]; } - - const value_type & - operator[](int i) const - { return (&a)[i]; } -}; - - -#if 1 -// Fast float implementation of a cubic bezier curve -template <> -class bezier_base : public std::unary_function -{ -public: - typedef float value_type; - typedef float time_type; -private: - affine_combo affine_func; - value_type a,b,c,d; - time_type r,s; - - value_type _coeff[4]; - time_type drs; // reciprocal of (s-r) -public: - bezier_base():r(0.0),s(1.0),drs(1.0) { } - bezier_base( - const value_type &a, const value_type &b, const value_type &c, const value_type &d, - const time_type &r=0.0, const time_type &s=1.0): - a(a),b(b),c(c),d(d),r(r),s(s),drs(1.0/(s-r)) { sync(); } - - void sync() - { -// drs=1.0/(s-r); - _coeff[0]= a; - _coeff[1]= b*3 - a*3; - _coeff[2]= c*3 - b*6 + a*3; - _coeff[3]= d - c*3 + b*3 - a; - } - - // Cost Summary: 4 products, 3 sums, and 1 difference. - inline value_type - operator()(time_type t)const - { t-=r; t*=drs; return _coeff[0]+(_coeff[1]+(_coeff[2]+(_coeff[3])*t)*t)*t; } - - void set_rs(time_type new_r, time_type new_s) { r=new_r; s=new_s; drs=1.0/(s-r); } - void set_r(time_type new_r) { r=new_r; drs=1.0/(s-r); } - void set_s(time_type new_s) { s=new_s; drs=1.0/(s-r); } - const time_type &get_r()const { return r; } - const time_type &get_s()const { return s; } - time_type get_dt()const { return s-r; } - - //! Bezier curve intersection function - /*! Calculates the time of intersection - ** for the calling curve. - */ - time_type intersect(const bezier_base &x, time_type t=0.0,int i=15)const - { - //BROKEN - the time values of the 2 curves should be independent - value_type system[4]; - system[0]=_coeff[0]-x._coeff[0]; - system[1]=_coeff[1]-x._coeff[1]; - system[2]=_coeff[2]-x._coeff[2]; - system[3]=_coeff[3]-x._coeff[3]; - - t-=r; - t*=drs; - - // Newton's method - // Inner loop cost summary: 7 products, 5 sums, 1 difference - for(;i;i--) - t-= (system[0]+(system[1]+(system[2]+(system[3])*t)*t)*t)/ - (system[1]+(system[2]*2+(system[3]*3)*t)*t); - - t*=(s-r); - t+=r; - - return t; - } - - value_type & - operator[](int i) - { return (&a)[i]; } - - const value_type & - operator[](int i) const - { return (&a)[i]; } -}; - - -// Fast double implementation of a cubic bezier curve -template <> -class bezier_base : public std::unary_function -{ -public: - typedef double value_type; - typedef float time_type; -private: - affine_combo affine_func; - value_type a,b,c,d; - time_type r,s; - - value_type _coeff[4]; - time_type drs; // reciprocal of (s-r) -public: - bezier_base():r(0.0),s(1.0),drs(1.0) { } - bezier_base( - const value_type &a, const value_type &b, const value_type &c, const value_type &d, - const time_type &r=0.0, const time_type &s=1.0): - a(a),b(b),c(c),d(d),r(r),s(s),drs(1.0/(s-r)) { sync(); } - - void sync() - { -// drs=1.0/(s-r); - _coeff[0]= a; - _coeff[1]= b*3 - a*3; - _coeff[2]= c*3 - b*6 + a*3; - _coeff[3]= d - c*3 + b*3 - a; - } - - // 4 products, 3 sums, and 1 difference. - inline value_type - operator()(time_type t)const - { t-=r; t*=drs; return _coeff[0]+(_coeff[1]+(_coeff[2]+(_coeff[3])*t)*t)*t; } - - void set_rs(time_type new_r, time_type new_s) { r=new_r; s=new_s; drs=1.0/(s-r); } - void set_r(time_type new_r) { r=new_r; drs=1.0/(s-r); } - void set_s(time_type new_s) { s=new_s; drs=1.0/(s-r); } - const time_type &get_r()const { return r; } - const time_type &get_s()const { return s; } - time_type get_dt()const { return s-r; } - - //! Bezier curve intersection function - /*! Calculates the time of intersection - ** for the calling curve. - */ - time_type intersect(const bezier_base &x, time_type t=0.0,int i=15)const - { - //BROKEN - the time values of the 2 curves should be independent - value_type system[4]; - system[0]=_coeff[0]-x._coeff[0]; - system[1]=_coeff[1]-x._coeff[1]; - system[2]=_coeff[2]-x._coeff[2]; - system[3]=_coeff[3]-x._coeff[3]; - - t-=r; - t*=drs; - - // Newton's method - // Inner loop: 7 products, 5 sums, 1 difference - for(;i;i--) - t-= (system[0]+(system[1]+(system[2]+(system[3])*t)*t)*t)/ - (system[1]+(system[2]*2+(system[3]*3)*t)*t); - - t*=(s-r); - t+=r; - - return t; - } - - value_type & - operator[](int i) - { return (&a)[i]; } - - const value_type & - operator[](int i) const - { return (&a)[i]; } -}; - -//#ifdef __FIXED__ - -// Fast double implementation of a cubic bezier curve -/* -template <> -template -class bezier_base > : std::unary_function,fixed_base > -{ -public: - typedef fixed_base value_type; - typedef fixed_base time_type; - -private: - affine_combo affine_func; - value_type a,b,c,d; - time_type r,s; - - value_type _coeff[4]; - time_type drs; // reciprocal of (s-r) -public: - bezier_base():r(0.0),s(1.0),drs(1.0) { } - bezier_base( - const value_type &a, const value_type &b, const value_type &c, const value_type &d, - const time_type &r=0, const time_type &s=1): - a(a),b(b),c(c),d(d),r(r),s(s),drs(1.0/(s-r)) { sync(); } - - void sync() - { - drs=time_type(1)/(s-r); - _coeff[0]= a; - _coeff[1]= b*3 - a*3; - _coeff[2]= c*3 - b*6 + a*3; - _coeff[3]= d - c*3 + b*3 - a; - } - - // 4 products, 3 sums, and 1 difference. - inline value_type - operator()(time_type t)const - { t-=r; t*=drs; return _coeff[0]+(_coeff[1]+(_coeff[2]+(_coeff[3])*t)*t)*t; } - - void set_rs(time_type new_r, time_type new_s) { r=new_r; s=new_s; drs=time_type(1)/(s-r); } - void set_r(time_type new_r) { r=new_r; drs=time_type(1)/(s-r); } - void set_s(time_type new_s) { s=new_s; drs=time_type(1)/(s-r); } - const time_type &get_r()const { return r; } - const time_type &get_s()const { return s; } - time_type get_dt()const { return s-r; } - - //! Bezier curve intersection function - //! Calculates the time of intersection - // for the calling curve. - // - time_type intersect(const bezier_base &x, time_type t=0,int i=15)const - { - value_type system[4]; - system[0]=_coeff[0]-x._coeff[0]; - system[1]=_coeff[1]-x._coeff[1]; - system[2]=_coeff[2]-x._coeff[2]; - system[3]=_coeff[3]-x._coeff[3]; - - t-=r; - t*=drs; - - // Newton's method - // Inner loop: 7 products, 5 sums, 1 difference - for(;i;i--) - t-=(time_type) ( (system[0]+(system[1]+(system[2]+(system[3])*t)*t)*t)/ - (system[1]+(system[2]*2+(system[3]*3)*t)*t) ); - - t*=(s-r); - t+=r; - - return t; - } - - value_type & - operator[](int i) - { return (&a)[i]; } - - const value_type & - operator[](int i) const - { return (&a)[i]; } -}; -*/ -//#endif - -#endif - - - -template -class bezier_iterator -{ -public: - - struct iterator_category {}; - typedef V value_type; - typedef T difference_type; - typedef V reference; - -private: - difference_type t; - difference_type dt; - bezier_base curve; - -public: - -/* - reference - operator*(void)const { return curve(t); } - const surface_iterator& - - operator++(void) - { t+=dt; return &this; } - - const surface_iterator& - operator++(int) - { hermite_iterator _tmp=*this; t+=dt; return _tmp; } - - const surface_iterator& - operator--(void) - { t-=dt; return &this; } - - const surface_iterator& - operator--(int) - { hermite_iterator _tmp=*this; t-=dt; return _tmp; } - - - surface_iterator - operator+(difference_type __n) const - { return surface_iterator(data+__n[0]+__n[1]*pitch,pitch); } - - surface_iterator - operator-(difference_type __n) const - { return surface_iterator(data-__n[0]-__n[1]*pitch,pitch); } -*/ - -}; - -template -class bezier : public bezier_base -{ -public: - typedef V value_type; - typedef T time_type; - typedef float distance_type; - typedef bezier_iterator iterator; - typedef bezier_iterator const_iterator; - - distance_func dist; - - using bezier_base::get_r; - using bezier_base::get_s; - using bezier_base::get_dt; - -public: - bezier() { } - bezier(const value_type &a, const value_type &b, const value_type &c, const value_type &d): - bezier_base(a,b,c,d) { } - - - const_iterator begin()const; - const_iterator end()const; - - time_type find_closest(bool fast, const value_type& x, int i=7)const - { - if (!fast) - { - value_type array[4] = { - bezier::operator[](0), - bezier::operator[](1), - bezier::operator[](2), - bezier::operator[](3)}; - return NearestPointOnCurve(x, array); - } - else - { - time_type r(0), s(1); - float t((r+s)*0.5); /* half way between r and s */ - - for(;i;i--) - { - // compare 33% of the way between r and s with 67% of the way between r and s - if(dist(operator()((s-r)*(1.0/3.0)+r), x) < - dist(operator()((s-r)*(2.0/3.0)+r), x)) - s=t; - else - r=t; - t=((r+s)*0.5); - } - return t; - } - } - - distance_type find_distance(time_type r, time_type s, int steps=7)const - { - const time_type inc((s-r)/steps); - if (!inc) return 0; - distance_type ret(0); - value_type last(operator()(r)); - - for(r+=inc;r 0.1 2 3 4 5.6 - */ - void subdivide(bezier *left, bezier *right, const time_type &time = (time_type)0.5) const - { - time_type t=(time-get_r())/get_dt(); - bezier lt,rt; - - value_type temp; - const value_type& a((*this)[0]); - const value_type& b((*this)[1]); - const value_type& c((*this)[2]); - const value_type& d((*this)[3]); - - //1st stage points to keep - lt[0] = a; - rt[3] = d; - - //2nd stage calc - lt[1] = affine_func(a,b,t); - temp = affine_func(b,c,t); - rt[2] = affine_func(c,d,t); - - //3rd stage calc - lt[2] = affine_func(lt[1],temp,t); - rt[1] = affine_func(temp,rt[2],t); - - //last stage calc - lt[3] = rt[0] = affine_func(lt[2],rt[1],t); - - //set the time range for l,r (the inside values should be 1, 0 respectively) - lt.set_r(get_r()); - rt.set_s(get_s()); - - lt.sync(); - rt.sync(); - - //give back the curves - if(left) *left = lt; - if(right) *right = rt; - } - - - void evaluate(time_type t, value_type &f, value_type &df) const - { - t=(t-get_r())/get_dt(); - - const value_type& a((*this)[0]); - const value_type& b((*this)[1]); - const value_type& c((*this)[2]); - const value_type& d((*this)[3]); - - const value_type p1 = affine_func( - affine_func(a,b,t), - affine_func(b,c,t) - ,t); - const value_type p2 = affine_func( - affine_func(b,c,t), - affine_func(c,d,t) - ,t); - - f = affine_func(p1,p2,t); - df = (p2-p1)*3; - } - -private: - /* - * Bezier : - * Evaluate a Bezier curve at a particular parameter value - * Fill in control points for resulting sub-curves if "Left" and - * "Right" are non-null. - * - * int degree; Degree of bezier curve - * value_type *VT; Control pts - * time_type t; Parameter value - * value_type *Left; RETURN left half ctl pts - * value_type *Right; RETURN right half ctl pts - */ - static value_type Bezier(value_type *VT, int degree, time_type t, value_type *Left, value_type *Right) - { - int i, j; /* Index variables */ - value_type Vtemp[W_DEGREE+1][W_DEGREE+1]; - - /* Copy control points */ - for (j = 0; j <= degree; j++) - Vtemp[0][j] = VT[j]; - - /* Triangle computation */ - for (i = 1; i <= degree; i++) - for (j =0 ; j <= degree - i; j++) - { - Vtemp[i][j][0] = (1.0 - t) * Vtemp[i-1][j][0] + t * Vtemp[i-1][j+1][0]; - Vtemp[i][j][1] = (1.0 - t) * Vtemp[i-1][j][1] + t * Vtemp[i-1][j+1][1]; - } - - if (Left != NULL) - for (j = 0; j <= degree; j++) - Left[j] = Vtemp[j][0]; - - if (Right != NULL) - for (j = 0; j <= degree; j++) - Right[j] = Vtemp[degree-j][j]; - - return (Vtemp[degree][0]); - } - - /* - * CrossingCount : - * Count the number of times a Bezier control polygon - * crosses the 0-axis. This number is >= the number of roots. - * - * value_type *VT; Control pts of Bezier curve - */ - static int CrossingCount(value_type *VT) - { - int i; - int n_crossings = 0; /* Number of zero-crossings */ - int sign, old_sign; /* Sign of coefficients */ - - sign = old_sign = SGN(VT[0][1]); - for (i = 1; i <= W_DEGREE; i++) - { - sign = SGN(VT[i][1]); - if (sign != old_sign) n_crossings++; - old_sign = sign; - } - - return n_crossings; - } - - /* - * ControlPolygonFlatEnough : - * Check if the control polygon of a Bezier curve is flat enough - * for recursive subdivision to bottom out. - * - * value_type *VT; Control points - */ - static int ControlPolygonFlatEnough(value_type *VT) - { - int i; /* Index variable */ - distance_type distance[W_DEGREE]; /* Distances from pts to line */ - distance_type max_distance_above; /* maximum of these */ - distance_type max_distance_below; - time_type intercept_1, intercept_2, left_intercept, right_intercept; - distance_type a, b, c; /* Coefficients of implicit */ - /* eqn for line from VT[0]-VT[deg] */ - /* Find the perpendicular distance */ - /* from each interior control point to */ - /* line connecting VT[0] and VT[W_DEGREE] */ - { - distance_type abSquared; - - /* Derive the implicit equation for line connecting first * - * and last control points */ - a = VT[0][1] - VT[W_DEGREE][1]; - b = VT[W_DEGREE][0] - VT[0][0]; - c = VT[0][0] * VT[W_DEGREE][1] - VT[W_DEGREE][0] * VT[0][1]; - - abSquared = (a * a) + (b * b); - - for (i = 1; i < W_DEGREE; i++) - { - /* Compute distance from each of the points to that line */ - distance[i] = a * VT[i][0] + b * VT[i][1] + c; - if (distance[i] > 0.0) distance[i] = (distance[i] * distance[i]) / abSquared; - if (distance[i] < 0.0) distance[i] = -(distance[i] * distance[i]) / abSquared; - } - } - - /* Find the largest distance */ - max_distance_above = max_distance_below = 0.0; - - for (i = 1; i < W_DEGREE; i++) - { - if (distance[i] < 0.0) max_distance_below = MIN(max_distance_below, distance[i]); - if (distance[i] > 0.0) max_distance_above = MAX(max_distance_above, distance[i]); - } - - /* Implicit equation for "above" line */ - intercept_1 = -(c + max_distance_above)/a; - - /* Implicit equation for "below" line */ - intercept_2 = -(c + max_distance_below)/a; - - /* Compute intercepts of bounding box */ - left_intercept = MIN(intercept_1, intercept_2); - right_intercept = MAX(intercept_1, intercept_2); - - return 0.5 * (right_intercept-left_intercept) < BEZIER_EPSILON ? 1 : 0; - } - - /* - * ComputeXIntercept : - * Compute intersection of chord from first control point to last - * with 0-axis. - * - * value_type *VT; Control points - */ - static time_type ComputeXIntercept(value_type *VT) - { - distance_type YNM = VT[W_DEGREE][1] - VT[0][1]; - return (YNM*VT[0][0] - (VT[W_DEGREE][0] - VT[0][0])*VT[0][1]) / YNM; - } - - /* - * FindRoots : - * Given a 5th-degree equation in Bernstein-Bezier form, find - * all of the roots in the interval [0, 1]. Return the number - * of roots found. - * - * value_type *w; The control points - * time_type *t; RETURN candidate t-values - * int depth; The depth of the recursion - */ - static int FindRoots(value_type *w, time_type *t, int depth) - { - int i; - value_type Left[W_DEGREE+1]; /* New left and right */ - value_type Right[W_DEGREE+1]; /* control polygons */ - int left_count; /* Solution count from */ - int right_count; /* children */ - time_type left_t[W_DEGREE+1]; /* Solutions from kids */ - time_type right_t[W_DEGREE+1]; - - switch (CrossingCount(w)) - { - case 0 : - { /* No solutions here */ - return 0; - } - case 1 : - { /* Unique solution */ - /* Stop recursion when the tree is deep enough */ - /* if deep enough, return 1 solution at midpoint */ - if (depth >= MAXDEPTH) - { - t[0] = (w[0][0] + w[W_DEGREE][0]) / 2.0; - return 1; - } - if (ControlPolygonFlatEnough(w)) - { - t[0] = ComputeXIntercept(w); - return 1; - } - break; - } - } - - /* Otherwise, solve recursively after */ - /* subdividing control polygon */ - Bezier(w, W_DEGREE, 0.5, Left, Right); - left_count = FindRoots(Left, left_t, depth+1); - right_count = FindRoots(Right, right_t, depth+1); - - /* Gather solutions together */ - for (i = 0; i < left_count; i++) t[i] = left_t[i]; - for (i = 0; i < right_count; i++) t[i+left_count] = right_t[i]; - - /* Send back total number of solutions */ - return (left_count+right_count); - } - - /* - * ConvertToBezierForm : - * Given a point and a Bezier curve, generate a 5th-degree - * Bezier-format equation whose solution finds the point on the - * curve nearest the user-defined point. - * - * value_type& P; The point to find t for - * value_type *VT; The control points - */ - static void ConvertToBezierForm(const value_type& P, value_type *VT, value_type w[W_DEGREE+1]) - { - int i, j, k, m, n, ub, lb; - int row, column; /* Table indices */ - value_type c[DEGREE+1]; /* VT(i)'s - P */ - value_type d[DEGREE]; /* VT(i+1) - VT(i) */ - distance_type cdTable[3][4]; /* Dot product of c, d */ - static distance_type z[3][4] = { /* Precomputed "z" for cubics */ - {1.0, 0.6, 0.3, 0.1}, - {0.4, 0.6, 0.6, 0.4}, - {0.1, 0.3, 0.6, 1.0}}; - - /* Determine the c's -- these are vectors created by subtracting */ - /* point P from each of the control points */ - for (i = 0; i <= DEGREE; i++) - c[i] = VT[i] - P; - - /* Determine the d's -- these are vectors created by subtracting */ - /* each control point from the next */ - for (i = 0; i <= DEGREE - 1; i++) - d[i] = (VT[i+1] - VT[i]) * 3.0; - - /* Create the c,d table -- this is a table of dot products of the */ - /* c's and d's */ - for (row = 0; row <= DEGREE - 1; row++) - for (column = 0; column <= DEGREE; column++) - cdTable[row][column] = d[row] * c[column]; - - /* Now, apply the z's to the dot products, on the skew diagonal */ - /* Also, set up the x-values, making these "points" */ - for (i = 0; i <= W_DEGREE; i++) - { - w[i][0] = (distance_type)(i) / W_DEGREE; - w[i][1] = 0.0; - } - - n = DEGREE; - m = DEGREE-1; - for (k = 0; k <= n + m; k++) - { - lb = MAX(0, k - m); - ub = MIN(k, n); - for (i = lb; i <= ub; i++) - { - j = k - i; - w[i+j][1] += cdTable[j][i] * z[j][i]; - } - } - } - - /* - * NearestPointOnCurve : - * Compute the parameter value of the point on a Bezier - * curve segment closest to some arbitrary, user-input point. - * Return the point on the curve at that parameter value. - * - * value_type& P; The user-supplied point - * value_type *VT; Control points of cubic Bezier - */ - static time_type NearestPointOnCurve(const value_type& P, value_type VT[4]) - { - value_type w[W_DEGREE+1]; /* Ctl pts of 5th-degree curve */ - time_type t_candidate[W_DEGREE]; /* Possible roots */ - int n_solutions; /* Number of roots found */ - time_type t; /* Parameter value of closest pt */ - - /* Convert problem to 5th-degree Bezier form */ - ConvertToBezierForm(P, VT, w); - - /* Find all possible roots of 5th-degree equation */ - n_solutions = FindRoots(w, t_candidate, 0); - - /* Compare distances of P to all candidates, and to t=0, and t=1 */ - { - distance_type dist, new_dist; - value_type p, v; - int i; - - /* Check distance to beginning of curve, where t = 0 */ - dist = (P - VT[0]).mag_squared(); - t = 0.0; - - /* Find distances for candidate points */ - for (i = 0; i < n_solutions; i++) - { - p = Bezier(VT, DEGREE, t_candidate[i], (value_type *)NULL, (value_type *)NULL); - new_dist = (P - p).mag_squared(); - if (new_dist < dist) - { - dist = new_dist; - t = t_candidate[i]; - } - } - - /* Finally, look at distance to end point, where t = 1.0 */ - new_dist = (P - VT[DEGREE]).mag_squared(); - if (new_dist < dist) - { - dist = new_dist; - t = 1.0; - } - } - - /* Return the point on the curve at parameter value t */ - return t; - } -}; - -_ETL_END_NAMESPACE - -/* === E X T E R N S ======================================================= */ - -/* === E N D =============================================================== */ - -#endif