浏览量:591

9.3. 包围球

本篇文章禁止用于任何商业目的,版权申明、版本说明等见《前言》。

PDF文档和源码下载地址:https://github.com/twinklingstar20/Programmers_Computational_Geometry

9.3. 包围球

9.3.1. 最小包围球的创建

问题描述:给定一个点集\(S = \left\{ {{P_0},{P_1}, \cdots ,{P_{n - 1}}} \right\}\),求包围该点集的最小包围球。

早在1982年,Megiddo[9]提出了第一个在三维空间下求最小包围球的线性算法,但是该算法实现难度较大。

Ritter[7](1990)提出了一种估计算法,该算法非常简单,但并不能得到最小包围球。主要思想是从点集\(S\)中,找到\(x\)值、\(y\)值、\(z\)值最小和最大的点,从这三组点对中找出距离最大的一组点对,取这组点对的中心作为球体的中心\(c\),取点对距离的一半作为球的半径\(r\)。接着,再次遍历点集\(S\),检查点\({{p}_{k}}\)与中心的距离\(d\),如果\(d\succ r\),说明点在球外,将半径更新为\(r'=(d+r)/2\),并将球心向该点移动\((d-r)/2\),即新的球心为\(c'=c+({{p}_{k}}-c)\cdot (d-r)/(2d)\),重复这个过程,直至点集\(S\)中所有的点都包含在球内。

Welzl[10](1991)提出了一种通过使用随机化的方法来解决问题,算法复杂度的期望是线性的,且算法相对简单、易实现。该算法不仅能得到二维、三维的最小包围球,而且能扩展到更高维的情况,通过“重要点前移”的优化步骤,能显著的提高算法的效率。Gärtner[11,12]详细地描述了该算法的实现,并提供了C++的代码实现,在CGAL[13]也有该算法的实现,本节将介绍该算法。

以三维空间下的最小包围球为例,定义\(\text{MB}\left( S \right)\)表示点集\(S\)的最小包围球,显然,\(\text{MB}\left( S \right)\)是唯一的。设\({{D}_{1}}\)\({{D}_{2}}\)是点集\(S\)的两个最小包围球,半径都是\(r\),球心分别是\({{c}_{1}}\)\({{c}_{2}}\)。有\(S\subset {{D}_{1}}\)\(S\subset {{D}_{2}}\),所以有\(S\subset {{D}_{1}}\bigcap {{D}_{2}}\)\({{D}_{1}}\bigcap {{D}_{2}}\)被包含在以\(({{c}_{1}}+{{c}_{2}})/2\)为球心,\(\sqrt{{{r}^{2}}-{{a}^{2}}}\)为半径的球内,其中\(a\)表示两个球心距离的一半。若\(a\ne 0\),则存在半径更小的、包含点集\(S\)的球,与假设条件矛盾,因此\(a=0\)\({{D}_{1}}\)\({{D}_{2}}\)是重合的。

\(P,R\)是点集\(S\)的两个子集,且满足\(P \cup R = S,P \cap R = \emptyset \),定义\(\text{MB}\left( P,R \right)\)表示以\(R\)中的点为球的边界且包含\(P\)的最小包围球,则有\({\rm{MB}}\left( S \right){\rm{ = MB}}\left( {S,\emptyset } \right)\),定义\(\overline {{\rm{MB}}} \left( R \right) = {\rm{MB}}\left( {\emptyset ,R} \right)\),即表示点集\(R\)可以确定的最小包围球。接下来,阐述三个与最小包围球相关的重要性质:

  • i.     若\(\text{MB}\left( P,R \right)\)存在,那么它是唯一的;
  • ii.    对于点\(p\in P\),若\(p\notin \text{MB}\left( P-\left\{ p \right\},R \right)\),说明\(p\)\(\text{MB}\left( P-\left\{ p \right\},R \right)\)的边界上,则有\(\text{MB}\left( P,R \right)=\text{MB}\left( P-\left\{ p \right\},R\bigcup \left\{ p \right\} \right)\)
  • iii.   若\(\text{MB}\left( P,R \right)\)存在,则\(P\)中最多存在\(\max \left\{ 0,4-\left| R \right| \right\}\)个点的子集\(B\),使得\(\text{MB}\left( P,R \right)=\)\(\text{MB}\left( B,R \right)=\overline{\text{MB}}\left( B\bigcup R \right)\),其中\(\left| R \right|\)表示点集\(R\)的个数。

性质i、ii,都较容易理解,对于性质iii,因为在三维空间下,最多4个点就可以确定一个以它们为边界的球,这4个点就称为球体的支撑点,这里暂不考虑4个点共面这种特殊性况。

现在利用上面的3个性质来计算\(\text{MB}\left( P,R \right)\),如果\(P = \emptyset \)或者\(\left| R \right|=3\),则直接计算出\(\overline{\text{MB}}\left( R \right)\)。否则,随机的选择一个点\(p\in P\)并计算\(D\leftarrow \text{MB}\left( P-\left\{ p \right\},R \right)\),如果\(D\)存在且\(p\in D\),则\(\text{MB}\left( P,R \right)=D\);相反,如果\(p\notin D\)\(\text{MB}\left( P,R \right)=\text{MB}\left( P-\left\{ p \right\},R\bigcup \left\{ p \right\} \right)\)。计算\(\text{MB}\left( P,R \right)\)的算法伪码如下所示:

