CAGD 学习笔记 | 有理曲线

在前几篇笔记中, 我们已经介绍了 Bézier 曲线和样条曲线, 它们分别是多项式曲线和分段多项式曲线, 其数学意义在于构造可以任意地逼近任何光滑曲线的参数曲线. 但是我们并不满足于是逼近, 我们希望精确地表示更多种类的曲线, 如圆锥曲线.

圆锥曲线的参数方程

以最简单的一类圆锥曲线为例: 二维平面中的单位圆, 其隐式表达式为 $x^2+y^2 = 1$, 很容易验证它不能被任何多项式精确表示. 但是可以用三角函数给出其参数方程

$$ x = \cos \theta,~ y = \sin \theta,~ \theta \in [0, 2\pi) $$

考虑三角函数的万能公式

$$ \cos \theta = \frac{1 - \tan ^2 \frac{\theta}{2}}{1 + \tan ^2 \frac{\theta}{2}},~ \sin \theta = \frac{2 \tan \frac{\theta}{2}}{1 + \tan ^2 \frac{\theta}{2}} $$

令 $t = \tan \frac{\theta}{2}$, 可以得到圆的另一个方程

$$ x = \frac{1-t^2}{1+t^2},~ y = \frac{2t}{1+t^2},~ t \in (-\infty, \infty) $$

形如这样的由多项式除法的商得到的函数被我们称为有理函数 (Rational Function), 有理曲线就是利用有理函数构造的曲线.

注 这里圆的参数方程缺少了一个点 $(-1,0)$, 我们可以通过分段的方式解决这一问题.

不难注意到, 其它的圆锥曲线也可以通过三角函数或双曲函数及对应的万能公式用有理函数来表示, 因此有理曲线包含了多项式曲线和圆锥曲线.

除法的意义

仍然以二维平面中的单位圆为例, 考虑三维空间中的多项式曲线

$$ x = 1-t^2,~ y = 2t,~ z = 1+t^2,~ t \in (-\infty, \infty) $$

以原点为相机, 在平面 $z=1$ 上做透视投影, 得到

$$ x = \frac{1-t^2}{1+t^2},~ y = \frac{2t}{1+t^2},~ z = 1,~ t \in (-\infty, \infty) $$

发现有理函数中多项式除法的意义可以解释为透视投影, 而单位圆就是三维空间中多项式曲线在二维平面上的投影. 类似地, 这也适用于任何其它的有理曲线.

齐次坐标

齐次坐标指的是一个 $n$ 维矢量的 $n+1$ 维表示, 它可以用于区分向量和点, 并将仿射变换和透视投影变换表示为线性变换.

对于 $d$ 维欧氏空间中的一个点, 若它的笛卡尔坐标 (Cartesian Coordinates) 为 $\boldsymbol{x} \in \mathbb{R}^d$, 则对应的齐次坐标 (Homogenous Coordinates) 是 $\bar{\boldsymbol{x}} = \left(\begin{array}{c} \boldsymbol{x} \\ 1 \end{array}\right)$ 或 $\left(\begin{array}{c} \omega\boldsymbol{x} \\ \omega \end{array}\right), \omega \neq 0$.

二维单位圆的参数方程在笛卡尔坐标下表示为

$$ \boldsymbol{p}(t) = \left( \begin{array}{c} (1-t^2)/(1+t^2) \\ 2t/(1+t^2) \end{array} \right) $$

则在齐次坐标下表示为

$$ \bar{\boldsymbol{p}}(t) = \left( \begin{array}{c} 1-t^2 \\ 2t \\ 1+t^2 \end{array} \right) $$

这样在齐次坐标下, 有理曲线就可以表示为多项式曲线的形式.

有理 Bézier 曲线

结合 Bézier 曲线与有理曲线, 我们得到有理 Bézier 曲线, 它的齐次坐标形式为

$$ \bar{\boldsymbol{p}}(t)=\sum_{k=0}^n B_k^{n}(t) \boldsymbol{q}_k,~ \boldsymbol{q}_k \in \mathbb{R}^{d+1},~ \boldsymbol{q}_k = \left(\begin{array}{c} \omega_k \boldsymbol{p}_k \\ \omega_k \end{array}\right),~ t \in [0,1] $$

