求解平面点集点间的点集凸包问题

求解平面点集点间的点集凸包问题

1简单有形凸包算法凸包问题是计算几何的基本问题。凸包算法广泛应用于计算机绘图、数据处理、自动设计、识别等领域。二维凸包问题可以分为基于多边形和基于平面点集两类凸包问题。简单多边形的凸包是指包含简单多边形的最小凸多边形,平面点集的凸包是指包含平面点集内所有点并且顶点属于平面点集的最小凸多边形。简单多边形凸包问题是平面点集凸包问题的一个特例。自1972年Sklansky最早提出具有线性时间复杂度的简单多边形凸包算法以来,涌现了很多凸包算法。Yao于1981年证明了Ω(nlogn)是所有平面点集凸包算法的时间复杂度下限,Kirkpatrick和Seidel于1986年证明了在同时考虑输入和输出的情况下,平面点集凸包算法的时间复杂度下限是Ω(nlogh),并提出了基于分枝定界法的算法复杂度为O(nlogh)的平面点集凸包算法。尽管Sklansky最先提出简单多边形的凸包算法,遗憾的是该算法是有缺陷的。McCallum于1979年提出了第一个正确的简单多边形凸包算法,而Melkman于1987年提出了一个实时的可以用于多边线的凸包算法。在国内,对于凸包算法的研究也很活跃。根据凸多边形的构成特性,本文先提出简单多边形凸包的新算法,进而扩展提出平面点集凸包的新算法。新的点集凸包算法的时间复杂度达到O(nlogh),其中,n为平面点集的点数,h为平面点集凸包的顶点数。2平面点集的确定本节给出算法中用到的基本概念。定义1设pi∈R2或R3,(i=1,2,…,n),设点集S={s|s=∑i=1nλipi‚∑i=1nλi≡1‚λ≥0‚λ∈R‚i=1‚2‚⋯‚n}(1)称包含S的最小凸多边形为S的凸包。定义2设Pi(xi,yi),i=1,2,…,n,是给定多边形的n个顶点,若对任意i,j,i≠j,i,j=1,2,…,n,线段PiPi+1与PjPj+1或是相邻且相交于一端点或是不相交,则称该多边形为简单多边形。定义3设A=(xa,ya)、B=(xb,yb)和C=(xc,yc)为XY平面内任意不同的三点,记L(AB)为A指向B的有向直线,用矢量叉积的方法判断点C与L(AB)的关系的侧向判别函数为S(A,B,C)=(xb-xa)×(yc-yb)-(xc-xb)×(yc-ya)(2)(1)如果S>0,点C在L(AB)的左侧;(2)如果S=0,点C在L(AB)上;(3)如果S<0,点C在L(AB)的右侧。定义4不论是简单多边形还是平面点集,称其凸包上的所有顶点构成的点列为凸包点列。沿一定的方向(顺时针或逆时针)连接而成的凸包点列是有序凸包点列。定义5设简单多边形的顶点列为Pi(xi,yi),i=1,2,…,n。设xmax=max(xi),xmin=min(xi),ymax=max(yi),ymin=min(yi),记PLD=(xmin,ymin),PLU=(xmin,ymax),PRD=(xmax,ymin),PRU=(xmax,xmax)为简单多边形的四个角点,称以PLD、PLU、PRU和PRD为顶点的矩形为简单多边形的矩形包围盒。定义简单多边形的8个极值点如下:(a)M1为x坐标等于xmin的点中y坐标最大的点;M8为x坐标等于xmin的点中y坐标最小的点;(b)M4为x坐标等于xmax的点中y坐标最大的点;M5为x坐标等于xmax的点中y坐标最小的点;(c)M6为y坐标等于ymin的点中x坐标最大的点;M7为y坐标等于ymin的点中x坐标最小的点;(d)M3为y坐标等于ymax的点中x坐标最大的点;M2为y坐标等于ymax的点中x坐标最小的点。以上四个角点、矩形包围盒和8个极值点的定义同样适用于平面点集。线段M1M2、M3M4、M5M6及M7M8将平面点集的矩形包围盒划分成五个子区域,见图1。落入子区域I、II、III和IV内的点可能成为凸包上的点。3凸包顶位置的划分为叙述方便,以下简称简单多边形为多边形。由点集凸包和极值点的定义知,Mi(i=1,2,…,8)必为多边形凸包的顶点,并且线段M8M1,M2M3,M4M5及M6M7为凸包边。因此,只需求解Mi与Mi+1(i=1,3,5,7)间的部分有序凸包点列即可。首先,约定记法并给出相关定义。若V为点,记V(x)、V(y)分别为点V的x坐标、y坐标;若A、B为相异两点,记R(AB)为由点A发出的指向点B的射线,记L(AB)为经过点A、B的有向直线(A指向B)。设Mi与Mi+1间的有序凸包点列(按顺时针方向连接)为H={H1,H2,…,Hr}(r≥2),H1=Mi,Hr=Mi+1。由定义4和定义5知,只有在MiMi+1连线(Mi指向Mi+1,i=1,3,5,7)左侧的点可能成为有序凸包点列中的点。又由于Mi与Mi+1之间的点列取自多边形,故只需考虑点列中在MiMi+1连线左侧的凸顶点,称满足以上条件的点为候选点。由于Mi与Mi+1(i=1,3,5,7)之间的点列的拓扑结构是一样的,故以求解M1到M2之间的有序凸包点列H为例进行说明。对于其他相邻极值点间的有序凸包点列可以类似求解。设M1与M2之间的点列为Q,若Q为空集,则H={M1,M2}。若Q非空,则令H={M1},r=1,再遍历Q中各点;若当前点为候选点,则进一步判断其能否成为凸包顶点。设V为当前候选点,若V是Q中的第一个候选点,则直接将V加到M1之后。称Q中第一个候选点以外的候选点为后续候选点,则判断后续候选点是否为凸包顶点以及如何更新H是算法的关键。下面,给出本文算法中用到的重要概念。设H={H1,H2,…,Hr}(r≥2)为当前有序凸包点列,则H中的点有如下有序性质:若m<n,则Hm(x)<Hn(x),Hm(y)<Hn(y)。称顺序连接H中的点及M2的连线为活动线,它将多边形矩形包围盒中对应极值点间的子区域IV划分为左右两部分,称子区域IV中在活动线左侧的部分(不含边界)为活动区。R(HjHj+1)(j=1,2,…,r-1)将活动区分割成许多与Hj(j=1,2,…,r)对应的活动分区。图2中的阴影部分即是与Hj对应的活动分区,其中的点在R(Hj-1Hj)的右侧,并在R(HjHj+1)的左侧或其上,称这样的活动分区为Hj的活动分区,记为AR(Hj)。称R(Hj-1Hj)、R(HjHj+1)分别为AR(Hj)的前、后边界。特别地,称在R(H1H2)左侧或其上的活动分区为H1的活动分区,记为AR(H1);称在R(Hr-1Hr)右侧并在R(HrM2)左侧的活动分区为Hr的活动分区,记为AR(Hr)。若AR(Hj)内的点不在其边界上,则该点具有如下性质:在其他分区两边界的同一侧,同在左侧或同在右侧。划分活动分区的目的就是为了便于更新当前有序凸包点列H。活动分区的定义及活动分区内点的性质为查找候选点所在的活动分区提供了依据。接下来,给出凸多边形的一个重要特性:沿一定方向(顺时针或逆时针)遍历凸多边形的边,凸多边形的内部位于遍历边的同一侧,我们称这一特性为凸多边形的边同侧特性。若顺时针遍历凸包的边,则凸包的内部在边的右侧。本文就是依据凸多边形的边同侧特性来判断后续候选点能否成为凸包顶点并更新H的。当后续候选点V位于活动区,即V位于活动线的左侧时,凸包的边同侧特性被破坏,说明V是新的凸包顶点,需要更新H。更新H的过程如下:先确定候选点V所在的活动分区,然后依据V所在的活动分区和候选点的位置信息,去除当前有序凸包点列中不再是凸包顶点的点,再把V放入到H中的合适位置使得凸包的边同侧特性得以满足。上述构造当前有序凸包点列的过程严格保证了凸包的边同侧特性,故该构造过程是正确的。以下是求解M1到M2之间的有序凸包点列H的算法。算法3.1SQRT_CHARRAY_SP(M1,M2,Q)Step1.若M1=M2,则H={M1},结束;否则,若Q=∅,则H={M1,M2},结束;若Q≠∅,设Q={Q1,…,Qs},s≥1,令H={M1},r=1,i=1,V=Q1,继续;Step2.若IsCandidate(V)=1,则DEAL_CANDIDATE_SP(V,H);否则,继续;Step3.i=i+1,若i≤s,则V=Qi,转Step2;否则,继续;Step4.r=r+1,Hr=M2。算法3.1中的IsCandidate(V)是判断点V是否为候选点的函数,设遍历多边形的方向为顺时针,V1和V2分别为V的前、后相邻顶点,依据候选点的定义:若S(M1,M2,V)>0且S(V1,V,V2)<0,则IsCandidate(V)=1;否则,IsCandidate(V)=0。DEAL_CANDIDATE_SP(V,H)是对候选点V进行处理并对H进行必要的更新的过程(见算法3.2)。算法3.2DEAL_CANDIDATE_SP(V,H)Step1.若r=1,则r=r+1,Hr=V,退出;否则,继续;Step2.若INAVR_SP(V,H)=0,则H不变,退出;否则,j=INAVR_SP(V,H),继续;Step3.若V(x)≥Hr(x),则r=j+1,Hr=V,退出;否则,继续;Step4.若S(V,M2,Hr)>0,则令W=Hr,r=j+1,Hr=V,r=r+1,Hr=W,退出;否则,r=j+1,Hr=V,退出。其中,INAVR_SP(V,H)过程判断V是否在活动区内,若是,则确定其所在的活动分区(见算法3.3)。依据V与活动线的位置关系进行判断,若V位于活动线的左侧,则V在活动区内。设j=INAVR_SP(V,H):若j=0,则点V不在活动区内;否则,AR(Hj)就是V所在活动分区。算法3.3INAVR_SP(V)Step1.若V(x)=Hr(x),转Step2;若Hr(x)<V(x)<M2(x),转Step3;否则,转Step5;Step2.若V(y)>Hr(y),转Step4;否则,返回0;Step3.若S(Hr-1,Hr,V)≥0,转Step4;若S(Hr-1,Hr,V)<0且S(Hr,M2,V)>0,返回r;否则,返回0;Step4.若S(H1,H2,V)≥0,则返回1;若S(Hr-1,Hr,V)<0,则返回r;否则,l=2,u=r-1,j=FIND_AR(V,l,u);Step5.若S(Hr-1,Hr,V)>0,转Step4;否则,返回0。算法3.3中的FIND_AR(V,l,u,H)是查找V所在活动分区的函数:在Hl与Hu之间查找Hj,使得S(Hj-1,Hj,V)<0且S(Hj,Hj+1,V)≥0,返回j,AR(Hj)就是V所在活动分区(特别地,当S(H1,H2,V)≥0时,j=1)。注意,当V(x)<Hr(x)且V在L(Hr-1Hr)之上或右侧时,由于点列Q取自简单多边形,V不可能在活动区内,否则,V与Hr的连线就与多边形的边相交,这与定义2矛盾,故在Step5中返回0。求解多边形凸包的完整算法如下:先找到多边形的极值点Mi(i=1,2,…,8),再求解得出Mi到Mi+1(i=1,3,5,7)之间的有序凸包点列,最后将Mi到Mi+1(i=1,3,5,7)之间的有序凸包点列顺序连接起来并去掉冗余的Mi点就得出多边形的凸包。下面,将新算法与已有算法进行比较。文中的算法是针对有序简单多边形的,若采用本文算法,依据活动区的定义,可以排除不在HrMi+1连线左侧的候选点,且只要用到算法3.2的Step1、Step2和Step3。而文中的算法在前瞻过程中不能排除这样的点,只能在回溯过程中加以去除,故本文算法可以减少不必要的判断。而且,依据活动分区内的点的性质,对于候选点所在活动分区的查找可以使用二分查找法,而文中的回溯过程是顺序进行的,假设当前有序凸包点列中有h个点,则文的回溯过程需要O(h)时间,而本文算法的回溯过程需要O(logh)时间。文与本文算法依据的原理相同,而本文算法要优于其算法。其一,本文算法先淘汰了非候选点的遍历点,而文的算法是求得可能成为凸包顶点的顶点列后,再淘汰不属于图1中子区域I、II、III和IV中的点,这样会增加对明显不能成为凸包顶点的点的多余判断,降低效率。其二,文的算法存在问题。其算法中的Step6提出,若当前顶点在已形成的凸包顶点栈中栈顶两元素连线之右,且其y坐标大于栈顶元素的y坐标,就将其加入凸包顶点栈内。如图3所示,当算法执行到点V,按照Step6,点V的y坐标大于栈顶元素Hr的y坐标,应加入点V到H中,再到Step7,发现点列中点遍历完毕,转Step9,加入M2,结束。从图3中可以看出,在加入M2之后,点V显然不再是凸包顶点。该算法的Step8中也存在类似问题,可能导致在构建凸包顶点列的过程中将明显不在凸包上的点(如图3中的V)加入到H中,影响对后续点的判断,且最终可能导致在生成的凸包点列中存在不是凸包顶点的点。4反向活动分区简单多边形是平面点集的特例,故扩展其凸包算法可以得到平面点集的凸包算法。平面点集与多边形不同之处在于:一是平面点集中的点没有凹、凸顶点之分;二是平面点集不能依据极值点自然进行点集的子集划分,三是平面点集中的点并不具有多边形顶点之间的某些特殊性质。因而,平面点集的凸包算法不是多边形凸包算法的简单扩展。定义平面点集的候选点为在极值点Mi到Mi+1(i=1,3,5,7)连线左侧的点。称落入图1中子区域I、II、III和IV(不含边界)中的候选点集分别为Mi到Mi+1(i=7,5,3,1)之间的候选点子集。候选点子集划分过程如下:先找到点集的极值点,然后遍历点集中的非极值点V,若V满足S(Mi,Mi+1,V)>0(i=1或3或5或7),则V为Mi到Mi+1之间的候选点子集中的点。同样,以求解M1到M2之间的有序凸包点列为例进行说明。设M1到M2之间的候选点子集为Q。以下是对给定Q进行求解M1到M2之间的有序凸包点列H的算法。算法4.1SQRT_CHARRAY_PPS(M1,M2,Q)Step1.若M1=M2,则H={M1},结束;否则,若Q=∅,则H={M1,M2},结束;若Q≠∅,设Q={Q1,…,Qs},s≥1,令H={M1},r=1,i=1,V=Q1,继续;Step2.DEAL_CANDIDATE_PPS(V),i=i+1;Step3.若i≤s,则V=Qi,转Step2;否则,继续;Step4.r=r+1,Hr=M2。算法4.1中的DEAL_CANDIDATE_PPS(V)(见算法4.2)是对Q中的点进行处理并对H进行必要的更新的过程。为说明该过程,需要先引入反向活动分区的概念。活动区定义仍然如前,称前面定义的关于Hj的活动分区为Hj的正向活动分区,仍然记为AR(Hj),FIND_AR(V,l,u)不变,是查找V所在正向活动分区的函数。图4中的阴影部分即是与Hj(j=2,…,r-1)所对应的反向活动分区,其中的点在R(HjHj-1)的右侧,并在R(Hj+1Hj)的左侧或其上,称这样的活动分区为Hj的反向活动分区,记为IAR(Hj)。称R(Hj+1Hj)、R(HjHj-1)分别为IAR(Hj)的前、后边界。于是,给定点V∈IAR(Hi)(l≤i≤u)和点Hj(l≤j≤u),有如下结论:当V在R(HjHj-1)右侧且在R(Hj+1Hj)左侧或其上,则V∈IAR(Hj);当V(y)≥Hj(y)或V(y)<Hj(y)且V在R(Hj+1Hj)右侧时,V∈IAR(Hi)(j+1≤i≤u);其他,V∈IAR(Hi)(l≤i≤j-1)。在此,定义查找V所在反向活动分区的函数为FIND_IAR(V,l,u,t,H):即在Hl与Hu之间查找Hk,使得S(Hk,Hk-1,V)<0且S(Hk+1,Hk,V)≥0,返回k,IAR(Hk)就是V所在的反向活动分区,若S(Hk+1,Hk,V)=0,即V在IAR(Hk)的前边界上,则t=1,否则,t=0。算法4.2DEAL_CANDIDATE_PPS(V)Step1.若r=1,则r=r+1,Hr=V,退出;否则,INAVR_PPS(V,m,n,t);Step2.若m=0,则H不变,退出;若n=0,转Step3;若n>0,转Step5;Step3.若V(x)≥Hr(x),则r=m+1,Hr=V,退出;否则,继续;Step4.若S(V,M2,Hr)>0,则令W=Hr,r=m+1,Hr=V,r=r+1,Hr=W,退出;否则,r=m+1,Hr=V,退出。Step5.若t=0,则删除H中在Hm之后并在Hn之前的点,令Hm+1=V,将原H中自Hn以后的点接续到Hm+1之后,再令r=r-n+m+2,退出;若t=1,则删除H中在Hm之后并在Hn+1之前的点,令Hm+1=V,将原H中自Hn+1以后的点接续到Hm+1之后,再令r=r-n+m+1,退出。算法4.2的Step1中的INAVR_PPS(V,m,n,t)过程判断点V是否在活动区内。若是,则确定V所在的活动分区(见算法4.3)。若m=0,则V不在活动区内;若m>0,则AR(Hm)就是V所在正向活动分区。若n>0,则IAR(Hn)就是V所在反向活动分区;若t=1,则V在IAR(Hn)的前边界上;若t=0,则V不在IAR(Hn)的前边界上。对于平面点集中的点,则要分两种情况来分别确定点V所在活动分区:(a)V(x)≥Hr(x)和V(x)<Hr(x)且V在L(Hr-1Hr)左侧;(b)V(x)<Hr(x)且V在L(Hr-1Hr)之上或右侧。对情况(a),若确定V在活动区内,则只要确定V所在的正向活动分区就可以更新H,具体确定过程与INAVR_SP(V)中的Step1到Step5一样。对情况(b),在简单多边形凸包算法中,由于利用了简单多边形的性质,在算法3.3确定点V所在活动分区的INAVR_SP(V)过程中的Step5排除V在活动区内的可能性。而对于平面点集则不能简单加以排除,需要先查找Hj,使得Hj(x)<V(x)≤Hj+1(x),若V在S(Hj-1Hj)之上或右侧,则V不在活动区,否则,V在活动区内。若确定V在活动区内,则要如算法4.3确定点V所在的正、反向活动分区以便更新H。算法4.3INAVR_PPS(V,m,n,t)Step1.若V(x)=Hr(x),令n=0,转Step2;若Hr(x)<V(x)<M2(x),令n=0,转Step3;否则,转Step5;Step2.若V(y)>Hr(y),转Step4;否则,则令m=0,退出;Step3.若S(Hr-1,Hr,V)≥0,转Step4;若S(Hr-1,Hr,V)<0且S(Hr,M2,V)>0,则令m=r,n=0,退出;其他,令m=0,退出;Step4.令l=1,u=r-1,m=FIND_AR(V,l,u,H),退出;Step5.若S(Hr-1,Hr,V)>0,令n=0,转Step4;否则,查找Hj,使得Hj(x)<V(x)≤Hj+1(x),若S(Hj,Hj+1,V)≤0,则令m=0,n=0,退出,若S(Hj,Hj+1,V)>0,则FIND_AVR_PPS(V,j,m,n,t,H);算法4.3中的FIND_AVR_PPS(V,j,m,n,t,H)是依据V所在的位置来确定V所在的正、反向活动分区的过程(见算法4.4)。假设点V在活动区内,则V在如图5中所示的阴影部分所示的活动分区中,该活动分区被两条直线Hj-1Hj和Hj+1Hj+2划分成子区域I、II、III和IV(不含直线Hj-1Hj和Hj+1Hj+2)。设直线Hj-1Hj与Hj+1Hj+2的交点为G。当j=r-2时,只有子区域I、III;当j=1时,若V在R(H3H2)左侧,则V在子区域I内,若V在R(H3H2)右侧,则V在子区域II内。表1是依据V所在的位置来确定m、n、t的取值表,其中有类似l=j+2,u=r-1的形式:若这样的形式出现在m取值列中,则m=FIND_AR(V,l,u,H);若这样的形式出现n取值列中,则n=FIND_IAR(V,l,u,t,H),同时可以得出相应的t值。m、n、t的取值确定下来之后,就可以对H进行必要的更新。算法4.4FIND_AVR_PPS(V,j,m,n,t,H)Step1.若j=1,令m=1,转Step2;若j=r-2,令n=r-1,若S(Hr-1,Hr,V)=0,则令t=1,若S(Hr-1,Hr,V)<0,则令t=0,转Step3;否则,转Step4;Step2.若S(H3,H2,V)=0,则令n=2,t=1,退出;若S(H3,H2,V)>0,则令n=2,t=0,退出;否则,设l=3,u=r-1,n=FIND_IAR(V,l,u,t,H),退出;Step3.若S(Hr-3,Hr-2,V)<0,则令m=r-2,退出;若S(Hr-3,Hr-2,V)=0,则令m=r-3,退出;否则,设l=1,u=r-2,m=FIND_AR(V,l,u,H),退出;Step4.依据表1来确定m、n、t。求解平面点集凸包的完整算法如下:先找到极值点Mi(i=1,2,…,8)并划分得到Mi到Mi+1(i=1,3,5,7)之间的候选点子集,再求解得出Mi到Mi+1(i=1,3,5,7)之间的有序凸包点列,最后将Mi到Mi+1(i=1,3,5,7)之间的有序凸包点列顺序连接起来并去掉冗余的Mi点就得出多边形的凸包。5以时间复杂度为on这里只分析平面点集凸包