计算\(\text{MB}\left( P,R \right)\),返回最小包围球\(D\)
1.    if \(P = = \emptyset \) or \(\left| R \right| = = 3\), then
2.           \(D\leftarrow \overline{\text{MB}}\left( R \right)\);
3.    else
4.           随机选择一个点\(p\in P\);
5.           \(D\leftarrow \text{MB}\left( P-\left\{ p \right\},R \right)\);
6.           if \(D\)存在 and \(p\notin D\), then
7.                  \(D\leftarrow \text{MB}\left( P-\left\{ p \right\},R\bigcup \left\{ p \right\} \right)\);
8.           end if;
9.    end if;

现在分析下算法的复杂度,设\({{t}_{j}}\left( n \right)\)表示计算\(\text{MB}\left( P,R \right)\)的步数,其中\(\left| P \right|=n,\left| R \right|=4-j\),则\(0\le j\le 4\),显然有\({{t}_{j}}\left( 0 \right)=1,{{t}_{0}}\left( n \right)=1\)。如果\(j\succ 0,n\succ 0\),算法中要么计算\(\text{MB}\left( P-\left\{ p \right\},R \right)\),要么计算\(\text{MB}\left( P-\left\{ p \right\},R\bigcup \left\{ p \right\} \right)\),由性质iii可知,后者发生的概率是\(j/n\),所以可以得到复杂度计算的递归不等式

\[{t_j}\left( n \right) \le {t_j}\left( {n - 1} \right) + 1 + \frac{j}{n}{t_{j - 1}}\left( {n - 1} \right)\]

其中\({{t}_{1}}\left( n \right)\prec n+1\)\({{t}_{2}}\left( n \right)\prec 3n\)\({{t}_{3}}\left( n \right)\prec 10n\)\({{t}_{4}}\left( n \right)\prec 41n\),“41”是一个常数,所以计算三维空间下的最小包围球,可以在复杂度期望为\(O\left( n \right)\)的时间内得到。

因为需要每次随机的选择\(P\)的一个点,我们可以采用的一个做法,就是提前对\(P\)中的点重排列,随机打乱它们的顺序,以后只需要顺序遍历点集\(P\)中的元素即可。

若序列的前4个点正好是点集\(S\)的最小包围球的支撑点,则算法的效率将大大的提高,相反,若序列最后的4个点是最小包围球的支撑点,则算法中就会进行大量不必要的计算。因此,启发式的优化步骤必需能将最小包围球的4个支撑点提到序列的前面。做法就是:若一个点\(p\)满足算法第6步的条件,则被认定为“重要的点”,把它移动到序列的最前面。

上面提供的伪码是为了使读者更好的理解算法并分析复杂度,接下来将完整地描述算法伪码,首先是\(\overline{\text{MB}}\left( R \right)\)的实现伪码:

