文章目录
免责
又来快乐地误人子弟了。 最近学习了一些计算机图形学中的曲线与曲面,感觉很多东西还是蛮有趣的,就打算写篇文章总结总结。
文章中如有疏漏,请联系 premierbob2015@gmail.com
进行更正,笔者感激不尽。
前言
在工业设计中,我们经常会遇到各种各样的绘图问题。随着计算机技术的发展,现代的工业设计正在逐步脱离纸质的工图,转而使用由计算机辅助绘图并设计的电子工图。这种电子版的设计图不仅便于存储与传输,还能更直接地转化为控制数控机床 / 3D打印机 的控制信号,实现工业生产的全自动化。
经过先前的学习,我们了解到在计算机中存储的图像大致可分为点阵图和矢量图两类。作为存储图像的不同格式,点阵图和矢量图各有千秋。而上文中所说的电子工图往往使用矢量图的格式进行存储。
优点 | 缺点 | |
---|---|---|
点阵图 | 易于从相机等设备直接获得,适合近似描述真实世界中的事物 | 使用像素矩阵存储,不能无限放大,放大会使得图像变得不清晰 |
矢量图 | 使用函数或参数方程存储图像,可以无限放大 | 难以从相机等设备直接获得 |
本文并不介绍这些图纸在计算机中具体以何种矢量图格式存储,只在数学层面介绍一些非常简单的曲线,而这些曲线在工业设计中有着广泛的应用。(当然,我对工业设计没啥兴趣,我只希望自己将来有一天,能够开发一款用于动漫的绘制辅助工具——属于是老二刺螈了。)
拉格朗日插值多项式
(不重合) 两点确定一条直线,(不共线) 三点确定一条二次函数,不难让我们联想到
n
n
n 个定点能确定一条的
n
−
1
n-1
n−1 次多项式曲线。形式化的来说给定函数上的
n
n
n 个横坐标互不相同的点:
P
1
=
(
x
1
,
y
1
)
,
P
2
=
(
x
2
,
y
2
)
,
⋯
,
P
n
=
(
x
n
,
y
n
)
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
=
a
n
−
1
x
n
−
1
+
a
n
−
2
x
n
−
2
+
⋯
+
a
1
x
+
a
0
y=a_{n-1}x^{n-1}+a_{n-2}x^{n-2}+\cdots+a_1x+a_0
y=an−1xn−1+an−2xn−2+⋯+a1x+a0
其中
a
0
,
a
1
,
⋯
,
a
n
−
1
a_{0}, a_1, \cdots, a_{n-1}
a0,a1,⋯,an−1 是函数中的参数,
x
x
x 是自变量,
y
y
y 是因变量。所谓确定就是指,将这
n
n
n 个点坐标带入函数中后,我们能够得到系数向量
a
0
,
a
1
,
⋯
,
a
n
−
1
a_0, a_1, \cdots, a_{n-1}
a0,a1,⋯,an−1 的至少一个解。在某些特定条件下,我们能够得到这些系数的一个唯一解,严谨的证明可以参考范德蒙德行列式,在此不再赘述。
对于拉格朗日插值而言,即使我们没有从数学上证明这种可解性按照下面给出的方法,我们也可以使用一种构造性的方式构造出一个恰好通过 P 1 , P 2 , ⋯ , P n P_1, P_2, \cdots, P_n P1,P2,⋯,Pn 这 n n n 个点的一个 n − 1 n-1 n−1 次函数。
第一步:得到一个基函数
既然我们不能很轻松地让我们的曲线经过所有点,那我们就先试着让我们的曲线经过其中的一个点而在其他点横坐标对应位置的值为零。令 x 1 , x 2 , ⋯ , x n x_1, x_2, \cdots, x_n x1,x2,⋯,xn 为互不相同的常数,我们能否构造一个 n − 1 n-1 n−1 次多项式函数 f 1 ( x ) f_1(x) f1(x) 使得 :
- f 1 ( x 1 ) = 1 f_1(x_1)=1 f1(x1)=1;
- f 1 ( x 2 ) = f 1 ( x 3 ) = ⋯ = f 1 ( x n ) = 0 f_1(x_2)=f_1(x_3)=\cdots=f_1(x_n)=0 f1(x2)=f1(x3)=⋯=f1(xn)=0 ;
很显然是可以的。首先看到这个函数
f
1
(
x
)
f_1(x)
f1(x) 有
n
−
1
n-1
n−1 个零点,分别是
x
2
,
x
3
,
⋯
,
x
n
x_2, x_3, \cdots, x_n
x2,x3,⋯,xn,我们很自然地想到设
g
1
(
x
)
=
(
x
−
x
2
)
(
x
−
x
3
)
⋯
(
x
−
x
n
)
g_1(x) = (x-x_2)(x-x_3)\cdots(x-x_n)
g1(x)=(x−x2)(x−x3)⋯(x−xn)。虽然这个函数已经满足了
g
1
(
x
2
)
=
g
1
(
x
3
)
=
⋯
=
g
1
(
x
n
)
=
0
g_1(x_2)=g_1(x_3)=\cdots=g_1(x_n)=0
g1(x2)=g1(x3)=⋯=g1(xn)=0 这一条件,但是我们并不能保证
g
1
(
x
1
)
=
1
g_1(x_1)=1
g1(x1)=1。由于
x
1
,
x
2
,
⋯
,
x
n
x_1, x_2, \cdots, x_n
x1,x2,⋯,xn 互不相同,所以
g
(
x
1
)
=
(
x
1
−
x
2
)
(
x
1
−
x
3
)
⋯
(
x
1
−
x
n
)
g(x_1)=(x_1-x_2)(x_1-x_3)\cdots(x_1-x_n)
g(x1)=(x1−x2)(x1−x3)⋯(x1−xn) 一定不等于零。不妨令:
f
1
(
x
)
=
g
1
(
x
)
g
1
(
x
1
)
f_1(x)=\frac{g_1(x)}{g_1(x_1)}
f1(x)=g1(x1)g1(x)
因为 x 1 x_1 x1 是常数,所以 g 1 ( x 1 ) g_1(x_1) g1(x1) 就是常数,而这样得到的 f 1 ( x ) f_1(x) f1(x) 是在 g 1 ( x ) g_1(x) g1(x) 基础之上进行了一次常数倍的缩放。此时我们所要求的两个条件就都被满足了。
第二步:得到所有基函数
在第一步中,我们定义了一个名为
f
i
f_i
fi 的函数。仿照第一步的思路,对于
i
=
1
,
2
,
⋯
,
n
i=1, 2, \cdots, n
i=1,2,⋯,n,我们定义:
f
i
(
x
)
=
g
i
(
x
)
g
i
(
x
i
)
f_i(x)=\frac{g_i(x)}{g_i(x_i)}
fi(x)=gi(xi)gi(x)
其中:
g
i
(
x
)
=
∏
j
=
1
,
2
,
⋯
n
且
j
≠
i
(
x
−
x
j
)
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)
我们就得到了 n n n 个函数 f 1 ( x ) , f 2 ( x ) , ⋯ , f n ( x ) f_1(x), f_2(x), \cdots, f_n(x) f1(x),f2(x),⋯,fn(x),根据第一步中的描述, f i ( x ) f_i(x) fi(x) 一定满足:
- f i ( x i ) = 1 f_i(x_i)=1 fi(xi)=1;
- f i ( x j ) = 0 , f_i(x_j)=0, fi(xj)=0, 对 j = 1 , 2 , ⋯ , n j=1, 2, \cdots, n j=1,2,⋯,n 且 j ≠ i j\neq i j=i;
这些函数如此的重要以至于我们要给这些函数起一个好听的名字——拉格朗日基函数。
第三步:对所有基函数进行线性组合
有了拉格朗日基函数后,我们令:
f ( x ) = y 1 ⋅ f 1 ( x ) + y 2 ⋅ f 2 ( x ) ⋯ + y n ⋅ f n ( 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) 一定经过点
P
1
(
x
1
,
y
1
)
,
P
2
(
x
2
,
y
2
)
,
⋯
P
n
(
x
n
,
y
n
)
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
=
x
1
x=x_1
x=x1 时,
f
1
(
x
)
=
f
1
(
x
1
)
=
1
f_1(x)=f_1(x_1)=1
f1(x)=f1(x1)=1 而
f
2
(
x
)
=
f
3
(
x
)
=
⋯
=
f
n
(
x
)
=
0
f_2(x)=f_3(x)=\cdots=f_n(x)=0
f2(x)=f3(x)=⋯=fn(x)=0。因此:
f
(
x
1
)
=
y
1
⋅
1
+
y
2
⋅
0
+
y
3
⋅
0
+
⋯
+
y
n
⋅
0
=
y
1
f(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 ( x 2 ) = y 2 , f ( x 3 ) = y 3 , ⋯ , f ( x n ) = y n f(x_2)=y_2, f(x_3)=y_3, \cdots, f(x_ n)=y_n f(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)。则可以计算出:
y 1 f 1 y_1f_1 y1f1 | y 2 f 2 y_2f_2 y2f2 | y 3 f 3 y_3f_3 y3f3 | y 4 f 4 y_4f_4 y4f4 | |
---|---|---|---|---|
y i f i ( x ) y_if_i(x) yifi(x) | 7 3 ( 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) | 5 3 ( x + 1 ) ( x − 0 ) ( x − 1 ) \frac 5 3 (x+1)(x-0)(x-1) 35(x+1)(x−0)(x−1) |
y i f i ( − 1 ) y_if_i(-1) yifi(−1) | − 14 -14 −14 | 0 0 0 | 0 0 0 | 0 0 0 |
y i f i ( 0 ) y_if_i(0) yifi(0) | 0 0 0 | − 4 -4 −4 | 0 0 0 | 0 0 0 |
y i f i ( 1 ) y_if_i(1) yifi(1) | 0 0 0 | 0 0 0 | − 4 -4 −4 | 0 0 0 |
y i f i ( 2 ) y_if_i(2) yifi(2) | 0 0 0 | 0 0 0 | 0 0 0 | 10 10 10 |
拉格朗日插值曲线绘制实践
限于篇幅我们把这部分内容放到了另一篇博客中,感兴趣的同学可以来看一看。
GGN_2015
拉格朗日插值曲线绘制实践
三次埃尔米特插值多项式
在工业设计中,有时我们不仅需要指定一条曲线经过哪些点,我们还需要指定这条曲线在某些位置的导数。
考虑一个最简单的情况:设 f ( x ) = a 3 ⋅ x 3 + a 2 ⋅ x 2 + a 1 ⋅ x + a 0 f(x)=a_3\cdot x^3+a_2\cdot x^2 + a_1 \cdot x + a_0 f(x)=a3⋅x3+a2⋅x2+a1⋅x+a0 是一个三次多项式。我们要找到一个 f ( x ) f(x) f(x) 使得它经过点 P 0 ( x 0 , y 0 ) P_0(x_0, y_0) P0(x0,y0) 和 P 1 ( x 1 , y 1 ) P_1(x_1, y_1) P1(x1,y1) 且在点 P 0 P_0 P0 处的导数为 y 0 ′ y_0' y0′,在点 P 1 P_1 P1 处的导数为 y 1 ′ y_1' y1′,其中 x 0 , x 1 , y 0 , y 1 , y 0 ′ , y 1 ′ x_0, x_1, y_0, y_1, y_0', y_1' x0,x1,y0,y1,y0′,y1′ 均为常数。
使用类似拉格朗日插值法中的构造方法,我们也能够得到这个 f ( x ) f(x) f(x)。
第一步:得到第一维基函数
设 x 0 ≠ x 1 x_0\neq x_1 x0=x1 且二者均为常数,令 p 0 ( x ) p_0(x) p0(x) 是一个三次多项式,它满足如下两个性质:
- p 0 ( x 0 ) = 1 p_0(x_0)=1 p0(x0)=1;
- p 0 ( x 1 ) = 0 , p 0 ′ ( x 1 ) = 0 p_0(x_1)=0, p_0'(x_1)=0 p0(x1)=0,p0′(x1)=0;
由于 p 0 ( x 1 ) = p 0 ′ ( x 1 ) = 0 p_0(x_1)=p_0'(x_1)=0 p0(x1)=p0′(x1)=0,我们很直观地感受到, x 1 x_1 x1 应该是函数 p 0 ( x ) p_0(x) p0(x) 的重根。不妨令 g 0 ( x ) = ( x − x 1 ) 2 g_0(x)=(x-x_1)^2 g0(x)=(x−x1)2,那么我们有 g 0 ′ ( x ) = 2 ( x − x 1 ) g_0'(x)=2(x-x_1) g0′(x)=2(x−x1),因此 g 0 ( x 1 ) = g 0 ′ ( x 1 ) = 0 g_0(x_1)=g_0'(x_1)=0 g0(x1)=g0′(x1)=0。
按照拉格朗日插值的思路,我们构造了一个函数:
p
0
(
x
)
=
g
0
(
x
)
g
0
(
x
0
)
p_0(x)=\frac{g_0(x)}{g_0(x_0)}
p0(x)=g0(x0)g0(x)
因为 x 0 x_0 x0 是常数,所以 g 0 ( x 0 ) g_0(x_0) g0(x0) 也是常数,我们可以看到:
p 0 ( x 0 ) = 1 ( 因 为 我 们 构 造 的 系 数 恰 好 是 g 0 ( x 0 ) 的 倒 数 ) p 0 ′ ( x 0 ) = g 0 ′ ( x 0 ) g 0 ( x 0 ) = 2 ( x 0 − x 1 ) ( x 0 − x 1 ) 2 = 2 x 0 − x 1 ≠ 0 p 0 ( x 1 ) = p 0 ′ ( x 1 ) = 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−x12=0=p0′(x1)=0
有了函数 p 0 p_0 p0,我们就可以在不改变 p 0 ( x 1 ) p_0(x_1) p0(x1) 和 p 0 ′ ( x 1 ) p_0'(x_1) p0′(x1) 的值的前提下随意控制 p 0 ( x 0 ) p_0(x_0) p0(x0) 的值了。但是现在我们还无法控制 p 0 ′ ( x 0 ) p_0'(x_0) p0′(x0) 的值。因此我们要在第二步构造一个函数,使得我们能够分离 p 0 ′ ( x 0 ) p_0'(x_0) p0′(x0) 和 p 0 ( x 0 ) p_0(x_0) p0(x0) 之间的倍数关系。
第二步:得到第二维基函数
类似第一步,我们还需要构造三次函数 q 0 ( x ) q_0(x) q0(x) 满足这样两条性质:
- q 0 ( x 0 ) = 0 , q 0 ′ ( x 0 ) = 1 q_0(x_0)=0, q_0'(x_0)=1 q0(x0)=0,q0′(x0)=1;
- q 0 ( x 1 ) = 0 , q 0 ′ ( x 1 ) = 0 q_0(x_1)=0, q_0'(x_1)=0 q0(x1)=0,q0′(x1)=0;
看起来也不难构造,零点又多了一个
x
0
x_0
x0。令
k
0
(
x
)
=
(
x
−
x
0
)
(
x
−
x
1
)
2
k_0(x)=(x-x_0)(x-x_1)^2
k0(x)=(x−x0)(x−x1)2,则:
k
0
′
(
x
)
=
(
x
−
x
1
)
2
+
2
(
x
−
x
0
)
(
x
−
x
1
)
=
(
x
−
x
1
)
(
3
x
−
2
x
0
−
x
1
)
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)
令:
q
0
(
x
)
=
k
0
(
x
)
k
0
′
(
x
0
)
q_0(x)=\frac{k_0(x)}{k_0'(x_0)}
q0(x)=k0′(x0)k0(x)
由于
x
0
x_0
x0 是常数因此
k
0
′
(
x
0
)
=
(
x
0
−
x
1
)
2
k_0'(x_0)=(x_0-x_1)^2
k0′(x0)=(x0−x1)2 也是常数,我们可以看到:
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
\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
第三步:对基函数进行线性组合
开门见山,令:
f
0
(
x
)
=
y
0
⋅
(
p
0
(
x
)
−
2
x
0
−
x
1
⋅
q
0
(
x
)
)
+
y
0
′
⋅
q
0
(
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−x12⋅q0(x))+y0′⋅q0(x)
为什么要这么构造呢? p 0 ( x ) p_0(x) p0(x) 在 x 0 x_0 x0 处的导数值为 2 x 0 − x 1 ≠ 0 \frac{2}{x_0-x_1}\neq 0 x0−x12=0,因此在它身上减去一个 2 x 0 − x 1 ⋅ q 0 ( x ) \frac{2}{x_0-x_1}\cdot q_0(x) x0−x12⋅q0(x) 就能使得它在 x 0 x_0 x0 处的导数为零,而由于 q 0 ( x ) q_0(x) q0(x) 在 x 1 x_1 x1 点的函数值和导数值都为零,因此无论我们如何将 p 0 ( x ) p_0(x) p0(x) 和 q 0 ( x ) q_0(x) q0(x) 进行线性组合,都不会影响函数在 x 1 x_1 x1 点处的函数值和导数值。
此时函数 f 0 ( x ) f_0(x) f0(x) 满足:
- f 0 ( x 0 ) = y 0 , f 0 ′ ( x 0 ) = y 0 ′ f_0(x_0)=y_0,f_0'(x_0)=y_0' f0(x0)=y0,f0′(x0)=y0′;
- f 0 ( x 1 ) = 0 , f 0 ′ ( x 1 ) = 0 f_0(x_1)=0, f_0'(x_1)=0 f0(x1)=0,f0′(x1)=0;
类似地,我们可以使用同样的方法得到一个函数 f 1 ( x ) f_1(x) f1(x),满足:
- f 1 ( x 0 ) = 0 , f 1 ′ ( x 0 ) = 0 f_1(x_0)=0,f_1'(x_0)=0 f1(x0)=0,f1′(x0)=0;
- f 1 ( x 1 ) = y 1 , f 1 ′ ( x 1 ) = y 1 ′ f_1(x_1)=y_1, f_1'(x_1)=y_1' f1(x1)=y1,f1′(x1)=y1′;
具体怎么做不再赘述,整理后得到的表达式我也不是很关心。我所关心的只有一件事那就是,当我们令 f ( x ) = f 0 ( x ) + f 1 ( x ) f(x)=f_0(x)+f_1(x) f(x)=f0(x)+f1(x),这样得到的 f ( x ) f(x) f(x) 一定经过 P 0 ( x 0 , y 0 ) P_0(x_0, y_0) P0(x0,y0) 和 P 1 ( x 1 , y 1 ) P_1(x_1, y_1) P1(x1,y1) 且在点 P 0 P_0 P0 处的导数为 y 0 ′ y_0' y0′,在点 P 1 P_1 P1 处的导数为 y 1 ′ y_1' y1′。
一些后话
如果我们想要找到一个函数使得他定义在 [ x 1 , x n ] [x_1, x_n] [x1,xn] 上且函数光滑(可导),且经过点 P 1 ( x 1 , y 1 ) , P 2 ( x 2 , y 2 ) , ⋯ , P n ( x n , y n ) 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),且在 P i P_i Pi 处导数值为 y i ′ y_i' yi′( i = 1 , 2 ⋯ , n i=1, 2\cdots, n i=1,2⋯,n),其中 x i , y i , y i ′ x_i, y_i, y_i' xi,yi,yi′ 均为常数( i = 1 , 2 , ⋯ , n i=1, 2, \cdots, n i=1,2,⋯,n),且满足 x 1 < x 2 < ⋯ < x n x_1<x_2<\cdots<x_n x1<x2<⋯<xn 。
我们可以在区间 [ x 1 , x 2 ] [x_1, x_2] [x1,x2] 上依据 y 1 , y 2 , y 1 ′ , y 2 ′ y_1, y_2, y_1', y_2' y1,y2,y1′,y2′ 构造一个三次埃尔米特插值函数,在 [ x 2 , x 3 ] [x_2, x_3] [x2,x3] 上依据 y 2 , y 3 , y 2 ′ , y 3 y_2, y_3, y_2', y_3 y2,y3,y2′,y3 构造插值函数, ⋯ \cdots ⋯,在 [ x n − 1 , x n ] [x_{n-1}, x_{n}] [xn−1,xn] 上根据 y n − 1 , y n , y n − 1 ′ , y n ′ y_{n-1}, y_{n}, y_{n-1}', y_n' yn−1,yn,yn−1′,yn′ 构造插值函数。最后再将每一段上的函数以分段函数的形式连接起来,由于连接点处的左右函数值、左右导数值均相等,因此函数在区间 ( x 1 , x n ) (x_1, x_n) (x1,xn) 上光滑。
规范化的三次埃尔米特多项式
如果我们的目的只是绘制一个函数,使用普通的埃尔米特多项式似乎已经实现了我们想要的。但是函数这东西,自变量 x x x 的每一个取值至多只有一个 y y y 值与之对应,如果我想绘制任意平面曲线甚至空间曲线,那上文给出的方法就不那么好用了。想到可以用参数方程的方法来实现。
- 下文中我们令 P 0 , P 1 , D 0 , D 1 P_0, P_1, D_0, D_1 P0,P1,D0,D1 是二维常向量,我们要构造参数方程 f ⃗ ( t ) \vec{f}(t) f(t) 使得 f ⃗ ( 0 ) = P 0 , f ⃗ ( 1 ) = P 1 , ∇ f ⃗ ( 0 ) = D 0 , ∇ f ⃗ ( 1 ) = D 1 \vec f (0)=P_0, \vec f(1) = P_1, \nabla \vec f(0) = D_0, \nabla \vec f(1)=D_1 f(0)=P0,f(1)=P1,∇f(0)=D0,∇f(1)=D1。
注:梯度算子 ∇ \nabla ∇,其实就是对函数值的 x x x, y y y 分量函数分别求导,例如当 f ⃗ ( t ) = [ f x ( t ) , f y ( t ) ] \vec f(t) = [f_x(t), f_y(t)] f(t)=[fx(t),fy(t)] 时 ∇ f ⃗ ( t ) = [ f x ′ ( t ) , f y ′ ( t ) ] \nabla \vec f(t)=[f_x'(t), f_y'(t)] ∇f(t)=[fx′(t),fy′(t)]。
对于一个标量函数 f ( t ) f(t) f(t) 而言,假如 P ∈ R 2 P\in \R^2 P∈R2 是一个二维常向量(即 P x , P y ∈ R P_x, P_y\in \R Px,Py∈R 是两个常数),那么 P f ( t ) = [ P x f ( t ) , P y f ( t ) ] Pf(t)=[P_xf(t), P_yf(t)] Pf(t)=[Pxf(t),Pyf(t)] 就是一个平面上的参数方程。对这个函数计算梯度,得到 ∇ ( P f ( t ) ) = [ P x f ′ ( t ) , P y f ′ ( t ) ] = P ⋅ f ′ ( t ) \nabla(Pf(t))=[P_xf'(t), P_yf'(t)]=P\cdot f'(t) ∇(Pf(t))=[Pxf′(t),Pyf′(t)]=P⋅f′(t)。此处的点运算 ‘ ⋅ \cdot ⋅’ 不表示向量的内积,而表示“数乘”,其中 P ∈ R 2 P\in \R^2 P∈R2 是 “向量”, f ( t ) ∈ R f(t)\in \R f(t)∈R 是定义在实数上的标量函数。
类似上文的方法,我们仍然可以构造出符合条件的向量函数 f ⃗ ( t ) \vec f(t) f(t)。
第一步:得到第一维基函数
在这一步中我们要构造标量函数 p 0 ( t ) p_0(t) p0(t) 和 p 1 ( t ) p_1(t) p1(t),其中 p 0 ( t ) p_0(t) p0(t) 满足:
- p 0 ( 0 ) = 1 p_0(0)=1 p0(0)=1;
- p 0 ( 1 ) = 0 , p 0 ′ ( 1 ) = 0 p_0(1)=0, p_0'(1)=0 p0(1)=0,p0′(1)=0;
同理, p 1 ( t ) p_1(t) p1(t) 满足:
- p 1 ( 1 ) = 1 p_1(1)=1 p1(1)=1;
- p 1 ( 0 ) = 0 , p 1 ′ ( 0 ) = 0 p_1(0)=0, p_1'(0)=0 p1(0)=0,p1′(0)=0;
为了方便表述,我们把 p 0 ( t ) p_0(t) p0(t) 称为 t = 0 t=0 t=0 处的常数控制基函数, p 1 ( t ) p_1(t) p1(t) 称为 t = 1 t=1 t=1 处的常数控制基函数,它们的存在就是为了用于将它们线性组合以调整函数在 t = 0 t=0 t=0 处以及 t = 1 t=1 t=1 处的函数值。
由于这两个函数都比较简单,稍微构造构造就能构造出来(拿眼睛看就能看出来,不信你试试):
p
0
(
t
)
=
(
t
−
1
)
2
p
1
(
t
)
=
t
2
\begin{aligned} p_0(t)&=(t-1)^2\\ p_1(t)&=t^2 \end{aligned}
p0(t)p1(t)=(t−1)2=t2
为了方便下一步中消去 p 0 ( t ) p_0(t) p0(t) 和 p 1 ( t ) p_1(t) p1(t) 对导数的影响,我们可以提前计算出 p 0 ′ ( 0 ) p_0'(0) p0′(0) 和 p 1 ′ ( 1 ) p_1'(1) p1′(1):
p 0 ′ ( t ) = 2 ( t − 1 ) ∴ p 0 ′ ( 0 ) = − 2 p 1 ′ ( t ) = 2 t ∴ p 1 ′ ( 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
第二步:得到第二维基函数
类似上文中对三次埃尔米特插值多项式的构造,我们要构造 q 0 ( t ) q_0(t) q0(t) 和 q 1 ( t ) q_1(t) q1(t),其中 q 0 ( t ) q_0(t) q0(t) 满足:
- q 0 ( 0 ) = 0 , q 0 ′ ( 0 ) = 1 q_0(0)=0, q_0'(0)=1 q0(0)=0,q0′(0)=1;
- q 0 ( 1 ) = 0 , q 0 ′ ( 1 ) = 0 q_0(1)=0, q_0'(1)=0 q0(1)=0,q0′(1)=0;
类似地, q 1 ( t ) q_1(t) q1(t) 满足:
- q 1 ( 0 ) = 0 , q 1 ′ ( 0 ) = 0 q_1(0) = 0, q_1'(0) = 0 q1(0)=0,q1′(0)=0;
- q 1 ( 1 ) = 0 , q 1 ′ ( 1 ) = 1 q_1(1) = 0, q_1'(1)=1 q1(1)=0,q1′(1)=1;
想解方程也不是不行,只不过确实可以直接看出结果:
q
0
(
t
)
=
t
(
t
−
1
)
2
q
1
(
t
)
=
t
2
(
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)
同样为了方便表述,我们把 q 0 ( t ) q_0(t) q0(t) 称为 t = 0 t=0 t=0 处的导数控制基函数, q 1 ( t ) q_1(t) q1(t) 称为 t = 1 t=1 t=1 处的导数控制基函数,它们有两重作用:一重是和 p 0 p_0 p0 或 p 1 p_1 p1 线性组合,从而消去常数控制基函数对导函数的影响,另外一重作用是设置导函数的值。接下来我们分别从这两个角度来阐述导数控制基函数的作用。
考虑,第一方面的作用,记:
g
0
(
t
)
=
p
0
(
t
)
−
p
0
′
(
0
)
⋅
q
0
(
t
)
=
(
t
−
1
)
2
+
2
t
(
t
−
1
)
2
=
(
t
−
1
)
2
(
2
t
+
1
)
g
1
(
t
)
=
p
1
(
t
)
−
p
1
′
(
1
)
⋅
q
1
(
t
)
=
t
2
−
2
t
2
(
t
−
1
)
=
t
2
(
3
−
2
t
)
\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)
这样我们得到的 g 0 ( t ) g_0(t) g0(t) 和 g 1 ( t ) g_1(t) g1(t) 就能够实现导数无关的常数控制,因为我们用 q 0 q_0 q0 抵消了 p 0 p_0 p0 在 x = 0 x=0 x=0 处对导数的影响,用 q 1 q_1 q1 抵消了 p 1 p_1 p1 在 x = 1 x=1 x=1 处的导数影响。为了和 p 0 p_0 p0 与 p 1 p_1 p1 加以区分,我们称 g 0 g_0 g0 和 g 1 g_1 g1 为导数无关的控制基函数。
顺便说一嘴, g 0 , g 1 , q 0 , q 1 g_0, g_1, q_0, q_1 g0,g1,q0,q1 这四个函数可以用待定系数法列方程求解,但是列方程多无聊啊,构造出来多有趣。
第三步:对基函数进行线性组合
f ⃗ ( t ) = P 0 ⋅ g 0 ( t ) + P 1 ⋅ g 1 ( t ) + D 0 ⋅ q 0 ( t ) + D 1 ⋅ q 1 ( 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]
不难证明这个东西就是我们想要的,把含
t
t
t 的解析式带入得到:
f
⃗
(
t
)
=
P
0
(
2
t
3
−
3
t
2
+
1
)
+
P
1
(
−
2
t
3
+
3
t
2
)
+
D
0
(
t
3
−
2
t
2
+
t
)
+
D
1
(
t
3
−
t
2
)
,
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 ∈ [ t 1 , t n ] \vec f(t),t\in[t_1, t_n] f(t),t∈[t1,tn] 使得 t = t 1 t=t_1 t=t1 时 f ⃗ ( t ) = P 1 \vec f(t) = P_1 f(t)=P1 且 ∇ f ⃗ ( t ) = D 1 \nabla \vec f(t) = D_1 ∇f(t)=D1, t = t 2 t=t_2 t=t2 时 f ⃗ ( t ) = P 2 \vec f(t) = P_2 f(t)=P2 且 ∇ f ⃗ ( t ) = D 2 \nabla \vec f(t) = D_2 ∇f(t)=D2, ⋯ \cdots ⋯, t = t n t=t_n t=tn 时 f ⃗ ( t ) = P n \vec f(t) = P_n f(t)=Pn 且 ∇ f ⃗ ( t ) = D n \nabla \vec f(t) = D_n ∇f(t)=Dn。分别对 [ t 1 , t 2 ] [t_1, t_2] [t1,t2], [ t 2 , t 3 ] [t_2, t_3] [t2,t3], ⋯ \cdots ⋯, [ t n − 1 , t n ] [t_{n-1}, t_n] [tn−1,tn] 建立一个规范化三次埃尔米特插值曲线后再用分段函数连接起来即可。
由于规范化三次埃尔米特插值曲线的自变量 t t t 的取值范围总是在 [ 0 , 1 ] [0,1] [0,1] 中,因此需要为每一段的参数进行换元。令 u i ( t ) = t − t i t i + 1 − t i u_i(t)=\frac{t-t_{i}}{t_{i+1}-t_i} ui(t)=ti+1−tit−ti 其中 i = 1 , 2 , ⋯ n − 1 i=1, 2, \cdots n-1 i=1,2,⋯n−1,则有当 t ∈ [ t i , t i + 1 ] t\in[t_i, t_{i+1}] t∈[ti,ti+1] 时, u i ( t ) u_i(t) ui(t) 从 0 0 0 线性地变化为 1 1 1,因此我们想要的分段函数即为:
f ⃗ ( t ) = { P 1 g 0 ( u 1 ( t ) ) + P 2 g 1 ( u 1 ( t ) ) + D 1 q 0 ( u 1 ( t ) ) + D 2 q 1 ( u 1 ( t ) ) , t ∈ [ t 1 , t 2 ) P 2 g 0 ( u 2 ( t ) ) + P 3 g 1 ( u 2 ( t ) ) + D 2 q 0 ( u 2 ( t ) ) + D 3 q 1 ( u 2 ( t ) ) , t ∈ [ t 2 , t 3 ) ⋯ P n − 1 g 0 ( u n − 1 ( t ) ) + P n g 1 ( u n − 1 ( t ) ) + D n − 1 q 0 ( u n − 1 ( t ) ) + D n q 1 ( u n − 1 ( t ) ) , t ∈ [ t n − 1 , t n ) \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)=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧P1g0(u1(t))+P2g1(u1(t))+D1q0(u1(t))+D2q1(u1(t)),P2g0(u2(t))+P3g1(u2(t))+D2q0(u2(t))+D3q1(u2(t)),⋯Pn−1g0(un−1(t))+Png1(un−1(t))+Dn−1q0(un−1(t))+Dnq1(un−1(t)),t∈[t1,t2)t∈[t2,t3)t∈[tn−1,tn)
在自变量 t t t 的变化过程中,这个分段函数每一段中的 u i u_i ui 都能从 0 0 0 开始到 1 1 1 结束——“规范化”想要表达的就是这重含义。表达式中的 g 0 , g 1 g_0, g_1 g0,g1 是导数无关的常数控制基函数, q 0 , q 1 q_0, q_1 q0,q1 是导数控制基函数。
具体而言,无论是描述空间中的曲线还是平面中的曲线,我们都可以使用这样的一种方式:先确定一系列的控制点(例如上文中的 P 0 , P 1 P_0, P_1 P0,P1),之后再将这些控制点与一组标量函数对应相乘,这组标量函数一般被称为基函数。
讲点儿时的回忆
贝塞尔曲线让我回忆起了我小时候经常玩的一个画图游戏:
- 在纸上画出两条登场的线段首位顺次相接,如下图中的 A B = B C AB=BC AB=BC 且使 A B AB AB 垂直于 B C BC BC;
- 在 A B AB AB 边上取一点 D D D,在 B C BC BC 边上取一点 E E E,使得 A D = B E AD=BE AD=BE;
- 在纸上借助直尺连接 D E DE DE;
- 随机更换 D D D 点的选取和 E E E 点的选取,但始终保持 A D = B E AD=BE AD=BE 不变,转步骤 3;
随着连接得到的线段越来越多,这些线段下方空白部分逐渐形成了一个曲线弧的形状(我上小学的时候画过不少这种东西,尽管那时的我不知道这个曲线是一条抛物线的一部分,一直以为它是一个圆弧)。
如何简单地说明这个曲线是一条抛物线呢,比较直观的方法就是先把焦点和准线找出来。对于这种绘制方式而言,想要证明绘制出的曲线是抛物线并不容易,但是我们可以反过来,去证明一个抛物线可以用这种方法进行绘制(虽然数学上这种证法并不严谨,但是我们做这个证明还有其他的用处)。
上图中红色的抛物线为
y
=
1
4
x
2
y=\frac 1 4 x^2
y=41x2,
M
M
M 是抛物线
A
C
AC
AC 段上的任意一点,过
M
M
M 做抛物线的切线
D
E
DE
DE,交
A
B
AB
AB、
B
C
BC
BC 于
D
,
E
D,E
D,E 两点,此时我们只需要证明
A
D
=
B
E
AD=BE
AD=BE 即可说明上文中的绘制方式能够绘制一条抛物线。
由于 y ′ = 1 2 x y'=\frac 1 2 x y′=21x,设 M M M 点横坐标为 x 0 x_0 x0,得到 M M M 点出的切线方程为 y = 1 2 x 0 ( x − x 0 ) + 1 4 x 0 2 y=\frac 1 2 x_0(x-x_0)+\frac 1 4 x_0^2 y=21x0(x−x0)+41x02,将切线方程分别与 A B AB AB 方程、 B C BC BC 方程联立可以解得 D , E D, E D,E 两点的坐标。最后我们可以验证 x D − x A = x E − x B x_D-x_A=x_E-x_B xD−xA=xE−xB。其中 x D = 1 2 x 0 − 1 , x E = 1 + 1 2 x 0 x_D=\frac 1 2 x_0-1, x_E=1+\frac 1 2 x_0 xD=21x0−1,xE=1+21x0 带入成立。
对此我们还可以验证得到 A D : D B = D M : M E AD:DB=DM:ME AD:DB=DM:ME,这也恰好为我们提供了一种二次贝塞尔曲线的绘制方式(几何作图法)。假设 A , B , C A, B, C A,B,C 是我们在平面上选取的三个点,我们根据这种比例绘制的方式可以得到 M M M 移动得到的的参数方程曲线。 设 t ∈ [ 0 , 1 ] t\in[0, 1] t∈[0,1],当 t = 0 t=0 t=0 时 M M M 与 A A A 点重合,当 t = 1 t=1 t=1 时 M M M 与 C C C 点重合。
首先我们能很容易的得到 D = ( 1 − t ) A + t B D=(1-t)A+tB D=(1−t)A+tB,换言之, t = 0 t=0 t=0 时 D D D 与 A A A 重合, t = 1 t=1 t=1 时 D D D 与 B B B 重合, t ∈ ( 0 , 1 ) t\in(0, 1) t∈(0,1) 时 D D D 在 A → B A\to B A→B 上做匀速运动。同理 E = ( 1 − t ) B + t C E=(1-t)B+tC E=(1−t)B+tC,而根据我们上文指出的比例关系,我们知道 D M : M E = A D : D B = t : ( 1 − t ) DM:ME=AD:DB=t:(1-t) DM:ME=AD:DB=t:(1−t),因此得到 M ( t ) = ( 1 − t ) D + t E = ( 1 − t ) 2 A + 2 t ( 1 − t ) B + t 2 C M(t)=(1-t)D+tE=(1-t)^2A+2t(1-t)B+t^2C M(t)=(1−t)D+tE=(1−t)2A+2t(1−t)B+t2C。当 t t t 取遍 [ 0 , 1 ] [0, 1] [0,1] 中的所有值时, M ( t ) M(t) M(t) 取遍抛物线在 A A A 到 C C C 之间的所有点。其中,我们称 ( 1 − t ) 2 (1-t)^2 (1−t)2, 2 t ( 1 − t ) 2t(1-t) 2t(1−t), t 2 t^2 t2 这三个函数为“二次伯恩斯坦基函数”。
三次贝塞尔曲线
不难发现其实“二次伯恩斯坦基函数”就是对
(
(
1
−
t
)
+
t
)
2
((1-t)+t)^2
((1−t)+t)2 进行牛顿二项式展开之后的每一项。类比地,我们可以得到四个三次伯恩斯坦基函数:
(
1
−
t
)
3
,
3
(
1
−
t
)
2
t
,
3
(
1
−
t
)
t
2
,
t
3
(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
,
D
A, B, C, D
A,B,C,D 时,我们可以确定一条三次贝塞尔函数:
M
(
t
)
=
(
1
−
t
)
3
A
+
3
(
1
−
t
)
2
t
B
+
3
(
1
−
t
)
t
2
C
+
t
3
D
,
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]
虽然有高次的贝塞尔函数,但是三次的贝塞尔曲线已经可以实现光滑光顺的拼接,所以我觉得不介绍高次贝塞尔曲线也罢。总而言之“
n
n
n 次的伯恩斯坦基函数”就是对
(
(
1
−
t
)
+
t
)
n
((1-t)+t)^n
((1−t)+t)n 进行牛顿二项式展开之后的每一项,其中的第
i
i
i 项为 (
i
=
0
,
1
,
⋯
,
n
i=0, 1, \cdots, n
i=0,1,⋯,n):
B
i
(
t
)
=
C
n
i
⋅
t
i
(
1
−
t
)
n
−
i
B_i(t)=C_n^i\cdot t^i(1-t)^{n-i}
Bi(t)=Cni⋅ti(1−t)n−i
我们将 n + 1 n+1 n+1 个 n n n 次伯恩斯坦基函数与 n n n 个控制点坐标对应相乘就得到了一条 n n n 次的贝塞尔曲线。美中不足的一点是,由于这 n + 1 n+1 n+1 个 n n n 次伯恩斯坦基函数的定义域都是全体实数,因此每当我们改变一个控制点的位置,整条曲线上的每一个点位置都会发生一定的改变,这不利于工业设计上的局部微调。
想象我们正在使用贝塞尔曲线绘制一位二次元少女的右脸,当我们觉得不满意时打算对这一段曲线进行微调,于是调整了右脸处的某个控制点,而这导致原本已经绘制得很不错的左脸发生了变形,这是一件多么悲伤的事。
三次贝塞尔曲线的光滑拼接
假设我们现在有两条贝塞尔曲线:
P
(
t
)
=
(
1
−
t
)
3
P
0
+
3
(
1
−
t
)
2
t
P
1
+
3
(
1
−
t
)
t
2
P
2
+
t
3
P
3
,
t
∈
[
0
,
1
]
Q
(
t
)
=
(
1
−
t
)
3
Q
0
+
3
(
1
−
t
)
2
t
Q
1
+
3
(
1
−
t
)
t
2
Q
2
+
t
3
Q
3
,
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]
其中 P 0 , P 1 , P 2 , P 3 , Q 0 , Q 1 , Q 2 , Q 3 P_0, P_1, P_2, P_3, Q_0, Q_1, Q_2, Q_3 P0,P1,P2,P3,Q0,Q1,Q2,Q3 是平面上的常数点。我们希望把曲线 Q Q Q 光滑地接在曲线 P P P 地后面,只需要保证:
- 连续: P ( 1 ) = Q ( 0 ) P(1)=Q(0) P(1)=Q(0);
- 梯度同向平行: ∇ 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
)
=
P
3
P(1)=P_3
P(1)=P3,
Q
(
0
)
=
Q
0
Q(0)=Q_0
Q(0)=Q0,所以想要保证连续只需要保证
P
3
=
Q
0
P_3=Q_0
P3=Q0。
∇
P
(
t
)
=
∇
(
(
1
−
t
)
(
(
1
−
t
)
2
P
0
+
3
(
1
−
t
)
t
P
1
+
3
t
2
P
2
)
+
t
3
P
3
)
=
−
(
(
1
−
t
)
2
P
0
+
3
(
1
−
t
)
t
P
1
+
3
t
2
P
2
)
+
(
1
−
t
)
(
⋯
)
+
3
t
2
P
3
\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
=
0
1-t=0
1−t=0 后也一定等于零),得到
∇
P
(
1
)
=
3
(
P
3
−
P
2
)
\nabla P(1)=3(P_3-P_2)
∇P(1)=3(P3−P2)。同理可以计算得到:
∇
Q
(
t
)
=
∇
(
(
1
−
t
)
3
Q
0
+
t
(
3
(
1
−
t
)
2
Q
1
+
3
(
1
−
t
)
t
Q
2
+
t
2
Q
3
)
)
=
−
3
(
1
−
t
)
2
Q
0
+
(
3
(
1
−
t
)
2
Q
1
+
3
(
1
−
t
)
t
Q
2
+
t
2
Q
3
)
+
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 ( Q 1 − Q 0 ) \nabla Q(0)=3(Q_1-Q_0) ∇Q(0)=3(Q1−Q0)。因此只要在 P 3 = Q 0 P_3=Q_0 P3=Q0 的基础上保证 P 2 P 3 ⟶ / / Q 0 Q 1 ⟶ \stackrel \longrightarrow {P_2P_3}//\stackrel \longrightarrow{Q_0Q_1} P2P3⟶//Q0Q1⟶ 即可保证(还要注意同向哦)。
三次伯恩斯坦基函数
我们看到三次的贝塞尔曲线有四个控制点:
M
(
t
)
=
(
1
−
t
)
3
A
+
3
(
1
−
t
)
2
t
B
+
3
(
1
−
t
)
t
2
C
+
t
3
D
,
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]
每当给定一个
t
t
t 时,曲线上的点
M
(
t
)
M(t)
M(t) 实际上就是对
A
,
B
,
C
,
D
A, B, C, D
A,B,C,D 四个控制点的坐标进行了一个加权求平均数的过程。当
t
→
0
t\to 0
t→0 时,
M
(
t
)
→
A
M(t)\to A
M(t)→A,当
t
→
1
t\to 1
t→1 时
M
(
t
)
→
D
M(t)\to D
M(t)→D。具体而言我们可以将
A
,
B
,
C
,
D
A, B, C, D
A,B,C,D 四个点分别对应的基函数绘制到一个坐标系中,让我们观察这四个点对整条曲线究竟起到了多大的作用。
可以看到,当
t
=
1
3
t=\frac 1 3
t=31 时,
B
B
B 前的系数变得最大,此时曲线最靠近
B
B
B 点;当
t
=
2
3
t=\frac 2 3
t=32 时,
C
C
C 前的系数变得最大,此时曲线最靠近
C
C
C 点。
换言之,这也为我们进一步研究曲线的表示提供了一个思路:给定 n n n 个控制点,我们只需要给出 n n n 个标量函数作为基函数和这 n n n 个点的坐标对应相乘,就能得到一条大趋势与这 n n n 个点的走势基本接近的曲线。这 n n n 个标量函数都是上凸单峰的,使得曲线在当前标量函数的极值点处最接近该控制点,从而使得曲线近似拟合了这 n n n 个控制点构成的折线。
由于 CSDN 对字数的限制,对贝塞尔曲线的绘制将在另一篇博客中进行讲解。
GGN_2015
计算机图形学中的曲线问题——贝塞尔曲线的绘制
三次均匀 B 样条曲线
由于三次贝塞尔曲线只能有四个控制点,贝塞尔曲线控制点越多次数越高。我们希望构造一种这样的曲线:
- 控制点个数可以无限增多,但函数次数保持为三次;
- 每个控制点只会影响自己附近的一小段曲线而对其他部分曲线没有任何影响;
- 曲线是光滑的。
局部基函数及其混淆
我们的目的非常明确,我们想要构造一种基函数,这种基函数只在很小的一段区间上有定义,而在此区间之外的部分均为零。例如这样的一个函数就是我们想要的一个函数:
N
1
(
t
)
=
{
1
,
0
≤
t
≤
1
0
,
elsewhere
N_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],因此我们称这个局部基函数的宽度为 1 1 1。类似的,如果一个定义在实数上的函数只在区间 [ 0 , L ] [0, L] [0,L] 上有非零值,我们可以称这个函数的“宽度”为 L L L 的局部基函数。
当我们用一组局部基函数与一组控制点对应相乘再求和从而得到目标曲线时,每一个控制点自然也只会影响曲线的一个局部(在局部基函数的宽度范围内会产生影响),而不会影响整个曲线上的绝大多数位置的取值。
假如我们已经有了一个宽度为 L L L 的局部基函数,我们可以使用一次混淆操作得到一个宽度为 L + 1 L+1 L+1 的局部基函数,这使得我们可以控制每个控制点在目标曲线上到底能够造成多大范围的影响,宽度越大影响范围自然也越大。
对于平台函数
N
1
N_1
N1 进行混淆,第一步我们要将
N
1
(
t
)
N_1(t)
N1(t) 右移一个单位长度得到
N
1
(
t
−
1
)
N_1(t-1)
N1(t−1),我们希望能够将
N
1
(
t
)
N_1(t)
N1(t) 与
N
1
(
t
−
1
)
N_1(t-1)
N1(t−1) 加权求平均,使得我们得到的新函数在
t
=
0
→
1
t=0\to1
t=0→1 段上函数值线性增加,在
t
=
1
→
2
t=1\to2
t=1→2 段上函数值线性减少。很直观的想法就是,我们可以在
N
1
(
t
)
N_1(t)
N1(t) 身上乘以一个单调递增的线性函数,并在
N
1
(
t
−
1
)
N_1(t-1)
N1(t−1) 身上乘以一个单调递减的线性函数:
N
2
(
t
)
=
t
⋅
N
1
(
t
)
+
(
2
−
t
)
⋅
N
1
(
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)
为什么我们在 N 1 ( t − 1 ) N_1(t-1) N1(t−1) 身上乘以的是一个 2 − t 2-t 2−t 而不是 1 − t 1-t 1−t 呢,这是因为 N 1 ( t − 1 ) N_1(t-1) N1(t−1) 的非零区间为 [ 1 , 2 ] [1, 2] [1,2],如果乘以 1 − t 1-t 1−t 会导致我们最终的结果出现负数,下图展示了 N 1 ( t ) , N 1 ( t − 1 ) N_1(t), N_1(t-1) N1(t),N1(t−1) 和 N 2 ( t ) N_2(t) N2(t):
其中绿色的虚线表示
N
1
(
t
)
N_1(t)
N1(t),蓝色的虚线代表
N
1
(
t
−
1
)
N_1(t-1)
N1(t−1),红色的虚线代表
N
2
(
t
)
N_2(t)
N2(t)。不难看出
N
2
(
t
)
N_2(t)
N2(t) 是一个宽度为
2
2
2 的局部基函数。这样得到的
N
2
(
t
)
N_2(t)
N2(t) 已经是一个最高次为
1
1
1 次的单峰函数,使用这种函数作为基函数,我们可以得到一段连续的折线。例如,我们可以令:
M
2
(
t
)
=
A
⋅
N
2
(
t
)
+
B
⋅
N
2
(
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]
这样得到的
M
2
(
t
)
M_2(t)
M2(t) 就是一个连接了
A
A
A 点和
B
B
B 点的线段。我们之所以要限制
t
∈
[
1
,
2
]
t\in[1, 2]
t∈[1,2] 是因为在
t
∈
[
0
,
1
]
t\in [0, 1]
t∈[0,1] 时只有
N
1
(
t
)
N_1(t)
N1(t) 非零,
N
1
(
t
−
1
)
N_1(t-1)
N1(t−1) 为零,我们潜在地在对
A
A
A 点和原点进行加权求平均,而我们显然不希望原点影响到我们的图像的形状。假如我们将定义域定义成
t
∈
[
0
,
3
]
t\in[0, 3]
t∈[0,3],会导致
M
2
(
t
)
M_2(t)
M2(t) 的图象变为了三段线段,第一段连接了原点和
A
A
A 点,第二段连接了
A
A
A 点和
B
B
B 点,第三段连接了
B
B
B 点和原点。我们还可以令:
M
3
(
t
)
=
A
⋅
N
2
(
t
)
+
B
⋅
N
2
(
t
−
1
)
+
C
⋅
N
2
(
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]
这样我们得到的
M
3
(
t
)
M_3(t)
M3(t) 就是两段线段的拼接,其中第一段连接了
A
A
A 点与
B
B
B 点,第二段连接了
B
B
B 点与
C
C
C 点。为什么能实现这个效果呢?是因为
N
2
(
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] 上
M
3
(
t
)
=
M
2
(
t
)
M_3(t)=M_2(t)
M3(t)=M2(t),因此为一段连接了
A
A
A 和
B
B
B 的线段。
更普遍地,设
N
L
N_L
NL 是一个宽度为
L
L
L 的局部基函数,对函数
N
L
N_L
NL 的混淆记为
M
i
x
[
N
L
]
Mix[N_L]
Mix[NL]:
N
L
+
1
(
t
)
=
M
i
x
[
N
L
]
(
t
)
=
t
L
⋅
N
L
(
t
)
+
(
L
+
1
)
−
t
(
L
+
1
)
−
1
⋅
N
L
(
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)
简而言之就是在 N L ( t ) N_L(t) NL(t) 身上乘以了一个他的非零区间内的线性递增函数,在 N L ( t − 1 ) N_L(t-1) NL(t−1) 身上乘以了一个他的非零区间内的线性递减函数。下图给出了 N 2 ( t ) , N 2 ( t − 1 ) N_2(t), N_2(t-1) N2(t),N2(t−1) 以及 N 3 ( t ) N_3(t) N3(t) 的函数图像:
其中绿色的曲线是
N
2
(
t
)
N_2(t)
N2(t),蓝色的曲线是
N
2
(
t
−
1
)
N_2(t-1)
N2(t−1),红色的曲线是
N
3
(
t
)
N_3(t)
N3(t)。使用
N
3
(
t
)
N_3(t)
N3(t) 作为基函数已经可以较好的定义出一条连续且光滑的曲线了,例如我们令:
M 3 ( t ) = A ⋅ N 3 ( t ) + B ⋅ N 3 ( t − 1 ) + C ⋅ N 3 ( t − 2 ) + D ⋅ N 3 ( 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] 也是为了避免原点位置影响曲线的绘制。而我们这里用来做基函数的 N 1 ( t ) , N 2 ( t ) , ⋯ N_1(t), N_2(t), \cdots N1(t),N2(t),⋯ 其实就是所谓的 B 样条基函数。具体而言,对于一个给定一个正整数 L L L,我们称 N L ( t ) N_L(t) NL(t) 为 L L L 阶( L − 1 L-1 L−1 次)B 样条基函数。
四阶三次均匀 B 样条曲线
对
N
3
(
t
)
N_3(t)
N3(t) 进行混淆可以得到
N
4
(
t
)
N_4(t)
N4(t),其中
N
4
(
t
)
N_4(t)
N4(t) 是一个宽度为
4
4
4 的局部基函数,我们称它为四阶三次 B 样条基函数,下图展示了它的大致形状:
二阶 B 样条至少要有两个控制点(不然线段怎么连),三阶 B 样条至少要有三个控制点(不然二次曲线怎么画),因此四阶 B 样条至少要有四个控制点。B 样条相对于贝塞尔曲线的有点在于,即使四阶三次 B 样条有多于四个控制点,它的次数仍然是三次,而且每个控制点的变化只会影响与他接近的一段曲线,而不会影响整个目标曲线(因为
N
4
(
t
)
N_4(t)
N4(t) 宽度为
4
4
4,这也就意味着与
N
4
(
t
)
N_4(t)
N4(t) 相乘的那个控制点不会在这个宽度之外造成任何影响)。
假设我们有
n
n
n 个常数点
P
0
,
P
1
,
⋯
,
P
n
−
1
P_0, P_1, \cdots, P_{n-1}
P0,P1,⋯,Pn−1,其中
n
≥
4
n\geq 4
n≥4,我们可以定义一条四阶三次均匀 B 样条曲线:
M
(
t
)
=
P
0
⋅
N
4
(
t
−
0
)
+
P
1
⋅
N
4
(
t
−
1
)
+
⋯
+
P
n
−
1
⋅
N
4
(
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 样条往往不会经过自己的第一个控制点 P 0 P_0 P0 和 P n − 1 P_{n-1} Pn−1,而贝塞尔曲线一定会经过这两个点。
B 样条绘制范例
随便画了一个有七个控制点的 B 样条,特别地,如果我们让
P
2
=
P
6
,
P
1
=
P
5
,
P
0
=
P
4
P_2=P_6, P_1=P_5, P_0=P_4
P2=P6,P1=P5,P0=P4,那么我们能够绘制一条闭合的 B 样条曲线。
换言之,任给一个多边形,我们都能找到这个多边形控制的一条闭合 B 样条曲线。如果这个多边形是一个凸多边形,感性理解一下不能得到这个闭合的 B 样条曲线一定完全位于整个多边形的内部。
上文中的所有图形都是使用几何画板绘制的感兴趣的同学可以联系我(邮箱见免责)获得这些几何画板文件,下图中给出了用于生成这两个 B 样条范例的配置:
未完待续 …
一个理念
不得不说我们的教材内容丰富充实,兼顾了严谨性与知识性。但我对我们的图形学教材仍有一个略有不满的地方,就是教材中没有贯彻一个目的性与实践性的假设:任何工业设计理论都带有其目的。当我们提出一个概念/公式时,指导我们的所谓“灵感”不是拍脑袋想出来的,更不是从天上掉下来的,而是来自真真正正的生产实践的需要。我们的图形学教材并不是没有这种理念,但是它没有贯彻这种理念,换言之,这种理念仅仅停留在每一章的简介部分,而没有真正深入到每一章的具体内容中。不要先给出一个复杂的公式,再利用这个公式解决问题,要先从最简单、最必要、最直观的概念入手逐渐让我们知道我们为什么要这么做(目的),而不是先告诉我们怎么做(方法)。否则学生很容易不知所云。