CG-三次样条函数

推导

参数连续性

给定2条曲线：
$$
\boldsymbol{F_1}(t), t \in [t_0, t_1]
$$

$$
\boldsymbol{F_2}(t), t \in [t_1, t_2]
$$

曲线$\boldsymbol{F_1}(t)$和$\boldsymbol{F_2}(t)$是$C^r$连续的，如果它们从$0$阶导数至$r$阶导数在$t_1$处完全相同

根据定义有：

$C^0$连续：曲线上的点连续，不存在尖点

$C^1$连续：曲线上点的一阶导数存在且连续

$C^2$连续：曲线上点的二阶导数存在且连续

几何连续性

由于参数连续性过于严格，在集合设计中不太直观，连续性依赖于参数的选择，同一条曲线，参数不同，连续阶也不同。但是，可以通过引入参数变换使其满足参数连续。

设$\varphi (t) (a \leq t \leq b)$ 是给定的曲线。若存在一个参数变换$ t = \rho (s)(a_1 \leq s \leq b_1)$, 使得$\varphi(\rho(s)) \in C^n[a_1, b_1]$, 且$\frac{d\varphi(\rho(s))}{ds} \neq 0$，则称曲线$\varphi (t) (a \leq t \leq b)$是$n$阶几何连续的曲线，记为：
$$
\varphi(t) \in GC^n[a, b] \ or \ \varphi(t) \in G^n[a,b]
$$
性质：

• 条件$\frac{d\varphi(\rho(s))}{ds} \neq 0$保证了曲线上没有奇点
• 几何连续性与参数选择无关，是曲线本身固有的性质
• $G^n$的条件比$C^n$要宽，曲线类型更多

根据定义有：

$G^0$连续：两条曲线有公共的端点

$G^1$连续：两条曲线在连接点处有公共的切线，即切线方向连续

$G^2$连续：两条曲线在连接点处有公共的曲率圆，即曲率连续

$C^2$连续

对于 $ R^{2} $ 内的一组点集$$ {P_i | i \in [1, n]} $$, 取参数集合为$${ t_i | t \in [1, n] }$$，假设经过点$$P_i = (x_i, y_i)$$ 和$$P_{i+1} = (x_{i+1}, y_{i+1})$$ 的曲线方程为$$\boldsymbol{F_i}(t) = (X_i(t), Y_i(t))$$, 则：
$$
\boldsymbol{F_i}(t) = \boldsymbol{a_i} + \boldsymbol{b_i}(t - t_i) + \boldsymbol{c_i}(t - t_i)^2 + \boldsymbol{d_i} (t - t_i)^3 \tag1
$$

$$
\boldsymbol{F^{‘}_i}(t) = \boldsymbol b_i + 2\boldsymbol{c_i}(t - t_i) + 3\boldsymbol{d_i}(t-t_i)^2 \tag2
$$

$$
\boldsymbol{F^{\‘\‘}_i}(t) = 2\boldsymbol{c_i} + 6\boldsymbol{d_i}(t-t_i) \tag3
$$

曲线经过点$P_i$, 则:
$$
\boldsymbol{a_i} = (x_i, y_i) = \boldsymbol{P_i} \tag4
$$
曲线经过点$$P_{i+1}$$且$C^0$连续，令$$h_{i} = t_{i+1}-t_{i}$$, 则：
$$
\boldsymbol{a_{i+1}} = \boldsymbol{a_i} + \boldsymbol{b_i}hi + \boldsymbol{c_i}h_i^2 + \boldsymbol{d_i}h_i^3 \tag5
$$
由于$C^1$连续，则：
$$
\boldsymbol{F^{‘}i}(t{i+1}) = \boldsymbol b_i + 2\boldsymbol{c_i}(t_{i+1} - t_i) + 3\boldsymbol{d_i}(t_{i+1}-t_i)^2 = \boldsymbol b_i + 2\boldsymbol{c_i}h_i + 3\boldsymbol{d_i}h_i^2 \tag6
$$

