文章目录

  • 免责
  • 前言
  • 拉格朗日插值多项式
    • 第一步:得到一个基函数
    • 第二步:得到所有基函数
    • 第三步:对所有基函数进行线性组合
    • 举例验证
    • 拉格朗日插值曲线绘制实践
  • 三次埃尔米特插值多项式
    • 第一步:得到第一维基函数
    • 第二步:得到第二维基函数
    • 第三步:对基函数进行线性组合
    • 一些后话
  • 规范化的三次埃尔米特多项式
    • 第一步:得到第一维基函数
    • 第二步:得到第二维基函数
    • 第三步:对基函数进行线性组合
    • 一些后话
  • 讲点儿时的回忆
  • 三次贝塞尔曲线
    • 三次贝塞尔曲线的光滑拼接
    • 三次伯恩斯坦基函数
  • 三次均匀 B 样条曲线
    • 局部基函数及其混淆
    • 四阶三次均匀 B 样条曲线
    • B 样条绘制范例
  • 未完待续 ...
  • 一个理念

免责

又来快乐地误人子弟了。 最近学习了一些计算机图形学中的曲线与曲面,感觉很多东西还是蛮有趣的,就打算写篇文章总结总结。
文章中如有疏漏,请联系 premierbob2015@gmail.com 进行更正,笔者感激不尽。

前言

在工业设计中,我们经常会遇到各种各样的绘图问题。随着计算机技术的发展,现代的工业设计正在逐步脱离纸质的工图,转而使用由计算机辅助绘图并设计的电子工图。这种电子版的设计图不仅便于存储与传输,还能更直接地转化为控制数控机床 / 3D打印机 的控制信号,实现工业生产的全自动化。

经过先前的学习,我们了解到在计算机中存储的图像大致可分为点阵图矢量图两类。作为存储图像的不同格式,点阵图和矢量图各有千秋。而上文中所说的电子工图往往使用矢量图的格式进行存储。

优点 缺点
点阵图 易于从相机等设备直接获得,适合近似描述真实世界中的事物 使用像素矩阵存储,不能无限放大,放大会使得图像变得不清晰
矢量图 使用函数或参数方程存储图像,可以无限放大 难以从相机等设备直接获得

本文并不介绍这些图纸在计算机中具体以何种矢量图格式存储,只在数学层面介绍一些非常简单的曲线,而这些曲线在工业设计中有着广泛的应用。(当然,我对工业设计没啥兴趣,我只希望自己将来有一天,能够开发一款用于动漫的绘制辅助工具——属于是老二刺螈了。)

拉格朗日插值多项式

(不重合) 两点确定一条直线,(不共线) 三点确定一条二次函数,不难让我们联想到 nnn 个定点能确定一条的 n−1n-1n−1 次多项式曲线。形式化的来说给定函数上的 nnn 个横坐标互不相同的点:
P1=(x1,y1),P2=(x2,y2),⋯,Pn=(xn,yn)P_1=(x_1, y_1), P_2=(x_2, y_2), \cdots, P_n=(x_n, y_n) P1​=(x1​,y1​),P2​=(x2​,y2​),⋯,Pn​=(xn​,yn​)
我们能够确定一个多项式函数:
y=an−1xn−1+an−2xn−2+⋯+a1x+a0y=a_{n-1}x^{n-1}+a_{n-2}x^{n-2}+\cdots+a_1x+a_0 y=an−1​xn−1+an−2​xn−2+⋯+a1​x+a0​
其中 a0,a1,⋯,an−1a_{0}, a_1, \cdots, a_{n-1}a0​,a1​,⋯,an−1​ 是函数中的参数,xxx 是自变量,yyy 是因变量。所谓确定就是指,将这 nnn 个点坐标带入函数中后,我们能够得到系数向量 a0,a1,⋯,an−1a_0, a_1, \cdots, a_{n-1}a0​,a1​,⋯,an−1​ 的至少一个解。在某些特定条件下,我们能够得到这些系数的一个唯一解,严谨的证明可以参考范德蒙德行列式,在此不再赘述。

对于拉格朗日插值而言,即使我们没有从数学上证明这种可解性按照下面给出的方法,我们也可以使用一种构造性的方式构造出一个恰好通过 P1,P2,⋯,PnP_1, P_2, \cdots, P_nP1​,P2​,⋯,Pn​ 这 nnn 个点的一个 n−1n-1n−1 次函数。

第一步:得到一个基函数

既然我们不能很轻松地让我们的曲线经过所有点,那我们就先试着让我们的曲线经过其中的一个点而在其他点横坐标对应位置的值为零。令 x1,x2,⋯,xnx_1, x_2, \cdots, x_nx1​,x2​,⋯,xn​ 为互不相同的常数,我们能否构造一个 n−1n-1n−1 次多项式函数 f1(x)f_1(x)f1​(x) 使得 :

  1. f1(x1)=1f_1(x_1)=1f1​(x1​)=1;
  2. f1(x2)=f1(x3)=⋯=f1(xn)=0f_1(x_2)=f_1(x_3)=\cdots=f_1(x_n)=0f1​(x2​)=f1​(x3​)=⋯=f1​(xn​)=0 ;

很显然是可以的。首先看到这个函数 f1(x)f_1(x)f1​(x) 有 n−1n-1n−1 个零点,分别是 x2,x3,⋯,xnx_2, x_3, \cdots, x_nx2​,x3​,⋯,xn​,我们很自然地想到设 g1(x)=(x−x2)(x−x3)⋯(x−xn)g_1(x) = (x-x_2)(x-x_3)\cdots(x-x_n)g1​(x)=(x−x2​)(x−x3​)⋯(x−xn​)。虽然这个函数已经满足了 g1(x2)=g1(x3)=⋯=g1(xn)=0g_1(x_2)=g_1(x_3)=\cdots=g_1(x_n)=0g1​(x2​)=g1​(x3​)=⋯=g1​(xn​)=0 这一条件,但是我们并不能保证 g1(x1)=1g_1(x_1)=1g1​(x1​)=1。由于 x1,x2,⋯,xnx_1, x_2, \cdots, x_nx1​,x2​,⋯,xn​ 互不相同,所以 g(x1)=(x1−x2)(x1−x3)⋯(x1−xn)g(x_1)=(x_1-x_2)(x_1-x_3)\cdots(x_1-x_n)g(x1​)=(x1​−x2​)(x1​−x3​)⋯(x1​−xn​) 一定不等于零。不妨令:
f1(x)=g1(x)g1(x1)f_1(x)=\frac{g_1(x)}{g_1(x_1)} f1​(x)=g1​(x1​)g1​(x)​

因为 x1x_1x1​ 是常数,所以 g1(x1)g_1(x_1)g1​(x1​) 就是常数,而这样得到的 f1(x)f_1(x)f1​(x) 是在 g1(x)g_1(x)g1​(x) 基础之上进行了一次常数倍的缩放。此时我们所要求的两个条件就都被满足了。

第二步:得到所有基函数

在第一步中,我们定义了一个名为 fif_ifi​ 的函数。仿照第一步的思路,对于 i=1,2,⋯,ni=1, 2, \cdots, ni=1,2,⋯,n,我们定义:
fi(x)=gi(x)gi(xi)f_i(x)=\frac{g_i(x)}{g_i(x_i)} fi​(x)=gi​(xi​)gi​(x)​

其中:
gi(x)=∏j=1,2,⋯n且j≠i(x−xj)g_i(x)=\prod_{ j=1,2,\cdots n 且 j\neq i }(x-x_j) gi​(x)=j=1,2,⋯n且j​=i∏​(x−xj​)

我们就得到了 nnn 个函数 f1(x),f2(x),⋯,fn(x)f_1(x), f_2(x), \cdots, f_n(x)f1​(x),f2​(x),⋯,fn​(x),根据第一步中的描述,fi(x)f_i(x)fi​(x) 一定满足:

  1. fi(xi)=1f_i(x_i)=1fi​(xi​)=1;
  2. fi(xj)=0,f_i(x_j)=0,fi​(xj​)=0, 对 j=1,2,⋯,nj=1, 2, \cdots, nj=1,2,⋯,n 且 j≠ij\neq ij​=i;

这些函数如此的重要以至于我们要给这些函数起一个好听的名字——拉格朗日基函数。

第三步:对所有基函数进行线性组合

有了拉格朗日基函数后,我们令:

