|
OpenSim
OpenSim 3.1
|
|
OpenSim
OpenSim 3.1
|
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...
#include <SegmentedQuinticBezierToolkit.h>
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 correesponds to the given x for a quintic Bezier curve. | |
| 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. | |
| 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. | |
| static double | calcQuinticBezierCurveDerivU (double u, const SimTK::Vector &pts, int order) |
| Calculates the value of a quintic Bezier derivative curve at value u. | |
| 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. | |
| 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. | |
| 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). | |
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 eachother 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 parcelled 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.
2. 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
3. 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 quantites 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.
4. 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.
5. The addition of additional Bezier control point design algorithms, to create 'S' shaped curves, and possibly do subdivision.
6. 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 usally possible to obtain superiour 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.
~9-25
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 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 matricies _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 respresented. |
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 (curviniess = 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 1-6 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 acheived 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 acheived by mulitplying 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 decribed 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 correesponds 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 correesponds 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:
1.8.1.2