C-3 样条曲线

C-3 样条曲线

插值三次样条函数

• 曲线的曲率$k(x)=\frac{y^{\prime \prime }}{(1+y^{\prime 2}{)}^{\frac{3}{2}}}$
• 三次样条函数的定义：设在区间[a,b]上给定的一个分割$\Delta :a=x_0<x_1...<x_{n-1}<x_n=b$，若满足以下条件，则称[a,b]上的一个函数$S(x)$是三次样条函数 。
  • 在每个小区间$[x_{i-1},x_i](i=1,2,3,...,n)$内$S(x)$是三次多项式
  • 在整个区间$[a,b]$上，$S(x)$为二阶连续可导的函数
  • $S(X_i)=y_i,(i=1,2,3,...,n)$则称$S(X)$为三次样条函数

关于样条函数可以这么想，想象一个木板上依次钉着一排钉子，其中钉子的纵坐标各不相同，横坐标依次排开(相当于木板是区间，钉子是分割)。然后将一个柔软的线依次按照钉子的奇偶上下绕过这些钉子。然后一用力拉紧这条线，这条线就会紧紧的绷在这些钉子之间，和每一个钉子也都有接触(相当于一个函数经过了这些点)，那么这根线就是样条函数。很容易想象，只是用柔软的线的话，该样条函数是C0连续的，也就是说在连接处只是函数值相同，导数并不相同。三次样条函数就是穿过这些钉子，但是保持着C2连续条件的函数。

• 样条函数的构造

  • 已经有了离散化的数据，如何去构造一个函数去连接它们，这一过程就是插值
  • 一般用插值的方式去构造样条函数
  • 常用的插值方法一般是线性插值方法，牛顿插值方法和Lagrange插值方法

• 线性插值方法

  • 就是用直线连接各个点，可以发现，上面举的直线和钉子的例子就是线性插值方法。

    当然线性插值有很多改进的版本，比如计算一下两个点的中点然后进行插值，还有一些其他更复杂的方法，数值计算方法这门课里有相关内容。线性插值是一种非常重要的插值方法，图形学很多地方都用的是线性插值，比如计算三角形内的颜色，双线性插值，三线性插值等等。因为很多情况下并不需要很高的连续性，或者在点足够密集的情况下，单纯将它们用直线连接起来反而是一种高效可行的办法。

• 牛顿插值方法

  • 其实这就是牛顿迭代法，其收敛速度是很快的，有一个刻画收敛速度的数学描述，我记得牛顿迭代法是$n^2$收敛，二阶收敛。
  • 计算差商： 首先，根据给定的数据点$(x_0,y_0),(x_1,y_1),...,(x_n,y_n)$，计算差商。差商是通过递归定义的，具体而言，$f[x_i,x_{i+1},...,x_{i+k}]$表示通过点 $x_i,x_{i+1},...,x_{i+k}$ 的差商，其中 $k$ 是差商的阶数。对于二阶差商，它的计算公式为：$f[x_i,x_{i+1}]=\frac{f(x_{i+1})-f(x_i)}{x_{i+1}-x_i}$，对于更高阶的差商，可以通过递归公式来计算。
  • 构建插值多项式： 一旦计算出所有的差商，可以使用这些差商构建插值多项式。牛顿插值多项式通常写成如下形式：$P(X)=f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+...+$$f[x_0,x_1,..x_n](x-x_0)(x-x_1)...(x-x_{n-1})$其中，$f[x_0,x_1,..x_i]$ 是通过数据点 $x_0,x_1,...,x_i$的差商，$f(x_0)$是插值的起点。

• Lagrange插值方法

  • 拉格朗日插值方法是Hermit插值方法的一种特例，直接采用拉格朗日基函数进行插值即可，给定的点还是$(x_0,y_0),(x_1,y_1),...,(x_n,y_n)$。

  • 基函数$L_i(x)={\prod }_{j=0,j\ne i}^n\frac{x-x_j}{x_i-x_j}$

  • 插值函数$P(x)={\sum }_i^n{y}_i\cdot L_i(x)$

  • 举例

    x	0	1	2
    y	5	3	9
    $L_i(x)$	$\frac{(x-1)(x-2)}{(0-1)(0-2)}=\frac{(x-1)(x-2)}{2}$	$\frac{(x-1)(x-2)}{(1-0)(1-2)}=-\frac{(x-0)(x-2)}{2}$	$\frac{(x-0)(x-1)}{(2-0)(2-1)}=\frac{(x-0)(x-1)}{2}$

​ $P(X)=5\frac{(x-1)(x-2)}{2}-3\frac{(x-0)(x-2)}{2}+9\frac{(x-0)(x-1)}{2}=4x^2-6x+5$

B样条函数

• 以贝塞尔基函数为插值函数的样条，因此称为B样条曲线

代码

• 用二次贝塞尔函数进行点的插值，这里没有做点移动后的C2条件的光滑处理
• 环境： Imgui、Implot、GLFW 、glad

void Node::updataContoral() {

	if (x.size() == 2&& x.size() > controlx.size() + 1) {
		controlx.push_back((x[0] + x[1]) * 1 / 2);
		controly.push_back((y[0] + y[1]) * 3 / 5);
	}
	//插入控制点
	if (x.size() >= 3&& x.size()>controlx.size() + 1) {
		controlx.push_back(2 * x[x.size() - 2] - controlx[controlx.size() - 1]);
		controly.push_back(2 * y[x.size() - 2] - controly[controly.size() - 1]);	
	}
}

void Node::createBesselNode(){

	updataContoral();
	BesselNode_x.clear();
	BesselNode_y.clear();
	//计算每一段的贝塞尔函数
	for (int i = 0; i < x.size() -1; i++) {
		Bessel(x[i], y[i], controlx[i], controly[i], x[i + 1], y[i + 1]);
	}
}

参考资料

• 某一本不知出处的《计算几何》 链接：https://pan.baidu.com/s/1eOHipTyQ2A_KmCkOa7687g 提取码：6666

• Games102的几何建模课程