它的笛卡尔坐标形式为

$$ \boldsymbol{p}(t)=\frac{\sum_{k=0}^n B_k^{n}(t)\omega_k\boldsymbol{p}_k}{\sum_{k=0}^n B_k^{n}(t)\omega_k} = \sum_{k=0}^n \frac{B_k^{n}(t)\omega_k}{\sum_{l=0}^n B_l^{n}(t)\omega_l} \boldsymbol{p}_k =: \sum_{k=0}^n R_k(t) \boldsymbol{p}_k,~ t \in [0,1] $$

注意到一些 Bézier 曲线的性质对于有理 Bézier 曲线仍然成立:

• 仿射无关性: 显然成立;
• 凸包性: 因为 $\sum_k^n R_k(t) = 1$, 只要 $\forall k, \omega_k>0$, 就有凸包性成立.
• 起点和终点处的导数与 $p_1-p_0$ 和 $p_n - p_{n-1}$ 平行: 对笛卡尔坐标形式的参数方程求导即可发现是容易证明的.

从上文中的讨论中我们知道, 有理 Bézier 曲线的本质是 Bézier 曲线的投影. 当 $\forall k, \omega_k = 1$ 时, 有理 Bézier 曲线退化为 Bézier 曲线.

我们可以将 $\omega_k$ 直观地理解为每一个控制点的权重, $\omega_k$ 增大时曲线会更靠近 $p_k$.

2 阶有理 Bézier 曲线

因为 2 阶有理 Bézier 曲线已经足以表示圆锥曲线, 我们可以可以从它们入手进行更深入的研究.

标准形式

考虑 2 阶有理 Bézier 曲线

$$ \boldsymbol{p}(t)= \frac{B_0^{2}(t) \omega_0 \boldsymbol{p}_0+B_1^{2}(t) \omega_1 \boldsymbol{p}_1+B_2^{2}(t) \omega_2 \boldsymbol{p}_2} {B_0^{2}(t) \omega_0+B_1^{2}(t) \omega_1+B_2^{2}(t) \omega_2},~ \boldsymbol{p}_k \in \mathbb{R}^d, \omega_k \in \mathbb{R} $$

作重参数化 (reparameterization)

$$ t = \frac{\tilde{t}}{\alpha(1-\tilde{t})+\tilde{t}} $$

$t$ 与 $\tilde{t}$ 的关系如下图所示. 则

$$ \begin{aligned} \boldsymbol{p}(t)&=\frac{\alpha^2(1-\tilde{t})^2 \omega_0 \boldsymbol{p}_0+2 \alpha(1-\tilde{t}) \tilde{t} \omega_1 \boldsymbol{p}_1+\tilde{t}^2 \omega_2 \boldsymbol{p}_2}{\alpha^2(1-\tilde{t})^2 \omega_0+2 \alpha(1-\tilde{t}) \tilde{t} \omega_1+\tilde{t}^2 \omega_2} \\ &=\frac{\alpha^2 B_0^{2}(\tilde{t}) \omega_0 \boldsymbol{p}_0+\alpha B_1^{2}(\tilde{t}) \omega_1 \boldsymbol{p}_1+B_2^{2}(\tilde{t}) \omega_2 \boldsymbol{p}_2}{\alpha^2 B_0^{2}(\tilde{t}) \omega_0+\alpha B_1^{2}(\tilde{t}) \omega_1+B_2^{2}(\tilde{t}) \omega_2} \end{aligned} $$

假设 $0 \leq \frac{\omega_2}{\omega_0}<\infty$, 令 $\alpha = \sqrt{\frac{\omega_2}{\omega_0}}$, 则有

$$ \boldsymbol{p}(t) = \frac{B_0^{2}(\tilde{t}) \boldsymbol{p}_0+B_1^{2}(\tilde{t}) \sqrt{\frac{1}{\omega_0 \omega_2}} \omega_1 \boldsymbol{p}_1+B_2^{2}(\tilde{t}) \boldsymbol{p}_2}{B_0^{2}(\tilde{t})+B_1^{2}(\tilde{t}) \sqrt{\frac{1}{\omega_0 \omega_2}} \omega_1+B_2^{2}(\tilde{t})} $$