f(x)=y1⋅f1(x)+y2⋅f2(x)⋯+yn⋅fn(x)f(x)=y_1\cdot f_1(x)+y_2\cdot f_2(x)\cdots +y_n\cdot f_n(x) f(x)=y1​⋅f1​(x)+y2​⋅f2​(x)⋯+yn​⋅fn​(x)

我们可以断言,函数 f(x)f(x)f(x) 一定经过点 P1(x1,y1),P2(x2,y2),⋯Pn(xn,yn)P_1(x_1, y_1), P_2(x_2, y_2), \cdots P_n(x_n, y_n)P1​(x1​,y1​),P2​(x2​,y2​),⋯Pn​(xn​,yn​)。例如,当 x=x1x=x_1x=x1​ 时,f1(x)=f1(x1)=1f_1(x)=f_1(x_1)=1f1​(x)=f1​(x1​)=1 而 f2(x)=f3(x)=⋯=fn(x)=0f_2(x)=f_3(x)=\cdots=f_n(x)=0f2​(x)=f3​(x)=⋯=fn​(x)=0。因此:
f(x1)=y1⋅1+y2⋅0+y3⋅0+⋯+yn⋅0=y1f(x_1)=y_1\cdot1+y_2\cdot0+y_3\cdot0+\cdots+y_n\cdot 0=y_1 f(x1​)=y1​⋅1+y2​⋅0+y3​⋅0+⋯+yn​⋅0=y1​

对于 f(x2)=y2,f(x3)=y3,⋯,f(xn)=ynf(x_2)=y_2, f(x_3)=y_3, \cdots, f(x_ n)=y_nf(x2​)=y2​,f(x3​)=y3​,⋯,f(xn​)=yn​ 我们也可以使用同样的方法进行验证。

举例验证

假设一条曲线经过点 (−1,−14),(0,−4),(1,−4),(2,10)(-1, -14), (0, -4), (1, -4), (2, 10)(−1,−14),(0,−4),(1,−4),(2,10)。则可以计算出:

y1f1y_1f_1y1​f1​ y2f2y_2f_2y2​f2​ y3f3y_3f_3y3​f3​ y4f4y_4f_4y4​f4​
yifi(x)y_if_i(x)yi​fi​(x) 73(x−0)(x−1)(x−2)\frac 7 3(x-0)(x-1)(x-2)37​(x−0)(x−1)(x−2) −2(x+1)(x−1)(x−2)-2(x+1)(x-1)(x-2)−2(x+1)(x−1)(x−2) 2(x+1)(x−0)(x−2)2(x+1)(x-0)(x-2)2(x+1)(x−0)(x−2) 53(x+1)(x−0)(x−1)\frac 5 3 (x+1)(x-0)(x-1)35​(x+1)(x−0)(x−1)
yifi(−1)y_if_i(-1)yi​fi​(−1) −14-14−14 000 000 000
yifi(0)y_if_i(0)yi​fi​(0) 000 −4-4−4 000 000
yifi(1)y_if_i(1)yi​fi​(1) 000 000 −4-4−4 000
yifi(2)y_if_i(2)yi​fi​(2) 000 000 000 101010

拉格朗日插值曲线绘制实践

限于篇幅我们把这部分内容放到了另一篇博客中,感兴趣的同学可以来看一看。

GGN_2015
拉格朗日插值曲线绘制实践

三次埃尔米特插值多项式

在工业设计中,有时我们不仅需要指定一条曲线经过哪些点,我们还需要指定这条曲线在某些位置的导数。

考虑一个最简单的情况:设 f(x)=a3⋅x3+a2⋅x2+a1⋅x+a0f(x)=a_3\cdot x^3+a_2\cdot x^2 + a_1 \cdot x + a_0f(x)=a3​⋅x3+a2​⋅x2+a1​⋅x+a0​ 是一个三次多项式。我们要找到一个 f(x)f(x)f(x) 使得它经过点 P0(x0,y0)P_0(x_0, y_0)P0​(x0​,y0​) 和 P1(x1,y1)P_1(x_1, y_1)P1​(x1​,y1​) 且在点 P0P_0P0​ 处的导数为 y0′y_0'y0′​,在点 P1P_1P1​ 处的导数为 y1′y_1'y1′​,其中 x0,x1,y0,y1,y0′,y1′x_0, x_1, y_0, y_1, y_0', y_1'x0​,x1​,y0​,y1​,y0′​,y1′​ 均为常数。

使用类似拉格朗日插值法中的构造方法,我们也能够得到这个 f(x)f(x)f(x)。

第一步:得到第一维基函数

设 x0≠x1x_0\neq x_1x0​​=x1​ 且二者均为常数,令 p0(x)p_0(x)p0​(x) 是一个三次多项式,它满足如下两个性质:

  1. p0(x0)=1p_0(x_0)=1p0​(x0​)=1;
  2. p0(x1)=0,p0′(x1)=0p_0(x_1)=0, p_0'(x_1)=0p0​(x1​)=0,p0′​(x1​)=0;

由于 p0(x1)=p0′(x1)=0p_0(x_1)=p_0'(x_1)=0p0​(x1​)=p0′​(x1​)=0,我们很直观地感受到,x1x_1x1​ 应该是函数 p0(x)p_0(x)p0​(x) 的重根。不妨令 g0(x)=(x−x1)2g_0(x)=(x-x_1)^2g0​(x)=(x−x1​)2,那么我们有 g0′(x)=2(x−x1)g_0'(x)=2(x-x_1)g0′​(x)=2(x−x1​),因此 g0(x1)=g0′(x1)=0g_0(x_1)=g_0'(x_1)=0g0​(x1​)=g0′​(x1​)=0。

按照拉格朗日插值的思路,我们构造了一个函数:
p0(x)=g0(x)g0(x0)p_0(x)=\frac{g_0(x)}{g_0(x_0)} p0​(x)=g0​(x0​)g0​(x)​

因为 x0x_0x0​ 是常数,所以 g0(x0)g_0(x_0)g0​(x0​) 也是常数,我们可以看到:

p0(x0)=1(因为我们构造的系数恰好是g0(x0)的倒数)p0′(x0)=g0′(x0)g0(x0)=2(x0−x1)(x0−x1)2=2x0−x1≠0p0(x1)=p0′(x1)=0\begin{aligned} p_0(x_0) &=1\;(因为我们构造的系数恰好是g_0(x_0)的倒数)\\ p_0'(x_0) &=\frac{g_0'(x_0)}{g_0(x_0)}=\frac{2(x_0-x_1)}{(x_0-x_1)^2}=\frac{2}{x_0-x_1}\neq 0\\ p_0(x_1) &=p_0'(x_1)=0 \end{aligned} p0​(x0​)p0′​(x0​)p0​(x1​)​=1(因为我们构造的系数恰好是g0​(x0​)的倒数)=g0​(x0​)g0′​(x0​)​=(x0​−x1​)22(x0​−x1​)​=x0​−x1​2​​=0=p0′​(x1​)=0​

有了函数 p0p_0p0​,我们就可以在不改变 p0(x1)p_0(x_1)p0​(x1​) 和 p0′(x1)p_0'(x_1)p0′​(x1​) 的值的前提下随意控制 p0(x0)p_0(x_0)p0​(x0​) 的值了。但是现在我们还无法控制 p0′(x0)p_0'(x_0)p0′​(x0​) 的值。因此我们要在第二步构造一个函数,使得我们能够分离 p0′(x0)p_0'(x_0)p0′​(x0​) 和 p0(x0)p_0(x_0)p0​(x0​) 之间的倍数关系。

第二步:得到第二维基函数

类似第一步,我们还需要构造三次函数 q0(x)q_0(x)q0​(x) 满足这样两条性质:

  1. q0(x0)=0,q0′(x0)=1q_0(x_0)=0, q_0'(x_0)=1q0​(x0​)=0,q0′​(x0​)=1;
  2. q0(x1)=0,q0′(x1)=0q_0(x_1)=0, q_0'(x_1)=0q0​(x1​)=0,q0′​(x1​)=0;

看起来也不难构造,零点又多了一个 x0x_0x0​。令 k0(x)=(x−x0)(x−x1)2k_0(x)=(x-x_0)(x-x_1)^2k0​(x)=(x−x0​)(x−x1​)2,则:
k0′(x)=(x−x1)2+2(x−x0)(x−x1)=(x−x1)(3x−2x0−x1)k_0'(x)=(x-x_1)^2+2(x-x_0)(x-x_1)=(x-x_1)(3x-2x_0-x_1) k0′​(x)=(x−x1​)2+2(x−x0​)(x−x1​)=(x−x1​)(3x−2x0​−x1​)