给定\(k\)个点的点集\(R=\left\{ {{R}_{0}},{{R}_{1}},\cdots ,{{R}_{k-1}} \right\}\),计算以该点集为边界的最小包围球,其中\(k\le 4\)\(R\)可以用数组来存储
createSphere(\(R\),\(k\))
1.    switch(\(k\))
2.           case 0:
3.                  \(center=\left( 0,0 \right)\);
4.                  \(radius=0\);
5.                  return;
6.           case 1:
7.                  \(center={{R}_{0}}\);
8.                  \(radius=0\);
9.                  return;
10.          case 2:
11.                 \(center=\left( {{R}_{0}}+{{R}_{1}} \right)/2\);
12.                 \(radius=\left\| {{R}_{1}}-{{R}_{0}} \right\|/2\);
13.                 return;
14.          case 3:
15.                 \(a=\left( {{R}_{1}}-{{R}_{0}} \right)\cdot \left( {{R}_{1}}-{{R}_{0}} \right),b=\left( {{R}_{1}}-{{R}_{0}} \right)\cdot \left( {{R}_{2}}-{{R}_{0}} \right),c=\left( {{R}_{2}}-{{R}_{0}} \right)\cdot \left( {{R}_{2}}-{{R}_{0}} \right)\);
16.                 \(d=ac-{{b}^{2}}\);
17.                 if \(d==0\), then
18.                        return;
19.                 end if;
20.                 \(s=\left( a-b \right)c/(2d),t=\left( c-b \right)a/(2d)\);
21.                 \(center={{R}_{0}}+s\left( {{R}_{1}}-{{R}_{0}} \right)+t\left( {{R}_{2}}-{{R}_{0}} \right)\);
22.                 \(radius=\left\| {{R}_{0}}-center \right\|\);
23.                 break;
24.          case 4:
25.                 \({{e}_{1}}=\left( {{X}_{1}},{{Y}_{1}},{{Z}_{1}} \right)={{R}_{1}}-{{R}_{0}},{{e}_{2}}=\left( {{X}_{2}},{{Y}_{2}},{{Z}_{2}} \right)={{R}_{2}}-{{R}_{0}},{{e}_{3}}=\left( {{X}_{3}},{{Y}_{3}},{{Z}_{3}} \right)={{R}_{3}}-{{R}_{0}}\);
26.                 \(V={{e}_{1}}\cdot \left( {{e}_{2}}\times {{e}_{3}} \right)\);
27.                 if \(V==0\), then
28.                        return;
29.                 end if;
30.                 \(V=2*V\);
31.                 \({{L}_{1}}={{e}_{1}}\cdot {{e}_{1}},{{L}_{2}}={{e}_{2}}\cdot {{e}_{2}},{{L}_{3}}={{e}_{3}}\cdot {{e}_{3}}\);
32.                 \(cente{{r}_{x}}={{R}_{0,x}}+\left[ \left( {{Y}_{2}}{{Z}_{3}}-{{Y}_{3}}{{Z}_{2}} \right){{L}_{1}}-\left( {{Y}_{1}}{{Z}_{3}}-{{Y}_{3}}{{Z}_{1}} \right){{L}_{2}}+\left( {{Y}_{1}}{{Z}_{2}}-{{Y}_{2}}{{Z}_{1}} \right){{L}_{3}} \right]/V\);
33.                 \(cente{{r}_{y}}={{R}_{0,y}}+\left[ -\left( {{X}_{2}}{{Z}_{3}}-{{X}_{3}}{{Z}_{2}} \right){{L}_{1}}+\left( {{X}_{1}}{{Z}_{3}}-{{X}_{3}}{{Z}_{1}} \right){{L}_{2}}-\left( {{X}_{1}}{{Z}_{2}}-{{X}_{2}}{{Z}_{1}} \right){{L}_{3}} \right]/V\);
34.                 \(cente{{r}_{z}}={{R}_{0,z}}+\left[ \left( {{X}_{2}}{{Y}_{3}}-{{X}_{3}}{{Y}_{2}} \right){{L}_{1}}-\left( {{X}_{1}}{{Y}_{3}}-{{X}_{3}}{{Y}_{1}} \right){{L}_{2}}+\left( {{X}_{1}}{{Y}_{2}}-{{X}_{2}}{{Y}_{1}} \right){{L}_{3}} \right]/V\);
35.                 \(radius=\left\| {{R}_{0}}-center \right\|\);
36.                 break;
37.   end switch;

2015-3-24 15-23-39

图9.10 经过不共线的3个点的最小包围球

在三维空间下,对于给定0个点、1个点、2个点的最小外接球的最小包围球的情况,比较容易理解。给定3个不共面的顶点,以它们为边界点的最小包围球,球心必定在3个点确定的平面上,如图9.10所示,设三个点分别是\({{R}_{0}},{{R}_{1}},{{R}_{2}}\),则球心\(C\)可以表示为

\[C = {R_0} + s\left( {{R_1} - {R_0}} \right) + t\left( {{R_2} - {R_0}} \right) \tag{9.23}\]

直线\(\overline{CM}\)是线段\(\overline{{{R}_{0}}{{R}_{2}}}\)的垂直平分线,因此有

\[{\left\| {\overline {{R_0}{R_2}} } \right\|^2} = {\left\| {{R_0} - {R_2}} \right\|^2} = 2\left( {C - {R_0}} \right) \cdot \left( {{R_2} - {R_0}} \right) \tag{9.24}\]

对于边\(\overline {{R_0}{R_1}} \)来说,同理,有

\[{\left\| {\overline {{R_0}{R_1}} } \right\|^2} = {\left\| {{R_0} - {R_1}} \right\|^2} = 2\left( {C - {R_0}} \right) \cdot \left( {{R_1} - {R_0}} \right) \tag{9.25}\]

把等式(9.23)代入等式(9.24)和(9.25),令\({{\left\| {{R}_{1}}-{{R}_{0}} \right\|}^{2}}=a,\left( {{R}_{1}}-{{R}_{0}} \right)\cdot \left( {{R}_{2}}-{{R}_{0}} \right)=b,\)\({{\left\| {{R}_{2}}-{{R}_{0}} \right\|}^{2}}=c\),则得到二元一次方程组