$$
\boldsymbol{F^{‘}{i+1}}(t{i+1}) = \boldsymbol b_{i+1} + 2\boldsymbol{c_{i+1}}(t_{i+1} - t_{i+1}) + 3\boldsymbol{d_{i+1}}(t_{i+1}-t_{i+1})^2 = \boldsymbol b_{i+1} \tag7
$$

联立(6)(7)式可得：
$$
\boldsymbol {b_{i+1}} = \boldsymbol b_i + 2\boldsymbol{c_i}h_i + 3\boldsymbol{d_i}h_i^2 \tag8
$$
由于$C^2$连续，则：
$$
2\boldsymbol{c_i} + 6\boldsymbol{d_i} h_i = 2\boldsymbol{c_{i+1}} \tag9
$$
设$\boldsymbol{m_i}$ 为插值点$ P_i $的二阶导数，则
$$
\boldsymbol{m_i} = 2\boldsymbol{c_i} \tag{10}
$$
根据(9)式得：
$$
\boldsymbol{m_i} + 6\boldsymbol{d_i} h_i = \boldsymbol{m_{i+1}} \tag{11}
$$
得到：
$$
\boldsymbol{d_i} = \frac{\boldsymbol {m_{i+1}} - \boldsymbol{m_i}}{6h_i} \tag{12}
$$
联立(5)(12)式可得：
$$
\boldsymbol {b_i} = \frac{\boldsymbol {P_{i+1}} - \boldsymbol {P_{i}} }{h_i} - \frac{\boldsymbol{m_i}h_i}{2} - \frac{(\boldsymbol{m_{i+1}} - \boldsymbol {m_{i}} ) h_i}{6} \tag{13}
$$
联立(4)(8)(10)(12)(13)可得：
$$
\boldsymbol{m_i} h_i + 2\boldsymbol {m_{i+1}}(h_i + h_{i+1}) + \boldsymbol {m_{i + 2}}h_{i + 1} = 6(\frac{\boldsymbol {P_{i+2}} - \boldsymbol{P_{i+1}} }{h_{i+1}} - \frac{\boldsymbol {P_{i+1}} - \boldsymbol{P_{i}}}{h_i} ), i \in [1, n-2] \tag{14}
$$
此时还缺少两个方程才能解出曲线的未知参数，有如下约束可以选择：

• 自由端：指定曲线在两个端点的二阶导数，即
  $$
  \boldsymbol {m_1} = \boldsymbol {m_n} = \boldsymbol {0} \tag{15}
  $$
  该样条为自然三次样条

• 夹持端：指定曲线在两个端点的一阶导数

• 其他：

  • 曲线两端为抛物线
  • 周期端
  • 混合边界条件

这里以自然边界条件为例，则线性方程组为：
$$
\begin{bmatrix}
1 & 0 & 0 & 0 & \cdots & 0 \
h_1 & 2(h_1 + h_2) & h_2 & 0 & \cdots & 0 \
0 & h_2 & 2(h_2 + h_3) & h_3 & \cdots & 0 \
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots \
0 & \cdots & 0 & h_{n-2} & 2(h_{n-2} + h_{n-1}) & h_{n-1} \
0 & \cdots & 0 & 0 & 0 & 1\
\end{bmatrix}

\begin{bmatrix}
m_1 \
m_2 \
m_3 \
\vdots \
m_{n-1} \
m_n \
\end{bmatrix}

6
\begin{bmatrix}
0 \
\frac{\boldsymbol {P_3} - \boldsymbol{P_2} }{h_2} - \frac{\boldsymbol {P_2} - \boldsymbol{P_1} }{h_1} \
\frac{\boldsymbol {P_4} - \boldsymbol{P_3} }{h_3} - \frac{\boldsymbol {P_3} - \boldsymbol{P_2} }{h_2} \
\vdots \
\frac{\boldsymbol {P_n} - \boldsymbol{P_{n-1}} }{h_{n-1}} - \frac{\boldsymbol {P_{n-1}} - \boldsymbol{P_{n-2}} }{h_{n-2}} \
0 \
\end{bmatrix}
$$
将该方程组的第1行第1列，最后1行最后1列去掉，满足三对角矩阵和对角线占优，可以使用追赶法求解方程组，具体方法见：