令:
q0(x)=k0(x)k0′(x0)q_0(x)=\frac{k_0(x)}{k_0'(x_0)}q0​(x)=k0′​(x0​)k0​(x)​

由于 x0x_0x0​ 是常数因此 k0′(x0)=(x0−x1)2k_0'(x_0)=(x_0-x_1)^2k0′​(x0​)=(x0​−x1​)2 也是常数,我们可以看到:
q0(x0)=0(因为k0(x0)=0)q0′(x0)=1(因为我们构造的系数恰好是k0′(x0)的倒数)q0(x1)=q0′(x1)=0\begin{aligned} q_0(x_0)&=0\;(因为k_0(x_0)=0)\\ q_0'(x_0)&=1\;(因为我们构造的系数恰好是k_0'(x_0)的倒数)\\ q_0(x_1)&=q_0'(x_1)=0 \end{aligned} q0​(x0​)q0′​(x0​)q0​(x1​)​=0(因为k0​(x0​)=0)=1(因为我们构造的系数恰好是k0′​(x0​)的倒数)=q0′​(x1​)=0​

第三步:对基函数进行线性组合

开门见山,令:
f0(x)=y0⋅(p0(x)−2x0−x1⋅q0(x))+y0′⋅q0(x)f_0(x)=y_0\cdot (p_0(x)-\frac{2}{x_0-x_1}\cdot q_0(x))+y_0'\cdot q_0(x) f0​(x)=y0​⋅(p0​(x)−x0​−x1​2​⋅q0​(x))+y0′​⋅q0​(x)

为什么要这么构造呢?p0(x)p_0(x)p0​(x) 在 x0x_0x0​ 处的导数值为 2x0−x1≠0\frac{2}{x_0-x_1}\neq 0x0​−x1​2​​=0,因此在它身上减去一个 2x0−x1⋅q0(x)\frac{2}{x_0-x_1}\cdot q_0(x)x0​−x1​2​⋅q0​(x) 就能使得它在 x0x_0x0​ 处的导数为零,而由于 q0(x)q_0(x)q0​(x) 在 x1x_1x1​ 点的函数值和导数值都为零,因此无论我们如何将 p0(x)p_0(x)p0​(x) 和 q0(x)q_0(x)q0​(x) 进行线性组合,都不会影响函数在 x1x_1x1​ 点处的函数值和导数值。

此时函数 f0(x)f_0(x)f0​(x) 满足:

  1. f0(x0)=y0,f0′(x0)=y0′f_0(x_0)=y_0,f_0'(x_0)=y_0'f0​(x0​)=y0​,f0′​(x0​)=y0′​;
  2. f0(x1)=0,f0′(x1)=0f_0(x_1)=0, f_0'(x_1)=0f0​(x1​)=0,f0′​(x1​)=0;

类似地,我们可以使用同样的方法得到一个函数 f1(x)f_1(x)f1​(x),满足:

  1. f1(x0)=0,f1′(x0)=0f_1(x_0)=0,f_1'(x_0)=0f1​(x0​)=0,f1′​(x0​)=0;
  2. f1(x1)=y1,f1′(x1)=y1′f_1(x_1)=y_1, f_1'(x_1)=y_1'f1​(x1​)=y1​,f1′​(x1​)=y1′​;

具体怎么做不再赘述,整理后得到的表达式我也不是很关心。我所关心的只有一件事那就是,当我们令 f(x)=f0(x)+f1(x)f(x)=f_0(x)+f_1(x)f(x)=f0​(x)+f1​(x),这样得到的 f(x)f(x)f(x) 一定经过 P0(x0,y0)P_0(x_0, y_0)P0​(x0​,y0​) 和 P1(x1,y1)P_1(x_1, y_1)P1​(x1​,y1​) 且在点 P0P_0P0​ 处的导数为 y0′y_0'y0′​,在点 P1P_1P1​ 处的导数为 y1′y_1'y1′​。

一些后话

如果我们想要找到一个函数使得他定义在 [x1,xn][x_1, x_n][x1​,xn​] 上且函数光滑(可导),且经过点 P1(x1,y1),P2(x2,y2),⋯,Pn(xn,yn)P_1(x_1, y_1), P_2(x_2, y_2), \cdots, P_n(x_n, y_n)P1​(x1​,y1​),P2​(x2​,y2​),⋯,Pn​(xn​,yn​),且在 PiP_iPi​ 处导数值为 yi′y_i'yi′​(i=1,2⋯,ni=1, 2\cdots, ni=1,2⋯,n),其中 xi,yi,yi′x_i, y_i, y_i'xi​,yi​,yi′​ 均为常数(i=1,2,⋯,ni=1, 2, \cdots, ni=1,2,⋯,n),且满足 x1<x2<⋯<xnx_1<x_2<\cdots<x_nx1​<x2​<⋯<xn​ 。

我们可以在区间 [x1,x2][x_1, x_2][x1​,x2​] 上依据 y1,y2,y1′,y2′y_1, y_2, y_1', y_2'y1​,y2​,y1′​,y2′​ 构造一个三次埃尔米特插值函数,在 [x2,x3][x_2, x_3][x2​,x3​] 上依据 y2,y3,y2′,y3y_2, y_3, y_2', y_3y2​,y3​,y2′​,y3​ 构造插值函数,⋯\cdots⋯,在 [xn−1,xn][x_{n-1}, x_{n}][xn−1​,xn​] 上根据 yn−1,yn,yn−1′,yn′y_{n-1}, y_{n}, y_{n-1}', y_n'yn−1​,yn​,yn−1′​,yn′​ 构造插值函数。最后再将每一段上的函数以分段函数的形式连接起来,由于连接点处的左右函数值、左右导数值均相等,因此函数在区间 (x1,xn)(x_1, x_n)(x1​,xn​) 上光滑。

规范化的三次埃尔米特多项式

如果我们的目的只是绘制一个函数,使用普通的埃尔米特多项式似乎已经实现了我们想要的。但是函数这东西,自变量 xxx 的每一个取值至多只有一个 yyy 值与之对应,如果我想绘制任意平面曲线甚至空间曲线,那上文给出的方法就不那么好用了。想到可以用参数方程的方法来实现。

  • 下文中我们令P0,P1,D0,D1P_0, P_1, D_0, D_1P0​,P1​,D0​,D1​ 是二维常向量,我们要构造参数方程 f⃗(t)\vec{f}(t)f​(t) 使得 f⃗(0)=P0,f⃗(1)=P1,∇f⃗(0)=D0,∇f⃗(1)=D1\vec f (0)=P_0, \vec f(1) = P_1, \nabla \vec f(0) = D_0, \nabla \vec f(1)=D_1f​(0)=P0​,f​(1)=P1​,∇f​(0)=D0​,∇f​(1)=D1​。

注:梯度算子 ∇\nabla∇,其实就是对函数值的 xxx,yyy 分量函数分别求导,例如当 f⃗(t)=[fx(t),fy(t)]\vec f(t) = [f_x(t), f_y(t)]f​(t)=[fx​(t),fy​(t)] 时 ∇f⃗(t)=[fx′(t),fy′(t)]\nabla \vec f(t)=[f_x'(t), f_y'(t)]∇f​(t)=[fx′​(t),fy′​(t)]。

对于一个标量函数 f(t)f(t)f(t) 而言,假如 P∈R2P\in \R^2P∈R2 是一个二维常向量(即 Px,Py∈RP_x, P_y\in \RPx​,Py​∈R 是两个常数),那么 Pf(t)=[Pxf(t),Pyf(t)]Pf(t)=[P_xf(t), P_yf(t)]Pf(t)=[Px​f(t),Py​f(t)] 就是一个平面上的参数方程。对这个函数计算梯度,得到 ∇(Pf(t))=[Pxf′(t),Pyf′(t)]=P⋅f′(t)\nabla(Pf(t))=[P_xf'(t), P_yf'(t)]=P\cdot f'(t)∇(Pf(t))=[Px​f′(t),Py​f′(t)]=P⋅f′(t)。此处的点运算 ‘⋅\cdot⋅’ 不表示向量的内积,而表示“数乘”,其中 P∈R2P\in \R^2P∈R2 是 “向量”,f(t)∈Rf(t)\in \Rf(t)∈R 是定义在实数上的标量函数。

