浏览量:547

8.5. 快速凸包算法

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

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

8.5. 快速凸包算法

8.5.1. 基本思想

快速凸包算法也可以看成是增量算法的一种变种,与随机增量算法相比,它们的不同就在于每次迭代从面的外部点集中选择的点不同。随机增量算法从外部点集中随机的选择一个点,但是快速凸包算法是选择距离最远的一个点,著名的开源代码Qhull[12]、 CGAL[13]都有快速凸包算法的C++实现。

在介绍快速三维凸包算法前,先介绍算法中会频繁使用到的两个基本的几何操作:

1.     给定3个三维空间上的点,确定一个平面。

平面可以用方程\(\vec{n}\cdot P+d=0\)表示,设3个点分别是\(\{{{v}_{0}},{{v}_{1}},{{v}_{2}}\}\),则有\(\vec{n}=({{v}_{1}}-{{v}_{0}})\times ({{v}_{2}}-{{v}_{0}})\)\(d=-\vec{n}\cdot {{v}_{0}}\),叉积的顺序影响法向量的方向,若\(\vec{n}\)是零向量,说明\({{v}_{0}},{{v}_{1}},{{v}_{2}}\)三点共线。

2.     计算三维空间上的点到平面的有符号距离。

设三维点为\(P\),平面为\(\Gamma :\vec{n}\cdot P+d=0\),那么点\(P\)到平面\(\Gamma \)的有符号距离表示为\(dist(P)=(\vec{n}\cdot P+d)/\left\| {\vec{n}} \right\|\)。若\(dist(P)\succ 0\),称点P在平面上方(Above the Plane),;若\(dist(P)\prec 0\),称点P在平面下方(Below the Plane);若\(dist(P)\text{=}0\),称点P在平面上(On the Plane)。

快速凸包算法是基于Beneath Beyond理论,增量算法也同样基于该理论,该理论如下所示,这里只考虑三维空间下的凸包:

\(H\)是一个\({{R}^{3}}\)空间上的凸包,\(p\)\({{R}^{3}}-H\)上的一个点。\(F\)是凸包\(conv(p\bigcup H)\)上的面,当且仅当

  1. \(F\)是凸包\(H\)的一个面且点\(p\)在面\(F\)的下方;
  2. \(F\)不是凸包\(H\)的一个面,\(F\)是由凸包\(H\)的边和点\(p\)构成,点\(p\)在该边相邻的一个面的上方,在该边相邻的另一个面的下方。

2015-3-24 11-43-45

图8.19 凸包\(conv(p\bigcup H)\)

若点\(p\)\(H\)内,则\(conv(p\bigcup H)\)\(H\)重合,显然,结论成立;若点\(p\)\(H\)外,分为两种情况,以图8.19所示为例,\(H\)是由极点\(\{{{p}_{0}},{{p}_{1}},{{p}_{2}},{{p}_{3}},{{p}_{4}},{{p}_{5}}\}\)构成的凸包,\(conv(p\bigcup H)\)是由极点\(\{p,{{p}_{0}},{{p}_{1}},{{p}_{2}},{{p}_{4}},{{p}_{5}}\}\)构成的凸包。对于凸包\(H\)来说,点\(p\)在三角形面片\(\Delta {{p}_{0}}{{p}_{4}}{{p}_{3}}\),\(\Delta {{p}_{0}}{{p}_{3}}{{p}_{5}}\),\(\Delta {{p}_{2}}{{p}_{3}}{{p}_{4}}\),\(\Delta {{p}_{2}}{{p}_{5}}{{p}_{3}}\)的上方,点\(p\)在三角形面片\(\Delta {{p}_{0}}{{p}_{1}}{{p}_{4}}\),\(\Delta {{p}_{0}}{{p}_{5}}{{p}_{1}}\),\(\Delta {{p}_{2}}{{p}_{4}}{{p}_{1}}\),\(\Delta {{p}_{2}}{{p}_{1}}{{p}_{5}}\)的下方,因此在边\(\overline{{{p}_{2}}{{p}_{4}}}\),\(\overline{{{p}_{4}}{{p}_{0}}}\),\(\overline{{{p}_{0}}{{p}_{5}}}\),\(\overline{{{p}_{5}}{{p}_{2}}}\)一侧的面片是在点\(p\)的上方,另一侧的面片在点\(p\)的下方,我们称这样的边为临界边。当一个面\(F\)\(conv(p\bigcup H)\)上时,若点\(p\)在面\(F\)的下方,则\(F\)\(H\)的一个面,否则,它是由点\(p\)与临界边构成的面片。

