B-样条曲线

作者：追风剑情发布于：2019-4-16 12:17 分类：Algorithms

以 Bernstein 基函数构造的 Bézier 曲线有许多优点，但也有不足：一是控制多边形的顶点个数决定了 Bézier 曲线的阶次，即 n+1 个顶点的控制多边形必然会产生 n 次 Bézier 曲线，并且当 n 较大时，控制多边形对曲线的控制将会减弱。二是 Bézier 曲线不能做局部修改，即改变某一个控制点的位置对整条曲线都有影响。其原因主要是 Bernstein 基函数在整个开区间 (0，1) 的范围内均不为 0，所以曲线在开区间内任何一点的值均要受到全部顶点的影响，改变其中某一顶点的位置对整个曲线均有影响。

B样条方法保留了Bézier方法的优点，克服其由于整体表示带来的不具备局部性质的缺点，具有表示与设计自由型曲线曲面的强大功能。另外，B样条方法目前已成为关于工业产品几何定义国标标准的有理B样条方法的基础。因此，B样条方法是形状数学描述的主流法之一。关于B样条的理论早在1946年就由Schoenberg 提出，但论文直到1967年才发表。1972年 de Boor 与 Cox 分别独立地给出关于样条计算的标准算法。但该方法作为在 CAGD中的一个形状数学描述的基本方法，是由 Gordon 和 Riesenfeld 于1974年在研究 Bézier 方法的基础上引人的。他们拓广了Bézier曲线，用B样条基代替Bernstein基，从而改进了 Bézier 控制多边形与 Bernstein 多项式次数有关和整体逼近的弱点。

B样条曲线的定义

给定n+1个控制点P0,P1,...Pn，它们所确定的k阶B样条曲线是

式中，基函数N_i,k(u)递归定义如下：

从定义的公式可知，高阶曲线由两个低一阶的相同曲线(形状相同，位置不同，位置相距1个节点区间)通过线性插值的方式混合而成。

因为u的下标最大值为n+k，所以i+l的最大值为n+k，即i+l=n+k，得i=n+k-l。又因为u的下标最小值为0，所以$0 \le i \le n+k-l$。

式中，u₀，u₁，...，u_n+k是一个非递减的序列，称为节点；(u₀，u₁，...，u_n+k)称为节点向量。定义中可能出现0/0，这时约定为0。

节点向量（u₀，u₁，...，u_n+k）所包含的n+k个区间并非都在该曲线的定义域内，其中两端各k-1个节点区间，不能作为B样条曲线的定义区间。这是因为n+1个顶点中最前的k个顶点pi，(i=0，1，...，k-1)定义了样条曲线的首段曲线，其定义区间为u∈[u_k-1，u_k)；随后的k个顶点P_i，(i=1，2，...，k)定义了第二段曲线，其定义区间为u∈[u_k，u_k+1)；...；最后的k个顶点P_i，(i=n-k+1，n-k+2，...，n)定义了末段曲线，其定义区间为u∈[u_n，u_n+1]。于是，得到k阶B样条曲线的定义域为u∈[u_k-1，u_n+1]。共含有n-k+2个节点区间(包括零长度的节点区间)。若其中不含重节点，则对应B样条曲线包含n-k+2段。也可看到，节点向量两侧各k-1个节点区间上的那些B样条基函数因其权性不成立，不能构成基函数。

由k阶B样条曲线的递归定义可以看出：
1）对n+1个控制点，曲线由n+1个混合函数所描述。
2）每个混合函数N_i,_k(u)定义在u取值范围的k个子区间，以节点向量值u_i为起点。
3）参数u的取值范围由n+k+1个给定节点向量值分成n+k个子区间。
4）节点向量(u₀，u₁，...，u_n+k)所生成的B样条曲线仅定义在从节点值u_k-1到节点值u_n+1的区间上。
5）任一控制点可以影响最多k个曲线段的形状。
6）P(u)是分段参数多项式，P(u)在每一区间u∈[u_i，u_i+1]，(k-1≤i≤n)上都是次数不高于k-1的多项式。