令 $\omega = \sqrt{\frac{1}{\omega_0 \omega_2}} \omega_1$, 则上式为

$$ \boldsymbol{p}(t) = \frac{B_0^{2}(\tilde{t}) \boldsymbol{p}_0+B_1^{2}(\tilde{t}) \omega \boldsymbol{p}_1+B_2^{2}(\tilde{t}) \boldsymbol{p}_2}{B_0^{2}(\tilde{t})+B_1^{2}(\tilde{t}) \omega + B_2^{2}(\tilde{t})} $$

发现 2 阶有理 Bézier 曲线的权重事实上只有一个自由度, W.L.O.G. 我们可以设 $\omega_0 = \omega_2 = 1$. 这一形式被称为 2 阶有理 Bézier 曲线的标准形式.

定理 对于标准形式的 2 阶有理 Bézier 曲线, 其形状取决于 $\omega$ 的值:

• $0 < \omega < 1$ 时, 曲线是椭圆;
• $\omega = 1$ 时, 曲线是抛物线;
• $\omega > 1$ 时, 曲线是双曲线.

证明思路是, $\boldsymbol{p}(t) = R_0(t) \boldsymbol{p}_0 + R_1(t) \boldsymbol{p}_1 + R_2(t) \boldsymbol{p}_2$, 其中有 $R_0(t) + R_1(t) + R_2(t) = 1$, 这可以理解为曲线上一点在三个控制点的三角形中的重心坐标. 任取其中两个值 $R_{i_1}, R_{i_2}$, 都可以仿射变换为平面上的坐标, 因此只需研究这两者满足的二次型即可. 可惜这里空间太小, 具体的证明写不下.

一般有理曲线转换为有理 Bézier 曲线

仍然以二维平面中的单位圆为例, 它在齐次坐标下的参数方程为

$$ \bar{\boldsymbol{p}}(t) = \left( \begin{array}{c} 1-t^2 \\ 2t \\ 1+t^2 \end{array} \right) $$

通过基的变换或是极形式的方法, 可以得到

$$ \bar{\boldsymbol{p}}(t) = B_0^{2}(t)\left(\begin{array}{l} 1 \\ 0 \\ 1 \end{array}\right)+B_1^{2}(t)\left(\begin{array}{l} 1 \\ 1 \\ 1 \end{array}\right) +B_2^{2}(t)\left(\begin{array}{l} 0 \\ 2 \\ 2 \end{array}\right) $$

通过上文中的推导, 我们知道可以作重参数化

$$ t = \frac{\tilde{t}}{\alpha(1-\tilde{t})+\tilde{t}},~ \alpha = \sqrt{\frac{\omega_2}{\omega_0}} = \sqrt{2} $$

以及

$$\omega = \sqrt{\frac{1}{\omega_0 \omega_2}} \omega_1 = \frac{\sqrt{2}}{2}$$

得到标准形式

$$ \boldsymbol{p}(t)=\frac{B_0^{2}(t)\left(\begin{array}{l} 1 \\ 0 \end{array}\right)+\frac{\sqrt{2}}{2} B_1^{2}(t)\left(\begin{array}{l} 1 \\ 1 \end{array}\right)+B_2^{2}(t)\left(\begin{array}{l} 0 \\ 1 \end{array}\right)}{B_0^{2}(t)+\frac{\sqrt{2}}{2} B_1^{2}(t)+B_2^{2}(t)} $$

一般情况下单位圆的有理 Bézier 曲线可以表示如下

对偶圆锥曲线

考虑标准形式的 2 阶有理 Bézier 曲线

$$ \boldsymbol{p}(t) = \frac{(1-t)^2 \boldsymbol{p}_0+ 2t(1-t) \omega \boldsymbol{p}_1+ t^2 \boldsymbol{p}_2}{(1-t)^2+ 2t(1-t) \omega + t^2},~ t \in [0,1] $$

