浏览量:366

8.4. 礼物包裹算法

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

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

8.4. 礼物包裹算法

8.4.1. 基本思想

礼物包裹算法最早由Chand&Kapur[3] (1970)提出的,它不仅可以实现二维、三维凸包,还可以实现更高维的凸包,算法复杂度是\(O(hn)\)\(h\)表示输出的面的数量,\(n\)表示点集的个数,算法复杂度跟输出面相关。直观上,三维的礼物包裹算法,可以看做是在点集的外面包围了一张纸,每次更新一个新的顶点,用纸覆盖住它,就像包裹礼物一样,直至覆盖整个点集。本小节只描述三维的礼物包裹算法,对于更高维的情况,可以参考Chand&Kapur[3] (1970)以及Preparata&Shamos[4](1985)的描述。

礼物包裹算法基于一个简单的理论:任意一个三维凸包的每一条边有且只有两个相邻面。三维凸包上的每个面用三角形面片表示,若一个凸包有\(h\)个面,则凸包共有\(3h/2\)条边,根据欧拉公式可知,顶点的个数是\(n=h/2+2\)

设三维凸包上的面用三角形面片表示,本节暂不考虑各种退化情况,三维凸包的礼物包裹算法的算法伪码如下所示:

三维空间下点集\(S=\{{{p}_{0}},{{p}_{1}},...,{{p}_{n-1}}\}\)的凸包
1.    面集合为Q,存储边的集合\(T\);
2.    \(Q \leftarrow \emptyset \);
3.    \(T \leftarrow \emptyset \);
4.    找到一个初始的凸包面\(F\);
5.    \(Q\leftarrow F\);
6.    \(T\leftarrow \)\(F\)的每条边信息和计数器;   //计数器初始值为1
7.    while( \(Q\)非空 )
8.           \(F\leftarrow Q\);
9.           for 面\(F\)的每一条边\(e\), then
10.                 根据边\(e\),从集合\(T\)中检索出相应的边结构\(E\);
11.                 if \(E.count\prec 2\), then
12.                        以边\(e\)为边界,找到新的面\(F'\);
13.                        \(Q\leftarrow F'\);
14                         更新边集合\(T\);   //若面\(F'\)的边存储在集合\(T\)中,对应边的计数器值加1;
//否则,创建新的边结构,存入集合\(T\)中,初始计数器的值为1;
15.                 end if;
16.          end for;
17.   end while;

实现礼物包裹算法时,必需解决算法中要处理的三个关键的子问题:

  1. 如何存储边结构的集合\(T\)
  2. 给定一个凸包上的面,如何确定与该面片上任意一条边相邻的三角形面片;
  3. 如何确定在凸包上的初始三角形面片;

子问题1:如何设计存储边的集合的\(T\)数据结构。

首先,确定两点大小关系的对比规则,设两点\({{P}_{0}}({{x}_{0}},{{y}_{0}},{{z}_{0}})\)\({{P}_{1}}({{x}_{1}},{{y}_{1}},{{z}_{1}})\),若\({{P}_{0}}<{{P}_{1}}\),当且仅当\(({{x}_{0}}<{{x}_{1}})\vee ({{x}_{0}}={{x}_{1}}\wedge {{y}_{0}}<{{y}_{1}})\vee ({{x}_{0}}={{x}_{1}}\wedge {{y}_{0}}={{y}_{1}}\wedge {{z}_{0}}<{{z}_{1}})\);若\({{P}_{0}}={{P}_{1}}\),当且仅当\({{x}_{0}}={{x}_{1}}\wedge \)\({{y}_{0}}={{y}_{1}}\wedge {{z}_{0}}={{z}_{1}}\);否则,\({{P}_{0}}>{{P}_{1}}\)。设两条边\({{e}_{0}}({{p}_{0}},{{q}_{0}})\)\({{e}_{1}}({{p}_{1}},{{q}_{1}})\),若\({{e}_{0}}<{{e}_{1}}\),当且仅当\(({{p}_{0}}<{{p}_{1}})\vee ({{p}_{0}}={{p}_{1}}\wedge {{q}_{0}}<{{q}_{1}})\);若\({{e}_{0}}={{e}_{1}}\),当且仅当\({{p}_{0}}={{p}_{1}}\wedge {{q}_{0}}={{q}_{1}}\);否则,\({{e}_{0}}>{{e}_{1}}\)

接着,考虑存储边的数据结构的设计,每个边结构表示为\(\{p,q,count\}\),其中\(p,q\)表示边的两个端点,\(count\)记录与边相邻面的个数,通过计数器可以判定与该边是否存在相邻面,因为凸包上的每一条边当且仅当有2个相邻面。

把边结构存储在一个平衡二叉树的数据结构中,边结构的查询、插入操作只需要\(O(\log m)\),其中\(m\)表示边的个数。

子问题2:给定一个凸包上的三角形面片,如何确定与该面片上任意一条边相邻的三角形面片。

2015-3-24 10-32-44

图8.10 三角形是\({{p}_{0}}{{p}_{1}}{{p}_{2}}\)凸包上的面,在平面\(\pi \)上,不同的点\({{p}_{i}}\)和边\({{p}_{0}}{{p}_{1}}\)确定的平面跟平面\(\pi \)有不同的夹角大小

设三个点\(\{{{p}_{0}},{{p}_{1}},{{p}_{2}}\}\)形成的三角形是凸包上的面片,三个点所在的平面是\(\pi \),从点集\(\{{{p}_{0}},{{p}_{1}},...,{{p}_{n}}\}\)中到找一个点\({{p}_{i}}\),使得由点\({{p}_{0}}\),\({{p}_{1}}\)\({{p}_{i}}\)形成的平面\({{\pi }_{i}}\)与平面\(\pi \)的夹角最大,那么由点\({{p}_{i}}\)与边\(\overline{{{p}_{0}}{{p}_{1}}}\)生成的三角形面片即是凸包的面。以图8.10所示为例,点\({{p}_{3}}\),\({{p}_{4}}\),\({{p}_{5}}\),\({{p}_{6}}\)与边\(\overline{{{p}_{0}}{{p}_{1}}}\)分别可以确定四个平面,这四个平面与平面\(\pi \)有一个夹角,选择夹角最大的平面,即由点\({{p}_{6}}\)与边\({{p}_{0}}{{p}_{1}}\)确定的平面,因此\(\{{{p}_{0}},{{p}_{1}},{{p}_{6}}\}\)生成的三角形面片是凸包上的面。

2015-3-24 10-33-02

图8.11 计算由点\({{p}_{i}}\),\({{p}_{0}}\)\({{p}_{1}}\)确定的面,与平面\(\pi \)的夹角

在介绍了采用的策略后,接下来介绍如何实施策略,即如何判断哪个点与三角形面片确定的平面\({{\pi }_{i}}\)与平面\(\pi \)之间的夹角最大。如图8.11所示,设平面\(\pi \)的法向量是\(\vec{n}\),若三角形点的顺序如该图所示,则法向量为\(\vec{n}=\left( -({{p}_{1}}-{{p}_{0}})\times ({{p}_{0}}-{{p}_{2}}) \right)/\left\| ({{p}_{1}}-{{p}_{0}})\times ({{p}_{0}}-{{p}_{2}}) \right\|\),即\(\left\| {\vec{n}} \right\|=1\),方向为垂直于平面\(\pi \)向上。因为\(\Delta {p_0}{p_1}{p_2}\)在凸包上,所以对于任意一个点\({{p}_{i}}\)\({{p}_{i}}\in \{{{p}_{0}},\cdots ,{{p}_{n}}\}\)),都在平面上(On the plane)或者平面的上方空间(Above the plane)。指定一个方向向量\(\vec{a}\),令\(\vec{a}\)同时垂直于法向量\(\vec{n}\)和边\(\overline{{{p}_{0}}{{p}_{1}}}\),则向量\(\vec{a}=\left( ({{p}_{1}}-{{p}_{0}})\times \vec{n} \right)/\left\| ({{p}_{1}}-{{p}_{0}})\times \vec{n} \right\|\),即有\(\left\| {\vec{a}} \right\|=1\)。点\({{p}_{2}}\)\(\vec{a}\)所在的方向上。连接点\({{p}_{i}}\)和点\({{p}_{0}}\),得到向量\({{v}_{i}}={{p}_{i}}-{{p}_{0}}\),计算\({{v}_{i}}\)在法向量\(\vec{n}\)和向量\(\vec{a}\)方向上的投影,分别是\(\vec{n}\cdot {{v}_{i}}\)\(\vec{a}\cdot {{v}_{i}}\),则对于不在平面\(\pi \)上的点\({{p}_{i}}\)的夹角可以用等式(8.7)表示:

\[{\rho _i} = - \frac{{\vec a \cdot {v_i}}}{{\vec n \cdot {v_i}}} \tag{8.7}\]

遍历点集上的每个点,选择\({{\rho }_{i}}\)值最大的点,就要所求的点,即:

\[\rho = \max {\rho _i} \tag{8.8}\]

 

2015-3-24 10-33-54

图8.12 计算经过点\({{p}_{i}}\)和边\({{p}_{0}}{{p}_{1}}\)的面的法向量

由点\({{p}_{0}}\)向点\({{p}_{1}}\)方向观察图8.11,可以得到横截面如图8.12所示,设\({{\vec{n}}_{i}}\)是由点\({{p}_{i}}\),\({{p}_{0}}\),\({{p}_{1}}\)确定的平面的法向量,向量\({{v}_{i}}'\)是向量\({{v}_{i}}\)在如9.12所示的横截面上的投影,易知向量\({{v}_{i}}'\)与法向量\({{\vec{n}}_{i}}\)垂直。设\({{\vec{n}}_{i}}\)\(\vec{a}\)方向上的投影为\({{\lambda }_{1}}\),在\(\vec{n}\)方向上的投影为\({{\lambda }_{2}}\),由相似三角形的性质可得:

\[\frac{{{\lambda _1}}}{{{\lambda _2}}} = \frac{{\vec n \cdot {v_i}}}{{\vec a \cdot {v_i}}} \tag{8.9}\]

那么法向量\({{\vec{n}}_{i}}\)的方向可以表示为

\[{\vec n_i} = (\vec n \cdot {v_i}) \cdot \vec a - (\vec a \cdot {v_i}) \cdot \vec n \tag{8.10}\]

可以注意到,我们只需要找到\({{\rho }_{i}}\)值最大的点,并不关心具体的值是多少,因此在计算向量\(\vec{n}\)\(\vec{a}\),可以不对它们进行归一化处理,即\(\vec{n}=-({{p}_{1}}-{{p}_{0}})\times ({{p}_{0}}-{{p}_{2}}),\vec{a}=({{p}_{1}}-{{p}_{0}})\times \vec{n}\)。未归一化的向量\(\vec{n}\)\(\vec{a}\)也不会对等式(8.10)中法向量\({{\vec{n}}_{i}}\)的方向造成影响,读者可以自行证明。确定\({{\rho }_{i}}\)值最大的点需要遍历点集,算法的复杂度为\(O(n)\)

子问题3:如何确定在凸包上的初始三角形面片。

算法的起始阶段,先排除点集中重复的点,找到初始的凸包面\(\Delta {{p}_{0}}{{p}_{1}}{{p}_{2}}\),如图8.13所示,从点集中找到最小的点\({{p}_{0}}\),即点的\(X\)分量最小,若\(X\)分量最小的点有多个,则取其中\(Y\)分量最小的点,若\(X,Y\)分量最小的点有多个,则取其中\(Z\)分量最小的点。最小的点\({{p}_{0}}\)是唯一的,经过点\({{p}_{0}}\)作一个平面\({{\pi }_{0}}\),平面的法向量表示为\({{\vec{n}}_{0}}=({{a}_{1}},0,0),{{a}_{1}}\ne 0\),任取一个向量\({{\vec{a}}_{0}}=(0,{{d}_{1}},0),{{d}_{1}}\ne 0\),遍历点集上的点,使用等式(8.7)计算得到\({{\rho }_{i}}\)最大的点,设该点是\({{p}_{1}}\)。采用等式(8.10),计算出平面\({{\pi }_{1}}\)的法向量\({{\vec{n}}_{1}}=({{\vec{n}}_{0}}\cdot {{v}_{1}})\cdot {{\vec{a}}_{0}}-({{\vec{a}}_{0}}\cdot {{v}_{1}})\cdot {{\vec{n}}_{0}}\),其中\({{v}_{1}}={{p}_{1}}-{{p}_{0}}\),易知法向量\({{\vec{n}}_{1}}\)\(({{b}_{1}},{{b}_{2}},0)\)的形式,第3个分量为0。接着,指定一个与法向量\({{\vec{n}}_{1}}\)和边\(\overline{{{p}_{0}}{{p}_{1}}}\)垂直的向量\({{\vec{a}}_{1}}=({{e}_{1}},{{e}_{2}},{{e}_{3}})\),可以采用等式\({{\vec{a}}_{1}}=({{p}_{1}}-{{p}_{0}})\times {{\vec{n}}_{1}}\)来计算。遍历点集上的点,使用等式(8.8)计算出\({{\rho }_{i}}\)最大的点,设该点是\({{p}_{2}}\)。至此,初始凸包上的三角形面片的三个点\({{p}_{0}},{{p}_{1}},{{p}_{2}}\)都计算出来了,算法需要三次遍历点集,算法规模是\(O(3n)\)

2015-3-24 10-34-09

图8.13 在凸包上的初始三角形面片\(\Delta {p_0}{p_1}{p_2}\)

综上,在礼物包裹算法中,确定初始的三角形面片,需要时间\(O(3n)\);每确定一个新的凸包面,需要时间\(O(logm)+O(n)\),其中\(m\)表示边,又\(m\le 3(n-2)\),即时间复杂度为\(O(n)\)。因此,三维礼物包裹算法的时间复杂度为\(O(nh)\)

8.4.2. 特殊情况

1.    多点共面

如果算法实现中,没有考虑多点共面的情况,会引起错误。如图8.14所示,求与平面\({{\pi }_{0}}\)上边\({{p}_{0}}{{p}_{1}}\)相邻在面,得到平面\({{\pi }_{1}}\)的夹角最大,但是在平面\({{\pi }_{1}}\)上,存在4个或者大于4个点,图中有7个点\(\{{{p}_{0}},{{p}_{1}},{{p}_{3}},{{p}_{4}},{{p}_{5}},{{p}_{6}},{{p}_{7}}\}\)在平面\({{\pi }_{1}}\)。若随机地选择一个点来更新平面F',可能会选择点\({{p}_{5}}\),也可能是点\({{p}_{7}}\),但是显然点\({{p}_{7}}\)不是凸包的极点,就引发了错误。若选择与边距离最远的点来更新平面F',就会选择点\({{p}_{4}}\),它一定是凸包上的极点,但可能会引发生成的面的边可能会发生交叉的情况。如图8.15所示,考虑边\(\overline{{{p}_{0}}{{p}_{1}}}\)的相邻面,最远点是\({{p}_{4}}\),连接3点得到凸包上的面\(\overline{{{p}_{0}}{{p}_{1}}{{p}_{4}}}\);接着考虑边\(\overline{{{p}_{1}}{{p}_{6}}}\)的相邻面,最远点是\({{p}_{3}}\),连接3点得到凸包上的另外一个三角形面片\(\overline{{{p}_{1}}{{p}_{6}}{{p}_{3}}}\),这就造成了边发生相交的情况。

2015-3-24 11-18-43

图8.14 多个点共面的情况

2015-3-24 11-18-51

图8.15 选择与边距离最远点更新新面,可能产生的问题

这里介绍一种方法,设三维空间上的点集\(R=\{{{q}_{0}},{{q}_{1}},...,{{q}_{m}}\},m\ge 3\)在同一个平面上,把点集\(R\)投影到二维空间上,得到新的点集\(R'=\{{{q}_{0}}',{{q}_{1}}',...,{{q}_{m}}'\}\)。三维点投影到二维空间上,可以采用的方法是:设点集所在的平面的法向量为\(\vec{m}\),移除法向量\(\vec{m}\)的绝对值最大的分量,即对于任意一个点\({{q}_{i}}=({{x}_{i}},{{y}_{i}},{{z}_{i}})\),投影后的点\({{q}_{i}}'\)如等式(8.11)所示。

\[{q_i}' = \left\{ {\begin{array}{*{20}{c}}{({y_i},{z_i}),\left\| {{m_x}} \right\| = \max \{ \left\| {{m_x}} \right\|,\left\| {{m_y}} \right\|,\left\| {{m_z}} \right\|\} }\\{({x_i},{z_i}),\left\| {{m_y}} \right\| = \max \{ \left\| {{m_x}} \right\|,\left\| {{m_y}} \right\|,\left\| {{m_z}} \right\|\} }\\{({x_i},{y_i}),\left\| {{m_z}} \right\| = \max \{ \left\| {{m_x}} \right\|,\left\| {{m_y}} \right\|,\left\| {{m_z}} \right\|\} }\end{array}} \right. \tag{8.11}\]

2015-3-24 11-19-00

图8.16 把凸包上的点,划分为不相交的三角形面

使用二维凸包算法,求出点集\(R'\)的凸包,设求出的二维凸包为\(CH'=\{{{p}_{0}}',{{p}_{1}}',\cdots ,{{p}_{k}}'\}\),那么三维点构成的对应的凸多边形为\(CH=\{{{p}_{0}},{{p}_{1}},\cdots ,{{p}_{k}}\}\),其中,点\({{p}_{i}}'\)是点\({{p}_{i}}\)在二维空间对应的投影点,\(0\le i\le k\)。可以把\(CH\)划分成\(\Delta {{p}_{0}}{{p}_{1}}{{p}_{2}}\),\(\Delta {{p}_{0}}{{p}_{2}}{{p}_{3}}\),...,\(\Delta {{p}_{0}}{{p}_{m-1}}{{p}_{m}}\)\(m-1\)个三角形面片,这种方法生成的三角形面片不会发生相交的情况。以图8.14为例,凸包为\({{p}_{0}}{{p}_{1}}{{p}_{6}}{{p}_{5}}{{p}_{4}}{{p}_{3}}\),就可以划分为\(\Delta {{p}_{0}}{{p}_{1}}{{p}_{6}}\),\(\Delta {{p}_{0}}{{p}_{6}}{{p}_{5}}\),\(\Delta {{p}_{0}}{{p}_{5}}{{p}_{4}}\),\(\Delta {{p}_{0}}{{p}_{4}}{{p}_{3}}\)共4个三角形面,如图8.16所示。

在确定凸包上的初始三角形面片时,也可能有多个点在同一个面上,此时不需将共面点集投影到二维平面上求出子凸包,只需选择与指定边距离最远的点即可,该点一定是凸包的极点,避免了不必要的算法操作。

2.     分母为零的情况

若所求的点\({{p}_{i}}\)在平面\(\pi \)上,则\({{v}_{i}}\)在法向量方向\(\vec{n}\)上的投影\(\vec{n}\cdot {{v}_{i}}=0\),显然,等式(8.7)的分母为零。此时,需要根据分子\(\vec{a}\cdot {{v}_{i}}\)的符号判断该点是否有效。若\(\vec{a}\cdot {{v}_{i}}\ge 0\),说明该点是无效点,直接忽略;若\(\vec{a}\cdot {{v}_{i}}\prec 0\),说明该点是有效的,则夹角最大的点在平面\(\pi \)上。以图8.17所示为例,点\({{p}_{i}},1\le i\le 9\)在平面\(\pi \)上,在求与三角形面\(\Delta {{p}_{0}}{{p}_{1}}{{p}_{2}}\)的边\(\overline{{{p}_{0}}{{p}_{1}}}\)相邻的凸包面时,点\({{p}_{6}},{{p}_{7}},{{p}_{8}},{{p}_{9}}\)是有效点,但是点\({{p}_{0}},{{p}_{1}},{{p}_{2}},{{p}_{3}},{{p}_{4}},{{p}_{5}}\)都是无效点。

2015-3-24 11-19-15

图8.17 分母为零的情况

8.4.3. 算法实现

存储点的数据结构可以表示为\(\{x,y,z,IsOn\}\),其中,\(x,y,z\)分别表示点的\(X,Y,Z\)方向上的坐标,\(IsOn\)是标志位,是一个布尔数,表示指定点是否是凸包上的极点,这个数据在输出结果时会起到作用,对于标志位为\(FALSE\)的点可以删除,保留下标志位为\(TRUE\)的点。起始时,对于任意一个点\(p\in S\)\(p.IsOn\)的值都为\(FALSE\)

前面介绍了边结构可以用\(\{p,q,count\}\)来表示,包括边的两个端点\(p,q\)和相邻面个数的记数器\(count\)。可以采用点的索引来记录边,减少存储空间,因此边结构中的\(p,q\)表示点的索引,用\(V\)表示全局点集,那么边结构\(e\)的端点信息就是\(V[e.p]\)\(V[e.q]\)。使用平衡二叉树来作为存储边结构的容器,用\(T\)表示,边结构的查询、插入操作只需要\(O(\log m)\)的时间,其中,\(m\)表示边的个数。对存储边结构的平衡二叉树\(T\)的查询和插入的表示方法如下所示:

  1. \(T.\operatorname{find}(\{p,q\})\),从\(T\)中检索出以\(V[p],V[q]\)为端点的边结构,并将它返回;
  2. \(T.insert(\{p,q\})\),若\(T\)中存在以\(V[p],V[q]\)为端点的边,则将相应的边结构的\(count\)值加1;否则,创建一个边结构\(\{p,q,1\}\),并把它插入\(T\)中。

三维凸包上的面用三角形面片表示,存储面的数据结构可以表示为\(\{{{v}_{0}},{{v}_{1}},{{v}_{2}}\}\),与边结构类似,为了减少存储内存,\({{v}_{0}},{{v}_{1}},{{v}_{2}}\)表示点的索引,如果给定一个三角形面片\(F\),则它的3个顶点信息就是\(V[F.{{v}_{0}}],V[F.{{v}_{1}}],V[F.{{v}_{2}}]\)。三角形面片的朝向是固定的,我们规定,通过等式\(\vec{n}=(V[F.{{v}_{1}}]-V[F.{{v}_{0}}])\times (V[F.{{v}_{2}}]-V[F.{{v}_{0}}])\),计算出的面片的法向量指向凸包外。

算法中还需要设计一种存储三角形面片的容器,由于我们只需要将三角形面片插入容器首部(尾部)元素、删除容器首部(尾部)元素、遍历操作,所以可以采用链表这种数据结构作为存储三角形面片的容器。它一方向可以用于存储候选三角形面片;另一方面用于存储求出的凸包面片,最终作为结果输出。对链表数据结构\(List\)主要的操作的表示方法如下所示:

  1. \(List.push(O)\),把元素\(O\)插入链表\(List\)内;
  2. \(List.pop()\),从链表\(List\)中弹出一个元素,并将它返回;

下面给出算法的伪码,考虑到了输入的点集存在重合或在同一个平面上,以及上一节描述的多点共面、分母为零的情况。

采用礼物包裹算法求三维空间下点集\(S=\{{{p}_{0}},{{p}_{1}},...,{{p}_{n-1}}\}\)的凸包
\(\text{GIFT }\!\!\_\!\!\text{ WRAPPING}(S)\)
1.    对点集\(S\)中的点从小到大排序;
2.    移除点集\(S\)中重复的点,得到新的点集\(V=\{{{q}_{0}},{{q}_{1}},\cdots ,{{q}_{m-1}}\}\);
3.    if \(\left\| V \right\|\le 3\) || 点集\(V\)共面, then
4.           return \(\emptyset \);
5.    end if;
6.    \(Z = \emptyset \);                                                          // 存储用于输出的面片,链表
7.    \(Q = \emptyset \);                                                          // 存储候选面片,链表
8.    \(T = \emptyset \);                                                          // 存储边结构,平衡二叉树
9.    \(F=\operatorname{INIT}\_FIRST\_FACET(V)\);
10.   \(T.\operatorname{insert}(\{p:F.{{v}_{0}},q:F.{{v}_{1}},count:1\})\);
11.   \(T.\operatorname{insert}(\{p:F.{{v}_{1}},q:F.{{v}_{2}},count:1\})\);
12.   \(T.\operatorname{insert}(\{p:F.{{v}_{1}},q:F.{{v}_{0}},count:1\})\);
13.   \(Q.push(F)\);
14.   \(Z.\operatorname{push}(F)\);
15.   while( \(Q\)非空 )
16.          \(F=Q.\text{pop}()\);
17.          \(\vec{n}=-(V[F.{{v}_{1}}]-V[F.{{v}_{0}}])\times (V[F.{{v}_{2}}]-V[F.{{v}_{0}}])\);
18.          for ( i=0 ; i<3 ; i+=1 )
19.                \(p = F.{v_i},q = F.{v_{(i + 1)\% 3}}\);
20.                 \(E=T.\operatorname{find}(\{p,q\})\);
21.                 if \(E.count\ge 2\), then
22.                        continue;
23.                 end if;
24.                 \(\vec{a}=\vec{n}\times (V[p]-V[q])\);
25.                 \(\{PS,{{k}_{max}}\}=\operatorname{GET}\_EXTREME(V,\vec{n},\vec{a},V[p])\);
26.                 \(\vec{d}=(V[p]-V[q])\times (V[{{k}_{\max }}]-V[q])\);
27.                 if \(\left\| PS \right\|\succ 1\), then
28.                        \(CH=SUB\_CH(\vec{d},E,PS)\);
29.                 else
30.                        \(CH=PS\)
31.                        \(CH.\operatorname{push}\_front(q)\);
32.                        \(CH.\operatorname{push}\_front(p)\);
33.                 end if;
34.                 for each \(k\in CH\)
35.                        \(V[k].IsOn=TRUE\);
36.                 end for;
37.                 \(m=\left\| CH \right\|\)                                                   //凸包\(CH\)的顶点个数
38.                 \({{h}_{0}}=CH[0],{{h}_{1}}=CH[1],{{h}_{2}}=CH[2],h=2\);
39.                 while(\(h\prec m\))
40.                        \(T.\operatorname{insert}(\{{{h}_{2}},{{h}_{1}}\})\);
41.                        \(F'=\{{{h}_{0}},{{h}_{2}},{{h}_{1}}\}\);
42.                        \(Z.\operatorname{push}(F')\);
43.                        \(Q.\operatorname{push}(F')\);
44.                        \(h+=1\);
45.                        if \(h==m-1\), then
46.                               \(T.\operatorname{insert}(\{{{h}_{0}},{{h}_{2}}\})\);
47.                        else
48.                               \({{h}_{1}}={{h}_{2}}\);
49.                               \({{h}_{2}}=CH[h]\);
50.                        end if;
51.                 end while;
52.          end for;
53.   end while;
54.   return \(Z\);

算法第1~5步,排除点集\(S\)中重复的点、以及点集\(S\)在同一个平面上的情况。第6~17步,初始化必要的存储结构,确定在凸包上的初始三角形面片。接着,进入while循环,从\(Q\)中取出一个候选面片\(F\),分别处理面\(F\)的3条边,从\(T\)中查找到相应存储有相应边的边结构\(E\)。如果\(E.count\ge 2\),说明该边已存在两个相邻面,不作任何处理,否则,计算出与其相邻面的\(F'\),将面\(F'\)的另外两条边插入到平衡二叉树\(T\)中,可以采用\(T.insert(\{p,q\})\)方法。

其中,函数\(\operatorname{GET}\_EXTREME(V,\vec{n},\vec{a},p)\)实现的功能是:给定一个凸包上的面\(F\),确定与面片\(F\)上的边\(E\)相邻的三角形面片,\(PS\)表示一个点的索引集合,因为有可能与\(E\)相邻的三角形面片所在的平面上有多个点,即前面叙述的多点共面的情况,\({{k}_{\max }}\)表示\(PS\)中与边\(E\)距离最远的点的索引。设面\(F\)所在的平面为\(\pi \),算法的伪码如下所示:

\(\operatorname{GET}\_EXTREME(V,\vec{n},\vec{a},p)\)
return (\(\{PS,{{k}_{max}}\}\))
1.    \(i=0\);
2.    \(PS = \emptyset \)
3.    while(\(i\prec \left\| V \right\|\))                                                     //确定\(vk,uk,d,{{k}_{\max }},PS\)的初值
4.           \(v=V[i]-p\);
5.           \(vk=\vec{a}\cdot v,uk=\vec{n}\cdot v\);
6.           \(i+=1\);
7.           if (\(uk==0\))&&(\(vk\ge 0\)), then
8.                  continue;
9.           else if \(uk\ne 0\), then
10.                 \(d=-vk/uk\);
11.          end if;
12.          \({{k}_{\max }}=k\);
13.          \(PS.\operatorname{push}(i)\);
14.          break;
15.   end while;
16.   while(\(i\prec \left\| V \right\|\))
17.          \(v=V[i]-p\);
18.          \(tpVk=\vec{a}\cdot v,tpUk=\vec{n}\cdot v\);
19.          \(i+=1\);
20.          if (\(tpVk==0\))&&(\(tpUk==0\)), then                     //待测点\(V[i]\)在边\(E\)
21.                 continue;
22.          else if(\(uk==0\)), then                              //遇到2个或2个以上分母为零的点
23.                 if (\(tpUk==0\))&&(\(tpVk\prec vk\)), then
24.                        \(vk=tpVk,{{k}_{\max }}=i\);
25.                 end if;
26.                 if \(tpVk\prec 0\), then                               //分母为零时,有效点的判定条件
27.                        \(PS.\operatorname{push}(i)\);
28.                 end if;
29.          else if (\(tpUk==0\))&&(\(tpVk\prec 0\)), then           //第1次遇到分母为零的点
30.                 \(vk=tpVk,uk=0\);
31.                 \({{k}_{\max }}=i\);
32.                 \(PS = \emptyset \);
33.                 \(PS.\operatorname{push}(i)\);
34.          else                                                          //分母不为零
35.                 \(tpD=-tpVk/tpUk\);                                 //参见等式(8.7)
36.                 if \(tpD\succ d\), then
37.                        \(vk=tpVk,uk=tpUk\);
38.                        \(d=tpD\);
39.                        \({{k}_{\max }}=i\);
40.                        \(PS = \emptyset \);
41.                        \(PS.\operatorname{push}(i)\);
42.                 else if \(tpD==d\), then                      //存在多点共面的情况
43.                        if \(tpUk\succ uk\), then
44.                               \(vk=tpVk,uk=tpUk\);
45.                               \({{k}_{\max }}=i\);
46.                        end if;
47.                 end if;
48.          end if;
49.   end while;
50.   return (\(\{PS,{{k}_{max}}\}\))

其中,函数\(\operatorname{INIT}\_FIRST\_FACET(V)\)实现的功能是:确定在凸包上的初始三角形面片,算法的伪码如下所示:

\(\operatorname{INIT}\_FIRST\_FACET(V)\)
1.    设\(0\le p\prec \left\| V \right\|\),且\(V[p]\)是点集\(V\)中的最小点的索引;
2.    \(\vec{n}=(1,0,0),\vec{a}=(0,1,0)\);
3.    \(\{PS,q\}=\operatorname{GET}\_EXTREME(V,\vec{n},\vec{a},V[p])\);
4.    \(v=V[q]-V[p]\);
5.    \(vk=\vec{a}\cdot v,uk=\vec{n}\cdot v\);
6.    \(\vec{n}=-vk\cdot \vec{n}+uk\cdot \vec{a}\);
7.    \(\vec{a}=v\times \vec{n}\);
8.    \(\{PS,r\}=\operatorname{GET}\_EXTREME(V,\vec{n},\vec{a},V[q])\);
9.    return \(\{p,r,q\}\);

其中,函数\(SUB\_CH(\vec{d},E,PS)\)实现的功能是:解决多点共面这种退化情况,把三维的点\(\{E.p,E.q\}\bigcup PS\)投影到二维空间上,再根据二维凸包算法计算出凸多边形,最后把结果返回,即\(CH=SUB\_CH(\vec{d},E,PS)\)。可以确定\(\{E.p,E.q\}\subseteq CH\),函数\(SUB\_CH(\vec{d},E,PS)\)需要确保凸多边形\(CH\)是以\(E.p,E.q,{{c}_{0}},{{c}_{1}},\cdots \)这样一个方向。三维点投影到二维,以及求二维凸包算法,在前面都已详细介绍,这里不再赘述算法。

2015-3-24 11-19-32

图8.18 礼物包裹算法过程

最后,以图8.18所示为例来演示礼物包裹算法\(\text{GIFT }\!\!\_\!\!\text{ WRAPPING}(S)\)的流程,可以结合伪码,模拟下算法的过程。图中的点集为\(S=\{{{p}_{0}},\cdots ,{{p}_{8}}\}\),分布在一个\(2\times 2\times 2\)的矩形方块上,点集的坐标如表8.4所示。礼物包裹算法会先确定一个初始的三角形面片,设为\(\Delta {{p}_{0}}{{p}_{4}}{{p}_{3}}\),把面片存入\(Q\)\(Z\)。接着,进入while循环,从\(Q\)中取出第1个面片\(\Delta {{p}_{0}}{{p}_{4}}{{p}_{3}}\),先处理面片上的边\(\overline{{{p}_{0}}{{p}_{4}}}\),产生新的面片\(\Delta {p_0}{p_2}{p_4}\);再处理边\(\overline{{{p}_{4}}{{p}_{3}}}\),此时出现多点共面的情况,点\({{p}_{6}},{{p}_{7}}\)与边\(\overline{{{p}_{4}}{{p}_{3}}}\)共面,且产生的夹角最大,随机启动“多点共面情况”的子算法,最终产生两个面片\(\Delta {{p}_{3}}{{p}_{4}}{{p}_{6}}\)\(\Delta {{p}_{4}}{{p}_{7}}{{p}_{6}}\);最后处理边\(\overline{{{p}_{3}}{{p}_{0}}}\),产生新的面片\(\Delta {p_3}{p_0}{p_1}\)。在处理完第1个三角形面片后,候选面片集\(Q=\{\Delta {{p}_{0}}{{p}_{2}}{{p}_{4}},\Delta {{p}_{3}}{{p}_{4}}{{p}_{6}},\Delta {{p}_{4}}{{p}_{7}}{{p}_{6}},\Delta {{p}_{3}}{{p}_{0}}{{p}_{1}}\}\)。由于\(Q\)非空,第2次进行while循环内的迭代操作,处理面片\(\Delta {p_0}{p_2}{p_4}\),先处理边\(\overline{{{p}_{0}}{{p}_{2}}}\),生成新的面片\(\Delta {p_0}{p_1}{p_2}\);再处理边\(\overline{{{p}_{2}}{{p}_{4}}}\),生成\(\Delta {{p}_{2}}{{p}_{7}}{{p}_{4}}\);最后处理边\(\overline{{{p}_{4}}{{p}_{0}}}\),由于该边的两邻面个数为2,已不存在其它相邻面,不作处理,此时候选面片集\(Q=\{\Delta {{p}_{3}}{{p}_{4}}{{p}_{6}},\Delta {{p}_{4}}{{p}_{7}}{{p}_{6}},\Delta {{p}_{3}}{{p}_{0}}{{p}_{1}},\Delta {{p}_{0}}{{p}_{1}}{{p}_{2}},\Delta {{p}_{2}}{{p}_{7}}{{p}_{4}}\}\)。由于\(Q\)非空,从中取出一个平面,重复上述操作,直至\(Q\)为空,输出三维凸多边体的面片,算法结束。如表1.4所示,由于篇幅限制,只记录了前4个面片的处理结果。

表8.4 点集坐标和算法循环体内的处理结果

2015-3-24 11-41-39

 

spacer

Leave a reply