API  4.5
For C++ developers
List of all members | Static Public Member Functions | Nested Classes
OpenSim::SegmentedQuinticBezierToolkit Class Reference


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::Vec6 &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::Vec6 > &bezierPtsX)
 
static double calcQuinticBezierCurveVal (double u, const SimTK::Vec6 &pts)
 
Calculates the value of a quintic Bezier curve at value u. More...
 
static double calcQuinticBezierCurveDerivU (double u, const SimTK::Vec6 &pts, int order)
 
Calculates the value of a quintic Bezier derivative curve at value u. More...
 
static double calcQuinticBezierCurveDerivDYDX (double u, const SimTK::Vec6 &xpts, const SimTK::Vec6 &ypts, int order)
 
Calculates the value of dydx of a quintic Bezier curve derivative at u. More...
 
static ControlPointsXY 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::Array_< SimTK::Vec6 > &ctrlPtsX, const SimTK::Array_< SimTK::Vec6 > &ctrlPtsY, const SimTK::Array_< SimTK::Spline > &aSplineUX, bool flag_intLeftToRight, const std::string &name)
 
This function numerically integrates the Bezier curve y(x). More...
 

Nested Classes

struct  ControlPointsXY
 

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

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

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

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

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

  1. The addition of additional Bezier control point design algorithms, to create 'S' shaped curves, and possibly do subdivision.
  2. 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::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
xA value that is interpolated by the set of Bezier curves
bezierPtsXA vector of 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);

◆ calcIndex() [2/2]

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

◆ calcNumIntBezierYfcnX()

static SimTK::Matrix OpenSim::SegmentedQuinticBezierToolkit::calcNumIntBezierYfcnX ( const SimTK::Vector &  vX,
double  ic0,
double  intAcc,
double  uTol,
int  uMaxIter,
const SimTK::Array_< SimTK::Vec6 > &  ctrlPtsX,
const SimTK::Array_< SimTK::Vec6 > &  ctrlPtsY,
const SimTK::Array_< SimTK::Spline > &  aSplineUX,
bool  flag_intLeftToRight,
const std::string &  name 
)
static


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

Parameters
vXValues of x to evaluate the integral of y(x)
ic0The initial value of the integral
intAccAccuracy of the integrated solution
uTolTolerance on the calculation of the intermediate u term
uMaxIterMaximum number of iterations allowed for u to reach its desired tolerance.
ctrlPtsXAn n-vector of 6 Bezier control points for x(u)
ctrlPtsYAn n-vector of 6 Bezier control points for y(u)
aSplineUXThe array of spline objects that approximate u(x) on each Bezier interval
flag_intLeftToRightSetting 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.
nameName 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, adding more control points 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
std::array<SimTK::Vec6> _ctrlPtsX { SimTK::Vec6(
1,
1.01164,
1.01164,
1.02364,
1.02364,
1.04) };
std::array<SimTK::Vec6> _ctrlPtsY { SimTK::Vec6(
0,
3.10862e-16,
3.10862e-16,
0.3,
0.3,
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),_ctrlPtsX[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,_ctrlPtsX, _ctrlPtsY,
_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 ControlPointsXY 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
x0First intercept x location
y0First intercept y location
dydx0First intercept slope
x1Second intercept x location
y1Second intercept y location
dydx1Second intercept slope
curvinessA 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 pair of Vec6 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;
SegmentedQuinticBezierToolkit::ControlPointsXY pXY = SegmentedQuinticBezierToolkit::
calcQuinticBezierCornerControlPoints(x0, y0, dydx0,x1,y1,dydx01,
c);

◆ calcQuinticBezierCurveDerivDYDX()

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


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

Parameters
uThe u value of interest. Note that u must be [0,1]
xptsThe 6 control points associated with the x axis
yptsThe 6 control points associated with the y axis
orderThe order of the derivative. Currently only orders from 1-6 can be evaluated
Exceptions
OpenSim::Exception-If u is outside [0,1] -If the order is less than 1 -If the order is greater than 6
Return values
Thevalue 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:

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

◆ calcQuinticBezierCurveDerivU()

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


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

Parameters
uThe independent variable of a Bezier curve, which ranges between 0.0 and 1.0.
ptsThe locations of the control points in 1 dimension.
orderThe desired order of the derivative. Order must be >= 1
Exceptions
OpenSim::Exceptionif -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::Vec6 vX(-2, 0, 0, 4, 4, 6);
double dxdu = calcQuinticBezierCurveDerivU(u,vX,1);

◆ calcQuinticBezierCurveVal()

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


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

Parameters
uThe independent variable of a Bezier curve, which ranges between 0.0 and 1.0.
ptsThe locations of the control points in 1 dimension.
Exceptions
OpenSim::Exception-If u is outside of [0,1]
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::Vec6 vX(-2, 0, 0, 4, 4, 6);
yVal = SegmentedQuinticBezierToolkit::
calcQuinticBezierCurveVal(u,vX);

◆ calcU()

static double OpenSim::SegmentedQuinticBezierToolkit::calcU ( double  ax,
const SimTK::Vec6 &  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
axThe x value
bezierPtsXThe 6 Bezier point values
splineUXThe spline for the approximate x(u) curve
tolThe desired tolerance on u.
maxIterThe 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::Vec6 vX(-2, 0, 0, 4, 4, 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: