等分3阶曲线的数值积分函数
按弧长等分
按 t 值线性分段
以上图形由 Mathematica 10 绘制, 分点坐标由下文的 C 代码中 PointOnCubicBezier 函数计算得到
Mathematica 10 代码
Needs["Splines`"]
pts = {{-31, -33}, {30, -36}, {-25, 28}, {33, 6}};
pps = {{-31.000000, -33.000000}, {-26.473223, -33.100813}, \
{-21.803112, -32.900533}, {-17.181859, -32.293437}, {-12.692476, \
-31.144602}, {-8.405203, -29.239359}, {-4.712460, -26.475325}, \
{-1.850744, -22.848809}, {0.068124, -18.613653}, {1.204749, \
-14.075212}, {1.833144, -9.500808}, {2.253196, -4.918588}, {2.785189, \
-0.247070}, {3.872591, 4.229112}, {6.350148, 8.141665}, {10.399117,
10.227684}, {15.030105, 10.588897}, {19.565813,
10.029625}, {24.141002, 8.952093}, {28.544624,
7.596929}, {33.000000, 6.000000}};
tps = {{-31.000000, -33.000000}, {-22.691375, -32.966625}, \
{-15.951000, -32.043000}, {-10.607125, -30.343875}, {-6.488000, \
-27.984000}, {-3.421875, -25.078125}, {-1.237000, -21.741000}, \
{0.238375, -18.087375}, {1.176000, -14.232000}, {1.747625, \
-10.289625}, {2.125000, -6.375000}, {2.479875, -2.602875}, {2.984000,
0.912000}, {3.809125, 4.054875}, {5.127000, 6.711000}, {7.109375,
8.765625}, {9.928000, 10.104000}, {13.754625,
10.611375}, {18.761000, 10.173000}, {25.118875,
8.674125}, {33.000000, 6.000000}};
Graphics[{Spline[pts, Bezier],
Point[pts], {RGBColor[1, 0, 0], Point[pps]}}]
Graphics[{Spline[pts, Bezier],
Point[pts], {RGBColor[0, 0, 1], Point[tps]}}]
函数
double get_t_of_BezierDegree3_Equal_parts(Point2D* cps, int i_of, int numOfParts)
获取曲线上 从 P0 起始, 到全长的 i_of / numOfParts 点处的 t 值. 仍可优化加速
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define SQR(x) ((x)*(x))
#define EPS (1e-3) /* EPS 积分步长, 一个小正数, 在运算范围内, 此值越小, 积分计算越准确, 但耗时也会越长 */
#define EPS_DIV_2 (EPS / 2)
#define ONE (999999e-6)
/*
3 阶曲线参数方程:
B(t) = P0 (1 - t)^3 + 3 P1 t (1 - t)^2 + 3 P2 t^2 (1 - t) + P3 t^3;
二维平面坐标分量参数方程:
x(t) = P0x (1 - t)^3 + 3 P1x t (1 - t)^2 + 3 P2x t^2 (1 - t) + P3x t^3;
y(t) = P0y (1 - t)^3 + 3 P1y t (1 - t)^2 + 3 P2y t^2 (1 - t) + P3y t^3;
分量微分 (x 分量微分 和 y 分量微分)
-3 P0x (1 - t)^2 + 3 P1x (1 - t)^2 - 6 P1x (1 - t) t + 6 P2x (1 - t) t - 3 P2x t^2 + 3 P3x t^2
-3 P0y (1 - t)^2 + 3 P1y (1 - t)^2 - 6 P1y (1 - t) t + 6 P2y (1 - t) t - 3 P2y t^2 + 3 P3y t^2
弧长微分 = Sqrt( x 分量微分 ^ 2 + y 分量微分 ^ 2 ) (原理为勾股定理, 如果是 3 维或更高维曲线, 则将所有 坐标分量微分的平方 求总和再开方)
3*Sqrt[(P0x*(-1 + t)^2 + t*(-2*P2x + 3*P2x*t - P3x*t) + P1x*(-1 + 4*t - 3*t^2))^2 + (P0y*(-1 + t)^2 + t*(-2*P2y + 3*P2y*t - P3y*t) + P1y*(-1 + 4*t - 3*t^2))^2]
*/
typedef struct {
double x;
double y;
}
Point2D;
// 3 阶曲线在 t 处的弧长微分
double BezierDegree3_dL(Point2D* cps, double t) {
// return 3*sqrt((cps[0].x * (-1 + t)^2 + t*(-2*cps[2].x + 3*cps[2].x * t - cps[3].x * t) + cps[1].x * (-1 + 4*t - 3*t^2))^2 +
// (cps[0].y * (-1 + t)^2 + t*(-2*cps[2].y + 3*cps[2].y * t - cps[3].y * t) + cps[1].y * (-1 + 4*t - 3*t^2))^2);
double v1, v2, v3, v4;
v1 = (-1 + t);
v2 = (-1 + 4 * t - 3 * t * t);
v3 = (cps[0].x * SQR(v1) + t * (-2 * cps[2].x + 3 * cps[2].x * t - cps[3].x * t) + cps[1].x * v2);
v4 = (cps[0].y * SQR(v1) + t * (-2 * cps[2].y + 3 * cps[2].y * t - cps[3].y * t) + cps[1].y * v2);
return 3 * sqrt(SQR(v3) + SQR(v4));
}
// 3 阶曲线的全长 (辛普森积分法)
// Integrate[f(t), {t, a, b}] ~= (f(a)+4*f(m)+f(b)) * (b-a) / 6 ; 其中 m = (a + b) / 2
// 此处积分步长 b-a = EPS
double BezierDegree3_Length(Point2D* cps) {
double a, m, b, bl;
bl = 0;
for (a = 0, m = a + EPS_DIV_2, b = a + EPS; a < ONE; a = b, b += EPS, m += EPS) {
bl += BezierDegree3_dL(cps, a) + 4 * BezierDegree3_dL(cps, m) + BezierDegree3_dL(cps, b);
}
return bl * EPS / 6;
}
// 把曲线等分成 numOfParts 段, 获取曲线上 从 P0 起始, 到全长的 i_of / numOfParts 点处的 t 值
// 例: 曲线被分成 numOfParts = 20 段, 从起点到终点一共 21 个分点;
// i_of = 3, 将返回第 3 个分点 Q (P0~~Q 的弧长 = 曲线全长 * 3 / 20) 处的 t 值
double get_t_of_BezierDegree3_Equal_parts(Point2D* cps, int i_of, int numOfParts) {
double a, m, b, bl, len;
if (i_of <= 0 || numOfParts <= 0) return 0;
if (i_of >= numOfParts) return 1;
// 要达到的长度
len = BezierDegree3_Length(cps) * i_of / numOfParts * 6 / EPS;
bl = 0;
for (a = 0, m = a + EPS_DIV_2, b = a + EPS; a < ONE; a = b, b += EPS, m += EPS) {
bl += BezierDegree3_dL(cps, a) + 4 * BezierDegree3_dL(cps, m) + BezierDegree3_dL(cps, b);
if (bl >= len) return a;
}
return 1;
}
// 返回由控制点集 cps 确定的曲线上, t 处的坐标, 结果存储在 pPoint 指向的二维坐标组中
void PointOnCubicBezier(Point2D* pPoint, Point2D* cps, double t) {
if (pPoint == NULL) return;
double v, vSQ, tSQ, c0, c1, c2, c3;
v = (1 - t);
vSQ = SQR(v);
tSQ = t * t;
c0 = vSQ * v; // (1 - t)^3
c1 = 3 * t * vSQ; // 3 t (1 - t)^2
c2 = 3 * tSQ * v; // 3 t^2 (1 - t)
c3 = tSQ * t; // t^3
// B(t) = P0 (1 - t)^3 + 3 P1 t (1 - t)^2 + 3 P2 t^2 (1 - t) + P3 t^3;
pPoint -> x = cps[0].x * c0 + cps[1].x * c1 + cps[2].x * c2 + cps[3].x * c3;
pPoint -> y = cps[0].y * c0 + cps[1].y * c1 + cps[2].y * c2 + cps[3].y * c3;
}
int main() {
// pts = {{-31, -33}, {30, -36}, {-25, 28}, {33, 6}};
Point2D cps[] = {
{-31, -33},
{30, -36},
{-25, 28},
{33, 6}
};
// 计算并输出曲线的全长
double bl;
bl = BezierDegree3_Length(cps);
printf("BezierDegree3_Length is %f\n", bl);
int i, numOfParts = 20;
double t;
Point2D pps[numOfParts + 1]; // 用于存储分点坐标
printf("弧长等分点\n");
for (i = 0; i <= numOfParts; i++) {
t = get_t_of_BezierDegree3_Equal_parts(cps, i, numOfParts);
PointOnCubicBezier(pps + i, cps, t);
printf("%d, %f, {%f, %f}\n", i, t, (pps + i)->x, (pps + i)->y);
}
printf("t 线性分点\n");
for (i = 0; i <= numOfParts; i++) {
t = (double) i / numOfParts;
PointOnCubicBezier(pps + i, cps, t);
printf("%d, %f, {%f, %f}\n", i, t, (pps + i)->x, (pps + i)->y);
}
return 0;
}
运行输出
BezierDegree3_Length is 92.896940
弧长等分点
0, 0.000000, {-31.000000, -33.000000}
1, 0.026000, {-26.473223, -33.100813}
2, 0.056000, {-21.803112, -32.900533}
3, 0.090000, {-17.181859, -32.293437}
4, 0.129000, {-12.692476, -31.144602}
5, 0.175000, {-8.405203, -29.239359}
6, 0.227000, {-4.712460, -26.475325}
7, 0.284000, {-1.850744, -22.848809}
8, 0.343000, {0.068124, -18.613653}
9, 0.402000, {1.204749, -14.075212}
10, 0.460000, {1.833144, -9.500808}
11, 0.519000, {2.253196, -4.918588}
12, 0.583000, {2.785189, -0.247070}
13, 0.653000, {3.872591, 4.229112}
14, 0.733000, {6.350148, 8.141665}
15, 0.807000, {10.399117, 10.227684}
16, 0.864000, {15.030105, 10.588897}
17, 0.907000, {19.565813, 10.029625}
18, 0.943000, {24.141002, 8.952093}
19, 0.973000, {28.544624, 7.596929}
20, 1.000000, {33.000000, 6.000000}
t 线性分点
0, 0.000000, {-31.000000, -33.000000}
1, 0.050000, {-22.691375, -32.966625}
2, 0.100000, {-15.951000, -32.043000}
3, 0.150000, {-10.607125, -30.343875}
4, 0.200000, {-6.488000, -27.984000}
5, 0.250000, {-3.421875, -25.078125}
6, 0.300000, {-1.237000, -21.741000}
7, 0.350000, {0.238375, -18.087375}
8, 0.400000, {1.176000, -14.232000}
9, 0.450000, {1.747625, -10.289625}
10, 0.500000, {2.125000, -6.375000}
11, 0.550000, {2.479875, -2.602875}
12, 0.600000, {2.984000, 0.912000}
13, 0.650000, {3.809125, 4.054875}
14, 0.700000, {5.127000, 6.711000}
15, 0.750000, {7.109375, 8.765625}
16, 0.800000, {9.928000, 10.104000}
17, 0.850000, {13.754625, 10.611375}
18, 0.900000, {18.761000, 10.173000}
19, 0.950000, {25.118875, 8.674125}
20, 1.000000, {33.000000, 6.000000}