MATLAB绘制B样条曲线 / 开普饭

1 B样条曲线

1.1 B样条曲线定义

B样条方法具有表示与设计自由型曲线曲面的强大功能，是形状数学描述的主流方法之一，另外B样条方法是目前工业产品几何定义国际标准——有理B样条方法(NURBS)的基础。B样条方法兼备了Bezier方法的一切优点，包括几何不变性，仿射不变性等等，同时克服了Bezier方法中由于整体表示带来不具有局部性质的缺点（移动一个控制顶点将会影响整个曲线）。B样条曲线方程可写为：

p(u)=∑i=0ndiNi,k(u)" role="presentation" style="position: relative;">p(u)=∑i=0ndiNi,k(u)p(u)=∑i=0ndiNi,k(u)

其中，di(i=0,1...n)" role="presentation" style="position: relative;">di(i=0,1...n)di(i=0,1...n)为控制顶点（坐标），Ni,k(i=0,1...n)" role="presentation" style="position: relative;">Ni,k(i=0,1...n)Ni,k(i=0,1...n)为k" role="presentation" style="position: relative;">kk次规范B样条基函数，最高次数是k" role="presentation" style="position: relative;">kk。基函数是由一个称为节点矢量的非递减参数u" role="presentation" style="position: relative;">uu的序列U" role="presentation" style="position: relative;">UU：u0≤u1≤...≤un+k+1" role="presentation" style="position: relative;">u0≤u1≤...≤un+k+1u0≤u1≤...≤un+k+1所决定的k" role="presentation" style="position: relative;">kk次分段多项式。
B样条的基函数通常采用Cox-deBoor递推公式：

