RT。求平面上点集的凸包。
1. GrahamScan算法，《算法导论》上的例子，先找到y最小的点O，以O建立极坐标，其它点按极角排序后再从头开始扫描（配合stack实现）。
2.BruteForce算法，依赖定理：如果一个点在平面上某三个点组成的三角形内，那么这个点不可能是凸包上的点。
所以暴力的思路是平面上的点每4个进行枚举，并判断是否满足定理，若满足，则删除这个点继续找；一直找到没有满足定理的点为止。。
3.DivideAndConquer思路：有很多种，这里我实现的只是严格意义上最后一步的分治，基于递归的分治暂时还没想明白；
思路即：首先找点集x-坐标的中位数(排序法O(nlogn)，或随机法更快O(n))，将每次将点集分为左右两部分。
分别递归求左右部分的凸包，这里我就是用Graham求的凸包。这样会过滤掉许多内部点，将左边凸包上的点序列称为left，右边凸包上的点序列称为right。
right按顺时针和逆时针分成两部分right1和right2（以y最大和y最小的点为界），得到这三个极角有序列表后，合并这三个有序表为一个list（merge的操作O(n)，原理类似归排的merge）
对新list直接进行GrahamScan，只需过滤掉极少的点即可得到大凸包。
ConvexHull.py:
#coding=utf-8
import math
import numpy
import pylab as pl
#画原始图
def drawGraph(x,y):
pl.figure(1)
pl.subplot(131) #1
pl.title("ConvexHull-GrahamScan")
pl.xlabel("x axis")
pl.ylabel("y axis")
pl.plot(x,y,'ro')
pl.subplot(132) #2
pl.title("ConvexHull-BruteForce")
pl.xlabel("x axis")
pl.ylabel("y axis")
pl.plot(x,y,'ro')
pl.subplot(133) #3
pl.title("ConvexHull-DivideConquer")
pl.xlabel("x axis")
pl.ylabel("y axis")
pl.plot(x,y,'ro')
#画凸包
def drawCH(CHQ):
x,y=[],[]
for p in CHQ:
x.append(p[0])
y.append(p[1])
pl.plot(x,y,color='blue',linewidth=2)
pl.plot(x[-1],y[-1],x[0],y[0])
lastx=[x[-1],x[0]]
lasty=[y[-1],y[0]]
pl.plot(lastx,lasty,color='blue',linewidth=2) #画最后一条封闭线
#存点的矩阵,每行一个点,列0->x坐标,列1->y坐标,列2->代表极角
def matrix(rows,cols):
cols=3
mat = [[0 for col in range (cols)]
for row in range(rows)]
return mat
#Graham用的叉积
def crossMut(stack,p3):
p2=stack[-1]
p1=stack[-2]
vx,vy=(p2[0]-p1[0],p2[1]-p1[1])
wx,wy=(p3[0]-p1[0],p3[1]-p1[1])
return (vx*wy-vy*wx)
#Graham扫描法O(nlogn),mat是经过极角排序的点集
def GrahamScan(x,y):
#print mat
max_num=len(x)
mat = matrix(max_num,3) #根据点数申请矩阵,mat[i][2]代表极角
minn = 300
for i in range(max_num): #存点集
mat[i][0],mat[i][1]=x[i],y[i]
if y[i]<minn: #(x[tmp],y[tmp])-->y轴最低坐标
minn=y[i]
tmp = i
d = {} #利用字典的排序
for i in range(max_num): #计算极角
if (mat[i][0],mat[i][1])==(x[tmp],y[tmp]):mat[i][2]=0
else:mat[i][2]=math.atan2((mat[i][1]-y[tmp]),(mat[i][0]-x[tmp]))
d[(mat[i][0],mat[i][1])]=mat[i][2]
lst=sorted(d.items(),key=lambda e:e[1]) #按极角由小到大排序
for i in range(max_num): #更新mat为排序后的矩阵
((x,y),eth0)=lst[i]
mat[i][0],mat[i][1],mat[i][2]=x,y,eth0
points=len(mat) #点数
stack=[]
stack.append((mat[0][0],mat[0][1])) #push p0
stack.append((mat[1][0],mat[1][1])) #push p1
stack.append((mat[2][0],mat[2][1])) #push p2
for i in range(3,points):
#print stack
p3=(mat[i][0],mat[i][1])
while crossMut(stack,p3)<0:stack.pop()
stack.append(p3)
return stack
#蛮力法判断叉积，返回ABxAC的向量中j的系数
def Product(A,B,C):
return (B[0]-A[0])*(C[1]-A[1])-(C[0]-A[0])*(B[1]-A[1])
#判断P是否在三角形ABC中
def isInTriangle(A,B,C,P):
if Product(A,B,P)>=0 and Product(B,C,P)>=0 and Product(C,A,P)>=0:
return 1
return 0
#凸包蛮力算法
def BruteForce(x,y):
max_num=len(x)
mat = matrix(max_num,3) #根据点数申请矩阵,mat[i][2]代表访问标记
for i in range(max_num): #存点集
mat[i][0],mat[i][1],mat[i][2]=x[i],y[i],1
#任选4个,即C(10,4)的功能
for a in range(0,max_num-3):
for b in range(a+1,max_num-2):
for c in range(b+1,max_num-1):
for d in range(c+1,max_num):
#如果在三角形中,则删除内部的点
#if 0 in (mat[a][2],mat[b][2],mat[c][2],mat[d][2]):continue
if isInTriangle(mat[b],mat[c],mat[d],mat[a]):mat[a][2]=0 #顺时针算一类
if isInTriangle(mat[b],mat[d],mat[c],mat[a]):mat[a][2]=0 #逆时针算另一类
if isInTriangle(mat[a],mat[c],mat[d],mat[b]):mat[b][2]=0
if isInTriangle(mat[a],mat[d],mat[c],mat[b]):mat[b][2]=0
if isInTriangle(mat[a],mat[b],mat[d],mat[c]):mat[c][2]=0
if isInTriangle(mat[a],mat[d],mat[b],mat[c]):mat[c][2]=0
if isInTriangle(mat[a],mat[b],mat[c],mat[d]):mat[d][2]=0
if isInTriangle(mat[a],mat[c],mat[b],mat[d]):mat[d][2]=0
#后处理,按极角排序,以便输出图形
#print mat
newmat = matrix(max_num,3) #newmat[i][2]是极角,和mat[i][2]不同
pos = 0 #记录newmat行数
minn = 300
for i in range(len(mat)):
if mat[i][2]==1:
if mat[i][1]<minn: #(mat[tmp][0],mat[tmp][1])-->y坐标最低的点
minn=mat[i][1]
tmp = i
newmat[pos][0],newmat[pos][1]=mat[i][0],mat[i][1]
pos+=1
d={} #排序字典
for i in range(pos): #计算极角
#print newmat[i][0],newmat[i][1],newmat[i][2]
if (newmat[i][0],newmat[i][1])==(mat[tmp][0],mat[tmp][1]):newmat[i][2]=0
else:newmat[i][2]=math.atan2((newmat[i][1]-mat[tmp][1]),(newmat[i][0]-mat[tmp][0]))
d[(newmat[i][0],newmat[i][1])]=newmat[i][2]
lst=sorted(d.items(),key=lambda e:e[1]) #按极角由小到大排序
#print lst
stack=[]
for i in range(pos): #更新mat为排序后的矩阵
((x,y),eth0)=lst[i]
stack.append((x,y))
return stack
#凸包分治算法
def DivideConquer(x,y):
max_num=len(x)
d={}
for i in range(max_num):
d[(x[i],y[i])]=x[i]
#print d
lst=sorted(d.items(),key=lambda e:e[1]) #首先将点集按x坐标排序
#for k in lst:print k
#print
left=len(lst)/2 #左边的个数
right=len(lst)-left #右边的个数
#print left,right
x,y=[],[] #画左半部分
for i in range(left):
((xi,yi),noUse)=lst[i]
x.append(xi)
y.append(yi)
CHQ_L=GrahamScan(x,y)
#print CHQ_L
x,y=[],[] #画右半部分
for i in range(right):
((xi,yi),noUse)=lst[left+i]
x.append(xi)
y.append(yi)
CHQ_R=GrahamScan(x,y)
for i in range(len(CHQ_R)): #找到y最高的位置high
if i==len(CHQ_R)-1 or CHQ_R[i][1]>CHQ_R[i+1][1]:
high = i
break
#将右半部分分成两个序列
sq2=[]
for i in range(high+1):
sq2.append((CHQ_R[i][0],CHQ_R[i][1]))
#print sq2
sq3=[]
for i in range(len(CHQ_R)-1,high,-1):
sq3.append((CHQ_R[i][0],CHQ_R[i][1]))
#print sq3
merge = CHQ_L+sq2+sq3 #合并操作，应该按顺序来merge
#print merge
x,y=[],[] #画左半部分
for i in range(len(merge)):
((xi,yi))=merge[i]
x.append(xi)
y.append(yi)
CHQ=GrahamScan(x,y)
return CHQ
#main
max_num = 30 #最大点数
x=100*numpy.random.random(max_num)#[0,100)
y=100*numpy.random.random(max_num)
drawGraph(x,y) #原始图
CHQ=GrahamScan(x,y)
pl.subplot(131)
drawCH(CHQ) #Graham凸包
CHQ=BruteForce(x,y)
pl.subplot(132)
drawCH(CHQ) #BruteForce凸包
CHQ=DivideConquer(x,y)
pl.subplot(133)
drawCH(CHQ) #DivideConquer凸包
pl.show()
ps:凸包的一个应用是求平面最远点对，源自定理：平面最远的点对一定均在凸包上。
时间： 2024-10-27 19:37:28