?這篇博客跟我上一篇博客《均勻三次B樣條曲線插值實現及MATLAB代碼》的內容有點像，只是在基函數的計算上不同，造成均勻/非均勻的區別。

參考資料：

[1]（這個PPT講得很通俗，但對于多插值點分段曲線的內容漏講了一個知識點）三次周期B樣條曲線的算法 - 百度文庫 (baidu.com)

[2]（這個介紹只有兩個插值點的三次B樣條曲線，是B樣條曲線最簡單的形式了吧~）(7條消息) 從B樣條的插值點反求控制點_cofd的專欄-CSDN博客

[3]（一本書，里面有講到整體參數和局部參數設置、節點矢量劃分等）《計算機輔助幾何設計與非均勻有理B樣條》

正文：

曲線插值一般指的是給定插值點，得出曲線的方程，曲線會經過所有的插值點。確定三次B樣條曲線的輸入量有兩種，一種是給出控制點和其它邊界條件，曲線一般不經過控制點；一種是給出插值點和其它邊界條件，曲線會經過所有插值點，顯然第二種輸入量使用更為廣泛，這里只介紹第二種情況。

①首先輸入插值點，這些插值點可以是二維(x,y)點，也可以是三維(x,y,z)點，甚至更多維都可以，處理方法都相似，以二維點為例，MATLAB代碼：

pointInput=[1 209.1;1.033 209.525;1.067 209.273;1.1 209.277;1.133 208.734;1.167 208.693;1.2 208.852;1.233 208.722;1.267 208.746;1.3 208.759];

這些點的分布如下圖的紅色拆線所示（不是藍色的*型點，藍色的是我最終擬合的曲線）。

??②計算N、n、k等值

N為插值點的個數，上面共有10個插值點，所以N=10。控制點個數為N+2個。n=N+1。k為次數，這里是三次B樣條曲線，所以k=3。MATLAB代碼：

N=length(pointInput); k=3; n=N-1;

③反算控制點

這里主要參考了資料[1]PPT中的內容，具體的請查看該PPT。三次B樣條曲線的表達式（1）為

?其中是第i段曲線的表達式，是局部參數，是第i個控制點。注意節點、節點矢量、整體參數、局部參數的區別（具體的要看資料[3]或其它相關的資料了，說明這些可能需要長篇大論哦~）。

反算控制點，我用了第三種邊界條件，如下圖

增加了邊界條件后，就能求如下圖的矩陣方程了，其實下面的矩陣方程是由表達式（1）得到的，具體可以參考資料[1]的PPT（注意表達式中的系數6不要看漏了）。

?

?對于二維點，我想對x和y坐標的B樣條曲線都求出來，所以用了兩組矩陣方程求解它們各自的控制點（對于三維點，還能用同樣的方法求z坐標的B樣條曲線）。MATLAB代碼：

%①先求x、y坐標，如果是三維點，還要求z b1=[0;pointInput(:,1);0]*6;%注意系數 px=Chase_method(A1,b1);%用追趕法求得所有控制點 b2=[0;pointInput(:,2);0]*6; py=Chase_method(A1,b2);

其中求三對角陣解用追趕法Chase_method（可以直接用常用的方法求解，但此處用追趕法時間復雜度更低），其MATLAB代碼：

function [ x ] = Chase_method( A, b ) %Chase method 追趕法求三對角矩陣的解 % A為三對角矩陣的系數，b為等式右端的常數項，返回值x即為最終的解 % 注：A盡量為方陣，b一定要為列向量 %% 求追趕法所需L及U T = A; for i = 2 : size(T,1)T(i,i-1) = T(i,i-1)/T(i-1,i-1);T(i,i) = T(i,i) - T(i-1,i) * T(i,i-1); end L = zeros(size(T)); L(logical(eye(size(T)))) = 1; %對角線賦值1 for i = 2:size(T,1)for j = i-1:size(T,1)L(i,j) = T(i,j);break;end end U = zeros(size(T)); U(logical(eye(size(T)))) = T(logical(eye(size(T)))); for i = 1:size(T,1)for j = i+1:size(T,1)U(i,j) = T(i,j);break;end end %% 利用matlab解矩陣方程的遍歷直接求解 y = L\b; x = U\y; end

④確定節點矢量

節點矢量是分段B樣條曲線進行分段的依據（請看下資料[3]或其它相應資料），我用哈德利-賈德方法確定節點矢量（節點矢量共有n+k+2個，注意這里是小寫n而不是大寫的N），即節點矢量中前k+1個節點取值0，后k+1個節點取值為1，重復度為k+1，然后中間的節點由哈德利-賈德方法計算得到。MATLAB代碼（注意MATLAB的數組和矩陣的下標是從1開始的，與公式中的從0開始有區別）：

%確定節點矢量 u=zeros(1,N+1+k+2);%0:1/(N+1+k+2-1):1; for i=1:k+1u(length(u)-i+1)=1; end l=zeros(1,N+1);%控制多邊形各邊長 for i=1:N+1%注意如果用到三維點，這里要作修改l(i)=sqrt((px(i+1)-px(i))^2+(py(i+1)-py(i))^2); end sum1=0;%用哈德利-賈德方法決定節點矢量 for g=k+1:N+1+1for j=g-k:g-1sum1=sum1+l(j);end end for i=k+2:N+2u(i)=sum(l(i-k:i-1))/sum1+u(i-1); end

計算得到了所有的控制點以及節點矢量，已經能知道該分段B樣條曲線的表達式了，下面就輸入一個點來測試。其實該分段B樣條曲線用參數方程表示，x、y坐標都以參數u來表示，同時以節點矢量中每個節點作為分段曲線分段的邊界，U中的元素的單調不減小的，后一個元素一不小于一個元素。下面MATLAB代碼使u取值為測試點uTest，然后先找到uTest屬于哪個分段，計算得到uTest屬于score分段：

uTest=0.5; score=find(uTest>=u,1,'last');%確定該參數方程的參數值所屬的分段 if uTest==1 score=find(u==1,1)-1;end score=score-k;%表示第score段樣條曲線

在用表達式(1)求出參數uTest對應的x、y坐標值之間，需要將整體變量uTest轉換為局部變量t。注意uTest在整體樣條曲線中的取值范圍是[0,1]，而t在第score分段曲線中的取值范圍是[0,1]。轉換公式為（為節點矢量中第i+1個元素，且滿足，這里的i可由score確定），從而，得到表達式(1)中的矩陣，MATLAB代碼：

uNode=1;%上面的uTest是整體參數，用來確定是哪段曲線，要將整體參數轉化為局部參數t t=(uTest-u(score+k))/(u(score+1+k)-u(score+k)); for i=1:kuNode=[t^i uNode]; end

對比我另一篇博客《均勻三次B樣條曲線插值實現及MATLAB代碼》，本篇博客的樣條曲線是非均勻的，即用的基函數并不是固定的值，而是由德布爾-考克斯遞推公式計算得到。先求出基函數，MATLAB代碼：

Nb=zeros(k+3,k+1); for j=1:k+1for i=1:n+1if j==1if i==score+kNb(i,j)=1; continue;elseNb(i,j)=0; continue;endelse%if u(i+j-1)-u(i)==0||u(i+j+1-1)-u(i+1)==0 continue;endif u(i+j-1)-u(i)==0&&u(i+j+1-1)-u(i+1)==0Nb(i,j)=0;elseif u(i+j-1)-u(i)==0Nb(i,j)=(u(i+j+1-1)-uTest)/(u(i+j+1-1)-u(i+1))*Nb(i+1,j-1);elseif u(i+j+1-1)-u(i+1)==0Nb(i,j)=(uTest-u(i))/(u(i+j-1)-u(i))*Nb(i,j-1);elseNb(i,j)=(uTest-u(i))/(u(i+j-1)-u(i))*Nb(i,j-1)+(u(i+j+1-1)-uTest)/(u(i+j+1-1)-u(i+1))*Nb(i+1,j-1);end%Nb(i,j)=(uTest-u(i))/(u(i+j-1)-u(i))*Nb(i,j-1)+(u(i+j+1-1)-uTest)/(u(i+j+1-1)-u(i+1))*Nb(i+1,j-1);endend end

最后由基函數和控制點計算得到夢寐以求的x、y坐標值，MATLAB代碼：

Qx=0;Qy=0; for i=score:score+k %注意非零基函數的個數Qx=Qx+px(i)*Nb(i,k+1);Qy=Qy+py(i)*Nb(i,k+1); end

附上全部MATLAB程序：

%三次B樣條插值，給定插值點及兩個邊界條件，反算出控制點，由控制點得到分段曲線 %三維點(x,y,z)和二維點(x,y)的插值算法類似，控制多邊形各邊長l要作修改 %pointInput插值點 pointInput=[1 209.1;1.033 209.525;1.067 209.273;1.1 209.277;1.133 208.734;1.167 208.693;1.2 208.852;1.233 208.722;1.267 208.746;1.3 208.759]; %pointInput=[0 10;5 8.66;10 0]; %用第三種邊界條件 https://wenku.baidu.com/view/2482396e011ca300a6c390b3.html N=length(pointInput); k=3; A1=eye(N+2,N+2)*4; A1(1,1)=6;A1(1,2)=-6;A1(N+2,N+1)=6;A1(N+2,N+2)=-6; for i=2:N+1A1(i,i-1)=1;A1(i,i+1)=1; end %①先求x、y坐標，如果是三維點，還要求z b1=[0;pointInput(:,1);0]*6;%注意系數 px=Chase_method(A1,b1);%用追趕法求得所有控制點 b2=[0;pointInput(:,2);0]*6; py=Chase_method(A1,b2); %確定節點矢量 u=zeros(1,N+1+k+2);%0:1/(N+1+k+2-1):1; for i=1:k+1u(length(u)-i+1)=1; end l=zeros(1,N+1);%控制多邊形各邊長 for i=1:N+1%注意如果用到三維點，這里要作修改l(i)=sqrt((px(i+1)-px(i))^2+(py(i+1)-py(i))^2); end sum1=0;%用哈德利-賈德方法決定節點矢量 for g=k+1:N+1+1for j=g-k:g-1sum1=sum1+l(j);end end for i=k+2:N+2u(i)=sum(l(i-k:i-1))/sum1+u(i-1); end %用基函數求值 %uTestArray=0:0.01:1; for uTest=0:0.01:1 %uTest=0.5; %u3Array=[uTest^3 uTest^2 uTest 1];%三次B樣條專屬 score=find(uTest>=u,1,'last');%確定該參數方程的參數值所屬的分段 if uTest==1 score=find(u==1,1)-1;end score=score-k;%表示第score段樣條曲線 %用德布爾-考克斯遞推公式計算基函數 Nb=zeros(k+3,k+1); for j=1:k+1for i=1:n+1if j==1if i==score+kNb(i,j)=1; continue;elseNb(i,j)=0; continue;endelse%if u(i+j-1)-u(i)==0||u(i+j+1-1)-u(i+1)==0 continue;endif u(i+j-1)-u(i)==0&&u(i+j+1-1)-u(i+1)==0Nb(i,j)=0;elseif u(i+j-1)-u(i)==0Nb(i,j)=(u(i+j+1-1)-uTest)/(u(i+j+1-1)-u(i+1))*Nb(i+1,j-1);elseif u(i+j+1-1)-u(i+1)==0Nb(i,j)=(uTest-u(i))/(u(i+j-1)-u(i))*Nb(i,j-1);elseNb(i,j)=(uTest-u(i))/(u(i+j-1)-u(i))*Nb(i,j-1)+(u(i+j+1-1)-uTest)/(u(i+j+1-1)-u(i+1))*Nb(i+1,j-1);end%Nb(i,j)=(uTest-u(i))/(u(i+j-1)-u(i))*Nb(i,j-1)+(u(i+j+1-1)-uTest)/(u(i+j+1-1)-u(i+1))*Nb(i+1,j-1);endend end Qx=0;Qy=0; for i=score:score+k %注意非零基函數的個數Qx=Qx+px(i)*Nb(i,k+1);Qy=Qy+py(i)*Nb(i,k+1); end plot(Qx,Qy,'*');hold on; end plot(pointInput(:,1),pointInput(:,2),'r'); function [ x ] = Chase_method( A, b ) %Chase method 追趕法求三對角矩陣的解 % A為三對角矩陣的系數，b為等式右端的常數項，返回值x即為最終的解 % 注：A盡量為方陣，b一定要為列向量 %% 求追趕法所需L及U T = A; for i = 2 : size(T,1)T(i,i-1) = T(i,i-1)/T(i-1,i-1);T(i,i) = T(i,i) - T(i-1,i) * T(i,i-1); end L = zeros(size(T)); L(logical(eye(size(T)))) = 1; %對角線賦值1 for i = 2:size(T,1)for j = i-1:size(T,1)L(i,j) = T(i,j);break;end end U = zeros(size(T)); U(logical(eye(size(T)))) = T(logical(eye(size(T)))); for i = 1:size(T,1)for j = i+1:size(T,1)U(i,j) = T(i,j);break;end end %% 利用matlab解矩陣方程的遍歷直接求解 y = L\b; x = U\y; end

效果圖，藍色點的是曲線擬合的點，紅色拆線是由插值點組成的：

總結

以上是生活随笔為你收集整理的非均匀三次B样条曲线插值实现及MATLAB代码的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。