B-spline反求控制頂點

B樣條曲線拟合時（https://blog.csdn.net/hanmingjunv5/article/details/106137002），曲線是不通過資料點的，這樣對于曲線的插補來說，是不合理的。是以，需要根據給定的資料點來求解控制頂點，使拟合的曲線通過全部資料點。

推導

曲線定義：

C ( u ) = ∑ i = 0 n N i , p ( u ) P i ( 1 ) C(u)=\sum^n_{i=0}N_{i,p}(u)P_i \quad \quad \quad (1) C(u)=i=0∑n​Ni,p​(u)Pi​(1)

B-Spline曲線存在n+1個位置的控制頂點。

将每個控制頂點的基函數寫成矩陣形式如下：

N = [ N 0 , p ( t 0 ) N 1 , p ( t 0 ) N 2 , p ( t 0 ) ⋯ N n , p ( t 0 ) N 0 , p ( t 1 ) N 1 , p ( t 1 ) N 2 , p ( t 1 ) ⋯ N n , p ( t 1 ) ⋯ N 0 , p ( t n ) N 1 , p ( t n ) N 2 , p ( t n ) ⋯ N n , p ( t n ) ] N=\left[ \begin{matrix} N_{0,p}(t_0) & N_{1,p}(t_0) & N_{2,p}(t_0) &\cdots & N_{n,p}(t_0)\\ N_{0,p}(t_1) & N_{1,p}(t_1) & N_{2,p}(t_1) &\cdots & N_{n,p}(t_1)\\ &\cdots\\ N_{0,p}(t_n) & N_{1,p}(t_n) & N_{2,p}(t_n) &\cdots & N_{n,p}(t_n)\\ \end{matrix} \right] N=⎣⎢⎢⎡​N0,p​(t0​)N0,p​(t1​)N0,p​(tn​)​N1,p​(t0​)N1,p​(t1​)⋯N1,p​(tn​)​N2,p​(t0​)N2,p​(t1​)N2,p​(tn​)​⋯⋯⋯​Nn,p​(t0​)Nn,p​(t1​)Nn,p​(tn​)​⎦⎥⎥⎤​

控制頂點寫成矩陣形式如下：

P = [ p ( 01 ) p ( 02 ) p ( 03 ) ⋯ p ( 0 s ) p ( 11 ) p ( 12 ) p ( 13 ) ⋯ p ( 1 s ) ⋯ p ( n 1 ) p ( n 2 ) p ( n 3 ) ⋯ p ( n s ) ] P=\left[ \begin{matrix} p(01) & p(02) & p(03) &\cdots & p(0s)\\ p(11) & p(12) & p(13) &\cdots & p(1s)\\ &\cdots\\ p(n1) & p(n2) & p(n3) &\cdots & p(ns)\\ \end{matrix} \right]\\ P=⎣⎢⎢⎡​p(01)p(11)p(n1)​p(02)p(12)⋯p(n2)​p(03)p(13)p(n3)​⋯⋯⋯​p(0s)p(1s)p(ns)​⎦⎥⎥⎤​

s為點位資料的次元，例如3D資料，s=3
           

資料點寫成矩陣形式如下：

C = [ c ( 01 ) c ( 02 ) c ( 03 ) ⋯ c ( 0 s ) c ( 11 ) c ( 12 ) c ( 13 ) ⋯ c ( 1 s ) ⋯ c ( n 1 ) c ( n 2 ) c ( n 3 ) ⋯ c ( n s ) ] C=\left[ \begin{matrix} c(01) & c(02) & c(03) &\cdots & c(0s)\\ c(11) & c(12) & c(13) &\cdots & c(1s)\\ &\cdots\\ c(n1) & c(n2) & c(n3) &\cdots & c(ns)\\ \end{matrix} \right]\\ C=⎣⎢⎢⎡​c(01)c(11)c(n1)​c(02)c(12)⋯c(n2)​c(03)c(13)c(n3)​⋯⋯⋯​c(0s)c(1s)c(ns)​⎦⎥⎥⎤​

則式子(1)可以寫成如下形式：

C = N ∗ P 因 此 ： P = N − 1 ∗ C \mathbf{C=N*P}\\ 是以：\\ \mathbf{P=N^{-1}*C} C=N∗P是以：P=N−1∗C

即可求解控制頂點P。

求解基函數系數矩陣

使用de-boor算法計算基函數系數矩陣，僞代碼如下。

Input: n, p, m, u, and m+1 clamped knots { u0, ..., um }
Output: Coefficients N0,p(u), N1,p(u), ..., Nn,p(u) in N[0], N[1], ..., N[n]
Algorithm:
	Initialize N[0..n] to 0; // initialization
	if u = u0 then // rule out special cases
		N[0] = 1.0;
		return
	else u = um then
		N[n] = 1.0
		return
	end if
	// now u is between u0 and um

	Let u be in knot span [uk,uk+1);
	N[k] := 1.0 // degree 0 coefficient
	for d :=1 to p do // degree d goes from 1 to p
		begin
			N[k-d] =  * N[(k-d)+1]; // right (south-west corner) term only
			for i := k-d+1 to k-1 do // compute internal terms
				N[i] :=  * N[i] +  * N[i+1];
			N[k] =  * N[k]; // let (north-west corner) term only
	end
	// array N[0..n] has the coefficients.
           

python結果仿真-3次B-Spline曲線

圖中，

紅色點位原始資料點；

黑色點為求解的控制點；

紅色曲線為由原始資料點拟合的B-Spline曲線；

黑色曲線為由控制頂點拟合的B-Spline曲線。

可知，黑色曲線通過全部資料點，算法驗證成功。

節點向量的選取

在計算控制頂點時，發現以下問題。

資料點總共n+1個，如果由此資料點構造節點向量，使用均勻向量法，得到n-3個非0非1控制節點。同樣，控制頂點的數量也為n+1個，加上首尾0、1的節點，總共n-1個節點。二計算基函數系數矩陣時，需要n+1個節點，是以，好少兩個節點。

這裡有兩個處理方法：

方法一：

仍然采用均勻向量法，在計算基函數系數矩陣時，添加在首尾添加兩個節點：

原 始 使 用 的 節 點 向 量 ： U = [ 0 , u 1 , u 2 , ⋯   , u n − 3 , 1 ] 增 加 後 的 節 點 向 量 ： U = [ 0 , t 1 , u 1 , u 2 , ⋯   , u n − 3 , t 2 , 1 ] t 1 = u 1 / 2 t 2 = ( n − 3 + 1 ) / 2 \begin{aligned} &原始使用的節點向量：\\ &\mathbf{U}=[0,u_1,u_2,\cdots,u_{n-3},1]\\ &增加後的節點向量：\\ &\mathbf{U}=[0,t_1,u_1,u_2,\cdots,u_{n-3},t_2,1]\\ &t_1=u_1/2\\ &t_2=({n-3}+1)/2\\ \end{aligned} ​原始使用的節點向量：U=[0,u1​,u2​,⋯,un−3​,1]增加後的節點向量：U=[0,t1​,u1​,u2​,⋯,un−3​,t2​,1]t1​=u1​/2t2​=(n−3+1)/2​

這樣計算後，即可按照僞代碼中的算法計算n+1個控制頂點。

方法二：

使用弧長法計算節點向量。

由n+1個資料點，得到n段線段長度:

L = [ l 1 , l 2 , l 3 , ⋯   , l n ] S = s u m ( L ) L = [l_1,l_2,l_3,\cdots,l_n]\\ S = sum(L) L=[l1​,l2​,l3​,⋯,ln​]S=sum(L)

則原始節點向量定義為：

k n o t s = [ 0 , 0 , 0 , 0 , l 1 + l 2 S , l 1 + l 2 + l 3 S , ⋯   , l 1 + l 2 + ⋯ + l n − 1 S , 1 , 1 , 1 , 1 ] knots = [0,0,0,0,\frac{l_1+l_2}{S},\frac{l_1+l_2+l_3}{S},\cdots,\frac{l_1+l_2+\cdots+l_{n-1}}{S},1,1,1,1] knots=[0,0,0,0,Sl1​+l2​​,Sl1​+l2​+l3​​,⋯,Sl1​+l2​+⋯+ln−1​​,1,1,1,1]

增加後的節點向量：

k n o t s ′ = [ 0 , 0 , 0 , 0 , l 1 S , l 1 + l 2 S , l 1 + l 2 + l 3 S , ⋯   , l 1 + l 2 + ⋯ + l n − 1 S , l 1 + l 2 + ⋯ + l n − 1 + l n S , 1 , 1 , 1 , 1 ] knots^{'} = [0,0,0,0,\mathbf{\frac{l_1}{S}},\frac{l_1+l_2}{S},\frac{l_1+l_2+l_3}{S},\cdots,\frac{l_1+l_2+\cdots+l_{n-1}}{S},\frac{l_1+l_2+\cdots+l_{n-1}+\mathbf{l_n}}{S},1,1,1,1] knots′=[0,0,0,0,Sl1​​,Sl1​+l2​​,Sl1​+l2​+l3​​,⋯,Sl1​+l2​+⋯+ln−1​​,Sl1​+l2​+⋯+ln−1​+ln​​,1,1,1,1]

按照僞代碼中的算法計算n+1個控制頂點。