三次样条插值

曲线曲面

发布日期: 2020-11-10

文章字数: 1.1k

阅读时长: 5 分

三次样条函数

对于给定\(\mathbb{R}^2\)域内一组有序点\(\{\mathbf{Q}_i\}_{i=0}^n\)，对应参数为\(\{t_i\}_{i=0}^n\)，找到一组分段三次多项式函数插值每个数据点。

考虑插值点处\(C^2\)连续（平滑顶点）

对于\(n+1\)个数据点，相邻两点之间用三次多项式表达，则有\(n\)段三次多项式，用待定系数的方式设第\(i(i=0,1,\dots,n)\)段向量型多项式函数为 \[\mathbf{P}_i(t)=a_{i,0}+a_{i,1}t+a_{i,2}t^2+a_{i,3}t^3\] 从上式可以看出总共\(4n\)个未知数，\(C^2\)连续需要满足约束条件：

每段多项式在端点处插值 \[\mathbf{P}_i(t_i)=\mathbf{Q}_i\qquad \mathbf{P}_i(t_{i+1})=\mathbf{Q}_{i+1}\qquad i = 0,1,\dots,n\]
中间点一阶导连续 \[\mathbf{P}'_{i-1}(t_i-0)=\mathbf{P}'_i(t_i+0)\qquad i=1,2,\dots,n-1\]
中间点二阶导连续 \[\mathbf{P}''_{i-1}(t_i-0)=\mathbf{P}''_i(t_i+0)\qquad i=1,2,\dots,n-1\]

这里总共\(4n-2\)个方程，还需要两个约束条件，这里我们选择首尾两点二阶导等于0，即自然边界条件： \[\mathbf{P}''_0(t_0)=\mathbf{P}''_{n-1}(t_n)=0\qquad i=1,2,\dots,n-1\] 这里我们引入中间变量弯矩，参考此文，可得每段多项式的表达式： \[\mathbf{P}_i(t)=\frac{M_i}{6h_i}(t_{i+1}-t)^3+\frac{M_{i+1}}{6h_i}(t-t_i)^3+\left(\frac{Q_{i+1}}{h_i}-\frac{M_{i+1}h_i}{6}\right)(t-t_i)+\left(\frac{\mathbf{Q}_i}{h_i}-\frac{M_ih_i}{6}\right)(t_{i+1}-t)\] 其中\(h_i=t_{i+1}-t_i\)，\(M\)是如下三对角线性系统的解，\(M_0=M_n=0\) \[\begin{equation} \left( \begin{array}{ccccc} u_1 & h_1 & & &\\ h_1 & u_2 & h_2 & &\\ &\ddots&\ddots&\ddots&\\ &&h_{n-3}&u_{n-2}&h_{n-2}\\ &&&h_{n-2}&u_{n-1} \end{array} \right)\left( \begin{array}{c} M_1\\M_2\\\vdots\\M_{n-1} \end{array} \right)=\left( \begin{array}{c} v_1\\v_2\\\vdots\\v_{n-1} \end{array} \right) \end{equation}\] 其中\(u_i=2(h_i+h_{i-1})\)，\(v_i=\frac{6}{h_i}(\mathbf{Q}_{i+1}-\mathbf{Q}_i)-\frac{6}{h_{i-1}}(\mathbf{Q}_i-\mathbf{Q}_{i-1})\) 平滑顶点绘制（全局C^2连续）

考虑插值点处\(C^0\)连续（角部顶点）

对于每一段给定插值点处的左导数\(f'_-\)和右导数\(f'_+\)，加上插值该段的两点，4个未知数4个方程可以求解出该段的三次多项式函数：

例如在\([t_i,t_{i+1}]\)区间： \[\begin{equation} \left\{ \begin{array}{ccccccccc} a_{i,0} & + & a_{i,1}t_i & + & a_{i,2}t_i^2 & + & a_{i,2}t_i^3 & = & \mathbf{Q}_i\\ a_{i,0} & + & a_{i,1}t_{i+1} & + & a_{i,2}t_{i+1}^2 & + & a_{i,2}t_{i+1}^3 & = & \mathbf{Q}_{i+1}\\ &&a_{i,1}&+&2a_{i,2}t_i&+&3a_{i,3}t_i^2&=&\mathbf{P}'_i(t_i+0)\\ &&a_{i,1}&+&2a_{i,2}t_{i+1}&+&3a_{i,3}t_{i+1}^2&=&\mathbf{P}'_{i+1}(t_i-0) \end{array} \right. \end{equation}\] 求解该式可得到每段的三次多项式。角部顶点绘制（型值点处为C^0连续，其他点C^2连续）

考虑插值点处\(G^1\)连续（直线点）

在角部顶点基础上保证一点的左导数与右导数方向在同一直线即可。即上面方程中\(\mathbf{P}'_i(t_i+0)=\lambda\mathbf{P}'_{i+1}(t_i-0)\) 角部顶点绘制（型值点处为G^1连续，其他点C^2连续）

主要代码

高斯-塞德尔迭代求解三对角：

void GaussSeidel(std::vector<float> h, std::vector<float> u, std::vector<float> v, std::vector<float>* M) {
    std::vector<float> _M0 = *M; 
    while (true) {
        for (int i=0; i<u.size(); ++i) { 
            (*M)[i+1]=(v[i]-h[i]*(*M)[i]-h[i+1]*_M0[i+2])/u[i];
        } 
        // 迭代停止条件
        if ((Eigen::VectorXf::Map((*M).data(),(*M).size())-Eigen::VectorXf::Map(_M0.data(),_M0.size())).norm()<EPSILON) {
            break;
        }
        _M0 = *M;
    }
}

求解给定导数后的分段样条：

if (validDerivative) {
    float x, y;
    auto A = [=](int i) -> Eigen::Matrix4f {
        Eigen::Matrix4f a;
        a << 1, t[i], t[i] * t[i], t[i] * t[i] * t[i],
            1, t[i+1], t[i+1]*t[i+1], t[i+1]*t[i+1]*t[i+1], 0, 1, 2 * t[i], 3 * t[i] * t[i],
        0, 1, 2 * t[i+1], 3*t[i+1]*t[i+1];
        return a; 
    };
    auto B = [=](int i, int xy) -> Eigen::Vector4f { 
        Eigen::Vector4f b;
        b << points[i][xy],
        points[i + 1][xy], 
        (*derivative)[i].second[xy],
        (*derivative)[i + 1].first[xy];
        return b; 
    };
    // 四维矩阵可用Eigen，不会影响交互延迟
    Eigen::Vector4f ax = A(segment_index).colPivHouseholderQr().solve(B(segment_index, 0));
    Eigen::Vector4f ay = A(segment_index).colPivHouseholderQr ().solve(B(segment_index, 1));
    x = ax[0]+ax[1]*t0+ax[2]*t0*t0+ax[3]*t0*t0*t0; 
    y = ay[0]+ay[1]*t0+ay[2]*t0*t0+ay[3]*t0*t0*t0; 
    return Ubpa::pointf2(x, y);