API
4.2
For C++ developers

This is a low level Quintic Bezier curve class that contains functions to design continuous sets of 'C' shaped Bezier curves, and to evaluate their values and derivatives. More...
Static Public Member Functions  
static double  calcU (double ax, const SimTK::Vector &bezierPtsX, const SimTK::Spline &splineUX, double tol, int maxIter) 
This function will compute the u value that corresponds to the given x for a quintic Bezier curve. More...  
static int  calcIndex (double x, const SimTK::Matrix &bezierPtsX) 
Given a set of Bezier curve control points, return the index of the set of control points that x lies within. More...  
static int  calcIndex (double x, const SimTK::Array_< SimTK::Vector > &bezierPtsX) 
static double  calcQuinticBezierCurveVal (double u, const SimTK::Vector &pts) 
Calculates the value of a quintic Bezier curve at value u. More...  
static double  calcQuinticBezierCurveDerivU (double u, const SimTK::Vector &pts, int order) 
Calculates the value of a quintic Bezier derivative curve at value u. More...  
static double  calcQuinticBezierCurveDerivDYDX (double u, const SimTK::Vector &xpts, const SimTK::Vector &ypts, int order) 
Calculates the value of dydx of a quintic Bezier curve derivative at u. More...  
static SimTK::Matrix  calcQuinticBezierCornerControlPoints (double x0, double y0, double dydx0, double x1, double y1, double dydx1, double curviness) 
Calculates the location of quintic Bezier curve control points to create a C shaped curve. More...  
static SimTK::Matrix  calcNumIntBezierYfcnX (const SimTK::Vector &vX, double ic0, double intAcc, double uTol, int uMaxIter, const SimTK::Matrix &mX, const SimTK::Matrix &mY, const SimTK::Array_< SimTK::Spline > &aSplineUX, bool flag_intLeftToRight, const std::string &name) 
This function numerically integrates the Bezier curve y(x). More...  
This is a low level Quintic Bezier curve class that contains functions to design continuous sets of 'C' shaped Bezier curves, and to evaluate their values and derivatives.
A set in this context is used to refer to 2 or more quintic Bezier curves that are continuously connected to each other to form one smooth curve, hence the name QuinticBezierSet.
In the special case when this class is being used to generate and evaluate 2D Bezier curves, that is x(u) and y(u), there are also functions to evaluate y(x), the first six derivatives of y(x), and also the first integral of y(x).
This class was not designed to be a stand alone Quintic Bezier class, but rather was developed out of necessity to model muscles. I required curves that, when linearly extrapolated, were C2 continuous, and by necessity I had to use quintic Bezier curves. In addition, the curves I was developing were functions in x,y space, allowing many of the methods (such as the evaluation of y(x) given that x(u) and y(u), the derivatives of y(x), and its first integral) to be developed, though in general this can't always be done.
I have parceled all of these tools into their own class so that others may more easily use and develop this starting point for their own means. I used the following text during the development of this code:
Mortenson, Michael E (2006). Geometric Modeling Third Edition. Industrial Press Inc., New York. Chapter 4 was quite helpful.
Future Upgrades
I think this is impossible because it is not possible, in general, to find the roots to a quintic polynomial, however, this fact may not preclude forming the inverse curve. The impossibility of finding the roots to a quintic polynomial was proven by Abel (Abel's Impossibility Theorem) and Galois
http://mathworld.wolfram.com/QuinticEquation.html
At the moment I am approximating the curve u(x) using cubic splines to return an approximate value for u(x), which I polish using Newton's method to the desired precision.
This is possible using the Divergence Theorem applied to 2D curves. A nice example application is shown in link 2 for computing the area of a closed cubic Bezier curve. While I have been able to get the simple examples to work, I have failed to successfully compute the area under a quintic Bezier curve correctly. I ran out of time trying to fix this problem, and so at the present time I am numerically computing the integral at a number of knot points and then evaluating the spline to compute the integral value.
a. http://en.wikipedia.org/wiki/Divergence_theorem b. http://objectmix.com/graphics/133553areaclosedbeziercurve.html
Currently the Bezier curve value and its derivative are computed separately to evaluate f and df in the Newton iteration to solve for u(x). Code optimized to compute both of these quantities at the same time could cut the cost of evaluating x(u) and dx/du in half. Since these quantities are evaluated in an iterative loop, this one change could yield substantial computational savings.
The function calcIndex computes which spline the contains the point of interest. This function has been implemented assuming a small number of Bezier curve sets, and so it simply linearly scans through the sets to determine the correct one to use. This function should be upgraded to use the bisection method if large quintic Bezier curve sets are desired.
I have exported all of the low level code as optimized code from Maple. Although the code produced using this means is reasonably fast, it is usually possible to obtain superior performance (and sometimes less round off error) by doing this work by hand.
Computational Cost Details All computational costs assume the following operation costs:
Operation Type : #flops +,,=,Boolean Op : 1 / : 10 sqrt: 20 trig: 40
These relative weightings will vary processor to processor, and so any of the quoted computational costs are approximate.

static 
Given a set of Bezier curve control points, return the index of the set of control points that x lies within.
x  A value that is interpolated by the set of Bezier curves 
bezierPtsX  A matrix of 6xn Bezier control points 
OpenSim::Exception  If the index is not located within this set of Bezier points 
Given a set of Bezier curve control points, return the index of the set of control points that x lies within. This function has been coded assuming a small number of Bezier curve sets (less than 10), and so, it simply scans through the Bezier curve sets until it finds the correct one.
Computational Costs Quoted for a Bezier curve set containing 1 to 5 curves.
~925
Example:

static 

static 
This function numerically integrates the Bezier curve y(x).
vX  Values of x to evaluate the integral of y(x) 
ic0  The initial value of the integral 
intAcc  Accuracy of the integrated solution 
uTol  Tolerance on the calculation of the intermediate u term 
uMaxIter  Maximum number of iterations allowed for u to reach its desired tolerance. 
mX  The 6xn matrix of Bezier control points for x(u) 
mY  The 6xn matrix of Bezier control points for y(u) 
aSplineUX  The array of spline objects that approximate u(x) on each Bezier interval 
flag_intLeftToRight  Setting this flag to true will evaluate the integral from the left most point to the right most point. Setting this flag to false will cause the integral to be evaluated from right to left. 
name  Name of caller. 
This function numerically integrates the Bezier curve y(x), when really both x and y are specified in terms of u. Evaluate the integral at the locations specified in vX and return the result.
Computational Costs
This the expense of this function depends on the number of points in vX, the points for which the integral y(x) must be calculated. The integral is evaluated using a RungeKutta 45 integrator, and so each point requires 6 function evaluations. (http://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method)
The cost of evaluating 1 Bezier curve y(x) scales with the number of points in xVal:
~1740 flops per point
The example below is quite involved, but just so it can show you an example of how to create the array of Spline objects that approximate the function u(x). Although the example has been created for only 1 Bezier curve set, simply changing the size and entries of the matrices _mX and _mY will allow multiple sets to be integrated.
Example:

static 
Calculates the location of quintic Bezier curve control points to create a C shaped curve.
x0  First intercept x location 
y0  First intercept y location 
dydx0  First intercept slope 
x1  Second intercept x location 
y1  Second intercept y location 
dydx1  Second intercept slope 
curviness  A parameter that ranges between 0 and 1 to denote a straight line or a curve 
OpenSim::Exception  If the curviness parameter is less than 0, or greater than 1; If the points and slopes are chosen so that an "S" shaped curve would be produced. This is tested by examining the points (x0,y0) and (x1,y1) together with the intersection (xC,yC) of the lines beginning at these points with slopes of dydx0 and dydx1 form a triangle. If the line segment from (x0,y0) to (x1,y1) is not the longest line segment, an exception is thrown. This is an overly conservative test as it prevents very deep 'V' shapes from being represented. 
Calculates the location of quintic Bezier curve control points to create a C shaped curve that intersects points 0 (x0, y0) and point 1 (x1, y1) with slopes dydx0 and dydx1 respectively, and a second derivative of 0. The curve that results can approximate a line (curviness = 0), or in a smooth C shaped curve (curviness = 1)
The current implementation of this function is not optimized in anyway and has the following costs:
Computational Costs
~55 flops
Example:

static 
Calculates the value of dydx of a quintic Bezier curve derivative at u.
u  The u value of interest. Note that u must be [0,1] 
xpts  The 6 control points associated with the x axis 
ypts  The 6 control points associated with the y axis 
order  The order of the derivative. Currently only orders from 16 can be evaluated 
OpenSim::Exception  If u is outside [0,1] If xpts is not 6 elements long If ypts is not 6 elements long If the order is less than 1 If the order is greater than 6 
The  value of (d^n y)/(dx^n) evaluated at u 
Calculates the value of dydx of a quintic Bezier curve derivative at u. Note that a 2D Bezier curve can have an infinite number of derivatives, because x and y are functions of u. Thus
dy/dx = (dy/du)/(dx/du) d^2y/dx^2 = d/du(dy/dx)*du/dx = [(d^2y/du^2)*(dx/du)  (dy/du)*(d^2x/du^2)]/(dx/du)^2 *(1/(dx/du)) etc.
Computational Costs
This obviously only functions when the Bezier curve in question has a finite derivative. Additionally, higher order derivatives are more numerically expensive to evaluate than lower order derivatives. For example, here are the number of operations required to compute the following derivatives
Name : flops dy/dx : ~102 d2y/dx2 : ~194 d3y/dx3 : ~321 d4y/dx4 : ~426 d5y/dx5 : ~564 d6y/dx6 : ~739
Example:

static 
Calculates the value of a quintic Bezier derivative curve at value u.
u  The independent variable of a Bezier curve, which ranges between 0.0 and 1.0. 
pts  The locations of the control points in 1 dimension. 
order  The desired order of the derivative. Order must be >= 1 
OpenSim::Exception  if u is outside [0,1] pts is not 6 elements long if order is less than 1 
Calculates the value of a quintic Bezier derivative curve at value u. This calculation is achieved by taking the derivative of the row vector uV and multiplying it by the 6x6 coefficient matrix associated with a quintic Bezier curve, by the vector of Bezier control points, pV, in a particular dimension.
Pseudo code for the first derivative (order == 1) would be
uV = [5*u^4 4*u^3 3*u^2 2u 1 0]; cM = [ 1 5 10 10 5 1; 5 20 30 20 5 0; 10 30 30 10 0 0; 10 20 10 0 0 0; 5 5 0 0 0 0; 1 0 0 0 0 0 ]; pV = [x1; x2; x3; x4; x5; x6]; dxdu = (uV*cM)*pV
Note that the derivative of uV only needed to be computed to compute dxdu. This process is continued for all 5 derivatives of x(u) until the sixth and all following derivatives, which are 0. Higher derivatives w.r.t. to U are less expensive to compute than lower derivatives.
Computational Costs
dy/dx : ~50 flops d2x/du2: ~43 flops d3x/du3: ~34 flops d4x/du4: ~26 flops d5x/du5: ~15 flops d6x/du6: ~1 flop
Example:

static 
Calculates the value of a quintic Bezier curve at value u.
u  The independent variable of a Bezier curve, which ranges between 0.0 and 1.0. 
pts  The locations of the control points in 1 dimension. 
OpenSim::Exception  If u is outside of [0,1] if pts has a length other than 6 
Calculates the value of a quintic Bezier curve at value u. This calculation is achieved by multiplying a row vector comprised of powers of u, by the 6x6 coefficient matrix associated with a quintic Bezier curve, by the vector of Bezier control points, pV, in a particular dimension. The code to compute the value of a quintic bezier curve has been optimized to have the following cost:
Computational Costs
~54 flops
The math this function executes is described in pseudo code as the following:
uV = [u^5 u^4 u^3 u^2 u 1]; cM = [ 1 5 10 10 5 1; 5 20 30 20 5 0; 10 30 30 10 0 0; 10 20 10 0 0 0; 5 5 0 0 0 0; 1 0 0 0 0 0 ]; pV = [x1; x2; x3; x4; x5; x6]; xB = (uV*cM)*pV
Example:

static 
This function will compute the u value that corresponds to the given x for a quintic Bezier curve.
ax  The x value 
bezierPtsX  The 6 Bezier point values 
splineUX  The spline for the approximate x(u) curve 
tol  The desired tolerance on u. 
maxIter  The maximum number of Newton iterations allowed 
OpenSim::Exception  if ax is outside the range defined in this Bezier spline section if the desired tolerance is not met if the derivative goes to 0 to machine precision 
This function will compute the u value that corresponds to the given x for a quintic Bezier curve. This is accomplished by using an approximate spline inverse of u(x) to get a very good initial guess, and then one or two Newton iterations to polish the answer to the desired tolerance.
Computational Costs
~219 flops
Example: