计算几何02_三次样条曲线

一、样条

样条（Spline）函数是由舍恩伯格于1946年提出的。样条是富有弹性的细木条或有机玻璃条，它的作用相当于“万能”曲线板。早期船舶、汽车、飞机放样时用铅压铁压住样条，使其通过一系列型值点，调整压铁达到设计要求后绘制其曲线，称为样条曲线。这样设计曲线的方法在20世纪六七十年代得到了广泛应用。

二、几何连续性

2.1 连续性条件

通常单一的曲线段或曲面片难以表达复杂的形状，必须将一些曲线段拼接成组合曲线，或将一些曲面片拼接成组合曲面，才能表达复杂的形状。为了保证在结合点处光滑过渡，需要满足连续性条件。连续性条件有两种：参数连续性（Parametric Continuity）与几何连续性（Geometric Continuity）。

2.2 参数连续性

• 零阶参数连续性，记为C0，指相邻两段曲线在结合点处具有相同的坐标。
  

• 一阶参数连续性，记为C1，指相邻两段曲线在结合点处具有相同的一阶导数
  

• 二阶参数连续性，记为C2，指相邻两段曲线在结合点处具有相同的一阶导数和二阶导数。
  

2.3 几何连续性

与参数连续性不同的是，几何连续性只要求参数成比例而非相等。

• 零阶几何连续性，记为G0，与零阶参数连续性相同，即相邻两段曲线在结合点处有相同的坐标。
• 一阶几何连续性，记为G1，指相邻两段曲线在结合点处的一阶导数成比例，但大小不一定相等。
• 二阶几何连续性，记为G2，指相邻两段曲线在结合点处的一阶导数和二阶导数成比例，即曲率一致，但大小不一定相等。

三、三次样条曲线

• 三次样条曲线是组合曲线，它在相邻的两型值点之间，使用三次函数进行插值。如果把n段三次曲线连接起来，使两相邻曲线在连接点（称为节点）的斜率和曲率相等，就获得n段三次函数组合成的曲线，即三次样条曲线。
• 三次样条曲线是插值曲线，它通过所有型值点。
• 样条函数的理论和运用是从三次样条函数发展起来的。
• 在计算几何中，应用得最早、研究得最详细的也是三次样条函数，因为 三次样条函数在节点处具有C2连续性，而二阶连续性是大多数工程问题所需要的。
• 样条函数也是放样工艺中绘制曲线用的细木条的数学模型的线性近似，符合传统的光顺要求。

3.1 三次样条函数的定义

已知n个型值点Pi(xi, yi), i = 1, 2,…, n且a = x1 < x2 <…< xn = b。若 y=s(x)满足下列条件：
（1）型值点Pi在函数 y=s(x) 上。
（2） s(x) 在整个区间[a, b]上二阶连续可导。
（3）在每个子区间[xi, xi+1], i = 1, 2,…, n-1上，分段函数 都是参数x的三次多项式。则称函数 是过型值点的三次样条函数，由三次样条函数构成的曲线称为三次样条曲线。

3.2 三次样条函数的表达式

第i段的 si(x) 表示为:

式中，ai，bi，ci，di为待定系数，i=1，2 … ，n-1。

第i段曲线的首端通过Pi（xi，yi），末端通过Pi+1（xi+1，yi+1）
Pi处的二阶导数为Mi。
其中各型值点的弯矩Mi的力学解释如下：

子区间示意图如下：

用型值点处的二阶导数Mi表示的三次B样条函数为:

3.3 求解思路

3.4 函数推导

3.5 边界条件

3.5 推导过程

见个人资源：三次样条插值函数求解过程.pdf

四、算法实现

4.1 思路

（1）读入n个型值点且满足其x坐标递增。
（2）根据实际情况确定三次样条曲线的边界条件。
（3）计算曲线的系数，将其表达为型值点二阶导数的函数。
（4）用追赶法求解三弯矩方程。
（5）将求解出的参数代入三次样条函数的表达式中，构造三次样条曲线。
（6）循环访问每个节点。在每个子区间内，按照精度要求，使用直线段连接各段内的若干等分点，即可绘制出三次样条曲线。

4.2 代码

4.2.1 读取型值点

首先定义数据结构p2.cpp与头文件(主要如下)：

CP2::CP2(double x, double y)
{
	this->x = x;
	this->y = y;
}
1
2
3
4
5

之后在主函数中定义读取函数：

// CGeometricfiguretestViewmessage handlers
void CGeometricfiguretestView::ReadPoint()
{
	P[1].x = -340, P[1].y = -200;
	P[2].x = -150, P[2].y =  0;
	P[3].x =  0,   P[3].y =  -50;
	P[4].x =  100, P[4].y =  -100;
	P[5].x =  250, P[5].y =  -100;
	P[6].x =  350, P[6].y =  -50;
}
1
2
3
4
5
6
7
8
9
10

4.2.2 绘制型值点