类似上文的方法,我们仍然可以构造出符合条件的向量函数 f⃗(t)\vec f(t)f​(t)。

第一步:得到第一维基函数

在这一步中我们要构造标量函数 p0(t)p_0(t)p0​(t) 和 p1(t)p_1(t)p1​(t),其中 p0(t)p_0(t)p0​(t) 满足:

  1. p0(0)=1p_0(0)=1p0​(0)=1;
  2. p0(1)=0,p0′(1)=0p_0(1)=0, p_0'(1)=0p0​(1)=0,p0′​(1)=0;

同理,p1(t)p_1(t)p1​(t) 满足:

  1. p1(1)=1p_1(1)=1p1​(1)=1;
  2. p1(0)=0,p1′(0)=0p_1(0)=0, p_1'(0)=0p1​(0)=0,p1′​(0)=0;

为了方便表述,我们把 p0(t)p_0(t)p0​(t) 称为 t=0t=0t=0 处的常数控制基函数,p1(t)p_1(t)p1​(t) 称为 t=1t=1t=1 处的常数控制基函数,它们的存在就是为了用于将它们线性组合以调整函数在 t=0t=0t=0 处以及 t=1t=1t=1 处的函数值。

由于这两个函数都比较简单,稍微构造构造就能构造出来(拿眼睛看就能看出来,不信你试试):
p0(t)=(t−1)2p1(t)=t2\begin{aligned} p_0(t)&=(t-1)^2\\ p_1(t)&=t^2 \end{aligned} p0​(t)p1​(t)​=(t−1)2=t2​

为了方便下一步中消去 p0(t)p_0(t)p0​(t) 和 p1(t)p_1(t)p1​(t) 对导数的影响,我们可以提前计算出 p0′(0)p_0'(0)p0′​(0) 和 p1′(1)p_1'(1)p1′​(1):

p0′(t)=2(t−1)∴p0′(0)=−2p1′(t)=2t∴p1′(1)=2\begin{aligned} p_0'(t)&=2(t-1)& \therefore p_0'(0)&=-2\\ p_1'(t)&=2t& \therefore p_1'(1)&=2 \end{aligned} p0′​(t)p1′​(t)​=2(t−1)=2t​∴p0′​(0)∴p1′​(1)​=−2=2​

第二步:得到第二维基函数

类似上文中对三次埃尔米特插值多项式的构造,我们要构造 q0(t)q_0(t)q0​(t) 和 q1(t)q_1(t)q1​(t),其中 q0(t)q_0(t)q0​(t) 满足:

  1. q0(0)=0,q0′(0)=1q_0(0)=0, q_0'(0)=1q0​(0)=0,q0′​(0)=1;
  2. q0(1)=0,q0′(1)=0q_0(1)=0, q_0'(1)=0q0​(1)=0,q0′​(1)=0;

类似地,q1(t)q_1(t)q1​(t) 满足:

  1. q1(0)=0,q1′(0)=0q_1(0) = 0, q_1'(0) = 0q1​(0)=0,q1′​(0)=0;
  2. q1(1)=0,q1′(1)=1q_1(1) = 0, q_1'(1)=1q1​(1)=0,q1′​(1)=1;

想解方程也不是不行,只不过确实可以直接看出结果:
q0(t)=t(t−1)2q1(t)=t2(t−1)\begin{aligned} q_0(t)&=t(t-1)^2\\ q_1(t)&=t^2(t-1)\\ \end{aligned} q0​(t)q1​(t)​=t(t−1)2=t2(t−1)​

同样为了方便表述,我们把 q0(t)q_0(t)q0​(t) 称为 t=0t=0t=0 处的导数控制基函数,q1(t)q_1(t)q1​(t) 称为 t=1t=1t=1 处的导数控制基函数,它们有两重作用:一重是和 p0p_0p0​ 或 p1p_1p1​ 线性组合,从而消去常数控制基函数对导函数的影响,另外一重作用是设置导函数的值。接下来我们分别从这两个角度来阐述导数控制基函数的作用。

考虑,第一方面的作用,记:
g0(t)=p0(t)−p0′(0)⋅q0(t)=(t−1)2+2t(t−1)2=(t−1)2(2t+1)g1(t)=p1(t)−p1′(1)⋅q1(t)=t2−2t2(t−1)=t2(3−2t)\begin{aligned} g_0(t)&=p_0(t)-p_0'(0)\cdot q_0(t)=&(t-1)^2+2t(t-1)^2&=(t-1)^2(2t+1)\\ g_1(t)&=p_1(t)-p_1'(1)\cdot q_1(t)=&t^2-2t^2(t-1)&=t^2(3-2t) \end{aligned} g0​(t)g1​(t)​=p0​(t)−p0′​(0)⋅q0​(t)==p1​(t)−p1′​(1)⋅q1​(t)=​(t−1)2+2t(t−1)2t2−2t2(t−1)​=(t−1)2(2t+1)=t2(3−2t)​

这样我们得到的 g0(t)g_0(t)g0​(t) 和 g1(t)g_1(t)g1​(t) 就能够实现导数无关的常数控制,因为我们用 q0q_0q0​ 抵消了 p0p_0p0​ 在 x=0x=0x=0 处对导数的影响,用 q1q_1q1​ 抵消了 p1p_1p1​ 在 x=1x=1x=1 处的导数影响。为了和 p0p_0p0​ 与 p1p_1p1​ 加以区分,我们称 g0g_0g0​ 和 g1g_1g1​ 为导数无关的控制基函数

顺便说一嘴,g0,g1,q0,q1g_0, g_1, q_0, q_1g0​,g1​,q0​,q1​ 这四个函数可以用待定系数法列方程求解,但是列方程多无聊啊,构造出来多有趣。

第三步:对基函数进行线性组合

f⃗(t)=P0⋅g0(t)+P1⋅g1(t)+D0⋅q0(t)+D1⋅q1(t),t∈[0,1]\vec f(t)=P_0 \cdot g_0(t) +P_1 \cdot g_1(t)+D_0\cdot q_0(t)+D_1\cdot q_1(t), t\in[0, 1] f​(t)=P0​⋅g0​(t)+P1​⋅g1​(t)+D0​⋅q0​(t)+D1​⋅q1​(t),t∈[0,1]

不难证明这个东西就是我们想要的,把含 ttt 的解析式带入得到:
f⃗(t)=P0(2t3−3t2+1)+P1(−2t3+3t2)+D0(t3−2t2+t)+D1(t3−t2),t∈[0,1]\vec f(t)=P_0(2t^3-3t^2+1)+P_1(-2t^3+3t^2)+D_0(t^3-2t^2+t)+D_1(t^3-t^2), t\in[0, 1] f​(t)=P0​(2t3−3t2+1)+P1​(−2t3+3t2)+D0​(t3−2t2+t)+D1​(t3−t2),t∈[0,1]

一些后话

上文中讨论的规范化三次埃尔米特插值曲线可以用类似一般的三次埃尔米特插值曲线一样逐段连接起来,形成一条经过多个控制点的光滑曲线段。例如,我们可能找到这样的一个 f⃗(t),t∈[t1,tn]\vec f(t),t\in[t_1, t_n]f​(t),t∈[t1​,tn​] 使得 t=t1t=t_1t=t1​ 时 f⃗(t)=P1\vec f(t) = P_1f​(t)=P1​ 且 ∇f⃗(t)=D1\nabla \vec f(t) = D_1∇f​(t)=D1​,t=t2t=t_2t=t2​ 时 f⃗(t)=P2\vec f(t) = P_2f​(t)=P2​ 且 ∇f⃗(t)=D2\nabla \vec f(t) = D_2∇f​(t)=D2​,⋯\cdots⋯,t=tnt=t_nt=tn​ 时 f⃗(t)=Pn\vec f(t) = P_nf​(t)=Pn​ 且 ∇f⃗(t)=Dn\nabla \vec f(t) = D_n∇f​(t)=Dn​。分别对 [t1,t2][t_1, t_2][t1​,t2​],[t2,t3][t_2, t_3][t2​,t3​],⋯\cdots⋯,[tn−1,tn][t_{n-1}, t_n][tn−1​,tn​] 建立一个规范化三次埃尔米特插值曲线后再用分段函数连接起来即可。