令 $\hat{t} = s(t) = \frac{t}{2t-1}$, 则 $t \in [0, \frac{1}{2})$ 时 $\hat{t} \in (-\infty, 0]$, $t \in (\frac{1}{2}$ 时 $\hat{t} \in [0, \infty)$. 进而有

$$ \begin{aligned} \boldsymbol{q}(t) := \boldsymbol{p}(\hat{t}) & = \boldsymbol{p}(s(t)) = \frac{(1-\hat{t})^2 \boldsymbol{p}_0+ 2\hat{t}(1-\hat{t}) \omega \boldsymbol{p}_1+ \hat{t}^2 \boldsymbol{p}_2}{(1-\hat{t})^2+ 2\hat{t}(1-\hat{t}) \omega + \hat{t}^2} \\ & = \frac{(1-t)^2 \boldsymbol{p}_0- 2t(1-t) \omega \boldsymbol{p}_1+ t^2 \boldsymbol{p}_2}{(1-t)^2- 2t(1-t) \omega + t^2},~ \end{aligned} $$

$\boldsymbol{q}(t)$ 称为 $\boldsymbol{p}(t)$ 的对偶圆锥曲线 (Dual Conic Section). 注意到, 求一段圆锥曲线的对偶曲线只需将 $\omega$ 换为 $-\omega$. 如下图所示, 对偶曲线是这一条圆锥曲线上剩余的部分.

有理曲线的性质

连续性

对于 2 阶有理曲线, 在它与它的对偶曲线的连接处只是 $C^1$ 的, 但它在几何上是光滑的.

考虑

$$ \begin{aligned} \boldsymbol{p}(t) &= \frac{(1-t)^2 \boldsymbol{p}_0+ 2t(1-t) \omega \boldsymbol{p}_1+ t^2 \boldsymbol{p}_2}{(1-t)^2+ 2t(1-t) \omega + t^2} \\ \boldsymbol{q}(t) &= \frac{(1-t)^2 \boldsymbol{p}_0- 2t(1-t) \omega \boldsymbol{p}_1+ t^2 \boldsymbol{p}_2}{(1-t)^2- 2t(1-t) \omega + t^2} \\ \end{aligned} $$

则

$$ \begin{aligned} \frac{\mathrm{d}}{\mathrm{d}t} \boldsymbol{p} (0) &= -\frac{\mathrm{d}}{\mathrm{d}t} \boldsymbol{q} (0) = 2 \omega (\boldsymbol{p}_1 - \boldsymbol{p}_0) \\ \frac{\mathrm{d}}{\mathrm{d}t} \boldsymbol{p} (1) &= -\frac{\mathrm{d}}{\mathrm{d}t} \boldsymbol{q} (1) = 2 \omega (\boldsymbol{p}_2 - \boldsymbol{p}_1) \\ \end{aligned} $$

参数化

弧长参数化是研究曲线的重要工具, 但是对于有理曲线, 要进行弧长参数化是很困难的, 我们有如下的定理.

定理 对任何的有理曲线 $\boldsymbol{p}(t) = \frac{\sum_{k=0}^n B_k^{n}(t)\omega_k\boldsymbol{p}_k}{\sum_{k=0}^n B_k^{n}(t)\omega_k}$, $t$ 是弧长参数 (或与弧长参数相差一个倍数) 当且仅当 $\boldsymbol{p}(t)$ 是一条直线.¹

注 尽管得到弧长参数很困难, 但对于任意一条有理曲线, 我们可以将其表示为弧长是时间的有理函数的有理曲线.²

NURBS

NURBS 是非均匀有理 B 样条曲线 (non-uniform rational B-spline) 的缩写, 它实际上就是用 B 样条来表示有理曲线. NURBS 是计算机图形学中最常用的曲线表示方法之一.

$$ \boldsymbol{p}(t) = \frac{\sum_{i=1}^n N_{i,d}(t) \omega_i \boldsymbol{p}_i}{\sum_{i=1}^n N_{i,d}(t) \omega_i} $$