Convex hull Algorithm-grahamscan+ violence + divided treatment

Rt. The convex hull of a point set on a plane.

1. Grahamscan algorithm, "Introduction to the algorithm" in the example, first find the y smallest point O, with O to establish the polar coordinates, other points by the angle of the order and then start from the beginning to scan (with stack implementation).

2.BruteForce algorithm, dependency theorem : If a point is within a triangle of three points on a plane, then this point cannot be a point on a convex hull.

So the idea of violence is that the dots on the plane are enumerated every 4, and whether the theorem is satisfied, if satisfied, then the point of deletion continues to be found; until the point where the theorem is not satisfied is found.

3.DivideAndConquer ideas: There are many kinds, here I achieve only the strict meaning of the last step of the division, based on recursive Division of the provisional not to understand;

The idea is: first find the point set X-coordinate of the median (sort O (nlogn), or Random method faster O (n)), each time the point set is divided into two parts.

The convex hull of the left and right part is recursively divided, and here I am the convex bag of Graham. This filters out many of the internal points, and the sequence of points on the left convex hull is referred to as leave, and the sequence of points on the right convex hull is referred to as.

Right is divided clockwise and counterclockwise into two parts right1 and right2 (bounded by y Max and y minimum points), and after obtaining the three ordered list of Poles, merge the three ordered tables into a list (the Operation O (n) of the merge, similar to the normalized merge)

To grahamscan the new list directly, you can get the large convex hull by filtering out the very few points.

convexhull.py:

#coding =utf-8import mathimport numpyimport Pylab as pl# draw original 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 (#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 DRAWC     H (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 ', Linewid th=2) #画最后一条封闭线 # Save point Matrix, one point per line, column 0->x coordinates, column 1->y coordinates, column 2-> for polar def Matrix (rows,cols): cols=3 mat = [[0 for Col in RA] Nge (cols)] for row in range (rows)] return Mat#graham with cross product def Crossmut (STACK,P3): P2=stack[-1] P1=sta Ck[-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 is a polar-sorted point set Def  Grahamscan (x, y): #print Mat Max_num=len (×) mat = Matrix (max_num,3) #根据点数申请矩阵, mat[i][2] for polar angle Minn = + 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 axis Lowest coordinate 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.appe nd ((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.appe nd (p3) return stack# brute force method to determine the cross product, return the coefficients of j in the vector of Abxac 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 (×) mat = Matrix (max_num,3) #根据点数申请矩阵, mat[i][2  ] represents the access token for I in range (max_num): #存点集 mat[i][0],mat[i][1],mat[i][2]=x[i],y[i],1 #任选4个, which is the function of 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): #如果在三角形中, delete the internal point #if 0 in (mat[a][2],mat[b][2],mat[c][2],m                    AT[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 IsIn Triangle (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 #后处理, sorted by polar angle for output graphics #print mat Newmat = Matrix (max_num,3) #newmat [i][2] is the polar angle, and mat[ I][2] Different pos = 0 #记录newmat行数 minn = + for I in range (Len (MAT)): If Mat[i][2]==1:if mat[i][1]& Lt;minn: # (mat[tmp][0],mat[tmp][1])-->y coordinates the lowest point minn=mat[i][1] tmp = i Newmat[po  S][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]=m Ath.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# convex packet divide algorithm def divideconquer (x, y): Max_num=len (×) 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=grah Amscan (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]>c HQ_R[I+1][1]: High = I broke #将右半部分分成两个序列 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 #合并操作 should be in order to 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 #mainmax_num = #最大点数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 (drawch) # Bruteforce Convex package Chq=divideconquer (x, y) Pl.subplot (133) Drawch (CHQ) #DivideConquer凸包pl. Show ()

PS: One application of convex hull is to find out the furthest point of plane, from the theorem: the point pair of the farthest plane must be on the convex package.

Convex hull Algorithm-grahamscan+ violence + divided treatment