Tridiagonal matrix algorithm

实现

github代码：

https://github.com/olleh-dlrow/Vortex/blob/master/Sandbox/src/CubicSplinesTest.h

数据结构

型值点具有如下信息：

• 位置
• 一阶左导数
• 一阶右导数
• 参数t
• 几何连续性类型

绘制自然三次样条

输入：型值点

输出：拟合后的三次样条曲线，型值点的参数值、导数值

根据公式计算每段曲线的未知参数，并初始化型值点的参数值和导数值，最后根据曲线计算若干离散点进行绘制

编辑型值点

选中型值点后，进入编辑模式，可以通过调整切线的控制点改变该点的切线值，根据点的几何连续性类型有如下情况：

$G^0$连续：切线无约束

$G^1$连续：左右切线向量的方向相同

$G^2 $连续：左右切线向量相同

移动型值点

在编辑模式下，拖动型值点即可移动该点，此时需要更新型值点两端的曲线，根据曲线两端点的函数值和一阶导数得到4个方程，可以解得曲线$$\boldsymbol{F_i} (t)$$的参数如下，其中$$h_{i} = t_{i+1} - t_{i}$$：
$$
\boldsymbol{a_i} = (x_i, y_i) = \boldsymbol{F_i}(t_i)
$$

$$
\boldsymbol{b_i} = \boldsymbol{F^{\‘}_{i+}}(t_i)
$$

$$
\boldsymbol{d_i} = \frac{2\boldsymbol{F_i}(t_i) + \boldsymbol{F^{\‘}{i+}}(t_i) h_i + \boldsymbol{F^{\‘}{i-}}(t_{i+1})h_i - 2\boldsymbol{F_i}(t_{i+1})}{h_i^3}
$$

$$
\boldsymbol{c_i} = \frac{\boldsymbol{F^{\‘}{i-}}(t{i+1}) - \boldsymbol{F^{\‘}{i+}}(t{i}) - 3\boldsymbol{d_i}h_i^2 }{2h_i}
$$

添加型值点

在完成初始化后，添加新的型值点不再需要保证$C^2$连续， 此时根据末端点$$P_{i}$$和新点$$P_{i+1}$$生成曲线$$\boldsymbol{F_i}(t)$$，需要4个约束方程来计算未知参数，一种方法是，使用两个点的函数值和$$P_i$$在$$\boldsymbol{F_{i-1}}(t)$$的一阶、二阶导数建立4个方程，该方法经过实现效果并不好，生成的新曲线容易出现极大值。这里介绍第二种方法，使用自然端约束，即$$P_{i+1}$$的二阶导数为$$0$$

根据条件建立4个方程，可以解得4个未知参数如下，其中$$h_i = t_{i+1} - t_i$$：
$$
\boldsymbol{a_i} = (x_i, y_i) = \boldsymbol{F_i}(t_i)
$$

$$
\boldsymbol{b_i} = \boldsymbol{F^{‘}{i+}}(t_i) = \boldsymbol{F^{‘}{i-}}(t_i)
$$

$$
\boldsymbol{d_i} = \frac{\boldsymbol{a_i} + \boldsymbol{b_i}h_i - \boldsymbol{F_i}(t_{i+1})}{2h_i^3}
$$

$$
\boldsymbol{c_i} = -3\boldsymbol{d_i}h_i
$$

该方法生成的切线向量长度适中，比较适合进一步编辑

删除型值点

鼠标悬停在某个型值点上后，按下E键即可删除该点，此时需要更新所有点的参数$t$，并重新计算曲线，方法同移动型值点处的计算方法