几何建模与处理之四 三次样条曲线

几何建模与处理之四 三次样条曲线

目录

• 几何建模与处理之四 三次样条曲线
  • 几何设计
    • 样条曲线
    • 自由曲线
  • 样条曲线的数学表达推导
    • 力学解释
    • 数学性质
    • 求解思路
      • 三次样条插值函数
      • 三弯矩方程
      • 简化技巧
  • 三次样条曲线
    • 样条函数局限性
    • 三次参数样条曲线
  • 曲线几何连续性
    • 参数连续性
    • 几何连续性
      • 几何连续性定义
      • 几何连续性性质
      • 几何连续性的具体形式
  • 作业
    • 追赶法 （Thomas Algorithm）
      • 推导一：
      • 推导二：
    • C2连续
    • 部分型值点G1连续
    • 部分型值点C0连续

几何设计

概念设计 >> 设计图纸 >> 数学（几何）表达

样条曲线

在早期没有计算机时，工业产品设计中工程师使用压铁使富有弹性的样条通过指定的样点，调整样条使它具有满意的形状，然后沿样条画出曲线，称为样条曲线。

• 样条：具有一定弹性的软木条

• 压铁：较重的铁块，用来固定样条所经过的点

  

自由曲线

在产品初始设计阶段，描述其外形的曲线或曲面常常只有大致的形状或只知道它通过的一些空间列（称为型值点），这类没有数学表达式的曲线或曲面称为自由曲线或自由曲面（Free from curve/surface）。

造型方法

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

样条曲线的数学表达推导

力学解释

曲线：弹性梁，弹性模量为E

压铁：载荷（力）

由材料力学中的贝努利-欧拉方程

[M(x)=EIk(x)=EIfrac{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.其次，假设曲线整体为C²连续，则相邻两端在拼接点要满足3个条件（C⁰连续、C¹连续，C²连续），共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) ag {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} ag {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}} ag {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}) ag {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)是对称的、三对角的和对角占优的，称为三弯矩方程组，它具有下列 形式：

[egin{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} egin{pmatrix} M_1\ M_2\ M_3\ vdots\ M_{n-2}\ M_{n-1} end{pmatrix}= egin{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} ]

​ 可用追赶法求解。

简化技巧

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)处完全相同。

• C⁰: position varies continuously

• C¹: First derivative is continuous across junction

  In other words: the velocity vector remains the same

• C² : Second derivative is continuous across junction

• The acceleration vector remains the same

参数连续性的不足：

参数连续性过于严格，在几何设计中不太直观。

连续性依赖于参数的选择，同一条曲线，参数不同，连续阶也不同。

几何连续性

几何连续性定义

设(psi(t)(a ≤t≤b))是给定的曲线。

若存在一个参数变换(t = ho(s)(a_1≤s ≤b_1))，使得(psi( ho(s))∈ C^n[a_1,b_1])，且(frac {dpsi( ho(s))}{ds}≠0)，

则称曲线(psi(t)(a ≤t≤b))是n阶几何连续的曲线,记为

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

几何连续性性质

• 条件(frac {dpsi( ho(s))}{ds}≠0)保证曲线上无奇点
• 几何连续性与参数选取无关，是曲线本身固有的几何性质
• (G^n)的条件比(C^n)的宽，曲线类型更多

几何连续性的具体形式

• (G^0)：表示两曲线有公共的连接端点，与(C^0)的条件一致
• (G^1)：两曲线在连接点处有公共的切线，即切线方向连续
• (G^2)：两曲线在连接点处有公共的曲率圆，即曲率连续

作业

模仿 PowerPoint 写一个曲线设计与编辑工具

• 输入有序点列（型值点），实时生成分段的三次样条曲线
• 可修改拖动型值点的位置（保持整条曲线 (C^2)）
• 可编辑型值点处的切线信息，成为 (G^1) 或 (G^0)

追赶法 （Thomas Algorithm）

形如：

[egin {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\ ldots\ 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} ag {1} ]

系数矩阵A为三对角矩阵：

[A=egin{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 end{pmatrix} ag {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 ]

得到方程组的简单表达：

[egin {cases} x_1+r_1x_2=y_1\ vdots\ x_k+r_kx_{k+1}=y_k\ vdots\ x_n=y_n end {cases} ag{3} ]

系数矩阵中的非零元素集中发布在主对角线和一条次主对角线上

[egin {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)：

[egin {cases} x_n=y_n\ x_k=y_k-r_kx_{k+1}&k=n-1,n-2,dots,1 end {cases} ]

推导二：

有方程组(Ax=d)

设矩阵(A)非奇异，(A)有Crout分解(A=LU)，记

[L=egin {pmatrix} eta_1\ alpha_2&eta_2\ &ddots&ddots\ &&alpha_{n-1}&eta_{n-1}\ &&&alpha_n&eta_n\ end {pmatrix} quad U=egin {pmatrix} 1&gamma_1\ &1&gamma_2\ &&ddots&ddots\ &&&1&gamma_{n-1}\ &&&&1 end {pmatrix} ]

由矩阵乘法，得到:

[eta_1=b_1\ eta_1*gamma_1=c_1 ightarrow gamma_1=frac{c_1}{eta_1}\ alpha_i=a_i(2le ile n)\ b_i=alpha_igamma_{i-1}+eta_i ightarrow eta_i=b_i-a_igamma_{i-1}(2le ile n)\ eta_i*gamma_i=c_i ightarrow gamma_i=frac{c_i}{eta_i}\ ]

引入中间变量(y)，令(y=Ux)，则有

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

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

[y_1=frac{d_1}{eta_1}\ y_i=frac{d_i-alpha_iy_{i-1}}{eta_i}\ (2le ile n) ]

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

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

综上，得到各量的递推公式：

[egin {cases} eta_1=b_1 quad gamma_1=frac{c_1}{eta_1} quad y_1= frac{d_1}{eta_1}\ alpha_i=a_i quad eta_i=b_i-a_igamma_{i-1} quad gamma_i=frac{c_i}{eta_i} quad y_i=frac{d_i-a_iy_{i-1}}{eta_i} &i=2,3,dots,n end {cases}\ egin {cases} x_n=y_n\ x_i=y_i-gamma_ix_{i+1}&i=n-1,n-2,dots,1 end {cases} ]

C²连续

先将曲线参数化，再使用追赶法求解(M)：

[egin{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}egin{pmatrix}M_1\M_2\M_3\vdots\M_{n-2}\M_{n-1}end{pmatrix}=egin{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) ]

部分型值点G¹连续

G¹连续：对于给定型值点处的左导数与右导数在同一直线上。

对给定型值点的左右两段曲线（若有），已知两个端点及在端点处的一阶导数，使用Hermite插值多项式进行求解。

Hermite插值多项式：

[h_0(x)=(1+2frac{x-x_0}{x_1-x_0})(frac{x-x_1}{x_0-x_1})^2\ h_1(x)=(1+2frac{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 ]

部分型值点C⁰连续

G⁰连续：表示两曲线有公共的连接端点。

同G¹，使用Hermite插值多项式进行求解。