快速凸包算法的伪码如下所示

用快速凸包算法求三维空间上的点集\(S=\{{{p}_{0}},{{p}_{1}},...,{{p}_{n-1}}\}\)的凸包
1.    初始化一个由四个点构成的四面体;
2.    for 四面体的每个面\(F\)
3.           for 点集中的每个未被分配的点\(p\)
4.                  if \(p\)在平面\(F\)的上方, then
5.                         把\(p\)分配到面\(F\)的外部点集中;
6.                  end if;
7.           end for;
8.    end for;
9.    把外部点集非空的面保存进待定面集\(Q\)中;
10.   for 待定面集\(Q\)中的每个面\(F\)
11.          从\(F\)的外部点集中找到距离面\(F\)最远的点\(p\);
12.          初始化可见面集\(V\),把面\(F\)保存进\(V\)中;
13.          for 可见面集\(V\)中每个面的未被访问过的邻居面\(N\)
14.                 if \(p\)在面\(N\)的上方, then
15.                        把\(N\)加到集合\(V\)中;
16.                 end if;
17.          end for;
18.          把集合\(V\)中每个面的外部点集汇总到一个点集\(L\)中;
19.          集合\(V\)中所有面的临界边,构成一个集合\(H\);
20.          for 集合\(H\)中的每个边界边\(R\)
21.                 连接点\(p\)和临界边\(R\),创建一个新的面,加入到集合新面集\(NS\);
22.                 更新新的面的相邻面;
23.          end for;
24.          for 新面集\(NS\)中的每一个新的面\(F'\)
25.                 for 点集\(L\)中每个未分配的点\(q\)
26.                        if \(q\)在平面\(F'\)的上方,then
27.                               把\(q\)分配到平面\(F'\)的外部点集中;
28.                        end if;
29.                 end for;
30.          end for
31.          若待定面集\(Q\)和可见面集\(V\)的交不为空,则从待定面集\(Q\)中移除它们的交集,并删除可见面集\(V\);
32.          把新面集\(NS\)中外部点集非空的新面\(F'\)添加到待定面集\(Q\)中;
33.   end for;

第1步,初始化一个由四个点构成的四面体。从点集中选择4个点生成初始的四面体,要保证这4个点生成的四面体是非退化的,尽可能选择跨度大的点。比如选择\(X\)分量最小的点和最大的点,\(Y\)最小的点和最大的点,\(Y\)最小的点和最大的点,这样可以使尽可能多的点包含在初始四面体内,达到排除尽可能多的非极点的目的。

第2~8步,存储每个面的数据结构包含面的3个顶点信息、它的外部点集和它的3个相邻面,面的外部点集指在面上方的、未分配的点构成的集合。“未分配”有两层意思:(1)该点不作为当前凸包的顶点,(2)其它面的外部点集不包含该点。若一个点在多个面的上方,则把它随机地分配到任意一个面的外部点集中。因此,若一个面的外部点集不为空,说明它一定不是的凸包\(conv(S)\)的面;相反,若一个面的外部点集为空,也不能保证它一定是凸包\(conv(S)\)的面。对于四面体的每个面\(F\),遍历点集,找到所有在面\(F\)上方的点,保存在面\(F\)的外部点集中,每个面结构都记录有一个外部点集,把外部点集非空的面保存在一个集合中,称这个集合为待定面集。

第9步,把外部点集非空的面保存进待定面集Q中。

第11步,从面\(F\)的外部点集中找到与面\(F\)距离最远的点\(p\),并且把点\(p\)从面\(F\)的外部点集中移除。

第12步,初始化可见面集\(V\),这里的可见面,是相对于点来说的,若点\(p\)在面的上方,则称该面是相对于点\(p\)的可见面。初始情况,可见面集\(V\)中只有一个面,就是面\(F\)

第13~17步,可见面集\(V\)中每个面中的未被访问过的邻居面\(N\),如果点\(p\)在面\(N\)的上方,把\(N\)加到集合\(V\)中。可见面集\(V\)的初始情况只有一个面\(F\),从面\(F\)向周围的面遍历,将所有对点\(p\)可见的面保存在可见集\(V\)中。以图8.19所示为例,设点\(p\)是面\(\Delta {p_0}{p_4}{p_3}\)的外部点集中距离最远的点,经过算法第13~17步,得到可见面集为\(V = \{ \Delta {p_0}{p_4}{p_3},\Delta {p_0}{p_3}{p_5},\Delta {p_2}{p_3}{p_4},\)\(\Delta {p_2}{p_5}{p_3}\} \)

第18步,把可见面集合\(V\)中每个面的外部点集汇总到一个点集\(L\)中,因为可见面集中的面需要从当前凸包中删除。

第19步,可见面集合\(V\)中所有可见面的临界边,构成一个集合\(V\),若一条边的两个相邻面中,一个面是可见的,另外一个面是不可见的,那么该边就是临界边。

第20~23步,连接点\(p\)和集合\(H\)中的边界边,创建出新的面,更新新的面的相邻面。

第24~30步,对于每个新的面\(F'\),遍历点集\(L\),如果对于点集\(L\)中未分配的点\(q\),它在在\(F'\)的上方,则把它添加到面\(F'\)的外部点集中。

第31步,若待定面集\(Q\)和可见面集\(V\)的交不为空,则从待定面集\(Q\)中移除它们的交集,用一个直观数学表示就是\(Q=Q-Q\bigcap V\),删除可见面集中保存的面。移除交集的目的是为了保证已经删除的可见面,不会在待定面集\(Q\)中,导致后面的重复处理。

第32~33步,对于每个新的面\(F'\),若它的外部点集非空,则把它添加到待定面集Q中,进行下次的迭代。

这里考虑一个问题,在快速凸包算法第2~8步和第24~30步中,若一个点在多个面的上方,则把它随机地分配到任意一个面的外部点集中,是否会造成错误,即点集\(S\)中存在一个极点最终不会是凸包\(conv(S)\)上的顶点。我们可以证明下这条结论:

若一个极点在多个面的上方,把它随机的分配给任意一个可见面,该点最终都会作为凸包的顶点。

我们用反证法来证明,若一个极点\(p\)未分配给任意一个面的外部点集,则它一定不会出现在凸包的顶点上,这就与它是极点相矛盾,因此它一定是初始化的四面体的某个面的外部点集。令\(p\)\(q\)在相同面的外部点集中,设\(p\)在面的上方但是不在由\(q\)与边生成的新的面的外部点集中,因此\(p\)在可见面上,但是在由\(q\)与临界边生成的新的面的下方,这表明\(p\)在凸包的内部,则它不是一个极点,发生矛盾。综上,我们可以得到结论。

快速凸包算法首先会初始化一个4个点的凸包,然后把点集分割成各个面的外部点集,挑选距离面最远的点进行处理。任意一个极点都不会因为分割为外部点集导致被丢弃。根据理论Beneath Beyond理论(1),一个面要么保持不变,要么根据理论Beneath Beyond理论(2)创建新的面,算法能产生已处理点的凸包,采用这种方法不断的加入新的点进行处理,直至点集上的点都处理结束为止。

算法的平均复杂度为\(O(n\log n)\),最坏情况下的复杂度为\(O({{n}^{2}})\)

8.5.2. 算法实现

为了方便后面的叙述,规定待定面集用\(Q\)表示,可见面集用\(V\)表示,临界边集用\(H\)表示。

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

存储临界边的数据结构可以表示为\(\{{{v}_{0}},{{v}_{1}},{{f}_{0}},{{f}_{1}}\}\),其中,\({{v}_{0}},{{v}_{1}}\)表示边的两个端点,我们规定临界边的方向为由\({{v}_{0}}\)指向\({{v}_{1}}\)\({{f}_{0}},{{f}_{1}}\)表示边的两个相邻面,因为它是临界边,所以两个相邻面中,必然有一个面是可见面,另一个面是不可见面,我们规定\({{f}_{0}}\)表示可见面,\({{f}_{1}}\)表示不可见面。算法中,最远点与它的临界边,构造新的面。

2015-3-24 11-44-01

图8.20 采用宽度优先搜索法遍历凸包

存储面的数据结构可以表示为\(\{{{v}_{0}},{{v}_{1}},{{v}_{2}},{{f}_{0}},{{f}_{1}},{{f}_{2}},OS,flag\}\),因为三维空间上凸包的每个面用三角形表示。其中,\({{v}_{0}},{{v}_{1}},{{v}_{2}}\)表示面片的三个顶点,顶点的顺序影响面的法向量的方向,我们规定采用等式\(\vec{n}=({{v}_{1}}-{{v}_{0}})\times ({{v}_{2}}-{{v}_{0}})\)计算出来的法向量指向凸包的外侧;\({{f}_{0}},{{f}_{1}},{{f}_{2}}\)表示与面片的三条边相邻的三个面,我们规定\({{f}_{0}}\)是边\(\overline{{{v}_{0}}{{v}_{1}}}\)的相邻面,\({{f}_{1}}\)是边\(\overline{{{v}_{1}}{{v}_{2}}}\)的相邻面, \({{f}_{2}}\)是边\(\overline{{{v}_{2}}{{v}_{0}}}\)的相邻面;\(OS\)表示面的外部点集;\(flag\)是标志位,通过标志位来判定面是否被访问过,以及若面被访问过,是否为可见面,通过标志位来确定哪些边属于临界边。采用这种面数据结构,给定一个面,就可以通过深度优先搜索或宽度优先搜索遍历凸包上的每一个面,如图8.20所示,采用宽度优先搜索遍历凸包,三角形上的数字就表示搜索的层数,这种结构对搜索某个点在凸包上的可见面的效率特别高。

存储点集和面集的数据结构,可以是链表,因为元素的插入、删除操作为\(O(1)\),算法中只需要对点集、面集进行元素的插入、删除、遍历操作。用\(List\)表示链表,对链表的操作可以表示为:

  1. \(List.\operatorname{empty}()\),若\(List\)为空,则返回\(TRUE\);否则,返回\(FALSE\)
  2. \(List.\operatorname{erase}(p)\),从\(List\)中删除\(p\)
  3. \(List.front()\)\(List.back()\)),返回\(List\)首部(尾部)的元素;
  4. \(List.pop\_front()\)\(List.pop\_back()\)),删除\(List\)首部(尾部)的元素;
  5. \(List.\operatorname{push}\_front(p)\)\(List.\operatorname{push}\_back(p)\)),把元素\(p\)插入到链表的首部(尾部);
  6. \(List.\operatorname{clear}()\),将链表中的元素清空;
  7. \(\operatorname{SPLICE}(dest,src)\),把链表\(src\)中的元素合并到链表\(dest\)的尾部,并将链表\(src\)清空;
  8. 给定链表的一个元素\(it\),它的上一个元素表示为\(it=List.\operatorname{last}(it)\),它的下一个元素表示为\(it=List.next(it)\)\(List.\operatorname{begin}()\)是链表起始的元素,但是\(List.end()\)是紧随着链表最后一个元素的标志位。链表的遍历可以表示为:

for (\(it=List.begin()\) ; \(it!=List.end()\) ; \(it=List.\operatorname{next}(it)\))       //顺序遍历
.................
end for;

for (\(it=List.rbegin()\) ; \(it!=List.rend()\) ; \(it=List.\operatorname{rnext}(it)\))   //逆序遍历
.................
end for;

每次迭代,所有的临界边会围成一圈,形成一个环,例如若只有三条临界边,它们一定是由三个点构成的环,设三个点是\(\{{{v}_{0}},{{v}_{1}},{{v}_{2}}\}\),则三条边可以表示为\(\overline{{{v}_{0}}{{v}_{1}}}\)\(\overline{{{v}_{1}}{{v}_{2}}}\)\(\overline{{{v}_{2}}{{v}_{0}}}\)。连接最远点\(p\)和临界边,构造新的平面后,需要更新新的平面的相邻面,所以临界边集合\(H\)采用的存储结构必需满足快速查询的要求。显然,链表不适合快速查找,可以采用平衡二叉树来存储临界边,使得元素的插入、删除和查询操作在\(O(logm)\)时间内完成,其中,\(m\)表示临界边的个数。插入平衡二叉树的元素由两部分组成——{关键字,数据},元素的关键字域作为二叉树中元素间对比用,数据域作为元素的数据部分。对于给定两个三维点\(p\)\(q\),若有\(p\prec q\),当且仅当\(({{p}_{x}}\prec {{q}_{x}})\vee ({{p}_{x}}={{q}_{x}}\wedge {{p}_{y}}\prec {{q}_{y}})\vee ({{p}_{x}}={{q}_{x}}\wedge {{p}_{y}}={{q}_{y}}\wedge {{p}_{z}}\prec {{q}_{z}})\);显然,当满足\((not\text{ }p\prec q)\wedge \)\((not\text{ }q\succ p)\)时,\(p=q\)

规定,存储在临界边集\(H\)中的元素的关键值域是一个三维坐标点,数据域是临界边结构。对平衡二叉树\(Tree\)的操作可以表示为:

  1. \(Tree.\operatorname{insert}(\{key,value\})\),将元素\(\{key,value\}\)插入到\(Tree\)中,若\(Tree\)中已经存在关键域为\(key\)的元素,则插入失败;
  2. \(Tree.\operatorname{empty}()\),若\(Tree\)为空,则返回\(TRUE\);否则,返回\(FALSE\)
  3. \(Tree.find(key)\),从\(Tree\)中检索出关键值域为\(key\)的元素,并返回它的数据域;
  4. \(Tree.erase(key)\),将关键值域为\(key\)的元素从\(Tree\)中删除。

快速凸包算法的伪码如下所示:

用快速凸包算法求三维空间上的点集\(S=\{{{p}_{0}},{{p}_{1}},...,{{p}_{n-1}}\}\)的凸包
1.    \(NS = \emptyset \);                                                               //链表结构,存储创建的新面
2.    \(Q = \emptyset \);                                                          //链表结构,存储待定面
3.    if  !\(\operatorname{INIT}\_TETRAHEDRON(S,NS)\), then //初始化四面体
4.           return \(FALSE\);
5.    end if;
6.    \(\operatorname{PARTITION}(Q,NS,S,R)\);                                  //\(R\)是一个面结构
7.    \(V = \emptyset \);                                                          //链表结构,存储可见面
8.    \(H = \emptyset \);                                                          //平衡二叉树结构,存储临界边
9.    \(L = \emptyset \);                                                          //链表结构,存储点集
10.   while( \(Q\)非空 )
11.          \(F=Q.\operatorname{front}()\);                               //此时,\(Q,V,H,L\)均为空
12.          \(Q.pop\_\operatorname{front}()\);
13.          \(p=\operatorname{FIND}\_FURTHEST(F)\);           //从面\(F\)的外部点集中找到与面\(F\)距离最远的点
14.          \((F.OS).\operatorname{erase}(p)\);
//从凸包上找到点\(p\)的可见面,存入\(V\)中,并将临界边插入到\(H\)
15.          \(F.flag=VISITED\);
16.          \(V.\operatorname{push}\_back(F)\);
17.          for (\(it=V.begin()\) ; \(it!=V.end()\) ; \(it=V.next(it)\))
18.                 for (\(i=0\) ; \(i\prec 3\) ; \(i=i+1\))
19.                        \(N=F.{{f}_{i}}\);
20.                        if \(N.flag==NOT\_VISITED\), then
21.                               \(N.flag==VISITED\);
22.                               \(\vec{n}=(N.{{v}_{1}}-N.{{v}_{0}})\times (N.{{v}_{2}}-N.{{v}_{0}}),d=-\vec{n}\cdot N.{{v}_{0}}\); //确定面\(N\)所在的平面
23.                               if \(\vec{n}\cdot p+d\succ 0\), then
24.                                      \(V.\operatorname{push}\_back(N)\);
25.                               else
26.                                      \(N.flag=BOUNDARY\);
27.                                      设\(E\)是一个临界边结构;
28.                                      \(E.{{f}_{0}}=it,E.{{f}_{1}}=N\);
29.                             \(E.{v_0} = it.{v_i},E.{v_1} = it.{v_{(i + 1)\% 3}}\);
30.                                      \(H.insert(\{E.{{v}_{0}},E\})\);
31.                               end if;
32.                        else if \(N.flag==BOUNDARY\), then
33.                               设\(E\)是一个临界边结构;
34.                               \(E.{{f}_{0}}=it,E.{{f}_{1}}=N\);
35.                        \(E.{v_0} = it.{v_i},E.{v_1} = it.{v_{(i + 1)\% 3}}\);
36.                               \(H.insert(\{E.{{v}_{0}},E\})\);
37.                        end if;
38.                 end for;
39.          end for;
//把集合\(V\)中每个面的外部点集汇总到一个点集\(L\)中,若待定面集\(Q\)和可见面集\(V\)的交不为空,则从待定面集\(Q\)中移除它们的交集
41.          for (\(it=V.begin()\) ; \(it!=V.end()\) ; \(it=V.next(it)\))
42.                 if not \((it.OS).{\mathop{\rm empty}\nolimits} ()\), then
43.                        \(\operatorname{SPLICE}(L,it.OS)\);
44.                 end if;
45.                 if \(it\in Q\), then
46.                        \(Q.\operatorname{erase}(it)\);
47.                 end if;
48.          end for;
//连接临界边和点\(p\),构造新的面,并将新的面存入\(NS\),更新受影响面的相邻面
49.          从\(H\)中获取一个元素\(E\);
50.          while(not \(H.\operatorname{empty}()\))
51.                 \((E.{{f}_{1}}).flag=NOT\_VISITED\);
52.                 \(F'.{{v}_{0}}=p,F'.{{v}_{1}}=E.{{v}_{0}},F'.{{v}_{2}}=E.{{v}_{1}}\);
53.                 \(F'.{{f}_{1}}=E.{{f}_{1}}\);
54.                 \(NS.\operatorname{push}\_back(F')\);
55.                 for( \(i=0\) ; \(i\prec 3\) ; \(i=i+1\) )
56.                        if \((E.{{f}_{1}}).{{f}_{i}}==E.{{f}_{0}}\), then
57.                               \((E.{{f}_{1}}).{{f}_{i}}==F'\);
58.                               break;
59.                        end if;
60.                 end for;
61.                 \(H.erase(E.{{v}_{0}})\);
62.                 \(E=H.find(E.{{v}_{1}})\)
63.          end while;
64. \(lit=NS.last(NS.end())\);
65.          \(it=V.\operatorname{begin}()\);
66.          while(\(it!=NS.end()\))
67.                 \(it.{{f}_{0}}=lit\);
68.                 \(lit.{{f}_{2}}=it\);
69.                 \(lit=it\);
70.                 \(it=NS.\operatorname{next}(it)\);
71.          end while;
72.          \(\operatorname{PARTITION}(Q,NS,L,R)\);
73.          \(V.\operatorname{clear}(),K.\operatorname{clear}(),L.\operatorname{clear}()\);
74.   end while;
75. 利用面结构\(R\),通过宽度优先搜索算法,得到三维凸包;
76.   return \(TRUE\);

2015-3-24 11-44-17

图8.21 初始四面体

上述伪码中,函数\(\operatorname{INIT}\_TETRAHEDRON(S,NS)\)用于初始化四面体,初始的四面体如图8.21所示,对于任意一个面\(f\),保证得到的法向量\(\vec{n}=(f.{{v}_{1}}-f.{{v}_{0}})\times (f.{{v}_{2}}-f.{{v}_{0}})\)的方向为垂直于面指向凸包外侧,算法伪码如下所示。

\(\operatorname{INIT}\_TETRAHEDRON(S,NS)\)          //\(S\)表示点集,\(NS\)表示面集
1.    if 点集\(S\)共线,then
2.           return \(FALSE\);
3.    end if;
4.    找到点集上\(X\)分量最小的和最大的点\({{p}_{0}},{{p}_{1}}\)\(Y\)分量最小的点\({{p}_{2}}\);
5.    if \({{p}_{0}},{{p}_{1}},{{p}_{2}}\)共线, then
6.           随机找到3个不共线的点\({{p}_{0}},{{p}_{1}},{{p}_{2}}\);
7.    end if;
8.    \(\vec{n}=({{p}_{1}}-{{p}_{0}})\times ({{p}_{2}}-{{p}_{0}}),d=-\vec{n}\cdot {{p}_{0}}\);
9.    \(i{{t}_{\min }}=i{{t}_{\max }}=it=S.\operatorname{begin}()\);
10.   \(dis{{t}_{\min }}=dis{{t}_{\max }}=\vec{n}\cdot it+d\);
11.   \(it=S.\operatorname{next}(S)\);
12.   while(\(it!=S.\operatorname{end}()\))
13.          \(dist=\vec{n}\cdot it+d\);
14.          if \(dist\prec dis{{t}_{\min }}\), then
15.                 \(dis{{t}_{\min }}=dist\);
16.                 \(i{{t}_{\min }}=it\);
17.          end if;
18.          if \(dist\succ dis{{t}_{\max }}\), then
19.                 \(dis{{t}_{\max }}=dist\);
20.                 \(i{{t}_{\max }}=it\);
21.          end if;
22.   end while;
23.   \({{p}_{3}}=i{{t}_{\max }}\);
24.   if \({{p}_{0}},{{p}_{1}},{{p}_{2}},{{p}_{3}}\)共面, then
25.          \(\operatorname{SWAP}({{p}_{0}},{{p}_{2}})\);           //交换点\({{p}_{0}}\)与点\({{p}_{2}}\)的值
26.          \({{p}_{3}}=i{{t}_{\min }}\);
27.          if \({{p}_{0}},{{p}_{1}},{{p}_{2}},{{p}_{3}}\)共面, then
28.                 return \(FALSE\);
29.          end if;
30.   end if;
31.   \({{f}_{0}}.{{v}_{0}}={{p}_{0}},{{f}_{0}}.{{v}_{1}}={{p}_{2}},{{f}_{0}}.{{v}_{2}}={{p}_{1}}\);    //4个面的顶点
32.   \({{f}_{1}}.{{v}_{0}}={{p}_{0}},{{f}_{1}}.{{v}_{1}}={{p}_{1}},{{f}_{1}}.{{v}_{2}}={{p}_{3}}\);
33.   \({{f}_{2}}.{{v}_{0}}={{p}_{0}},{{f}_{2}}.{{v}_{1}}={{p}_{3}},{{f}_{2}}.{{v}_{2}}={{p}_{2}}\);
34.   \({{f}_{3}}.{{v}_{0}}={{p}_{1}},{{f}_{2}}.{{v}_{1}}={{p}_{2}},{{f}_{2}}.{{v}_{2}}={{p}_{3}}\);
35.   \({{f}_{0}}.{{f}_{0}}={{f}_{2}},{{f}_{0}}.{{f}_{1}}={{f}_{3}},{{f}_{0}}.{{f}_{2}}={{f}_{1}}\); //4个面的相邻面
36.   \({{f}_{1}}.{{f}_{0}}={{f}_{0}},{{f}_{1}}.{{f}_{1}}={{f}_{3}},{{f}_{1}}.{{f}_{2}}={{f}_{2}}\);
37.   \({{f}_{2}}.{{f}_{0}}={{f}_{1}},{{f}_{2}}.{{f}_{1}}={{f}_{3}},{{f}_{2}}.{{f}_{2}}={{f}_{0}}\);
38.   \({{f}_{3}}.{{f}_{0}}={{f}_{0}},{{f}_{3}}.{{f}_{1}}={{f}_{2}},{{f}_{3}}.{{f}_{2}}={{f}_{1}}\);
39.   \(S.erase({{p}_{0}})\);\(S.erase({{p}_{1}})\);
40.   \(S.erase({{p}_{2}})\);\(S.erase({{p}_{3}})\);
41.   \(NS.\operatorname{push}\_back({{f}_{0}})\);\(NS.\operatorname{push}\_back({{f}_{1}})\);
42.   \(NS.\operatorname{push}\_back({{f}_{2}})\);\(NS.\operatorname{push}\_back({{f}_{3}})\);

在函数\(\operatorname{PARTITION}(Q,NS,L,R)\)中,\(Q\)表示待定面集,\(NS\)表示新创建的面集,\(L\)表示集点,\(R\)是一个面结构,通过它可以获取凸包上的所有面片,函数实现的功能是:若对于每个点\(p\in L\),存在一个面\(F\in NS\),则把点\(p\)加入到面\(F\)的外部点集中,若不存在这样的面,说明点\(p\)在当前凸包内部。

\(\operatorname{PARTITION}(Q,NS,L,R)\)
1.    for (\(fit=NS.begin()\) ; \(fit!=NS.end()\) ; \(fit=NS.next(fit)\))
2.           \(\vec{n}=(fit.{{p}_{1}}-fit.{{p}_{0}})\times (fit.{{p}_{2}}-fit.{{p}_{0}}),d=-\vec{n}\cdot fit.{{p}_{0}}\);
3.           for (\(vit=L.begin()\) ; \(vit!=L.end()\) ; \(vit=L.next(fit)\))
4.                  if \(\vec{n}\cdot vit+d\succ 0\), then
5.                         \(fit.OS.\operatorname{push}\_back(vit)\);
6.                         \(L.\operatorname{erase}(vit)\);
7.                  end if;
8.           end for;
9.    end for;
10.   for (\(fit=NS.begin()\) ; \(fit!=NS.end()\) ; \(fit=NS.next(fit)\))
11.          if not \(fit.OS.\operatorname{empty}()\), then
12.                 \(Q.\operatorname{push}\_back(fit)\);
13.          else
14.                 \(R=fit\);
15.          end if;
16.   end for;

算法结束。

spacer

Leave a reply