{Ni,0(u)={1, if ui≤u≤ui+10, othersNi,k=u−uiui+k−uiNi,k−1(u)+ui+k+1−uui+k+1−ui+1Ni+1,k−1(u)define 00=0" role="presentation" style="position: relative;">⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪Ni,0(u)={1, if ui≤u≤ui+10, othersNi,k=u−uiui+k−uiNi,k−1(u)+ui+k+1−uui+k+1−ui+1Ni+1,k−1(u)define 00=0{Ni,0(u)={1, if ui≤u≤ui+10, othersNi,k=u−uiui+k−uiNi,k−1(u)+ui+k+1−uui+k+1−ui+1Ni+1,k−1(u)define 00=0

式中 i" role="presentation" style="position: relative;">ii 为节点序号， k" role="presentation" style="position: relative;">kk 是基函数的次数，共有n+1" role="presentation" style="position: relative;">n+1n+1个控制顶点。注意区分节点和控制顶点，节点是在节点矢量U" role="presentation" style="position: relative;">UU中取得，控制顶点则是坐标点，决定B样条的控制多边形。Cox-deBoor递推公式是B样条曲线的定义的核心，该部分在程序中实现可采用递归的方式：

% BaseFunction.m文件function Nik_u = BaseFunction(i, k , u, NodeVector)% 计算基函数Ni,k(u),NodeVector为节点向量if k == 0       % 0次B样条    if (u >= NodeVector(i+1)) && (u < NodeVector(i+2))        Nik_u = 1.0;    else        Nik_u = 0.0;    endelse    Length1 = NodeVector(i+k+1) - NodeVector(i+1);    Length2 = NodeVector(i+k+2) - NodeVector(i+2);      % 支撑区间的长度    if Length1 == 0.0       % 规定0/0 = 0        Length1 = 1.0;    end    if Length2 == 0.0        Length2 = 1.0;    end    Nik_u = (u - NodeVector(i+1)) / Length1 * BaseFunction(i, k-1, u, NodeVector) ...        + (NodeVector(i+k+2) - u) / Length2 * BaseFunction(i+1, k-1, u, NodeVector);end

所给程序可用于计算基函数Ni,k(u)" role="presentation" style="position: relative;">Ni,k(u)Ni,k(u)的值，程序中对不同类型的B样条曲线区别在于节点矢量 NodeVector 的取值不同。

1.2 B样条曲线的分类

根据节点矢量中节点的分布情况不同，可以划分4中类型的B样条曲线：
1. 均匀B样条曲线
节点矢量中节点为沿参数轴均匀或等距分布。
2. 准均匀B样条曲线
其节点矢量中两端节点具有重复度k+1" role="presentation" style="position: relative;">k+1k+1，即u0=u1=...=uk" role="presentation" style="position: relative;">u0=u1=...=uku0=u1=...=uk，un+1=un+2=...=un+k+1" role="presentation" style="position: relative;">un+1=un+2=...=un+k+1un+1=un+2=...=un+k+1，所有的内节点均匀分布，具有重复度1。
3. 分段Bezier曲线
其节点矢量中两端节点的重复度与类型2相同，为k+1" role="presentation" style="position: relative;">k+1k+1。不同的是内节点重复度为k" role="presentation" style="position: relative;">kk。该类型有限制条件，控制顶点数减1必须等于次数的正整数倍，即nk=正整数" role="presentation" style="position: relative;">nk=正整数nk=正整数。
4. 一般非均匀B样条曲线
对任意分布的节点矢量U=[u0,u1...un+k+1]" role="presentation" style="position: relative;">U=[u0,u1...un+k+1]U=[u0,u1...un+k+1]，只要在数学上成立都可选取。

2 B样条曲线的绘制

2.1 节点矢量的确定

不同类型的B样条曲线区别主要在于节点矢量，对于具有n+1" role="presentation" style="position: relative;">n+1n+1个控制顶点(P0,P1,...,Pn)" role="presentation" style="position: relative;">(P0,P1,...,Pn)(P0,P1,...,Pn)的 k" role="presentation" style="position: relative;">kk 次B样条曲线，无论是哪种类型都具有n+k+2" role="presentation" style="position: relative;">n+k+2n+k+2个节点([u0,u1...un+k+1])" role="presentation" style="position: relative;">([u0,u1...un+k+1])([u0,u1...un+k+1])。

根据图示，三种类型的B样条曲线对应的节点矢量分别为：

[01727374757671]" role="presentation" style="position: relative;">[01727374757671][01727374757671]

[0 0 013231 1 1]" role="presentation" style="position: relative;">[0 0 013231 1 1][0 0 013231 1 1]

[0 0 012121 1 1]" role="presentation" style="position: relative;">[0 0 012121 1 1][0 0 012121 1 1]

需要注意的是分段Bezier曲线必须满足nk=正整数" role="presentation" style="position: relative;">nk=正整数nk=正整数。

这里给出准均匀B样条和分段Bezier曲线的生成节点矢量的代码，均匀B样条的很简单就不列出了。假设共n+1个控制顶点，k次B样条，输入参数为 n, k ，输出节点矢量NodeVector。

准均匀B样条曲线的节点矢量生成：

% U_quasi_uniform.m文件function NodeVector = U_quasi_uniform(n, k)% 准均匀B样条的节点向量计算，共n+1个控制顶点，k次B样条NodeVector = zeros(1, n+k+2);piecewise = n - k + 1;       % 曲线的段数if piecewise == 1       % 只有一段曲线时，n = k    for i = n+2 : n+k+2        NodeVector(1, i) = 1;    endelse    flag = 1;       % 不止一段曲线时    while flag ~= piecewise        NodeVector(1, k+1+flag) = NodeVector(1, k + flag) + 1/piecewise;        flag = flag + 1;    end    NodeVector(1, n+2 : n+k+2) = 1;end

分段Bezier曲线的节点矢量生成：

% U_piecewise_Bezier.m文件function NodeVector = U_piecewise_Bezier(n, k)% 分段Bezier曲线的节点向量计算，共n+1个控制顶点，k次B样条% 分段Bezier端节点重复度为k+1，内间节点重复度为k,且满足n/k为正整数if ~mod(n, k) && (~mod(k, 1) && k>=1)   % 满足n是k的整数倍且k为正整数    NodeVector = zeros(1, n+k+2);   % 节点矢量长度为n+k+2    NodeVector(1, n+2 : n+k+2) = ones(1, k+1);  % 右端节点置1    piecewise = n / k;      % 设定内节点的值    Flg = 0;    if piecewise > 1        for i = 2 : piecewise            for j = 1 : k                NodeVector(1, k+1 + Flg*k+j) = (i-1)/piecewise;            end            Flg = Flg + 1;        end    endelse    fprintf('error!\n');end

2.2 B样条曲线的绘制

根据B样条曲线的定义公式，曲线上任一点坐标值是参数变量u" role="presentation" style="position: relative;">uu的函数，用矩阵形式表示

p(u)=(d0d1…dn)(N0,k(u)N1,k(u)⋮Nn,k(u))" role="presentation" style="position: relative;">p(u)=(d0d1…dn)⎛⎝⎜⎜⎜⎜⎜N0,k(u)N1,k(u)⋮Nn,k(u)⎞⎠⎟⎟⎟⎟⎟p(u)=(d0d1…dn)(N0,k(u)N1,k(u)⋮Nn,k(u))

只需要确定控制顶点di" role="presentation" style="position: relative;">didi、曲线的次数k" role="presentation" style="position: relative;">kk 以及基函数Ni,k(u)" role="presentation" style="position: relative;">Ni,k(u)Ni,k(u)，就完全确定了曲线。

B样条曲线的绘制函数：

% DrawSpline.m文件function DrawSpline(n, k, P, NodeVector)% B样条的绘图函数% 已知n+1个控制顶点P(i), k次B样条，P是2*(n+1)矩阵存控制顶点坐标, 节点向量NodeVectorplot(P(1, 1:n+1), P(2, 1:n+1),...                    'o','LineWidth',1,...                    'MarkerEdgeColor','k',...                    'MarkerFaceColor','g',...                    'MarkerSize',6);line(P(1, 1:n+1), P(2, 1:n+1));Nik = zeros(n+1, 1);for u = 0 : 0.005 : 1-0.005    for i = 0 : 1 : n        Nik(i+1, 1) = BaseFunction(i, k , u, NodeVector);    end    p_u = P * Nik;    if u == 0        tempx = p_u(1,1);        tempy = p_u(2,1);        line([tempx p_u(1,1)], [tempy p_u(2,1)],...            'Marker','.','LineStyle','-', 'Color',[.3 .6 .9], 'LineWidth',3);    else        line([tempx p_u(1,1)], [tempy p_u(2,1)],...            'Marker','.','LineStyle','-', 'Color',[.3 .6 .9], 'LineWidth',3);        tempx = p_u(1,1);        tempy = p_u(2,1);    endend

调用 DrawSpline(n, k, P, NodeVector) 函数就能绘制曲线，注意输入变量要正确。

下面给出绘制三种不同B样条曲线的命令流，可以参考比较每种类型之间的区别。

% 绘制三种类型的B样条曲线，需要前面所给的所有.m文件clear all;%控制顶点P = [9.036145, 21.084337, 37.607573, 51.893287, 61.187608;    51.779661, 70.084746, 50.254237, 69.745763, 49.576271];n = 4; k = 2;flag = 2;% flag = 1，绘制均匀B样条曲线% flag = 2, 绘制准均匀B样条曲线% flag = 3, 绘制分段Bezier曲线switch flag    case 1        NodeVector = linspace(0, 1, n+k+2); % 均匀B样条的节点矢量        % 绘制样条曲线        plot(P(1, 1:n+1), P(2, 1:n+1),...                        'o','LineWidth',1,...                        'MarkerEdgeColor','k',...                        'MarkerFaceColor','g',...                        'MarkerSize',6);        line(P(1, 1:n+1), P(2, 1:n+1));        Nik = zeros(n+1, 1);        for u = k/(n+k+1) : 0.001 : (n+1)/(n+k+1)            % for u = 0 : 0.005 : 1            for i = 0 : 1 : n                Nik(i+1, 1) = BaseFunction(i, k , u, NodeVector);            end        p_u = P * Nik;        line(p_u(1,1), p_u(2,1), 'Marker','.','LineStyle','-', 'Color',[.3 .6 .9]);        end    case 2        NodeVector = U_quasi_uniform(n, k); % 准均匀B样条的节点矢量        DrawSpline(n, k, P, NodeVector);    case 3        NodeVector = U_piecewise_Bezier(n, k);  % 分段Bezier曲线的节点矢量        DrawSpline(n, k, P, NodeVector);    otherwise        fprintf('error!\n');end

三种类型的B样条曲线：
1. 均匀B样条曲线

2. 准均匀B样条曲线

3. 分段Bezier曲线

参考文献：

[1] 施法中. 计算机辅助几何设计与非均匀有理B样条(修订版)[M]. 北京: 高等教育出版社, 2013 : 217-248.

MATLAB绘制B样条曲线