+++ /dev/null
-/*! ========================================================================
-** 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 <cmath> // for ldexp
-// #include <ETL/fixed> // 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<typename V,typename T> 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 <typename V,typename T=float>
-class bezier_base : public std::unary_function<T,V>
-{
-public:
- typedef V value_type;
- typedef T time_type;
-
-private:
- value_type a,b,c,d;
- time_type r,s;
-
-protected:
- affine_combo<value_type,time_type> 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<value_type,time_type> &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<value_type,time_type> &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<float,float> : public std::unary_function<float,float>
-{
-public:
- typedef float value_type;
- typedef float time_type;
-private:
- affine_combo<value_type,time_type> 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<value_type,time_type> &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<double,float> : public std::unary_function<float,double>
-{
-public:
- typedef double value_type;
- typedef float time_type;
-private:
- affine_combo<value_type,time_type> 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<value_type,time_type> &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 T,unsigned int FIXED_BITS>
-class bezier_base<fixed_base<T,FIXED_BITS> > : std::unary_function<fixed_base<T,FIXED_BITS>,fixed_base<T,FIXED_BITS> >
-{
-public:
- typedef fixed_base<T,FIXED_BITS> value_type;
- typedef fixed_base<T,FIXED_BITS> time_type;
-
-private:
- affine_combo<value_type,time_type> 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<value_type,time_type> &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 <typename V, typename T>
-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<V,T> 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 <typename V,typename T=float>
-class bezier : public bezier_base<V,T>
-{
-public:
- typedef V value_type;
- typedef T time_type;
- typedef float distance_type;
- typedef bezier_iterator<V,T> iterator;
- typedef bezier_iterator<V,T> const_iterator;
-
- distance_func<value_type> dist;
-
- using bezier_base<V,T>::get_r;
- using bezier_base<V,T>::get_s;
- using bezier_base<V,T>::get_dt;
-
-public:
- bezier() { }
- bezier(const value_type &a, const value_type &b, const value_type &c, const value_type &d):
- bezier_base<V,T>(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<V,T>::operator[](0),
- bezier<V,T>::operator[](1),
- bezier<V,T>::operator[](2),
- bezier<V,T>::operator[](3)};
- float t = NearestPointOnCurve(x, array);
- return t > 0.999999 ? 0.999999 : t < 0.000001 ? 0.000001 : t;
- }
- 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);
- distance_type ret(0);
- value_type last(operator()(r));
-
- for(r+=inc;r<s;r+=inc)
- {
- const value_type n(operator()(r));
- ret+=dist.uncook(dist(last,n));
- last=n;
- }
- ret+=dist.uncook(dist(last,operator()(r)))*(s-(r-inc))/inc;
-
- return ret;
- }
-
- distance_type length()const { return find_distance(get_r(),get_s()); }
-
- /* 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 *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
projects / synfig.git / blobdiff