接上文,本篇文章讲解GJK算法论文的第四、第五部分前半部分,这是整个算法最为核心的部分。第四部分阐述了算法的核心思想,而第五部分则给出了可行的算法流程。不废话,直接进入正题。
The Theoretical Algorithm
内容详解
这一节的主要内容就是介绍寻找
v
(
K
)
v(K)
v(K)的算法,其中
K
K
K是凸的并且是紧的。如果
K
K
K是多面体,那么这个算法能在有限次迭代后找到
v
(
K
)
v(K)
v(K)。
这个算法的主要思想如下:在
K
K
K中生成一系列的多面体
c
o
V
k
co\space V_k
co Vk,并且序列
v
(
c
o
V
k
)
v(co\space V_k)
v(co Vk)会收敛于
v
(
K
)
v(K)
v(K)。为了表述这个算法,需要引进以下定理:
定理1:设
K
⊂
R
m
K \subset R^m
K⊂Rm,
K
K
K是闭集且是凸的,然后定义函数
g
K
:
R
m
→
R
g_K:R^m\to R
gK:Rm→R为:
g
K
(
x
)
=
∣
x
∣
2
+
h
K
(
−
x
)
,
x
∈
K
(
25
)
g_K(x) = |x|^2 + h_K(-x), x \in K \kern3em \lparen25\rparen
gK(x)=∣x∣2+hK(−x),x∈K(25)那么有以下三个结论:
1)若
g
K
(
x
)
>
0
g_K(x) \gt 0
gK(x)>0,那么在线段
{
x
,
s
K
(
−
x
)
}
\lbrace x,s_K(-x) \rbrace
{x,sK(−x)}间必定存在一个点
z
z
z,使得
z
z
z满足
∣
z
∣
<
∣
x
∣
|z| \lt |x|
∣z∣<∣x∣;
2)只有当
g
K
(
x
)
=
0
g_K(x) = 0
gK(x)=0时,
x
=
v
(
K
)
x = v(K)
x=v(K);
3)
∣
x
−
v
(
K
)
∣
2
⩽
g
K
(
x
)
|x - v(K)|^2 \leqslant g_K(x)
∣x−v(K)∣2⩽gK(x)
为了不把思路带偏,定理1的证明这里就不再陈述(因为不了解对理解整个算法没有太大的区别,个人观点),其证明的主要思想是三角不等式,论文的后面也有证明的过程。
接下来阐述距离算法:
1)设
V
0
=
{
y
1
,
y
2
.
.
.
y
v
}
V_0 = \lbrace y_1,y_2...y_v \rbrace
V0={y1,y2...yv},其中
y
1
,
y
2
.
.
.
y
v
∈
K
y_1,y_2...y_v \in K
y1,y2...yv∈K,并且
1
⩽
v
⩽
m
+
1
1 \leqslant v \leqslant m+1
1⩽v⩽m+1,这一步就是说一开始要在集合中随机选取v个点,当然这个随机是可以伴随一定的策略的。这里再次提醒,
K
K
K必须是多面体,如果是像球这样的物体,这个算法是求不了的,这是因为球会让这个算法无限迭代下去(不过mesh的话也不会有真正意义上的球)。
2)求出
v
k
=
v
(
c
o
V
k
)
v_k = v(co \space V_k)
vk=v(co Vk),就是求出$
y
1
,
y
2
.
.
.
y
v
y_1,y_2...y_v
y1,y2...yv这些点构成的凸包到原点的最近的点是什么,而这个点代表的向量会代入到(15)式中进行计算。这个怎么找呢,论文的下一节才给出了算法,所以这里先不提。
3)如果
g
K
(
v
k
)
=
0
g_K(v_k) = 0
gK(vk)=0(要注意里面的大小写,大写代表集合,小写代表标号),那么这个点就是我们要找的点,这个时候就可以退出循环。否则进入下一步
4)设
V
k
+
1
=
V
k
^
⋃
{
s
K
(
−
v
k
)
}
V_{k+1} = \hat{V_k} \bigcup \lbrace s_K(-v_k)\rbrace
Vk+1=Vk^⋃{sK(−vk)},其中
V
k
^
\hat{V_k}
Vk^指的是当
v
k
∈
c
o
V
k
^
v_k \in co \space \hat{V_k}
vk∈co Vk^时,
V
k
^
\hat{V_k}
Vk^中的顶点数最少,举个例子,在二维平面下一个三角形到原点最近的点有三种情况,一是原点落在三角形内部,那么原点需要用三个顶点来进行表示((12)式);二是最近的点落在三角形边上,那么这个点只需要两个顶点就可以表示;还有就是离原点最近的点就是其中一个顶点,那么这个点只用一个点就可以进行表示。然后后面的
{
s
K
(
−
v
k
)
}
\lbrace s_K(-v_k)\rbrace
{sK(−vk)}表示选择一个尽量靠近原点(向量的长度短),并且离
v
k
v_k
vk远的点A,点A必定会更靠近原点。然后将点A和
V
k
^
\hat{V_k}
Vk^的并集输入到2)中继续执行。
结合上面对算法的阐述以及论文的Fig.3,这个算法应该是比较清晰的了,而且论文中也有计算过程的描述,所以这里就不多说明了。下面是论文中的图,这里的
V
0
=
{
z
1
,
z
2
,
z
3
}
V_0=\lbrace z_1,z_2,z_3 \rbrace
V0={z1,z2,z3}
接下来是收敛性证明,因为会把思路带偏,这里也不再做过多阐述,有兴趣的可以自己去看
初探The Distance Subalgorithm
内容详解
这一节的内容是对上一节的距离算法的一些细节进行补充,主要是步骤2),这一个小节的论述是在是太精彩了,喜欢数学以及图形的你绝对不能错过。
首先还是设定条件,设
Y
=
{
y
1
,
y
2
.
.
.
y
v
}
⊂
R
m
Y = \lbrace y_1,y_2...y_v \rbrace \subset R^m
Y={y1,y2...yv}⊂Rm,即
Y
=
V
k
Y = V_k
Y=Vk,这个
V
k
V_k
Vk就是上一节提到的集合。那么从这个算法最终可以将
v
(
c
o
Y
)
v(co\space Y)
v(co Y)表示成如下形式:
v
(
c
o
Y
)
=
∑
i
∈
I
S
λ
i
y
i
v(co\space Y) = \sum_{i \in I_S} \lambda^i y_i
v(co Y)=i∈IS∑λiyi
∑
i
∈
I
S
λ
i
=
1
λ
i
>
0
i
∈
I
S
⊂
{
1
,
2
,
.
.
.
,
v
}
\sum_{i \in I_S} \lambda^i = 1 \lambda^i > 0\space i \in I_S \subset \lbrace 1, 2,..., v \rbrace
i∈IS∑λi=1λi>0 i∈IS⊂{1,2,...,v}
Y
S
=
{
y
i
,
i
∈
I
S
}
i
s
a
f
f
i
n
e
l
y
i
n
d
e
p
e
n
d
e
n
t
(
29
)
Y_S = \lbrace y_i, i \in I_S \rbrace \space is \space affinely \space independent \kern3em \lparen29\rparen
YS={yi,i∈IS} is affinely independent(29)其中
S
S
S表示集合
Y
Y
Y的所有非空索引集中的其中一个。可以发现,上面的式子其实就是公式(14)。这里需要注意的是
Y
S
Y_S
YS就是距离算法中4)的
V
k
^
\hat{V_k}
Vk^,之前提到说
V
k
^
\hat{V_k}
Vk^是包含最少顶点的集合,这样能提高迭代的效率。
当
v
v
v的值很小的时候,即便我们对所有顶点的所有组合进行最近点的计算,效率也是很高的,要计算的次数可以通过组合的公式得到:
σ
=
∑
k
=
1
v
[
v
!
/
(
k
!
(
v
−
k
)
!
)
]
\sigma = \sum^v_{k = 1}[v! / (k!(v - k)!)]
σ=k=1∑v[v!/(k!(v−k)!)]例如当
v
=
4
v = 4
v=4时,那么
σ
=
15
\sigma = 15
σ=15。回到距离算法,基于(12)式,可以知道,在三维空间中,
V
k
^
\hat{V_k}
Vk^最多拥有4个点。如果还是不太清楚为什么最多是4,我这里展开解释一下。可能会有人觉得,5个点不可以吗?我们假设可以,那么要表示凸包中的点A的时候,我们会用5个点来表示这个点A,可是在三维空间里面,我们真的需要用5个点来表示一个点吗?因为用4个点的时候我们已经能构成一个四面体了,这个四面体已经能够包含空间内的点了,再多一个是可以,但是没必要,因为这样一来我们要处理的数据就更多了,而且得到的结果说不定还更复杂,而
V
k
^
\hat{V_k}
Vk^仅仅需要最少的点(再次强调)。然后我们提到
Y
S
=
V
k
^
Y_S = \hat{V_k}
YS=Vk^,所以下文中的
Y
i
Y_i
Yi就代表
V
i
^
\hat{V_i}
Vi^.
(30)式使用了一堆奇怪的符号,定理3也使用了一堆奇怪的符号,一开始看得我晕头转向,然后发现在附录二中有详细的说明,所以我们先跳到附录二,然后再回来看定理3
Appendix Ⅱ
涉及的概念
在附录2中将会使用到下面几个概念,需要对其有较为清晰的认识才能保证后面不会晕,有些概念不是三言两语能解释清楚的,所以原谅我直接放上连链接。
- 代数余子式(cofactor):百度百科
- 克莱姆法则(Cramer’s Rule):百度百科
- 多元函数的极值与局部极值:多元函数的极值指的是当 x = x 0 x = x_0 x=x0时,对于 ∀ x ′ ∈ E \forall x' \in E ∀x′∈E都有 f ( x ′ ) ⩽ f ( x 0 ) f(x') \leqslant f(x_0) f(x′)⩽f(x0),那么 f ( x 0 ) f(x_0) f(x0)就是函数 f f f在集合 E E E中的极值。至于局部极值,简单来说,就是在这个点 x l x_l xl附近, f ( x l ) f(x_l) f(xl)最大或最小,然后这个“附近”不一定是全局,它可以是半径为1的球内(三维中),也可以是半径为10的球内,只要在某个半径内满足即可。如果点 f ( x l ) f(x_l) f(xl)是局部极值,且对 x l x_l xl内每个变量的偏导数都存在( x l x_l xl是向量),那么这些偏导数都为0。我这里的描述并不严瑾,主要是要展开的话又会涉及一系列的概念。这里记住一点就是局部极值的偏导数为0(当然前提是存在,不过针对论文描述的问题,偏导数必定存在)就够了。
内容详解
附录2主要内容是要证明定理3,但现在我们抛开这个,把附录2当成是独立的篇章,那么附录2的主要内容就是确定(29)式中
λ
i
\lambda^i
λi的值
首先,我们考虑如何确定
v
(
a
f
f
Y
S
)
v(aff\space Y_S)
v(aff YS)的值,其中
Y
S
∈
Y
Y_S \in Y
YS∈Y。如果
Y
S
Y_S
YS中只有1个点,那么直接求距离即可,没有太多讨论的必要,所以这里假设
Y
S
Y_S
YS至少含有2个点,且后面提到的
r
r
r也至少为2。
v
(
a
f
f
Y
S
)
v(aff\space Y_S)
v(aff YS)的形式如下:
v
(
a
f
f
Y
S
)
=
∑
i
=
2
r
λ
i
x
i
(
a
)
v(aff\space Y_S) = \sum^r_{i = 2} \lambda^i x_i \kern3em \lparen a\rparen
v(aff YS)=i=2∑rλixi(a)根据仿射集的定义有:
λ
1
=
1
−
∑
i
=
2
r
λ
i
(
b
)
\lambda^1 = 1 - \sum^r_{i = 2} \lambda^i \kern3em \lparen b\rparen
λ1=1−i=2∑rλi(b)那么点
v
(
a
f
f
Y
S
)
v(aff\space Y_S)
v(aff YS)到原点的距离可以表示成
λ
2
,
λ
3
.
.
.
λ
r
\lambda^2, \lambda^3...\lambda^r
λ2,λ3...λr的函数:
f
(
λ
2
,
λ
3
.
.
.
λ
r
)
=
∣
x
1
+
∑
i
=
2
r
λ
i
(
x
i
−
x
1
)
∣
2
(
c
)
f(\lambda^2, \lambda^3...\lambda^r) = |x_1 + \sum^r_{i = 2}\lambda^i(x_i - x_1)|^2 \kern3em \lparen c\rparen
f(λ2,λ3...λr)=∣x1+i=2∑rλi(xi−x1)∣2(c)上面那个式子可以把求和的项拆开,就会得到一般的距离公式。然后关键点来了,
v
(
a
f
f
Y
S
)
v(aff\space Y_S)
v(aff YS)是离原点最近的点,那代表什么?没错,就是
f
(
λ
2
,
λ
3
.
.
.
λ
r
)
f(\lambda^2, \lambda^3...\lambda^r)
f(λ2,λ3...λr)取最小值,之前我们提到了局部最值的概念,这里就正好可以用上。首先,我们要求
f
(
λ
2
,
λ
3
.
.
.
λ
r
)
f(\lambda^2, \lambda^3...\lambda^r)
f(λ2,λ3...λr)的偏导数,简单起见,这里设r = 3,那么对©式求偏导数可得:
∂
f
∂
λ
2
=
2
λ
2
(
x
2
−
x
1
)
⋅
(
x
1
+
λ
2
(
x
2
−
x
1
)
+
λ
3
(
x
3
−
x
1
)
)
=
λ
2
(
x
2
−
x
1
)
⋅
(
λ
1
x
1
+
λ
2
x
2
+
λ
3
x
3
)
=
0
\frac{\partial f}{\partial \lambda^2} = 2\lambda^2(x_2 - x_1)·(x_1 + \lambda^2(x_2 - x_1) + \lambda^3(x_3 - x_1)) = \lambda^2(x_2 - x_1)·(\lambda^1x_1 + \lambda^2x_2 + \lambda^3x_3) = 0
∂λ2∂f=2λ2(x2−x1)⋅(x1+λ2(x2−x1)+λ3(x3−x1))=λ2(x2−x1)⋅(λ1x1+λ2x2+λ3x3)=0
∂
f
∂
λ
3
=
2
λ
3
(
x
3
−
x
1
)
⋅
(
x
1
+
λ
2
(
x
2
−
x
1
)
+
λ
3
(
x
3
−
x
1
)
)
=
λ
3
(
x
3
−
x
1
)
⋅
(
λ
1
x
1
+
λ
2
x
2
+
λ
3
x
3
)
=
0
\frac{\partial f}{\partial \lambda^3} = 2\lambda^3(x_3 - x_1)·(x_1 + \lambda^2(x_2 - x_1) + \lambda^3(x_3 - x_1)) = \lambda^3(x_3 - x_1)·(\lambda^1x_1 + \lambda^2x_2 + \lambda^3x_3) = 0
∂λ3∂f=2λ3(x3−x1)⋅(x1+λ2(x2−x1)+λ3(x3−x1))=λ3(x3−x1)⋅(λ1x1+λ2x2+λ3x3)=0然后结合(b)式,再将这几个式子写成矩阵相乘的形式,那么我们可以得到:
∣
1
1
1
(
x
2
−
x
1
)
⋅
x
1
(
x
2
−
x
1
)
⋅
x
2
(
x
2
−
x
1
)
⋅
x
3
(
x
3
−
x
1
)
⋅
x
1
(
x
3
−
x
1
)
⋅
x
2
(
x
3
−
x
1
)
⋅
x
3
∣
∣
λ
1
λ
2
λ
3
∣
=
∣
1
0
0
∣
(
d
)
\begin{vmatrix} 1 & 1 & 1 \\ (x_2 - x_1)·x_1 & (x_2 - x_1)·x_2 & (x_2 - x_1)·x_3 \\ (x_3 - x_1)·x_1 & (x_3 - x_1)·x_2 & (x_3 - x_1)·x_3 \end{vmatrix} \begin{vmatrix} \lambda^1 \\ \lambda^2 \\ \lambda^3 \end{vmatrix} = \begin{vmatrix} 1 \\ 0 \\ 0 \end{vmatrix} \kern3em \lparen d\rparen
∣∣∣∣∣∣1(x2−x1)⋅x1(x3−x1)⋅x11(x2−x1)⋅x2(x3−x1)⋅x21(x2−x1)⋅x3(x3−x1)⋅x3∣∣∣∣∣∣∣∣∣∣∣∣λ1λ2λ3∣∣∣∣∣∣=∣∣∣∣∣∣100∣∣∣∣∣∣(d)这就是论文中的(41)式(系数2约掉了)!接下来我们要做的事情就是通过这个式子求出
λ
1
,
λ
2
,
λ
3
\lambda^1, \lambda^2 ,\lambda^3
λ1,λ2,λ3。求解方程组有挺多方法的,这里用的是克莱姆法则(Cramer’s Rule),因为向量(1, 0, 0)很特殊,在用余子式求解行列式的时候,除了第一行向量(1, 0, 0)所在的那一列的余子式,其余余子式的行列式都为0,这样我们就知道了实际上符号
Δ
i
(
Y
S
)
\Delta_i(Y_S)
Δi(YS)代表的是克莱姆法则中向量(1, 0, 0)的1所在的行列的代数余子式, 而
Δ
(
Y
S
)
\Delta(Y_S)
Δ(YS)代表的是系数矩阵的行列式。拿(d)式做例子,求
λ
2
\lambda^2
λ2,那么根据克莱姆法则有:
Δ
2
(
Y
S
)
=
∣
1
1
1
(
x
2
−
x
1
)
⋅
x
1
0
(
x
2
−
x
1
)
⋅
x
3
(
x
3
−
x
1
)
⋅
x
1
0
(
x
3
−
x
1
)
⋅
x
3
∣
=
(
−
1
)
1
+
2
∣
(
x
2
−
x
1
)
⋅
x
1
(
x
2
−
x
1
)
⋅
x
3
(
x
3
−
x
1
)
⋅
x
1
(
x
3
−
x
1
)
⋅
x
3
∣
\Delta_2(Y_S) = \begin{vmatrix} 1 & 1 & 1 \\ (x_2 - x_1)·x_1 & 0 & (x_2 - x_1)·x_3 \\ (x_3 - x_1)·x_1 & 0 & (x_3 - x_1)·x_3 \end{vmatrix} = (-1)^{1 + 2}\begin{vmatrix} (x_2 - x_1)·x_1 & (x_2 - x_1)·x_3 \\ (x_3 - x_1)·x_1 & (x_3 - x_1)·x_3 \end{vmatrix}
Δ2(YS)=∣∣∣∣∣∣1(x2−x1)⋅x1(x3−x1)⋅x11001(x2−x1)⋅x3(x3−x1)⋅x3∣∣∣∣∣∣=(−1)1+2∣∣∣∣(x2−x1)⋅x1(x3−x1)⋅x1(x2−x1)⋅x3(x3−x1)⋅x3∣∣∣∣那么
λ
2
=
Δ
2
(
Y
S
)
/
Δ
(
Y
S
)
\lambda^2 = \Delta_2(Y_S) / \Delta(Y_S)
λ2=Δ2(YS)/Δ(YS)。至此,我们给出了仿射集系数的解的形式。
下一篇文章我们会继续证明定理3,并且会回到第五章,继续论文里的算法,敬请期待。