\[\left\{ \begin{array}{l}sa + tb = \frac{a}{2}\\sb + tc = \frac{c}{2}\end{array} \right. \tag{9.26}\]

解方程组,得

\[\left\{ \begin{array}{l}s = \frac{{\left( {a - b} \right)c}}{{2\left( {ac - {b^2}} \right)}}\\t = \frac{{\left( {c - b} \right)a}}{{2\left( {ac - {b^2}} \right)}}\end{array} \right. \tag{9.27}\]

对于不共线的三个点,一定有\(ac-{{b}^{2}}\ne 0\)

对于经过不共面的4个点的最小包围球,就是求这4个点的外接球。现在考虑计算\(\text{MB}\left( P,R \right)\)的实现伪码:

给定\(k\)个点的点集\(R=\left\{ {{R}_{0}},{{R}_{1}},\cdots ,{{R}_{k-1}} \right\}\)\(n\)个点的点集\(P=\left\{ {{P}_{0}},{{P}_{1}},\cdots ,{{P}_{n-1}} \right\}\)\(R\)可用数组来存储,\(P\)用链表来存储,\(E\)表示链表中的终止位置,计算\(R\)为边界且包含\(P\)内第1个元素始至\(E\)为止的最小包围球
smallestEnclosingSphere(\(R\),\(k\),\(P\),\(E\))
1.    createSphere(\(R\),\(k\));
2.    if \(k==4\), then
3.           return;
4.    end if;
5.    for ( \(p\)=\(P\).begin() ; \(address(p)\)!=\(E\) ; )
6.           if \(\left\| p-center \right\|\succ radius\), then
7.                  \({{R}_{k}}=p\);
8.                  smallestEnclosingSphere(\(R\),\(k+1\),\(P\),\(address(p)\));
9.                  \(q=next\left( P,address(p) \right)\);
10.                 \(P\).erase(\(address(p)\));
11.                 \(P\).push_front(\(p\));
12.                 \(p\)=\(q\);
13.          else
14.                 \(p=next\left( P,address(p) \right)\)
15.          end if;
16.   end for;

在有了两个主要的函数实现后,求最小包围球的算法伪码如下所示:

给定一个点集\(S=\left\{ {{P}_{0}},{{P}_{1}},\cdots ,{{P}_{n-1}} \right\}\),求包围该点集的最小包围球。
1.    随机打乱点集\(S\),设新的点集的序列表示为\(S=\left\{ {{Q}_{0}},{{Q}_{1}},\cdots ,{{Q}_{n-1}} \right\}\);
2.    \(R = \emptyset \);
3.    smallestEnclosingSphere(\(R\),0, \(S\),\(S.end()\));

9.3.2. 外接球

单形的外接球,就是经过单形所有顶点的球,满足限定条件:球心与所有的顶点之间的距离相同,设单形的顶点为\({{V}_{i}}=\left( {{x}_{i}},{{y}_{i}} \right)\)(不考虑顶点的顺序),球心为\(C\),半径为\(R\),则有

\[{\left\| {{V_i} - C} \right\|^2} = {R^2} \tag{9.28}\]

1.    二维外接圆

对于二维空间下,由不共线的3个顶点就可以确定一个外接圆,如图9.11所示,由等式(9.28)可以得到

2015-3-24 15-23-53

图9.11 在二维空间下,三角形的外接圆

\[\left\{ \begin{array}{l}{\left\| {{V_1} - C} \right\|^2} - {\left\| {{V_0} - C} \right\|^2} = 0\\{\left\| {{V_2} - C} \right\|^2} - {\left\| {{V_0} - C} \right\|^2} = 0\end{array} \right. \tag{9.29}\]

等式可以化简为

\[\left\{ \begin{array}{l}\left( {{V_1} - {V_0}} \right) \cdot \left( {C - {V_0}} \right) = \frac{1}{2}{\left\| {{V_1} - {V_0}} \right\|^2}\\\left( {{V_2} - {V_0}} \right) \cdot \left( {C - {V_0}} \right) = \frac{1}{2}{\left\| {{V_2} - {V_0}} \right\|^2}\end{array} \right. \tag{9.30}\]

等式(9.30)是类似于\(AX=a\)的方程形式,是一个方程组,根据克拉默法则,求解方程组,设\({{X}_{i}}={{x}_{i}}-{{x}_{0}},{{Y}_{i}}={{y}_{i}}-{{y}_{0}},{{L}_{1}}=\left\| {{V}_{1}}-{{V}_{0}} \right\|,{{L}_{2}}=\left\| {{V}_{2}}-{{V}_{0}} \right\|\),解得

\[\left\{ \begin{array}{l}x = {x_0} + \frac{1}{{4A}}\left( {{Y_2} \cdot L_1^2 - {Y_1} \cdot L_2^2} \right)\\y = {y_0} + \frac{1}{{4A}}\left( {{X_1} \cdot L_2^2 - {X_2} \cdot L_1^2} \right)\end{array} \right. \tag{9.31}\]

其中,\(A\)表示由3个点确定的三角形的有符号面积,由于3个点不共线,所以一定可以保证\(A\ne 0\)\(A\)可以表示为

\[A = \frac{1}{2}\left| {\begin{array}{*{20}{c}}{{X_1}}&{{Y_1}}\\{{X_2}}&{{Y_2}}\end{array}} \right| = \frac{1}{2}\left( {{V_1} - {V_0}} \right) \times \left( {{V_2} - {V_0}} \right) \tag{9.32}\]

那么,外接圆的半径为

\[R = \sqrt {{{\left( {x - {x_0}} \right)}^2} + {{\left( {y - {y_0}} \right)}^2}} \tag{9.33}\]

2.    三维外接圆

2015-3-24 15-24-06

图9.12 三维空间下,四面体的外接球

对于三维空间下,由不共面的4个顶点就可以确定一个外接球,如图9.12所示,与二维的情况类似,可以得到方程组

\[\left\{ \begin{array}{l}\left( {{V_1} - {V_0}} \right) \cdot \left( {C - {V_0}} \right) = \frac{1}{2}{\left\| {{V_1} - {V_0}} \right\|^2}\\\left( {{V_2} - {V_0}} \right) \cdot \left( {C - {V_0}} \right) = \frac{1}{2}{\left\| {{V_2} - {V_0}} \right\|^2}\\\left( {{V_3} - {V_0}} \right) \cdot \left( {C - {V_0}} \right) = \frac{1}{2}{\left\| {{V_3} - {V_0}} \right\|^2}\end{array} \right. \tag{9.34}\]

根据克拉默法则,求解方程组的解,设\({{X}_{i}}={{x}_{i}}-{{x}_{0}},{{Y}_{i}}={{y}_{i}}-{{y}_{0}},{{Z}_{i}}={{z}_{i}}-{{z}_{0}},{{L}_{1}}=\left\| {{V}_{1}}-{{V}_{0}} \right\|\)\({{L}_{2}}=\left\| {{V}_{2}}-{{V}_{0}} \right\|,{{L}_{3}}=\left\| {{V}_{3}}-{{V}_{0}} \right\|\),得

\[\left\{ \begin{array}{l}x = {x_0} + \frac{1}{{12V}}\left[ { + \left( {{Y_2}{Z_3} - {Y_3}{Z_2}} \right)L_1^2 - \left( {{Y_1}{Z_3} - {Y_3}{Z_1}} \right)L_2^2 + \left( {{Y_1}{Z_2} - {Y_2}{Z_1}} \right)L_3^2} \right]\\y = {y_0} + \frac{1}{{12V}}\left[ { - \left( {{X_2}{Z_3} - {X_3}{Z_2}} \right)L_1^2 + \left( {{X_1}{Z_3} - {X_3}{Z_1}} \right)L_2^2 - \left( {{X_1}{Z_2} - {X_2}{Z_1}} \right)L_3^2} \right]\\z = {z_0} + \frac{1}{{12V}}\left[ { + \left( {{X_2}{Y_3} - {X_3}{Y_2}} \right)L_1^2 - \left( {{X_1}{Y_3} - {X_3}{Y_1}} \right)L_2^2 + \left( {{X_1}{Y_2} - {X_2}{Y_1}} \right)L_3^2} \right]\end{array} \right. \tag{9.35}\]

其中,\(V\)表示由4个点确定的四面体的有符号体积,由于4个点不共面,所以一定可以保证\(V\ne 0\)\(V\)可以表示为

\[V = \frac{1}{6}\left| {\begin{array}{*{20}{c}}{{X_1}}&{{Y_1}}&{{Z_1}}\\{{X_2}}&{{Y_2}}&{{Z_2}}\\{{X_3}}&{{Y_3}}&{{Z_3}}\end{array}} \right| = \left( {{V_1} - {V_0}} \right) \cdot \left[ {\left( {{V_2} - {V_0}} \right) \times \left( {{V_3} - {V_0}} \right)} \right] \tag{9.36}\]

那么,外接球的半径为

\[R = \sqrt {{{\left( {x - {x_0}} \right)}^2} + {{\left( {y - {y_0}} \right)}^2} + {{\left( {z - {z_0}} \right)}^2}} \tag{9.37}\]

9.3.3. 内切球

单形的内切球就是包含于单形内的体积最大的球,这个球必须与这个单形的所有面相切,球心与所有面之间的距离都相同,设单形的顶点为\({{V}_{i}}=\left( {{x}_{i}},{{y}_{i}} \right)\)(不考虑顶点的顺序),球心为\(C\),半径为\(R\),若单形的一个面上指向内的单位长度法向量是\({{\vec{n}}_{i}}\),则有

\[{\vec n_i} \cdot \left( {C - {V_i}} \right) = R \tag{9.38}\]

1.    二维内切圆

2015-3-24 15-24-20

图9.13 二维空间下,三角形的内切圆

在二维空间下,单形顶点的个数为3,即计算三角形的内切圆,如图9.13所示。设\({{V}_{3}}={{V}_{0}}\),边的长度表示为\({{L}_{i}}=\left\| {{V}_{i+1}}-{{V}_{i}} \right\|\),边的方向可以表示为\({{\vec{d}}_{i}}=\left( {{V}_{i+1}}-{{V}_{i}} \right)/{{L}_{i}}\)\({{\vec{d}}_{i}}\)是单位长度,边上指向内的单位长度法线为\({{\vec{n}}_{i}}\),其中\(0\le i\le 2\)

圆心可以用3个顶点的参数方程表示\(C={{t}_{0}}{{V}_{0}}+{{t}_{1}}{{V}_{1}}+{{t}_{2}}{{V}_{2}}\),其中\({{t}_{0}}+{{t}_{1}}+{{t}_{2}}=1\)。把圆心的参数方程代入等式(9.38)中,得到等式

\[{t_2}{L_2}\left( {{{\vec d}_0} \cdot {{\vec n}_2}} \right) = {t_0}{L_0}\left( {{{\vec d}_1} \cdot {{\vec n}_0}} \right) = {t_1}{L_1}\left( {{{\vec d}_2} \cdot {{\vec n}_1}} \right) = R \tag{9.39}\]

对于三角形的面积\(A\)\(A\succ 0\)),有

\[2A = {L_0}{L_2}\left( {{{\vec d}_0} \cdot {{\vec n}_2}} \right) = {L_1}{L_0}\left( {{{\vec d}_1} \cdot {{\vec n}_0}} \right) = {L_2}{L_1}\left( {{{\vec d}_2} \cdot {{\vec n}_1}} \right) \tag{9.40}\]

把等式(9.40)代入等式(9.39)中,消去\(\left( {{{\vec d}_0} \cdot {{\vec n}_2}} \right),\left( {{{\vec d}_1} \cdot {{\vec n}_0}} \right),\left( {{{\vec d}_2} \cdot {{\vec n}_1}} \right)\),解得

\[\left\{ \begin{array}{l}{t_0} = \frac{{R{L_1}}}{{2A}}\\{t_1} = \frac{{R{L_2}}}{{2A}}\\{t_2} = \frac{{R{L_0}}}{{2A}}\end{array} \right. \tag{9.41}\]

其中,\(A=\left\| \left( {{V}_{1}}-{{V}_{0}} \right)\times \left( {{V}_{2}}-{{V}_{0}} \right) \right\|/2\),把等式(9.41)代入\({{t}_{0}}+{{t}_{1}}+{{t}_{2}}=1\),令\({{L}_{0}}+{{L}_{1}}+{{L}_{2}}=L\),计算出内切圆的半径为

\[R = \frac{{2A}}{L} \tag{9.42}\]

把等式(9.42)代入(9.41),计算出\({{t}_{0}},{{t}_{1}},{{t}_{2}}\)的值后,得到内切圆的圆心为

\[C = \frac{1}{L}\left( {{L_1}{V_0} + {L_2}{V_1} + {L_0}{V_2}} \right) \tag{9.43}\]

2.    三维内切球

2015-3-24 15-24-33

图9.14 三维空间下四面体的内切球

在三维空间下,单形顶点的个数为4,即计算四面体的的内切圆,如图9.14所示。设\({{A}_{0}}\)表示\(\Delta {{V}_{0}}{{V}_{1}}{{V}_{2}}\)的面积,该面的单位法向量表示为\({{\vec{n}}_{0}}=\left( {{V}_{1}}-{{V}_{0}} \right)\times \left( {{V}_{2}}-{{V}_{0}} \right)/\left( 2{{A}_{0}} \right)\)\({{A}_{1}}\)表示\(\Delta {{V}_{1}}{{V}_{2}}{{V}_{3}}\)的面积,该面的单位法向量表示为\({{\vec{n}}_{1}}=\left( {{V}_{3}}-{{V}_{1}} \right)\times \left( {{V}_{2}}-{{V}_{1}} \right)/\left( 2{{A}_{1}} \right)\)\({{A}_{2}}\)表示\(\Delta {{V}_{0}}{{V}_{2}}{{V}_{3}}\)的面积,该面的单位法向量表示为\({{\vec{n}}_{2}}=\left( {{V}_{3}}-{{V}_{2}} \right)\times \left( {{V}_{0}}-{{V}_{2}} \right)/\left( 2{{A}_{2}} \right)\)\({{A}_{3}}\)表示\(\Delta {{V}_{0}}{{V}_{1}}{{V}_{3}}\)的面积,该面的单位法向量表示为\({{\vec{n}}_{3}}=\left( {{V}_{0}}-{{V}_{3}} \right)\times \left( {{V}_{1}}-{{V}_{3}} \right)/\left( 2{{A}_{3}} \right)\)。对于三角形\(\Delta pqr\)的面积\(A\)的计算,可以采用等式

\[{A_{\Delta pqr}} = \frac{1}{2}\vec n \cdot \left[ {p \times q + q \times r + r \times p} \right] \tag{9.44}\]

其中, \(\vec{n}=\left[ \left( q-p \right)\times \left( r-p \right) \right]/\left\| \left( q-p \right)\times \left( r-p \right) \right\|\),可以保证\({{A}_{\Delta pqr}}\ge 0\)

与二维类似,球心可以用参数方程表示为\(C=\sum\nolimits_{i=0}^{3}{{{t}_{i}}{{V}_{i}}}\),其中\(\sum\nolimits_{i=0}^{3}{{{t}_{i}}}=1\)。内切球的半径可以表示为

\[R = \left\| {{{\vec n}_i} \cdot \left( {C - {V_i}} \right)} \right\|,0 \le i \le 3 \tag{9.45}\]

把等式\(\sum\nolimits_{i = 0}^3 {{t_i}} = 1\)代入(9.45)中,则有

\[{t_3}\left\| {{{\vec n}_0} \cdot \left( {{V_3} - {V_0}} \right)} \right\| = {t_0}\left\| {{{\vec n}_1} \cdot \left( {{V_0} - {V_1}} \right)} \right\| = {t_1}\left\| {{{\vec n}_2} \cdot \left( {{V_1} - {V_2}} \right)} \right\| = {t_2}\left\| {{{\vec n}_3} \cdot \left( {{V_2} - {V_3}} \right)} \right\| = R \tag{9.46}\]

四面体的体积\(V\)可以表示为

\[\begin{array}{l}{\rm{ }}\left\| {\left[ {\left( {{V_1} - {V_0}} \right) \times \left( {{V_2} - {V_0}} \right)} \right] \cdot \left( {{V_3} - {V_0}} \right)} \right\| = \left\| {\left[ {\left( {{V_3} - {V_1}} \right) \times \left( {{V_2} - {V_1}} \right)} \right] \cdot \left( {{V_0} - {V_1}} \right)} \right\|\\ = \left\| {\left[ {\left( {{V_3} - {V_2}} \right) \times \left( {{V_0} - {V_2}} \right)} \right] \cdot \left( {{V_1} - {V_2}} \right)} \right\| = \left\| {\left[ {\left( {{V_0} - {V_3}} \right) \times \left( {{V_1} - {V_3}} \right)} \right] \cdot \left( {{V_2} - {V_3}} \right)} \right\|\\ = 6V\end{array} \tag{9.47}\]

\({{\vec{n}}_{0}},{{\vec{n}}_{1}},{{\vec{n}}_{2}},{{\vec{n}}_{3}}\)的值和等式(9.47)代入(9.46)中,得到

\[\left\{ \begin{array}{l}{t_0} = \frac{{R{A_1}}}{{3V}}\\{t_1} = \frac{{R{A_2}}}{{3V}}\\{t_2} = \frac{{R{A_3}}}{{3V}}\\{t_3} = \frac{{R{A_0}}}{{3V}}\end{array} \right. \tag{9.48}\]

把等式(9.48)代入等式\(\sum\nolimits_{i=0}^{3}{{{t}_{i}}}=1\)中,令\(A={{A}_{0}}+{{A}_{1}}+{{A}_{2}}+{{A}_{3}}\),得内切球的半径为

\[R = \frac{{3V}}{A} \tag{9.49}\]

再把内切球的半径,代入等式(9.48)中,计算出\({{t}_{i}}\)的值后,代入球心的参数方程,得

\[C = \frac{1}{A}\left( {{A_1}{V_0} + {A_2}{V_1} + {A_3}{V_2} + {A_0}{V_3}} \right) \tag{9.50}\]

其中,四面体的体积\(V\),三角形的面积\({{A}_{i}}\)的计算,如(9.47)和等式(9.44)所示。

9.3.4. 线性对象与球的相交测试

问题描述:线性对象\(L(t)=P+t\cdot \vec{d}\)与以\(C\)为球心、\(R\)为半径的球体的相交测试,不需要保证\(\vec{d}\)是单位长度。

本节将描述两种线性对象与球体的相交测试方法,方法一是采用参数方程求解的方法,可以计算出直线、射线、线段与球的交点,方法二在Akenine-Möller等[14](2002)中有描述,只涉及射线与球的相交测试,采用几何关系来优化测试。此外,对二维空间下,线性对象与圆的相交测试,也可以采用本节介绍的两种方法,由于维度更小,所以更加简单,这里就不再赘述。

方法一. 参数方程法

球面可以用参数方程表示为

\[f(X) = {\left\| {X - C} \right\|^2} = {R^2} \tag{9.51}\]

现在考虑直线与球面的相交,把直线的参数方程代入等式(9.51)中,得

\[\begin{array}{l}{\left\| {\left( {P + t \cdot \vec d} \right) - C} \right\|^2} = {R^2}\\ \Leftrightarrow {\left\| {\vec d} \right\|^2}{t^2} + 2t\left[ {\vec d \cdot \left( {P - C} \right)} \right] + \left( {P - C} \right) \cdot \left( {P - C} \right) - {R^2} = 0\end{array} \tag{9.52}\]

这是类似于\(a{t^2} + 2bt + c = 0\)的一元二次方程,直接求解该方程的解,得

\[t = \frac{{ - \vec d \cdot \left( {P - C} \right) \pm \sqrt {{{\left[ {\vec d \cdot \left( {P - C} \right)} \right]}^2} - \left( {P - C} \right) \cdot \left( {P - C} \right) + {R^2}} }}{{{{\left\| {\vec d} \right\|}^2}}} \tag{9.53}\]

\(\Delta ={{\left[ \vec{d}\cdot \left( P-C \right) \right]}^{2}}-\left( P-C \right)\cdot \left( P-C \right)+{{R}^{2}}\),则解可以分为三种情况:

  • i.     若\(\Delta \prec 0\),方程无解,直线与球体不相交;
  • ii.    若\(\Delta =0\),方程有一个解,直线与球体相切;
  • iii.   若\(\Delta \succ 0\),方程有两个解,直线与球体相交。

2015-3-24 15-24-46

图9.15 直线与球相交的三种情况

直线与球体相交的三种情况,如图9.15所示。图(a)就是第i种情况,两者不相交;图(b)就是第ii种情况,相切于一点;图(c)就是第iii种情况,相交于两点。算法的伪码如下所示:

直线\(L(t)=P+t\cdot \vec{d}\)与以\(C\)为球心、\(R\)为半径的球体的相交测试,\(\vec{d}\)可以是非单位长度,若相交,求计算出所有交点
return (交点的个数)
1.    \(e=P-C\);
2.    \(b=\vec{d}\cdot e,c=e\cdot e-R\cdot R\);
3.    \(k=b\cdot b-c\);
4.    if \(k\prec 0\), then
5.           return 0;
6.    else if \(k==0\), then
7.           \(a=\vec{d}\cdot \vec{d}\);
8.           \(t=-b/a\);
9.           \(Q=P+t\cdot \vec{d}\);
10.          return 1;
11.   else
12.          \(k'=\sqrt{k}\);
13.          \(a=\vec{d}\cdot \vec{d}\);
14.          \({{t}_{1}}=\left( -b-k' \right)/a,{{t}_{2}}=\left( -b+k' \right)/a\);
15.          \({{Q}_{1}}=P+{{t}_{1}}\vec{d},{{Q}_{2}}=P+{{t}_{2}}\vec{d}\);
16.          return 2;
17.   end if;

对于射线与球体的相交测试中,需要保证\(t\)的值限定的区间范围\(\left[ 0,+\infty \right]\)内;对于线段与球体的相交测试中,需要保证\(t\)的值限定的区间范围\(\left[ 0,1 \right]\)内。

方法二. 射线与球体相交测试的优化法

设射线的参数方程为\(L(t)=P+t\cdot \vec{d}\),其中\(\vec{d}\)是单位长度。如图9.16(a)所示,计算出从射线原点出发到球体中心的向量\(\vec{a}=C-P\),计算出该向量的长度\({{a}^{2}}=\vec{a}\cdot \vec{a}\),若\({{a}^{2}}\prec {{R}^{2}}\),那么说明射线的原点在球体内部,一定与球体相交,如图9.16(c)所示,若只是想判断它们是否相交,算法至此就可以结束了。否则,计算向量\(\vec{a}\)沿着\(\vec{d}\)方向的投影:\(l=\vec{d}\cdot \vec{a}\),若\(l\prec 0\),说明球体在射线原点的后面,且射线的原点在球体外,因此可以判定射线与球体一定不相交,完成第一次排除测试,如图9.16(b)所示。否则,利用勾股定理,计算出球心与投影之间距离的平方:\({{m}^{2}}={{a}^{2}}-{{l}^{2}}\),若\({{m}^{2}}\succ {{R}^{2}}\),则可以判定射线与球体一定不相交,完成第二次排除测试。如果射线和球体经过两次排除测试,则就可以判定它们一定相交。如果只希望判定它们是否相交,算法至此就可以结束了。

2015-3-24 15-24-58

图9.16 射线与球体相交测试的优化方法

接下来,计算射线与球体的交点,计算出距离的平方,\({{q}^{2}}={{R}^{2}}-{{m}^{2}}\),因为\({{m}^{2}}\le {{R}^{2}}\),所以计算出\(q\)的值。相交点与射线原点之间的距离为\(t=l\pm q\),求第一个相交点,若\({{l}^{2}}\succ {{R}^{2}}\),说明射线原点在球体外,取\(t=l-q\),否则取\(t=l+q\),最后把\(t\)值代入射线的参数方程中,就求出第一个相交点。

射线与球体相交测试并计算交点的算法伪码如下所示:

射线\(L(t)=P+t\cdot \vec{d}\)与以\(C\)为球心、\(R\)为半径的球体的相交测试,\(\vec{d}\)是单位长度,若相交,计算出第一个交点
return ({DISJOINT, INTERSECTING})
1.  \(\vec{a}=C-P\);
2.  \(l=\vec{d}\cdot \vec{a}\);
3.  \({{a}^{2}}=\vec{a}\cdot \vec{a}\)
4.  if \({{a}^{2}}\succ {{R}^{2}}\) and \(l\prec 0\), then
5.      return DISJOINT;
6.  end if;
7.  \({{m}^{2}}={{a}^{2}}-{{l}^{2}}\);
8.  if \({{m}^{2}}\succ {{R}^{2}}\), then
9.      return DISJOINT;
10. end if;
11. \(q=\sqrt{{{R}^{2}}-{{m}^{2}}}\);
12. if \({{a}^{2}}\succ {{R}^{2}}\), thn
13.     \(t=l-q\);
14. else
15.     \(t=l+q\);
16. end if;
17. \(Q=P+t\cdot \vec{d}\);
18. return INTERSECTING;

若不需要计算射线与球体的交点,只需要判断它们是否相交,那么伪码如下所示:

射线\(L(t)=P+t\cdot \vec{d}\)与以\(C\)为球心、\(R\)为半径的球体的相交测试,\(\vec{d}\)是单位长度,不需要计算交点
return ({DISJOINT, INTERSECTING})
1.  \(\vec{a}=C-P\);
2.  \({{a}^{2}}=\vec{a}\cdot \vec{a}\);
3.  if \({{a}^{2}}\le {{R}^{2}}\), then
4.      return INTERSECTING;
5.  end if;
6.  \(l=\vec{d}\cdot \vec{a}\);
7.  if \(l\prec 0\), then
8.      return DISJOINT;
9.  end if;
10. \({{m}^{2}}={{a}^{2}}-{{l}^{2}}\);
11. if \({{m}^{2}}\succ {{R}^{2}}\), then
12.     return DISJOINT;
13. end if;
14. return INTERSECTING;

与参数方程法相比,具有相同的运算次数,但最大的不同是该方法较早地进行筛除计算,从而使这个算法总的开销要低一些。

spacer

Leave a reply