三次样条插值


三次样条函数

对于给定\(\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);


评论
 上一篇
简单细分曲线 简单细分曲线
本文介绍了简单细分曲线,特别关注于Chaikin方法生成的二次B-spline曲线细分。通过点分裂成边的拓扑规则,文章阐述了如何从一组有序点构成的多边形中逼近生成光滑曲线。
2020-12-23
下一篇 
多维放缩MDS算法 多维放缩MDS算法
在线性代数中我们接触到embedding(嵌入)这一术语,表示从高维到低维空间。考虑一个高维空间中的样本集在低维空间中的形态。 “embedding”也应用在计算机图形学的参数化领域中。
2020-11-09