void CGeometricfiguretestView::DrawDataPoint(CDC* pDC)//绘制型值点
{
	CBrush NewBrush, *pOldBrush;
	NewBrush.CreateSolidBrush(RGB(0, 0, 0));
	pOldBrush = pDC->SelectObject(&NewBrush);
	for( int i = 1; i < 7; i++)
		pDC->Ellipse(ROUND(P[i].x - 5), ROUND(P[i].y - 5), ROUND(P[i].x + 5), ROUND(P[i].y + 5));
	pDC->SelectObject(pOldBrush);
}
1
2
3
4
5
6
7
8
9

4.2.3 三次样条曲线绘制

具体过程见注释：

void CGeometricfiguretestView::DrawCubicSpline(CDC* pDC)//三次样条曲线
{
	int n = 6;
	const int dim = 7;//二维数组维数
	double b1 = 10, bn = 10;//边界条件:"夹持端",给出起点和终点的一阶导数
	double h[dim], lambda[dim], mu[dim], D[dim];
    double l[dim], m[dim], u[dim]; 
    double M[dim], K[dim];
    double a[dim], b[dim], c[dim], d[dim];
	for(int i = 1; i < n; i++)//计算hi＝xi+1-xi
    	h[i] = P[i+1].x - P[i].x;
	for(int i = 2; i < n; i++)
	{
		lambda[i] = h[i-1]/(h[i-1]+h[i]);//计算λ
		mu[i] = h[i]/(h[i-1]+h[i]);//计算μ
		D[i] = 6/(h[i-1]+h[i])*((P[i+1].y-P[i].y)/h[i]-(P[i].y-P[i-1].y)/h[i-1]);//计算D
    }
	D[1]=6*((P[2].y-P[1].y)/h[1]-b1)/h[1];//夹持端的D[1]
    D[n]=6*(bn-(P[n].y-P[n-1].y)/h[n-1])/h[n-1];//夹持端的D[n]
    mu[1]=1;//夹持端的μ[1],
	lambda[n]=1;//夹持端的λ[n]
	//追赶法求解三弯矩方程
    l[1]=2;
	u[1]=mu[1]/l[1];
    for(int i=2; i <= n; i++)
    {
		m[i]=lambda[i];
        l[i]=2-m[i]*u[i-1];
        u[i]=mu[i]/l[i];
    }
    K[1] = D[1]/l[1];//解LK=D
    for(int i = 2; i <= n;i++)
    {
		K[i]=(D[i]-m[i]*K[i-1])/l[i];
	}
	M[n] = K[n];//解UM=K
	for(int i = n-1; i >= 1;i--)
	{
		M[i]=K[i]-u[i]*M[i+1];
	}
	//计算三次样条函数的系数
	for(int i = 1; i < n; i++)
    {
		a[i] = P[i].y;
        b[i] = (P[i+1].y-P[i].y)/h[i] - h[i]*(M[i]/3+M[i+1]/6);
        c[i] = M[i]/2;
        d[i] = (M[i+1]-M[i])/(6*h[i]);
    }
	pDC->MoveTo(ROUND(P[1].x), ROUND(P[1].y));
	double xStep = 0.5;//x的步长
	double x, y;//当前点
	for(int i = 1; i < n; i++)//循环访问每个节点
    {
		for(x = P[i].x; x < P[i+1].x; x += xStep)//循环访问每个节点
		{
			y=a[i]+b[i]*(x-P[i].x)+c[i]*(x-P[i].x)*(x-P[i].x)+d[i]*(x-P[i].x)*(x-P[i].x)*(x-P[i].x);
			pDC->LineTo(ROUND(x), ROUND(y));//绘制样条曲线
        }
    }
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

4.2.4 主函数调用

void CGeometricfiguretestView::OnDraw(CDC* pDC)
{
	CTestDoc* pDoc = GetDocument();
	ASSERT_VALID(pDoc);
	if (!pDoc)
		return;
	// TODO: add draw code for native data here
	CRect rect;//定义客户区矩形
	GetClientRect(&rect);//获得客户区矩形的信息
	pDC->SetMapMode(MM_ANISOTROPIC);//自定义二维坐标系
	pDC->SetWindowExt(rect.Width(), rect.Height());//设置窗口范围
	pDC->SetViewportExt(rect.Width(), -rect.Height());//设置视区范围，且x轴水平向右为正，y轴垂直向上为正
	pDC->SetViewportOrg(rect.Width() / 2, rect.Height() / 2);//设置客户区中心为二维坐标系原点
	rect.OffsetRect(-rect.Width() / 2, -rect.Height() / 2);//rect矩形与客户区重合
	ReadPoint();
	DrawDataPoint(pDC);
	DrawCubicSpline(pDC);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

编译运行，可见如下：

我们也可以改变其中一个边界约束条件：

double b1 = 10, bn = -5;//边界条件:"夹持端",给出起点和终点的一阶导数
1

运行可以看到如下：