MATLAB繪制B樣條曲線

本文轉載自查看原文 2022-04-09 17:14 1819

1 B樣條曲線

1.1 B樣條曲線定義

B樣條方法具有表示與設計自由型曲線曲面的強大功能，是形狀數學描述的主流方法之一，另外B樣條方法是目前工業產品幾何定義國際標准——有理B樣條方法(NURBS)的基礎。B樣條方法兼備了Bezier方法的一切優點，包括幾何不變性，仿射不變性等等，同時克服了Bezier方法中由於整體表示帶來不具有局部性質的缺點（移動一個控制頂點將會影響整個曲線）。B樣條曲線方程可寫為：

其中，d_i(i=0,1...n)為控制頂點（坐標）， N_i,k(i=0,1...n)k次規范B樣條基函數，最高次數是k。基函數是由一個稱為節點矢量的非遞減參數u的序列U： u₀≤u₁≤...≤u_n+k+1所決定的k次分段多項式。

B樣條的基函數通常采用Cox-deBoor（德布爾—考克斯）遞推公式：

式中 i 為節點序號， k是基函數的次數，共有 n+1個控制頂點。
注意區分節點和控制頂點，節點是在節點矢量U中取得，控制頂點則是坐標點，決定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; end else 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

所給程序可用於計算基函數 $N_{i, k} (u)$

1.2 B樣條曲線的分類

根據節點矢量中節點的分布情況不同，可以划分4中類型的B樣條曲線：

1. 均勻B樣條曲線

節點矢量中節點為沿參數軸均勻或等距分布。

2. 准均勻B樣條曲線

其節點矢量中兩端節點具有重復度k+1，即u₀=u₁=...=u_k，u_n+1=u_n+2=...=u_n+k+1，所有的內節點均勻分布，具有重復度1。

3. 分段Bezier曲線

其節點矢量中兩端節點的重復度與類型2相同，為 $k + 1$

$k + 1$

對任意分布的節點矢量 $U = [u_{0}, u_{1} . . . u_{n + k + 1}]$

2 B樣條曲線的繪制

2.1 節點矢量的確定

不同類型的B樣條曲線區別主要在於節點矢量，對於具有n+1個控制頂點(P₀,P₁,...,P_n)的 k 次B樣條曲線，無論是哪種類型都具有n+k+2個節點([u₀,u₁...u_n+k+1])。

根據圖示，三種類型的B樣條曲線對應的節點矢量分別為：

需要注意的是分段Bezier曲線必須滿足 $\frac{n}{k} = 正整数$

這里給出准均勻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; end else 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 end else fprintf('error!\n'); end

2.2 B樣條曲線的繪制

只需要確定控制頂點 $d_{i}$

B樣條曲線的繪制函數：

% DrawSpline.m文件
 function DrawSpline(n, k, P, NodeVector) % B樣條的繪圖函數 % 已知n+1個控制頂點P(i), k次B樣條，P是2*(n+1)矩陣存控制頂點坐標, 節點向量NodeVector 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 = 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); end end

調用 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