本文描述曲线参数化的具体内容，由于篇幅限制，会分几篇叙述。曲线参数化的目的是将曲线使用参数u表示，如直线表示为：
C(u)=P+u⋅D,u∈(−∞,∞)
P为起始点，D为方向向量，u为参数。然而很多曲线的参数化表示并不像这么简单，一般分为如下三种：
1. 点集插值
意味着获取参数对{ti,Pi}，ti对应点集中的唯一点ti,ti即是点Pi的“abscissa”。此方法很多，分为均匀和非均匀，没有具体区别。
2. 插值方法
主要分为，Language法、Hermite法，多项式法。
2.1 Lagrange
有一组控制点P0，P1…….，Pn，计算其参数坐标ti(取值范围0到1)，使其满足下式：
此步称之为“Parameterization of a set of points”,意为将点可参数坐标t对应起来，但是如何求取？此处的ti可以任意规整求得，“只要能对应即可”。如以下采取定间距（uniform interval），
其插值公式如下，主要任务是求解Pi前的系数。有两种形式，基本式，迭代形式（recursive form）和矩阵形式(matrix form)。
Language插值推导：其基本思想是，假设过每个点的曲线是二次曲线，此二次曲线在本点取系数为1（过此点），在其余点取值为0，即上式中的Φ，称之为某个点的权值。然后由过所有点的二次曲线*权值，累加，表示这条曲线本身。如果有n+1个点，权值系数就是n次多项式，称为n次插值。所有的插值方法都有共同特点，高阶插值计算不稳定且计算量大。故实际使用时，使用低阶差值，然后做曲线拼接。
Language插值方法的矩阵形式：由插值公式可以看出，对于特定的n次Language插值，其系数是确定的。系数的具体形式取决于参数坐标ti的选取，如下述2阶（3个控制点）Language曲线插值公式：
故可以使用如下形式的矩阵相乘表示，
对三个点使用uniform Interpolation取得的参数坐标ti为{0,0.5,1},故曲线表示为，
故其M矩阵为，实际使用中，使用矩阵形式这个是很方便的。
以下采用按基本式实现Language插值，
static void ParaPoint(int n, double* T) { //uniform 计算点参数 double tmin = 0.0; double tmax = 1.0;
T[0] = tmin; T[n] = tmax; for (int i = 1; i < n; i++) { T[i] = tmin + (double(i) / double(n)) * (tmax - tmin); }}static double calculateCoeffi(int i, int n, double *T, double u) { double val = 1.0; for (int j = 0; j <= n; j++) { if (j != i) val *= (u - T[j]) / (T[i] - T[j]); }
return val;}//Language插值//返回n个控制点，参数为u的曲线上的点，u[0,1]static MyPoint interpolateCurve(int n, MyPoint* points,double u) { //计算点的参数化：首先要将离散点参数化，使ti对应于Pi。此处使用最简单的uniform parameterization double* T = new double[n+1]; ParaPoint(n,T); //计算权值 double* weights = new double[n+1]; for (int i = 0; i <= n; i++) { weights[i] = calculateCoeffi(i,n,T,u); } //计算最终的点 double x = 0.0; double y = 0.0; double z=0.0; for (int i = 0; i <= n; i++) { x += points[i].getX()*weights[i]; y += points[i].getY()*weights[i]; z += points[i].getZ()*weights[i]; } return MyPoint(x, y, z);}//计算过5个点的曲线const int n = 4;MyPoint points[n+1] = { MyPoint(0,0,0),MyPoint(10,-10,0),MyPoint(20,-5,0),MyPoint(30,20,0),MyPoint(60,20,0) };//获得此曲线上的n个点，并显示const int pnts = 100000;std::vector<MyPoint> outPoints;for (int i = 0; i <= pnts; i++) {//MyPoint p1 = interpolateCurve(2, points, i / 100);outPoints.push_back(interpolateCurve(n, points, double(i)/ pnts));}
2.2 Hermite
同时考虑考虑控制点和控制点处的斜率。每个控制点处的斜率难以控制，故先不描述。
2.3 控制多边形
此类方法的特点是，曲线不通过控制点，曲线的形状由控制点连线组成的多边形，及每个点的权系数决定。其要求为曲线“像”控制多边形，且光滑。
Bezier曲线，描述不变，由n+1个控制点，形成n阶Bezier曲线。其形式如下，
点Pi前的系数，称之为Bernstein多项式，表示为，
插值的主要任务是求解Bernstein多项式，可以直接求解。
de Casteljau Algorithm：使用上述Bernstein多项式多项式直接求解的缺点是，其不够稳定。实际中多使用de Casteljau Algorithm递归算法，此算法的几何意义更明确，且可以手工作图。其基本思想是，2点之间多次线性插值。如2个控制点，则就是直线本身。3个控制点，3次线性插值。如，
需要说明的是，de Casteljau与Bernstein形式其实是统一的。看如下3阶Bezier曲线的两种形式推导，
可以将3阶Bezier曲线的插值系数表示为矩阵形式，
使用Bernstein形式，实现Bezier曲线，
//阶乘double Factorial(int n){ if (n <= 1) return 1; return n*Factorial(n-1);}
//Bezier曲线static MyPoint interpolateBezierCurve(int n, MyPoint* points, double u) { //Bezier不需要将点参数化， //计算点的参数化：首先要将离散点参数化，使ti对应于Pi。此处使用最简单的uniform parameterization //double* T = new double[n + 1]; //ParaPoint(n, T); //计算权值 double* B = new double[n + 1]; for (int i = 0; i <= n; i++) { B[i] = Factorial(n) / (Factorial(i)*Factorial(n - i))*std::pow(u, i)*std::pow(1.0f - u, n - i); } //计算最终的点 double x = 0.0; double y = 0.0; double z = 0.0; for (int i = 0; i <= n; i++) { x += points[i].getX()*B[i]; y += points[i].getY()*B[i]; z += points[i].getZ()*B[i]; }
return MyPoint(x, y, z);}计算过5个点的曲线MyPoint points[n+1] = { MyPoint(0,0,0),MyPoint(10,-10,0),MyPoint(20,-5,0),MyPoint(30,20,0),MyPoint(60,20,0) };//获得此曲线上的n个点，并显示const int pnts = 100000;std::vector<MyPoint> outPoints;for (int i = 0; i <= pnts; i++) {//MyPoint p1 = interpolateCurve(2, points, i / 100);outPoints.push_back(interpolateBezierCurve(n, points, double(i)/ pnts));}
篇幅问题，后续文章继续。
路漫漫其修远兮，共勉。
参考：
Mesh generation application to Finite Element method
联系客服
