幾何模組化與處理之四 三次樣條曲線

幾何模組化與處理之四 三次樣條曲線

目錄

• ​幾何模組化與處理之四 三次樣條曲線​

• ​幾何設計​

• ​​樣條曲線​​
• ​​自由曲線​​

• ​樣條曲線的數學表達推導​

• ​​力學解釋​​
• ​​數學性質​​
• ​求解思路​

• ​​三次樣條插值函數​​
• ​​三彎矩方程​​
• ​​簡化技巧​​

• ​三次樣條曲線​

• ​​樣條函數局限性​​
• ​​三次參數樣條曲線​​

• ​曲線幾何連續性​

• ​​參數連續性​​
• ​幾何連續性​

• ​​幾何連續性定義​​
• ​​幾何連續性性質​​
• ​​幾何連續性的具體形式​​

• ​作業​

• ​追趕法 （Thomas Algorithm）​

• ​​推導一：​​
• ​​推導二：​​

• ​​C2連續​​
• ​​部分型值點G1連續​​
• ​​部分型值點C0連續​​

幾何設計

概念設計 >> 設計圖紙 >> 數學（幾何）表達

樣條曲線

在早期沒有計算機時，工業産品設計中工程師使用壓鐵使富有彈性的樣條通過指定的樣點，調整樣條使它具有滿意的形狀，然後沿樣條畫出曲線，稱為樣條曲線。

• 樣條：具有一定彈性的軟木條
• 壓鐵：較重的鐵塊，用來固定樣條所經過的點

  

自由曲線

在産品初始設計階段，描述其外形的曲線或曲面常常隻有大緻的形狀或隻知道它通過的一些空間列（稱為型值點），這類沒有數學表達式的曲線或曲面稱為自由曲線或自由曲面（Free from curve/surface）。

造型方法

• 拟合（fitting）
• 插值（interpolation）
• 逼近（approximation）

樣條曲線的數學表達推導

力學解釋

曲線：彈性梁，彈性模量為E

壓鐵：載荷（力）

由材料力學中的貝努利-歐拉方程

\[M(x)=EIk(x)=EI\frac{y''(x)}{(1+(y'(x))^2)^{\frac32}}\]

其中，M(x)：曲線彎矩 E：木樣條的彈性模量 I：幾何慣性矩 k(x)：曲線的曲率

小擾度假設：y'(x)<<1（彎角不大于45°）

\[M(x)\approx EIy''(x)\]

因兩壓鐵之間無外力，M(x)為一次

\[M(x)=ax+b\]

兩壓貼間y(x)為三次函數，即樣條曲線為分段三次函數。

數學性質

分段3次多項式的好處

• 2次多項式無法表達拐點，不夠自由
• 高次（4次及以上）多項式拐點多，次數若較高計算易出現較大誤差

求解思路

每段為3次多項式，有4個變量：

\[y_i(x)=a_i+b_ix+c_ix^2+d_ix^3\]

假設有n+1個型值點（n段），則總共有4n個變量。

1.首先，曲線要插值型值點，有n+1個限制條件；

2.其次，假設曲線整體為C2連續，則相鄰兩端在拼接點要滿足3個條件（C0連續、C1連續，C2連續），共3n-3個條件。

共4n-2個方程，再加2個額外條件（邊界條件），即可唯一整條曲線。

三次樣條插值函數

定義： 函數\(S(x) \in C^2 [a, b]\)，且在每個小區間 \([x_i , x_{i+1} ]\) 上是三次多項式，其中 $a = x_0 < x_1 < . . . < x_n = b $是給定節點，則稱 \(S(x)\) 是節點 \(x_0, x_1, . . . , x_n\) 上的三次樣條函數。 若在節點\(x_i\)上給定函數值$ y_i = f(x_i)(i = 0, 1, . . . , n)$，并成立：

\[S(x_i) = y_i (i = 0, 1, . . . , n)\]

則稱$ S(x)$ 為三次樣條插值函數。

根據 \(S(x)\) 在 \([a, b]\) 上二階導數連續，在節點 \(x_i(i = 1, 2, . . . , n − 1)\) 處應滿足連續性條件：

\[S(x_i − 0) = S(x_i + 0),\\S ′ (x_i − 0) = S ′ (x_i + 0), \\

S ′′(x_i − 0) = S ′′(x_i + 0)

\]

共有 3n − 3 個條件，再加上 $S(x) $滿足插值條件 \(S(x_i) = y_i (i = 0, 1, . . . , n)\)，共 4n − 2 個條件，是以還需 要兩個條件才能确定\(S(x)\)。

通常可在區間 \([a, b]\) 端點 \(a = x_0, b = x_n\) 處各加一個條件（稱為邊界條件），可根據實際問題的要求給定。常見的有以下三種：

1. 已知兩端的一階導數值，即：

   \[S ′ (x_0) = f ′_0 , S′ (x_n) = f ′_n\]

2. 已知兩端的二階導數值，即：

   \[S ′′(x_0) = f ′′_0 , S′′(x_n) = f ′′_n\]

   其特殊情況：

   \[S ′′(x_0) = S ′′(x_n) = 0\]

   稱為自然邊界條件

3. 當 \(f(x)\) 是以 $x_n − x_0 $為周期的周期函數時，則要求 \(S(x)\) 也是周期函數。這時邊界條件應滿足：

   \[S(x_0 + 0) = S(x_n − 0), \\S ′ (x_0 + 0) = S ′ (x_n − 0),\\

   S ′′(x_0 + 0) = S ′′(x_n − 0)

   而此時$ y_0 = y_n$。這樣确定的樣條函數 \(S(x)\)，稱為周期樣條函數。

三彎矩方程

​ 将在區間$ [x_i , x_{i+1}] \(上表示\) S(x)$的多項式記為 \(S_i(x)\)，現在來推導區間 \([x_i , x_{i+1}]\) 上$ S_i(x)$ 的 表達式。

首先，定義一組數 \(M_i = S ′′(x_i)\)。由于\(S_i(x)\) 是$[x_i , x_{i+1}] $上的三次多項式，是以 $S ′′_i (x) \(是滿 足\) S ′′ _i (x_i) = M_i $和 \(S ′′_i (x_{i+1}) = M_{i+1}\) 的線性函數，是以說 \(S ′′ _i (x)\) 也是 \(M_i\) 和 \(M_{i+1}\) 之間的直線：

\[S ′′ _i (x) = \frac{M_i}{h_i} (x_{i+1} − x) + \frac{M_{i+1}}{h_i} (x − x_i)\]

其中 \(h_i = x_{i+1} − x_i\)，把這個函數積分兩次，其結果就是 \(S_i\) :

\[S_i(x) = \frac{M_i}{6h_i} (x_{i+1} − x)^3 + \frac{M_{i+1}}{6h_i} (x − x_i)^3 + C(x − x_i) + D(x_{i+1} − x)\]

其中 \(C\) 和 \(D\) 是積分常數, 将插值條件$S_i(x_i) = y_i $和 \(S_i(x_{i+1}) = y_{i+1}\) 作用在$ Si $上就可以确定 \(C\) 和 \(D\)，結 果為：

\[S_i(x) = \frac{M_i}{6h_i} (x_{i+1} − x)^3 + \frac{M_{i+1}}{6h_i} (x − x_i)^3 + (\frac{y_{i+1}}{hi} − \frac{M_{i+1}h_i}6 )(x − x_i) + ( \frac{y_i}{h_i} − \frac{M_ih_i} 6 )(x_{i+1} − x)\tag {1}\]

上式稱為三次樣條插值函數的三彎矩方程。是以一旦确定了 $M_0, M_1, . . . , M_n \(的值以後，可以用上式算出\) S(x)$ 在區間 \([x_0, x_n]\) 内任意一點 \(x\) 處的值。

​ 下面，用$ S ′ (x)$ 的連續性條件來确定 $M_0, M_1, . . . , M_n $, 在内節點 \(x_i\) 上，一定有 \(S ′ _{i−1} (x_i) = S ′ _i (x_i)\)。 對(1)式求微分可得 \(S ′ _i (x)\)，然後做替換 $x = x_i $并化簡得：

\[S ′ _i (x_i) = − \frac{h_i} 3 M_i − \frac{h_i} 6 M_{i+1} − \frac{y_i} {h_i} + \frac{y_{i+1}} {h_i} \tag {2}\]

同理，可由(1)式得到$ S ′ _{i−1} (x)$，有：

\[S ′ _{i-1} (x_i) = \frac{h_{i-1}} 6 M_{i-1} + \frac{h_{i-1}} 3 M_i − \frac{y_{i-1}} {h_{i-1}} + \frac{y_i} {h_{i-1}} \tag {3}\]

當(2)和(3)右端項建立一個等式時，其結果可寫為：

\[h_{i−1}M_{i−1} + 2(h_i + h_{i−1})M_i + h_iM_{i+1} = \frac 6 h_i (y_{i+1} − y_i) − \frac 6 {h_{i−1}} (y_i − y_{i−1}) \tag {4}\]

這個等式僅僅對 \(i = 1, 2, . . . , n − 1\) 成立，進而給出了一個含有 n + 1 個未知量 $M_0, M_1, . . . , M_n \(的 n − 1 次方程的線性方程組。這時，取自然邊界條件，即\) M_0 = M_n = 0$，由此得到的樣條函數稱為自然三次樣條。

​ 對于 \(1 ≤ i ≤ n − 1\)，以及$ M_0 = M_n = 0$，方程組(4)是對稱的、三對角的和對角占優的，稱為三彎矩方程組，它具有下列 形式：

\[\begin{pmatrix}u_1&h_1\\

h_1&u_2&h_2\\

&h_2&u_3&h_3\\

&&\ddots&\ddots&\ddots\\

&&&h_{n-3}&u_{n-2}&h_{n-2}\\

&&&&h_{n-2}&u_{n-1}

\end{pmatrix}

\begin{pmatrix}

M_1\\

M_2\\

M_3\\

\vdots\\

M_{n-2}\\

M_{n-1}

\end{pmatrix}=

v_1\\

v_2\\

v_3\\

v_{n-2}\\

v_{n-1}

其中：

\[h_i = x_{i+1} − x_i\\ u_i = 2(h_i + h_{i−1})\\

b_i = \frac 6 {h_i} (y_{i+1} − y_i)\\

v_i = b_i − b_{i−1}

​ 可用追趕法求解。

簡化技巧

Hermite型插值多項式

• 兩點及其一階導數（切線）

Lidstone型插值多項式

• 兩點及二階導數（曲率）

三次基樣條

三次樣條曲線

樣條函數局限性

• 須有小擾度假定
• 不能處理多值問題
• 不能很好表達空間曲線
• 不具有坐标不變性（幾何不變）

三次參數樣條曲線

• 3個坐标分量 x，y，z 分别是參數 t 的三次樣條函數
• 對型值點做參數化
• 對3個坐标分量分别處理

曲線幾何連續性

參數連續性

在數學分析/高等數學中,我們所說的“連續性”(光滑性）是指“參數連續性”:

給定⒉條曲線 \(x_1(t)定義在[t_o,t_1],x_2(t)定義在[t_1,t_2]\)，曲線\(x_1\)和\(x_2\),在\(t_1\)稱為\(C^r\)連續的，如果它們的從\(0^{th}\)(0階)至\(r^{th}\)(r階）的導數向量在\(t_1\)處完全相同。

• C0: position varies continuously

• C1: First derivative is continuous across junction

  In other words: the velocity vector remains the same

• C2 : Second derivative is continuous across junction
• The acceleration vector remains the same

參數連續性的不足：

參數連續性過于嚴格，在幾何設計中不太直覺。

連續性依賴于參數的選擇，同一條曲線，參數不同，連續階也不同。

幾何連續性

幾何連續性定義

設\(\psi(t)(a ≤t≤b)\)是給定的曲線。

若存在一個參數變換\(t = \rho(s)(a_1≤s ≤b_1)\)，使得\(\psi(\rho(s))∈ C^n[a_1,b_1]\)，且\(\frac {d\psi(\rho(s))}{ds}≠0\)，

則稱曲線\(\psi(t)(a ≤t≤b)\)是n階幾何連續的曲線,記為

\[\psi(t)\in GC^n[a,b]或\psi(t)\in G^n[a,b]\]

幾何連續性性質

• 條件\(\frac {d\psi(\rho(s))}{ds}≠0\)保證曲線上無奇點
• 幾何連續性與參數選取無關，是曲線本身固有的幾何性質
• \(G^n\)的條件比\(C^n\)的寬，曲線類型更多

幾何連續性的具體形式

• \(G^0\)：表示兩曲線有公共的連接配接端點，與\(C^0\)的條件一緻
• \(G^1\)：兩曲線在連接配接點處有公共的切線，即切線方向連續
• \(G^2\)：兩曲線在連接配接點處有公共的曲率圓，即曲率連續

作業

模仿 PowerPoint 寫一個曲線設計與編輯工具

• 輸入有序點列（型值點），實時生成分段的三次樣條曲線
• 可修改拖動型值點的位置（保持整條曲線 \(C^2\)）
• 可編輯型值點處的切線資訊，成為 \(G^1\) 或 \(G^0\)

追趕法 （Thomas Algorithm）

形如：

\[\begin {cases}b_1x_1+c_1x_2 &=d_1\\

a_2x_1+b_2x_2+c_2x_3 &=d_2\\

\ldots\\

a_kx_{k-1}+b_kx_k+c_kx_{k+1}&=d_k\\

a_{n-1}x_{k-2}+b_{n-1}x_{n-1}+c_{n-1}x_{n}&=d_{n-1}\\

a_nx_{n-1}+b_nx_n &=d_n\\

\end {cases}\tag {1}

系數矩陣A為三對角矩陣：

\[A=\begin{pmatrix}b_1&c_1\\

a_2&b_2&c_2\\

&\ddots&\ddots&\ddots\\

&&a_k&b_k&c_k\\

&&&\ddots&\ddots&\ddots\\

&&&&a_{n-1}&b_{n-1}&c_{n-1}\\

&&&&&a_{n}&b_n

\tag {2}

追趕法是高斯消去法的一種簡化形式，分消元與回代兩個過程。

推導一：

先将（1）第一個方程中\(x_1\)的系數化為1

\[x_1+\frac{c_1}{b_1}x_2=\frac{d_1}{b_1}\\[4ex]x_1+r_1x_2=y_1\\[2ex]

r_1=\frac{c_1}{b_1},y_1=\frac{d_1}{b_1}

對第二個方程消元\(x_1\)

\[x_2+\frac{c_2}{b_2-r_1a_2}x_3=\frac{d_2-y_1a_2}{b_2-r_1a_2}\\[4ex]\]

一步一步地順序消元，得第k個方程為

\[x_k+\frac{c_k}{b_k-r_{k-1}a_k}x_{k+1}=\frac{d_k-y_{k-1}a_k}{b_k-r_{k-1}a_k}\quad k=2,3,...,n-1\]

記

\[r_k=\frac{c_k}{b_k-r_{k-1}a_k} \quad y_k=\frac{d_k-y_{k-1}a_k}{b_k-r_{k-1}a_k}

得到

\[x_k+r_kx_{k+1}=y_k\]

則第n-1個方程為\(x_{n-1}+r_{n-1}x_{n}=y_{n-1}\)，聯立第n個方程的原方程：\(a_nx_{n-1}+b_nx_n =d_n\)，得出

\[x_n=y_n\]

得到方程組的簡單表達：

\[\begin {cases}x_1+r_1x_2=y_1\\

x_k+r_kx_{k+1}=y_k\\

x_n=y_n

\end {cases}

\tag{3}

系數矩陣中的非零元素集中釋出在主對角線和一條次主對角線上

\[\begin {pmatrix}1&r_1\\

&1&r_2\\

&&\ddots&\ddots\\

&&&1&r_{n-1}\\

&&&&1

\end {pmatrix}

對方程組(3)自下而上逐漸回代，依次求出\(x_n,x_{n-1},\dots,x_1\)：

\[\begin {cases}x_n=y_n\\

x_k=y_k-r_kx_{k+1}&k=n-1,n-2,\dots,1

推導二：

有方程組\(Ax=d\)

設矩陣\(A\)非奇異，\(A\)有Crout分解\(A=LU\)，記

\[L=\begin {pmatrix}\beta_1\\

\alpha_2&\beta_2\\

&\ddots&\ddots\\

&&\alpha_{n-1}&\beta_{n-1}\\

&&&\alpha_n&\beta_n\\

\end {pmatrix} \quad

U=\begin {pmatrix}

1&\gamma_1\\

&1&\gamma_2\\

&&&1&\gamma_{n-1}\\

由矩陣乘法，得到:

\[\beta_1=b_1\\\beta_1*\gamma_1=c_1 \rightarrow \gamma_1=\frac{c_1}{\beta_1}\\

\alpha_i=a_i(2\le i\le n)\\

b_i=\alpha_i\gamma_{i-1}+\beta_i\rightarrow \beta_i=b_i-a_i\gamma_{i-1}(2\le i\le n)\\

\beta_i*\gamma_i=c_i \rightarrow \gamma_i=\frac{c_i}{\beta_i}\\

引入中間變量\(y\)，令\(y=Ux\)，則有

\[Ax=LUx=Ly=d\\Ly=d

已知L和d，由矩陣乘法，得到:

\[y_1=\frac{d_1}{\beta_1}\\y_i=\frac{d_i-\alpha_iy_{i-1}}{\beta_i}\\ (2\le i\le n)

由\(y=Ux\)，得\(X\)的表達式：

\[x_n=y_n\\x_i=y_i-\gamma_ix_{i+1}\\(1 \le i \le n-1)

綜上，得到各量的遞推公式：

\[\begin {cases}\beta_1=b_1 \quad \gamma_1=\frac{c_1}{\beta_1} \quad y_1= \frac{d_1}{\beta_1}\\

\alpha_i=a_i \quad \beta_i=b_i-a_i\gamma_{i-1} \quad \gamma_i=\frac{c_i}{\beta_i} \quad y_i=\frac{d_i-a_iy_{i-1}}{\beta_i} &i=2,3,\dots,n

\end {cases}\\

\begin {cases}

x_n=y_n\\

x_i=y_i-\gamma_ix_{i+1}&i=n-1,n-2,\dots,1

C2連續

先将曲線參數化，再使用追趕法求解\(M\)：

\[\begin{pmatrix}u_1&h_1\\h_1&u_2&h_2\\&h_2&u_3&h_3\\&&\ddots&\ddots&\ddots\\&&&h_{n-3}&u_{n-2}&h_{n-2}\\&&&&h_{n-2}&u_{n-1}\end{pmatrix}\begin{pmatrix}M_1\\M_2\\M_3\\\vdots\\M_{n-2}\\M_{n-1}\end{pmatrix}=\begin{pmatrix}v_1\\v_2\\v_3\\\vdots\\v_{n-2}\\v_{n-1}\end{pmatrix}\]

\[h_i = x_{i+1} − x_i\\ u_i = 2(h_i + h_{i−1})\\ b_i = \frac 6 {h_i} (y_{i+1} − y_i)\\ v_i = b_i − b_{i−1}\]

然後代入以下到

\[S_i(x) = \frac{M_i}{6h_i} (x_{i+1} − x)^3 + \frac{M_{i+1}}{6h_i} (x − x_i)^3 + (\frac{y_{i+1}}{hi} − \frac{M_{i+1}h_i}6 )(x − x_i) + ( \frac{y_i}{h_i} − \frac{M_ih_i} 6 )(x_{i+1} − x)\]

部分型值點G1連續

G1連續：對于給定型值點處的左導數與右導數在同一直線上。

對給定型值點的左右兩段曲線（若有），已知兩個端點及在端點處的一階導數，使用Hermite插值多項式進行求解。

Hermite插值多項式：

\[h_0(x)=(1+2\frac{x-x_0}{x_1-x_0})(\frac{x-x_1}{x_0-x_1})^2\\h_1(x)=(1+2\frac{x-x_1}{x_0-x_1})(\frac{x-x_0}{x_1-x_0})^2\\

H_0(x)=(x-x_0)(\frac{x-x_1}{x_0-x_1})^2\\

H_1(x)=(x-x_1)(\frac{x-x_0}{x_1-x_0})^2\\

H(x)=h_0y_0+h_1y_1+H_0y'_0+H_1y'_1

部分型值點C0連續

G0連續：表示兩曲線有公共的連接配接端點。

同G1，使用Hermite插值多項式進行求解。