由于规范化三次埃尔米特插值曲线的自变量 ttt 的取值范围总是在 [0,1][0,1][0,1] 中,因此需要为每一段的参数进行换元。令 ui(t)=t−titi+1−tiu_i(t)=\frac{t-t_{i}}{t_{i+1}-t_i}ui​(t)=ti+1​−ti​t−ti​​ 其中 i=1,2,⋯n−1i=1, 2, \cdots n-1i=1,2,⋯n−1,则有当 t∈[ti,ti+1]t\in[t_i, t_{i+1}]t∈[ti​,ti+1​] 时,ui(t)u_i(t)ui​(t) 从 000 线性地变化为 111,因此我们想要的分段函数即为:

f⃗(t)={P1g0(u1(t))+P2g1(u1(t))+D1q0(u1(t))+D2q1(u1(t)),t∈[t1,t2)P2g0(u2(t))+P3g1(u2(t))+D2q0(u2(t))+D3q1(u2(t)),t∈[t2,t3)⋯Pn−1g0(un−1(t))+Png1(un−1(t))+Dn−1q0(un−1(t))+Dnq1(un−1(t)),t∈[tn−1,tn)\vec f(t) = \left\{\begin{aligned} P_1g_0(u_1(t))+P_2g_1(u_1(t))+D_1q_0(u_1(t))+D_2q_1(u_1(t)),&t\in[t_1, t_2)\\ P_2g_0(u_2(t))+P_3g_1(u_2(t))+D_2q_0(u_2(t))+D_3q_1(u_2(t)),&t\in[t_2, t_3)\\ \cdots\\ P_{n-1}g_0(u_{n-1}(t))+P_ng_1(u_{n-1}(t))+D_{n-1}q_0(u_{n-1}(t))+D_nq_1(u_{n-1}(t)),&t\in[t_{n-1}, t_n)\\ \end{aligned}\right. f​(t)=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​P1​g0​(u1​(t))+P2​g1​(u1​(t))+D1​q0​(u1​(t))+D2​q1​(u1​(t)),P2​g0​(u2​(t))+P3​g1​(u2​(t))+D2​q0​(u2​(t))+D3​q1​(u2​(t)),⋯Pn−1​g0​(un−1​(t))+Pn​g1​(un−1​(t))+Dn−1​q0​(un−1​(t))+Dn​q1​(un−1​(t)),​t∈[t1​,t2​)t∈[t2​,t3​)t∈[tn−1​,tn​)​

在自变量 ttt 的变化过程中,这个分段函数每一段中的 uiu_iui​ 都能从 000 开始到 111 结束——“规范化”想要表达的就是这重含义。表达式中的 g0,g1g_0, g_1g0​,g1​ 是导数无关的常数控制基函数,q0,q1q_0, q_1q0​,q1​ 是导数控制基函数

具体而言,无论是描述空间中的曲线还是平面中的曲线,我们都可以使用这样的一种方式:先确定一系列的控制点(例如上文中的 P0,P1P_0, P_1P0​,P1​),之后再将这些控制点与一组标量函数对应相乘,这组标量函数一般被称为基函数

讲点儿时的回忆

贝塞尔曲线让我回忆起了我小时候经常玩的一个画图游戏:

  1. 在纸上画出两条登场的线段首位顺次相接,如下图中的 AB=BCAB=BCAB=BC 且使 ABABAB 垂直于 BCBCBC;
  2. 在 ABABAB 边上取一点 DDD,在 BCBCBC 边上取一点 EEE,使得 AD=BEAD=BEAD=BE;
  3. 在纸上借助直尺连接 DEDEDE;
  4. 随机更换 DDD 点的选取和 EEE 点的选取,但始终保持 AD=BEAD=BEAD=BE 不变,转步骤 3;

随着连接得到的线段越来越多,这些线段下方空白部分逐渐形成了一个曲线弧的形状(我上小学的时候画过不少这种东西,尽管那时的我不知道这个曲线是一条抛物线的一部分,一直以为它是一个圆弧)。

如何简单地说明这个曲线是一条抛物线呢,比较直观的方法就是先把焦点和准线找出来。对于这种绘制方式而言,想要证明绘制出的曲线是抛物线并不容易,但是我们可以反过来,去证明一个抛物线可以用这种方法进行绘制(虽然数学上这种证法并不严谨,但是我们做这个证明还有其他的用处)。


上图中红色的抛物线为 y=14x2y=\frac 1 4 x^2y=41​x2,MMM 是抛物线ACACAC 段上的任意一点,过 MMM 做抛物线的切线 DEDEDE,交 ABABAB、BCBCBC 于 D,ED,ED,E 两点,此时我们只需要证明 AD=BEAD=BEAD=BE 即可说明上文中的绘制方式能够绘制一条抛物线。

由于 y′=12xy'=\frac 1 2 xy′=21​x,设 MMM 点横坐标为 x0x_0x0​,得到 MMM 点出的切线方程为 y=12x0(x−x0)+14x02y=\frac 1 2 x_0(x-x_0)+\frac 1 4 x_0^2y=21​x0​(x−x0​)+41​x02​,将切线方程分别与 ABABAB 方程、BCBCBC 方程联立可以解得 D,ED, ED,E 两点的坐标。最后我们可以验证 xD−xA=xE−xBx_D-x_A=x_E-x_BxD​−xA​=xE​−xB​。其中 xD=12x0−1,xE=1+12x0x_D=\frac 1 2 x_0-1, x_E=1+\frac 1 2 x_0xD​=21​x0​−1,xE​=1+21​x0​ 带入成立。

对此我们还可以验证得到 AD:DB=DM:MEAD:DB=DM:MEAD:DB=DM:ME,这也恰好为我们提供了一种二次贝塞尔曲线的绘制方式(几何作图法)。假设 A,B,CA, B, CA,B,C 是我们在平面上选取的三个点,我们根据这种比例绘制的方式可以得到 MMM 移动得到的的参数方程曲线。 设 t∈[0,1]t\in[0, 1]t∈[0,1],当 t=0t=0t=0 时 MMM 与 AAA 点重合,当 t=1t=1t=1 时 MMM 与 CCC 点重合。

首先我们能很容易的得到 D=(1−t)A+tBD=(1-t)A+tBD=(1−t)A+tB,换言之,t=0t=0t=0 时 DDD 与 AAA 重合,t=1t=1t=1 时 DDD 与 BBB 重合,t∈(0,1)t\in(0, 1)t∈(0,1) 时 DDD 在 A→BA\to BA→B 上做匀速运动。同理 E=(1−t)B+tCE=(1-t)B+tCE=(1−t)B+tC,而根据我们上文指出的比例关系,我们知道 DM:ME=AD:DB=t:(1−t)DM:ME=AD:DB=t:(1-t)DM:ME=AD:DB=t:(1−t),因此得到 M(t)=(1−t)D+tE=(1−t)2A+2t(1−t)B+t2CM(t)=(1-t)D+tE=(1-t)^2A+2t(1-t)B+t^2CM(t)=(1−t)D+tE=(1−t)2A+2t(1−t)B+t2C。当 ttt 取遍 [0,1][0, 1][0,1] 中的所有值时,M(t)M(t)M(t) 取遍抛物线在 AAA 到 CCC 之间的所有点。其中,我们称 (1−t)2(1-t)^2(1−t)2,2t(1−t)2t(1-t)2t(1−t),t2t^2t2 这三个函数为“二次伯恩斯坦基函数”。

三次贝塞尔曲线

不难发现其实“二次伯恩斯坦基函数”就是对 ((1−t)+t)2((1-t)+t)^2((1−t)+t)2 进行牛顿二项式展开之后的每一项。类比地,我们可以得到四个三次伯恩斯坦基函数
(1−t)3,3(1−t)2t,3(1−t)t2,t3(1-t)^3,3(1-t)^2t,3(1-t)t^2,t^3 (1−t)3,3(1−t)2t,3(1−t)t2,t3

当我们有四个平面上的常数点 A,B,C,DA, B, C, DA,B,C,D 时,我们可以确定一条三次贝塞尔函数:
M(t)=(1−t)3A+3(1−t)2tB+3(1−t)t2C+t3D,t∈[0,1]M(t)=(1-t)^3A+3(1-t)^2tB+3(1-t)t^2C+t^3D,t\in[0, 1] M(t)=(1−t)3A+3(1−t)2tB+3(1−t)t2C+t3D,t∈[0,1]

虽然有高次的贝塞尔函数,但是三次的贝塞尔曲线已经可以实现光滑光顺的拼接,所以我觉得不介绍高次贝塞尔曲线也罢。总而言之“nnn 次的伯恩斯坦基函数”就是对 ((1−t)+t)n((1-t)+t)^n((1−t)+t)n 进行牛顿二项式展开之后的每一项,其中的第 iii 项为 (i=0,1,⋯,ni=0, 1, \cdots, ni=0,1,⋯,n):
Bi(t)=Cni⋅ti(1−t)n−iB_i(t)=C_n^i\cdot t^i(1-t)^{n-i} Bi​(t)=Cni​⋅ti(1−t)n−i

我们将 n+1n+1n+1 个 nnn 次伯恩斯坦基函数与 nnn 个控制点坐标对应相乘就得到了一条 nnn 次的贝塞尔曲线。美中不足的一点是,由于这 n+1n+1n+1 个 nnn 次伯恩斯坦基函数的定义域都是全体实数,因此每当我们改变一个控制点的位置,整条曲线上的每一个点位置都会发生一定的改变,这不利于工业设计上的局部微调。

想象我们正在使用贝塞尔曲线绘制一位二次元少女的右脸,当我们觉得不满意时打算对这一段曲线进行微调,于是调整了右脸处的某个控制点,而这导致原本已经绘制得很不错的左脸发生了变形,这是一件多么悲伤的事。

三次贝塞尔曲线的光滑拼接

假设我们现在有两条贝塞尔曲线:
P(t)=(1−t)3P0+3(1−t)2tP1+3(1−t)t2P2+t3P3,t∈[0,1]Q(t)=(1−t)3Q0+3(1−t)2tQ1+3(1−t)t2Q2+t3Q3,t∈[0,1]\begin{aligned} P(t)&=(1-t)^3P_0+3(1-t)^2tP_1+3(1-t)t^2P_2+t^3P_3&,t\in[0, 1]\\ Q(t)&=(1-t)^3Q_0+3(1-t)^2tQ_1+3(1-t)t^2Q_2+t^3Q_3&,t\in[0, 1] \end{aligned} P(t)Q(t)​=(1−t)3P0​+3(1−t)2tP1​+3(1−t)t2P2​+t3P3​=(1−t)3Q0​+3(1−t)2tQ1​+3(1−t)t2Q2​+t3Q3​​,t∈[0,1],t∈[0,1]​

其中 P0,P1,P2,P3,Q0,Q1,Q2,Q3P_0, P_1, P_2, P_3, Q_0, Q_1, Q_2, Q_3P0​,P1​,P2​,P3​,Q0​,Q1​,Q2​,Q3​ 是平面上的常数点。我们希望把曲线 QQQ 光滑地接在曲线 PPP 地后面,只需要保证:

  1. 连续: P(1)=Q(0)P(1)=Q(0)P(1)=Q(0);
  2. 梯度同向平行:∇P(1)⋅∇Q(0)=∣∇P(1)∣⋅∣∇Q(0)∣\nabla P(1) \cdot \nabla Q(0)= |\nabla P(1)|\cdot|\nabla Q(0)|∇P(1)⋅∇Q(0)=∣∇P(1)∣⋅∣∇Q(0)∣ (其中第一个点运算是向量内积,第二个是标量乘法);

由于 P(1)=P3P(1)=P_3P(1)=P3​,Q(0)=Q0Q(0)=Q_0Q(0)=Q0​,所以想要保证连续只需要保证 P3=Q0P_3=Q_0P3​=Q0​。
∇P(t)=∇((1−t)((1−t)2P0+3(1−t)tP1+3t2P2)+t3P3)=−((1−t)2P0+3(1−t)tP1+3t2P2)+(1−t)(⋯)+3t2P3\begin{aligned} \nabla P(t)&=\nabla ((1-t)((1-t)^2P_0+3(1-t)tP_1+3t^2P_2)+t^3P_3)\\ &=-((1-t)^2P_0+3(1-t)tP_1+3t^2P_2)+(1-t)(\cdots)+3t^2P_3 \end{aligned} ∇P(t)​=∇((1−t)((1−t)2P0​+3(1−t)tP1​+3t2P2​)+t3P3​)=−((1−t)2P0​+3(1−t)tP1​+3t2P2​)+(1−t)(⋯)+3t2P3​​

由于我们想计算的是 ∇P(1)\nabla P(1)∇P(1) 所以 ⋯\cdots⋯ 的部分我不想算了(反正乘上一个 1−t=01-t=01−t=0 后也一定等于零),得到 ∇P(1)=3(P3−P2)\nabla P(1)=3(P_3-P_2)∇P(1)=3(P3​−P2​)。同理可以计算得到:
∇Q(t)=∇((1−t)3Q0+t(3(1−t)2Q1+3(1−t)tQ2+t2Q3))=−3(1−t)2Q0+(3(1−t)2Q1+3(1−t)tQ2+t2Q3)+t(⋯)\begin{aligned} \nabla Q(t)&=\nabla ((1-t)^3Q_0 + t(3(1-t)^2Q_1+3(1-t)tQ_2+t^2Q_3))\\ &=-3(1-t)^2Q_0+(3(1-t)^2Q_1+3(1-t)tQ_2+t^2Q_3)+t(\cdots) \end{aligned} ∇Q(t)​=∇((1−t)3Q0​+t(3(1−t)2Q1​+3(1−t)tQ2​+t2Q3​))=−3(1−t)2Q0​+(3(1−t)2Q1​+3(1−t)tQ2​+t2Q3​)+t(⋯)​

得到:∇Q(0)=3(Q1−Q0)\nabla Q(0)=3(Q_1-Q_0)∇Q(0)=3(Q1​−Q0​)。因此只要在 P3=Q0P_3=Q_0P3​=Q0​ 的基础上保证 P2P3⟶//Q0Q1⟶\stackrel \longrightarrow {P_2P_3}//\stackrel \longrightarrow{Q_0Q_1}P2​P3​⟶​//Q0​Q1​⟶​ 即可保证(还要注意同向哦)。

三次伯恩斯坦基函数

我们看到三次的贝塞尔曲线有四个控制点:
M(t)=(1−t)3A+3(1−t)2tB+3(1−t)t2C+t3D,t∈[0,1]M(t)=(1-t)^3A+3(1-t)^2tB+3(1-t)t^2C+t^3D,t\in[0, 1] M(t)=(1−t)3A+3(1−t)2tB+3(1−t)t2C+t3D,t∈[0,1]

每当给定一个 ttt 时,曲线上的点 M(t)M(t)M(t) 实际上就是对 A,B,C,DA, B, C, DA,B,C,D 四个控制点的坐标进行了一个加权求平均数的过程。当 t→0t\to 0t→0 时,M(t)→AM(t)\to AM(t)→A,当 t→1t\to 1t→1 时 M(t)→DM(t)\to DM(t)→D。具体而言我们可以将 A,B,C,DA, B, C, DA,B,C,D 四个点分别对应的基函数绘制到一个坐标系中,让我们观察这四个点对整条曲线究竟起到了多大的作用。

可以看到,当 t=13t=\frac 1 3t=31​ 时,BBB 前的系数变得最大,此时曲线最靠近 BBB 点;当 t=23t=\frac 2 3t=32​ 时,CCC 前的系数变得最大,此时曲线最靠近 CCC 点。

换言之,这也为我们进一步研究曲线的表示提供了一个思路:给定 nnn 个控制点,我们只需要给出 nnn 个标量函数作为基函数和这 nnn 个点的坐标对应相乘,就能得到一条大趋势与这 nnn 个点的走势基本接近的曲线。这 nnn 个标量函数都是上凸单峰的,使得曲线在当前标量函数的极值点处最接近该控制点,从而使得曲线近似拟合了这 nnn 个控制点构成的折线。

由于 CSDN 对字数的限制,对贝塞尔曲线的绘制将在另一篇博客中进行讲解。

GGN_2015
计算机图形学中的曲线问题——贝塞尔曲线的绘制

三次均匀 B 样条曲线

由于三次贝塞尔曲线只能有四个控制点,贝塞尔曲线控制点越多次数越高。我们希望构造一种这样的曲线:

  1. 控制点个数可以无限增多,但函数次数保持为三次;
  2. 每个控制点只会影响自己附近的一小段曲线而对其他部分曲线没有任何影响;
  3. 曲线是光滑的。

局部基函数及其混淆

我们的目的非常明确,我们想要构造一种基函数,这种基函数只在很小的一段区间上有定义,而在此区间之外的部分均为零。例如这样的一个函数就是我们想要的一个函数:
N1(t)={1,0≤t≤10,elsewhereN_1(t)=\left\{\begin{aligned}1, &0\leq t \leq 1\\0,&\text{elsewhere} \end{aligned}\right. N1​(t)={1,0,​0≤t≤1elsewhere​

这个函数又被称作“平台函数”,因为它在非零区间 [0,1][0, 1][0,1] 上的取值始终是一。由于这个函数的非零区间为 [0,1][0, 1][0,1],因此我们称这个局部基函数的宽度为 111。类似的,如果一个定义在实数上的函数只在区间 [0,L][0, L][0,L] 上有非零值,我们可以称这个函数的“宽度”为 LLL 的局部基函数

当我们用一组局部基函数与一组控制点对应相乘再求和从而得到目标曲线时,每一个控制点自然也只会影响曲线的一个局部(在局部基函数的宽度范围内会产生影响),而不会影响整个曲线上的绝大多数位置的取值。

假如我们已经有了一个宽度为 LLL 的局部基函数,我们可以使用一次混淆操作得到一个宽度为 L+1L+1L+1 的局部基函数,这使得我们可以控制每个控制点在目标曲线上到底能够造成多大范围的影响,宽度越大影响范围自然也越大。

对于平台函数 N1N_1N1​ 进行混淆,第一步我们要将 N1(t)N_1(t)N1​(t) 右移一个单位长度得到 N1(t−1)N_1(t-1)N1​(t−1),我们希望能够将 N1(t)N_1(t)N1​(t) 与 N1(t−1)N_1(t-1)N1​(t−1) 加权求平均,使得我们得到的新函数在 t=0→1t=0\to1t=0→1 段上函数值线性增加,在 t=1→2t=1\to2t=1→2 段上函数值线性减少。很直观的想法就是,我们可以在 N1(t)N_1(t)N1​(t) 身上乘以一个单调递增的线性函数,并在 N1(t−1)N_1(t-1)N1​(t−1) 身上乘以一个单调递减的线性函数:
N2(t)=t⋅N1(t)+(2−t)⋅N1(t−1)N_2(t)=t\cdot N_1(t)+(2-t)\cdot N_1(t-1)N2​(t)=t⋅N1​(t)+(2−t)⋅N1​(t−1)

为什么我们在 N1(t−1)N_1(t-1)N1​(t−1) 身上乘以的是一个 2−t2-t2−t 而不是 1−t1-t1−t 呢,这是因为 N1(t−1)N_1(t-1)N1​(t−1) 的非零区间为 [1,2][1, 2][1,2],如果乘以 1−t1-t1−t 会导致我们最终的结果出现负数,下图展示了 N1(t),N1(t−1)N_1(t), N_1(t-1)N1​(t),N1​(t−1) 和 N2(t)N_2(t)N2​(t):


其中绿色的虚线表示 N1(t)N_1(t)N1​(t),蓝色的虚线代表 N1(t−1)N_1(t-1)N1​(t−1),红色的虚线代表 N2(t)N_2(t)N2​(t)。不难看出 N2(t)N_2(t)N2​(t) 是一个宽度为 222 的局部基函数。这样得到的 N2(t)N_2(t)N2​(t) 已经是一个最高次为 111 次的单峰函数,使用这种函数作为基函数,我们可以得到一段连续的折线。例如,我们可以令:
M2(t)=A⋅N2(t)+B⋅N2(t−1),t∈[1,2]M_2(t)=A\cdot N_2(t)+B\cdot N_2(t-1), t\in[1, 2] M2​(t)=A⋅N2​(t)+B⋅N2​(t−1),t∈[1,2]
这样得到的 M2(t)M_2(t)M2​(t) 就是一个连接了 AAA 点和 BBB 点的线段。我们之所以要限制 t∈[1,2]t\in[1, 2]t∈[1,2] 是因为在 t∈[0,1]t\in [0, 1]t∈[0,1] 时只有 N1(t)N_1(t)N1​(t) 非零,N1(t−1)N_1(t-1)N1​(t−1) 为零,我们潜在地在对 AAA 点和原点进行加权求平均,而我们显然不希望原点影响到我们的图像的形状。假如我们将定义域定义成 t∈[0,3]t\in[0, 3]t∈[0,3],会导致 M2(t)M_2(t)M2​(t) 的图象变为了三段线段,第一段连接了原点和 AAA 点,第二段连接了 AAA 点和 BBB 点,第三段连接了 BBB 点和原点。我们还可以令:
M3(t)=A⋅N2(t)+B⋅N2(t−1)+C⋅N2(t−2),t∈[1,3]M_3(t)=A\cdot N_2(t)+B\cdot N_2(t-1)+C\cdot N_2(t-2), t\in[1, 3] M3​(t)=A⋅N2​(t)+B⋅N2​(t−1)+C⋅N2​(t−2),t∈[1,3]
这样我们得到的 M3(t)M_3(t)M3​(t) 就是两段线段的拼接,其中第一段连接了 AAA 点与 BBB 点,第二段连接了 BBB 点与 CCC 点。为什么能实现这个效果呢?是因为 N2(t−2)N_2(t-2)N2​(t−2) 只在 t∈[2,4]t\in[2, 4]t∈[2,4] 区间内有非零值,因此在区间 t∈[1,2]t\in[1, 2]t∈[1,2] 上M3(t)=M2(t)M_3(t)=M_2(t)M3​(t)=M2​(t),因此为一段连接了 AAA 和 BBB 的线段。

更普遍地,设 NLN_LNL​ 是一个宽度为 LLL 的局部基函数,对函数 NLN_LNL​ 的混淆记为 Mix[NL]Mix[N_L]Mix[NL​]:
NL+1(t)=Mix[NL](t)=tL⋅NL(t)+(L+1)−t(L+1)−1⋅NL(t−1)N_{L+1}(t)=Mix[N_L](t)=\frac{t}{L}\cdot N_L(t)+\frac{(L+1) - t}{(L+1) - 1} \cdot N_L(t-1) NL+1​(t)=Mix[NL​](t)=Lt​⋅NL​(t)+(L+1)−1(L+1)−t​⋅NL​(t−1)

简而言之就是在 NL(t)N_L(t)NL​(t) 身上乘以了一个他的非零区间内的线性递增函数,在 NL(t−1)N_L(t-1)NL​(t−1) 身上乘以了一个他的非零区间内的线性递减函数。下图给出了 N2(t),N2(t−1)N_2(t), N_2(t-1)N2​(t),N2​(t−1) 以及 N3(t)N_3(t)N3​(t) 的函数图像:


其中绿色的曲线是 N2(t)N_2(t)N2​(t),蓝色的曲线是 N2(t−1)N_2(t-1)N2​(t−1),红色的曲线是 N3(t)N_3(t)N3​(t)。使用 N3(t)N_3(t)N3​(t) 作为基函数已经可以较好的定义出一条连续且光滑的曲线了,例如我们令:

M3(t)=A⋅N3(t)+B⋅N3(t−1)+C⋅N3(t−2)+D⋅N3(t−3),t∈[2,4]M_3(t)=A\cdot N_3(t)+B\cdot N_3(t-1)+C\cdot N_3(t-2)+D\cdot N_3(t-3), t\in[2, 4] M3​(t)=A⋅N3​(t)+B⋅N3​(t−1)+C⋅N3​(t−2)+D⋅N3​(t−3),t∈[2,4]

和先前的例子类似,在这里我们需要将定义域限制到 [2,4][2, 4][2,4] 也是为了避免原点位置影响曲线的绘制。而我们这里用来做基函数的 N1(t),N2(t),⋯N_1(t), N_2(t), \cdotsN1​(t),N2​(t),⋯ 其实就是所谓的 B 样条基函数。具体而言,对于一个给定一个正整数 LLL,我们称 NL(t)N_L(t)NL​(t) 为 LLL 阶(L−1L-1L−1 次)B 样条基函数。

四阶三次均匀 B 样条曲线

对 N3(t)N_3(t)N3​(t) 进行混淆可以得到 N4(t)N_4(t)N4​(t),其中 N4(t)N_4(t)N4​(t) 是一个宽度为 444 的局部基函数,我们称它为四阶三次 B 样条基函数,下图展示了它的大致形状:

二阶 B 样条至少要有两个控制点(不然线段怎么连),三阶 B 样条至少要有三个控制点(不然二次曲线怎么画),因此四阶 B 样条至少要有四个控制点。B 样条相对于贝塞尔曲线的有点在于,即使四阶三次 B 样条有多于四个控制点,它的次数仍然是三次,而且每个控制点的变化只会影响与他接近的一段曲线,而不会影响整个目标曲线(因为 N4(t)N_4(t)N4​(t) 宽度为 444,这也就意味着与 N4(t)N_4(t)N4​(t) 相乘的那个控制点不会在这个宽度之外造成任何影响)。

假设我们有 nnn 个常数点 P0,P1,⋯,Pn−1P_0, P_1, \cdots, P_{n-1}P0​,P1​,⋯,Pn−1​,其中 n≥4n\geq 4n≥4,我们可以定义一条四阶三次均匀 B 样条曲线:
M(t)=P0⋅N4(t−0)+P1⋅N4(t−1)+⋯+Pn−1⋅N4(t−(n−1)),t∈[3,n]M(t)=P_0\cdot N_4(t-0)+P_1\cdot N_4(t-1)+\cdots+P_{n-1}\cdot N_4(t-(n-1)),t\in[3, n] M(t)=P0​⋅N4​(t−0)+P1​⋅N4​(t−1)+⋯+Pn−1​⋅N4​(t−(n−1)),t∈[3,n]

和贝塞尔曲线相比,B 样条往往不会经过自己的第一个控制点 P0P_0P0​ 和 Pn−1P_{n-1}Pn−1​,而贝塞尔曲线一定会经过这两个点。

B 样条绘制范例

随便画了一个有七个控制点的 B 样条,特别地,如果我们让 P2=P6,P1=P5,P0=P4P_2=P_6, P_1=P_5, P_0=P_4P2​=P6​,P1​=P5​,P0​=P4​,那么我们能够绘制一条闭合的 B 样条曲线。

换言之,任给一个多边形,我们都能找到这个多边形控制的一条闭合 B 样条曲线。如果这个多边形是一个凸多边形,感性理解一下不能得到这个闭合的 B 样条曲线一定完全位于整个多边形的内部。


上文中的所有图形都是使用几何画板绘制的感兴趣的同学可以联系我(邮箱见免责)获得这些几何画板文件,下图中给出了用于生成这两个 B 样条范例的配置:

未完待续 …

一个理念

不得不说我们的教材内容丰富充实,兼顾了严谨性与知识性。但我对我们的图形学教材仍有一个略有不满的地方,就是教材中没有贯彻一个目的性与实践性的假设:任何工业设计理论都带有其目的。当我们提出一个概念/公式时,指导我们的所谓“灵感”不是拍脑袋想出来的,更不是从天上掉下来的,而是来自真真正正的生产实践的需要。我们的图形学教材并不是没有这种理念,但是它没有贯彻这种理念,换言之,这种理念仅仅停留在每一章的简介部分,而没有真正深入到每一章的具体内容中。不要先给出一个复杂的公式,再利用这个公式解决问题,要先从最简单、最必要、最直观的概念入手逐渐让我们知道我们为什么要这么做(目的),而不是先告诉我们怎么做(方法)。否则学生很容易不知所云。

计算机图形学中的曲线问题相关推荐

  1. 计算机图形学中的曲线问题——贝塞尔曲线的绘制

    贝塞尔曲线的绘制 由于 CSDN 的博客修改字数的限制,我们不得不将这一部分放到一个新的博客中.原文详见: GGN_2015 计算机图形学中的曲线问题 贝塞尔曲线的几何作图法 在上面介绍儿时的回忆中, ...

  2. 数学系列:数学在计算机图形学中的应用

    宇宙的琴弦 博客园 首页 新随笔 联系 订阅 管理 随笔 - 60   文章 - 0   评论 - 0 数学系列:数学在计算机图形学中的应用 Copyright © 1900-2016, NORYES ...

  3. 数学在计算机图形学中的应用

    数学在计算机图形学中的应用 刘利刚 中国科技大学 "学习计算机图形学需要多少的数学?"这是初学者最经常问的问题. 狭义的计算机图形学指的是传统的三维建模,绘制,动画等,而广义的计算 ...

  4. 《 线性代数及其应用 (原书第4版)》—— 2.7 计算机图形学中的应用

    本节书摘来自华章出版社< 线性代数及其应用 (原书第4版)>一书中的第2章,第2.7节,作者:(美)戴维C. 雷(David C. Lay)马里兰大学帕克学院 著刘深泉 张万芹 陈玉珍 包 ...

  5. 计算机图形学中的四元数(Quaternions)

    计算机图形学中的四元数 计算机图形学中的四元数(Quaternions) 前言 欧拉角 欧拉角的万向节死锁(Gimbal Lock)与插值问题 四元数 参考 计算机图形学中的四元数(Quaternio ...

  6. 欧拉公式在计算机图形学中的,计算机图形学 第九章课件.ppt

    <计算机图形学 第九章课件.ppt>由会员分享,提供在线免费全文阅读可下载,此文档格式为ppt,更多相关<计算机图形学 第九章课件.ppt>文档请在天天文库搜索. 1.甘朝华第 ...

  7. 计算机图形学中向量点乘和叉乘的用途_图形学笔记(一):基础知识

    从这便文章开始整理学习到的计算机图像学相关知识,原则是只写我没在网上找到清楚解释的内容,如果有很好的文章介绍相关内容,我会直接把链接贴上. 首先弄清 Computer Graphics和 Comput ...

  8. c语言 连通域算法 递归,VC++ 6.0编写计算机图形学中的种子填充算法,想用递归的八向连通域,求助!...

    VC++ 6.0编写计算机图形学中的种子填充算法,想用递归的八向连通域,求助!0 填充函数代码如下: void CComputerGraphicsView::PolygonFill2()//区域填充函 ...

  9. 计算机图形学中需要掌握的数学基础知识有哪些?

    计算机图形学中使用了大量数学知识,尤其是矩阵和线性代数.虽然我们倾向于认为3D图形编程是紧跟最新技术的领域之一(它在很多方面确实是),但它用到的很多技术实际上可以追溯到上百年前,其中一些甚至是由文艺复 ...

最新文章

  1. jvm:虚方法与非虚方法
  2. Windows Phone笔记索引(总)
  3. 企业五年后卓越或者死亡,数据战略是关键!
  4. 让网络更轻盈——网络功能虚拟化技术的现状和未来(中兴通讯)
  5. centos磁盘满了,查找大文件并清理
  6. 如何在计算机命令内转换操作盘,如何在命令行窗口中从驱动器C切换到驱动器D...
  7. Luogu6186 [NOI Online #1 提高组] 冒泡排序
  8. 中科院计算机和理论物理双硕士白,中科院研究生理论物理怎么不学相对论?
  9. 华为光猫设置及拨号连接下开启移动热点
  10. 学校机房计算机类型,学校计算机机房的管理和维护建议原稿(备份存档)
  11. 笔记本不小心网络重置后,不能上网,网络适配器存在感叹号
  12. flowchart.js使用总结
  13. 微信小程序解密出来是乱码的问题
  14. 华为的服务器虚拟化软件,服务器虚拟化软件
  15. OpenSSL ssL_read: Connection was aborted,errno 10053 报错
  16. C++的atof()
  17. 天天基金数据接口的处理
  18. paddlehub创意赛《王者四大美女--红昭愿》
  19. python 编程题 埃及金字塔罐子倒水
  20. CNCC 2017大会第一天,邱成桐,梅宏,沈向洋,李飞飞,汤道生,马维英都讲了什么?...

热门文章

  1. libc and libc++
  2. 屏蔽无法验证发行者,你确实要运行此软件吗?的提示
  3. 智能硬件产品经理需要哪些技术基础?
  4. Windows下查看mysql是否启动
  5. Python 实现 周志华 《机器学习》 BP算法(高级版)
  6. 支付宝报错: invalid-signature 错误原因: 验签出错,建议检查签名字符串或签名私钥与应用公钥是否匹配,网关生成的验签字符串为:xxx
  7. 基于PHP的篮球宝篮球娱乐网站
  8. 科学史上最伟大的十位单身科学家
  9. jsp项目发布到服务器
  10. 5码默认版块_用字节码解释try、catch、finally、i++、++i的执行结果?