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...

Detailed Description

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

Analytical Inverse to x(u).

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.

Analytical Integral of y(x)

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/133553-area-closed-bezier-curve.html

calcU

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.

calcIndex

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.

The addition of additional Bezier control point design algorithms, to create 'S' shaped curves, and possibly do subdivision.
Low level Code Optimization

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.

Author: Matt Millard

Version: 0.0

Member Function Documentation

◆ calcIndex() [1/2]

static int OpenSim::SegmentedQuinticBezierToolkit::calcIndex	(	double	x,
		const SimTK::Array_< SimTK::Vector > &	bezierPtsX
	)

static

◆ calcIndex() [2/2]

static int OpenSim::SegmentedQuinticBezierToolkit::calcIndex	(	double	x,
		const SimTK::Matrix &	bezierPtsX
	)

static

Given a set of Bezier curve control points, return the index of the set of control points that x lies within.

Parameters

x	A value that is interpolated by the set of Bezier curves
bezierPtsX	A matrix of 6xn Bezier control points

Exceptions

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.

    ~9-25

Example:

SimTK::Matrix mX(6,2);
 
//The first set of spline points
    mX(0,0) = -2;
    mX(1,0) = -1;
    mX(2,0) = -1;
    mX(3,0) = 1;
    mX(4,0) = 1;
    mX(5,0) = 2;
//The second set of spline points
    mX(0,1) = 2;
    mX(1,1) = 3;
    mX(2,1) = 3;
    mX(3,1) = 5;
    mX(4,1) = 5;
    mX(5,1) = 6;
 
//The value of x for which we want the index for
double xVal = 1.75;
int idx  = SegmentedQuinticBezierToolkit::calcIndex(xVal,mX);

◆ calcNumIntBezierYfcnX()

static SimTK::Matrix OpenSim::SegmentedQuinticBezierToolkit::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
	)

static

This function numerically integrates the Bezier curve y(x).

Parameters

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.

Returns: SimTK::Matrix Col 0: X vector, Col 1: int(y(x))

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 Runge-Kutta 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:

//Integrator and u tolerance settings
double INTTOL = 1e-12;
double UTOL   = 1e-14;
int    MAXITER= 10;
 
//Make up a Bezier curve - these happen to be the control points
//for a tendon curve
SimTK::Matrix _mX(6,1), _mY(6,1);
_mX(0)= 1;
_mX(1)= 1.01164;
_mX(2)= 1.01164;
_mX(3)= 1.02364;
_mX(4)= 1.02364;
_mX(5)= 1.04;
 
_mY(0) = 0;
_mY(1) = 3.10862e-16;
_mY(2) = 3.10862e-16;
_mY(3) = 0.3;
_mY(4) = 0.3;
_mY(5) = 1;
 
_numBezierSections = 1;
bool _intx0x1 = true; //Says we're integrating from x0 to x1
//Generate the set of splines that approximate u(x)
SimTK::Vector u(NUM_SAMPLE_PTS); //Used for the approximate inverse
SimTK::Vector x(NUM_SAMPLE_PTS); //Used for the approximate inverse
 
//Used to generate the set of knot points of the integral of y(x)    
SimTK::Vector xALL(NUM_SAMPLE_PTS*_numBezierSections-(_numBezierSections-1));
_arraySplineUX.resize(_numBezierSections);
int xidx = 0;
 
for(int s=0; s < _numBezierSections; s++){
    //Sample the local set for u and x
    for(int i=0;i<NUM_SAMPLE_PTS;i++){
        u(i) = ( (double)i )/( (double)(NUM_SAMPLE_PTS-1) );
        x(i) = SegmentedQuinticBezierToolkit::
            calcQuinticBezierCurveVal(u(i),_mX(s),_name);            
        if(_numBezierSections > 1){
           //Skip the last point of a set that has another set of points
           //after it. Why? The last point and the starting point of the
           //next set are identical in value.
            if(i<(NUM_SAMPLE_PTS-1) || s == (_numBezierSections-1)){
                xALL(xidx) = x(i);
                xidx++;
            }
        }else{
            xALL(xidx) = x(i);                
            xidx++;
        }
    }
    //Create the array of approximate inverses for u(x)    
    _arraySplineUX[s] = SimTK::SplineFitter<Real>::
        fitForSmoothingParameter(3,x,u,0).getSpline();
}
 
 
//Compute the integral of y(x) and spline the result    
 
SimTK::Vector yInt =  SegmentedQuinticBezierToolkit::
    calcNumIntBezierYfcnX(xALL,0,INTTOL, UTOL, MAXITER,_mX, _mY,
    _arraySplineUX,_name);
 
if(_intx0x1==false){
    yInt = yInt*-1;
    yInt = yInt - yInt(yInt.nelt()-1);
}
 
_splineYintX = SimTK::SplineFitter<Real>::
        fitForSmoothingParameter(3,xALL,yInt,0).getSpline();

◆ calcQuinticBezierCornerControlPoints()

static SimTK::Matrix OpenSim::SegmentedQuinticBezierToolkit::calcQuinticBezierCornerControlPoints	(	double	x0,
		double	y0,
		double	dydx0,
		double	x1,
		double	y1,
		double	dydx1,
		double	curviness
	)