从B样条曲线的这个递归定义可以看出，曲线与给定的阶数k及节点向量都有关系。就是说，即使k相同，选择不同的节点向量，也能得到不同的曲线。

任意的一阶B样条曲线就是控制点本身，可以看作是零次多项式。例如，n=2,k=1，控制点是P₀，P₁，P₂，这样应选择参数节点n+k+1=4个，设节点向量是(u₀，u₁，u₂，u₃)，按照定义，可写出三个基函数：

由上式可知所定义的B样条曲线是 $$ \begin{flalign} P(u)&=N_{0,1}(u)P_0+N_{1,1}(u)P_1+N_{2,1}(u)P_2 &\\ &= \left\{ \begin{aligned} P0 \quad u_0 \le u < u_1 \\ P1 \quad u_1 \le u < u_2 \\ P2 \quad u_2 \le u \le u_3 \\ \end{aligned} \right. \end{flalign} $$

二阶B样条$N_{i,2}(u)$由两个一阶B样条$N_{i,1}(u)$与$N_{i+1,1}(u)$递归推得，是它们的凸线性组合，即 $$ N_{i,2}(u)=\frac{u-u_i}{u_{i+1}-u_i}N_{i,1}(u)+\frac{u_{i+2}-u}{u_{i+2}-u_{i+1}}N_{i+1,1}(u) $$ 式中 $$ N_{i,1}(u)= \left\{ \begin{aligned} 1 \quad u \in [u_i, u_{i+1}) \\ 0 \quad u \notin [u_i, u_{i+1}) \\ \end{aligned} \right. $$ $$ N_{i+1,1}(u)= \left\{ \begin{aligned} 1 \quad u \in [u_{i+1}, u_{i+2}) \\ 0 \quad u \notin [u_{i+1}, u_{i+2}) \\ \end{aligned} \right. $$ 它们就像开关那样发生作用，即1表示接通，0表示断开。于是得到 $$ N_{i,2}(u)= \left\{ \begin{aligned} &\frac{u-u_i}{u_{i+1}-u_i} & u \in [u_{i}, u_{i+1}) \\ &\frac{u_{i+2}-u}{u_{i+2}-u_{i+1}} & u \in [u_{i+1}, u_{i+2}) \\ &0 & u \notin [u_{i}, u_{i+2}) \\ \end{aligned} \right. $$

节点向量为(0, 1, 2)的二阶B样条基函数如图4-23所示。如此继续下去，可由B样条$N_{i+1,1}(u)$与$N_{i+2,1}(u)$递归推得B样条$N_{i+1,2}(u)$。再由两个二阶B样条$N_{i,2}(u)$与$N_{i+1,2}(u)$进一步递归推得三阶B样条$N_{i,3}(u)$，如此继续下去可计算出其他的三阶B样条。四阶B样条的计算也以此类推。图4-24和图4-25分别给出了节点向量为(0,1,2,3)和(0,1,2,3,4)的三阶和四阶B样条。

递推公式表明，k阶B样条$N_{i,k}(u)$可由两个k-1阶B样条$N_{i,k-1}(u)$与$N_{i+1,k-1}(u)$递推得到。其凸线性组合的系数分别为$\frac{u-u_i}{u_{i+k-1}-u_i}$与$\frac{u_{i+k}-u}{u_{i+k}-u_{i+1}}$，两个系数的分母恰好是两个k-1阶B样条的支承区间，分子恰好是参数u把第i个k阶B样条$N_{i,k}(u)$的支承区间$[u_i, u_{i+k}]$划分成两部分的长度。

B样条基函数有下列性质：

（1）正性和局部性 $$ N_{i,k}(u) \left\{ \begin{aligned} &>0 & u_i < u < u_{i+k} \\ &=0 & \text{其它} \\ \end{aligned} \right. $$ 即$N_{i,k}(u)$在区间$(u_i, u_{i+k})$中为正，在其他地方$N_{i,k}(u)$为0。

（2）规范性 $$ 0 \le N_{i,k}(u) \le 1 $$

（3）权性
对从节点值$u_{k-1}$到$u_{n+1}$区间上的任一值u，全体基函数之和为1。 $$ \displaystyle\sum_{i=0}^{n}N_{i,k}(u)=1 $$

（4）递推性
由定义式表明。

借用网上的一张图片来理解控制点

P0，P1，P2，P3为控制点

借用网上的一张图来理解支撑区间与节点向量

[0, 2]、[1, 3]为两条二次样条曲线的支撑区间。(0, 1, 2, 3)为节点向量。

示例

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;

namespace WindowsFormsAppB
{
    public partial class Form1 : Form
    {
        public Form1()
        {
            InitializeComponent();
        }

        public void PaintControlPoint(PaintEventArgs e, BSpline.Point[] CP)
        {
            Graphics g = e.Graphics;
            //绘制控制点
            Pen dotPen = new Pen(Color.Red, 2);
            for (int i = 0; i < CP.Length; i++) {
                g.DrawEllipse(dotPen, (float)CP[i].x, (float)CP[i].y, 3, 3);
            }

            //绘制控制点连线
            Pen linePen = new Pen(Color.Red, 1);
            for (int i = 1; i < CP.Length; i++) {
                g.DrawLine(linePen, (float)CP[i].x, (float)CP[i].y, (float)CP[i-1].x, (float)CP[i-1].y);
            }
        }

        public void PaintCurvePoints(PaintEventArgs e, BSpline.Point[] pts)
        {
            Graphics g = e.Graphics;
            Pen myPen = new Pen(Color.Blue, 2);
            for (int i = 0; i < pts.Length; i++) {
                g.DrawEllipse(myPen, (float)pts[i].x, (float)pts[i].y, 1, 1);
            }
        }

        public void PaintCurveLerp(PaintEventArgs e, BSpline.Point[] CP, double[] knot, double u)
        {
            BSpline.Point tp = BSpline.Lerp(CP, knot, u);
            Graphics g = e.Graphics;
            Pen myPen = new Pen(Color.Green, 2);
            g.DrawEllipse(myPen, (float)tp.x, (float)tp.y, 1, 1);
        }

        private void Form1_Paint(object sender, PaintEventArgs e)
        {
            int k = 3;

            //控制点
            BSpline.Point[] CP = {
                //加个重复节点,让曲线经过起始控制点
                new BSpline.Point{ x=0, y=100},
                new BSpline.Point{ x=0, y=100},
                new BSpline.Point{ x=100, y=0},
                new BSpline.Point{ x=200, y=100},
                new BSpline.Point{ x=300, y=0},
                new BSpline.Point{ x=400, y=100},
                //加个重复节点,让曲线经过终止控制点
                new BSpline.Point{ x=400, y=100}
            };

            //节点向量
            double[] knot =
            {
                0, 100, 200, 300, 400, 500, 600, 700, 800, 900
            };

            //一次性计算出所有插值点
            BSpline.Point[] pts = BSpline.SplinePoints(CP, CP.Length - 1, knot, 100);
            PaintCurvePoints(e, pts);

            //控制点
            BSpline.Point[] CP1 = {
                new BSpline.Point{ x=0, y=300},
                new BSpline.Point{ x=0, y=300},
                new BSpline.Point{ x=100, y=200},
                new BSpline.Point{ x=200, y=300},
                new BSpline.Point{ x=300, y=200},
                new BSpline.Point{ x=400, y=300},
                new BSpline.Point{ x=400, y=300}
            };
            //------end

            //用插值方式绘制曲线
            double u0 = knot[k - 1];//起始插值点
            for (double u=u0; u<=700; u+=5)
            {
                PaintCurveLerp(e, CP1, knot, u);
            }
            //------end

            //绘制控制点
            PaintControlPoint(e, CP);
            PaintControlPoint(e, CP1);
        }
    }

    /// <summary>
    /// B-样条曲线
    /// </summary>
    public class BSpline
    {
        private const int k = 3; //三次B-样条曲线

        /// <summary>
        /// 插值方法
        /// </summary>
        /// <param name="CP">控制点坐标</param>
        /// <param name="knot">曲线结点向量</param>
        /// <param name="u">插值点</param>
        /// <returns></returns>
        public static Point Lerp(Point[] CP, double[] knot, double u)
        {
            int n = CP.Length;
            int i = 0;
            while ((i < n) && (u > knot[i + 1])) i++;
            return Deboor(CP, i, k, knot, u);
        }

        /// <summary>
        /// 计算曲线离散点序列
        /// </summary>
        /// <param name="CP">控制点坐标</param>
        /// <param name="n">控制点个数</param>
        /// <param name="k">曲线的阶数</param>
        /// <param name="knot">曲线结点向量</param>
        /// <param name="npoints">要计算出的离散点个数</param>
        /// <returns>采用德布尔(de Boor)算法生成的曲线上的离散点序列pts</returns>
        public static Point[] SplinePoints(Point[] CP, int n, double[] knot, int npoints)
        {
            Point[] pts = new Point[npoints];

            double u, delt;
            int i, j;
            //在每个节点区间，将参数t变化区间进行npoints等分
            delt = (knot[n + 1] - knot[k - 1]) / (double)npoints;
            i = k - 1;
            u = knot[k-1];
            for(j=0; j<npoints; j++)
            {
                //确定参数u所在的节点区间[ui, ui+1)
                while ((i < n) && (u > knot[i + 1])) i++;
                //在每个节点区间，分别求出npoints个离散点pts的坐标
                pts[j] = Deboor(CP, i, k, knot, u);
                u += delt;
            }
            return pts;
        }

        /// <summary>
        /// 用德布尔(de Boor)算法计算出插值点u的坐标
        /// 结点数(m)=控制点数(n)+次数(p)+1
        /// 举例: 14个控制点定义的6次B-样条曲线，其节点的数目是21=14+6+1
        /// </summary>
        /// <param name="CP">控制点坐标</param>
        /// <param name="i">第i个曲线段, i∈[0, n+k-l]</param>
        /// <param name="k">曲线的阶数</param>
        /// <param name="knot">曲线结点向量</param>
        /// <param name="u">变化范围为[ui, ui+1)</param>
        /// <returns>曲线在参数为t的坐标值</returns>
        public static Point Deboor(Point[] CP, int i, int k, double[] knot, double u)
        {
            double denom, alpha;
            Point[] p = new Point[k];
            const double epsilon = 0.0005;
            int index = 0;
            //p[]存放要参与计算的控制点
            for (int l = 0; l < k; l++)
            {
                index = i - k + l + 1;
                if (index < 0)
                    p[l] = CP[0].Clone();
                else if (index >= CP.Length)
                    p[l] = CP[CP.Length - 1].Clone();
                else
                    p[l] = CP[index].Clone();
            }
                

            //进行k-1次循环，即进行k-1级递推
            for(int r=1; r<k; r++)
            {
                //在每一级递推中，按照递减的顺序对控制顶点进行更新
                //按递减顺序更新，是为了确保已更新的控制顶点
                //不会对未更新的控制顶点的计算产生影响
                for(int m=k-1; m>=r; m--)
                {
                    int j = m + i - k + 1;
                    denom = knot[j + k - r] - knot[j];//分母为前一阶曲线的支撑区间
                    if (Math.Abs(denom) < epsilon)
                        alpha = 0;
                    else
                        alpha = (u - knot[j]) / denom;
                    //(1 - 比例因子) * 控制点坐标 + 比例因子 * 控制点坐标
                    p[m].x = (1 - alpha) * p[m - 1].x + alpha * p[m].x;
                    p[m].y = (1 - alpha) * p[m - 1].y + alpha * p[m].y;
                }
            }
            return p[k-1];
        }

        #region Point
        public class Point
        {
            public double x;
            public double y;
            public double z;

            public Point Clone()
            {
                Point p = new Point();
                p.x = x;
                p.y = y;
                p.z = z;
                return p;
            }
        }
        #endregion
    }
}

运行测试