static

Calculates the location of quintic Bezier curve control points to create a C shaped curve.

Parameters

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

Exceptions

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.

Returns: a SimTK::Matrix of 6 points Matrix(6,2) that correspond to the X, and Y control points for a quintic Bezier curve that has the above properties

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:

double x0 = 1;
double y0 = 0;
double dydx0 = 0;
double x1 = 1.04;
double y1 = 1;
double dydx1 = 43;
double c = 0.75;
 
SimTK::Matrix p0 = SegmentedQuinticBezierToolkit::
   calcQuinticBezierCornerControlPoints(x0, y0, dydx0,x1,y1,dydx01,
                                                         c);

◆ calcQuinticBezierCurveDerivDYDX()

static double OpenSim::SegmentedQuinticBezierToolkit::calcQuinticBezierCurveDerivDYDX	(	double	u,
		const SimTK::Vector &	xpts,
		const SimTK::Vector &	ypts,
		int	order
	)

static

Calculates the value of dydx of a quintic Bezier curve derivative at u.

Parameters

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 1-6 can be evaluated

Exceptions

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

Return values

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:

SimTK::Vector vX(6), vY(6);
 
double u = 0.5;
 
vX(0) = 1;
vX(1) = 1.01164;
vX(2) = 1.01164;
vX(3) = 1.02364;
vX(4) = 1.02364;
vY(5) = 1.04;
 
vY(0) = 0;
vY(1) = 3e-16;
vY(2) = 3e-16;
vY(3) = 0.3;
vY(4) = 0.3;
vY(5) = 1;
 
 
d2ydx2 = SegmentedQuinticBezierToolkit::calcQuinticBezierCurveDerivDYDX(
         u,vX, vY, 2);

◆ calcQuinticBezierCurveDerivU()

static double OpenSim::SegmentedQuinticBezierToolkit::calcQuinticBezierCurveDerivU	(	double	u,
		const SimTK::Vector &	pts,
		int	order
	)

static

Calculates the value of a quintic Bezier derivative curve at value u.

Parameters

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

Exceptions

OpenSim::Exception if -u is outside [0,1] -pts is not 6 elements long -if order is less than 1

Returns: The value of du/dx of Bezier curve located at u.

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:

double u = 0.5;
 
//Choose the control points
SimTK::Vector vX(6);
vX(0) = -2;
vX(1) = 0;
vX(2) = 0;
vX(3) = 4;
vX(4) = 4;
vX(5) = 6;
 
double dxdu  =calcQuinticBezierCurveDerivU(u,vX,1);

◆ calcQuinticBezierCurveVal()

static double OpenSim::SegmentedQuinticBezierToolkit::calcQuinticBezierCurveVal	(	double	u,
		const SimTK::Vector &	pts
	)

static

Calculates the value of a quintic Bezier curve at value u.

Parameters

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.

Exceptions

OpenSim::Exception -If u is outside of [0,1] -if pts has a length other than 6

Returns: The value of the Bezier curve located at u.

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:

double u = 0.5;
 
//Choose the control points
SimTK::Vector vX(6);
vX(0) = -2;
vX(1) = 0;
vX(2) = 0;
vX(3) = 4;
vX(4) = 4;
vX(5) = 6;
 
yVal = SegmentedQuinticBezierToolkit::
         calcQuinticBezierCurveVal(u,vX);

◆ calcU()

static double OpenSim::SegmentedQuinticBezierToolkit::calcU	(	double	ax,
		const SimTK::Vector &	bezierPtsX,
		const SimTK::Spline &	splineUX,
		double	tol,
		int	maxIter
	)

static

This function will compute the u value that corresponds to the given x for a quintic Bezier curve.

Parameters

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

Exceptions

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:

double xVal = 2;
 
//Choose the control points
SimTK::Vector vX(6);
vX(0) = -2;
vX(1) = 0;
vX(2) = 0;
vX(3) = 4;
vX(4) = 4;
vX(5) = 6;
 
SimTK::Vector x(100);
SimTK::Vector u(100);
 
//Create the splined approximate inverse of u(x)
for(int i=0; i<100; i++){
    u(i) = ( (double)i )/( (double)(100-1) );
    x(i) = SegmentedQuinticBezierToolkit::
        calcQuinticBezierCurveVal(u(i),vX); 
}
SimTK::Spline splineUX = SimTK::SplineFitter<Real>::
    fitForSmoothingParameter(3,x,u,0).getSpline();
 
//Now evaluate u at the given xVal
double u = SegmentedQuinticBezierToolkit::
         calcU(xVal,vX, splineUX, 1e-12,20);

The documentation for this class was generated from the following file:

OpenSim/Common/SegmentedQuinticBezierToolkit.h

Static Public Member Functions

Detailed Description

Member Function Documentation

◆ calcIndex() [1/2]

◆ calcIndex() [2/2]

◆ calcNumIntBezierYfcnX()

◆ calcQuinticBezierCornerControlPoints()

◆ calcQuinticBezierCurveDerivDYDX()

◆ calcQuinticBezierCurveDerivU()

◆ calcQuinticBezierCurveVal()

◆ calcU()