凸包求解三种实现

枚举方法实现凸包

定理依赖：如果一个点在平面上某三个点组成的三角形内，那么这个点不可能是凸包上的点

• 需要的基本知识

叉积的计算：
a(x1,y1) b(x2,y2)=> axb=x1y2-x2y1
转换到代码中即为
ABxAC=(B[0]-A[0])*(C[1]-A[1])-(C[0]-A[0])*(B[1]-A[1])

利用叉乘判断M是否在三角形ABC中：
若M在ABC中，则MA X MB,MB X MC,MC X MA
此三组的向量叉乘的结果都是同号的（或都正，或都负），即方向相同的，则说明点M在三角形每条边的同侧，即内部。否则必在外部！

枚举代码实现：

import math
import numpy
import pylab as pl
from matplotlib import pyplot as plt
import time

'''
依赖定理：如果一个点在平面上某三个点组成的三角形内，那么这个点不可能是凸包上的点
'''
"""
实现思路：
利用numpy来实现随机的点
构造一个n*3大小的矩阵，每一列代表的意义即为x y 极坐标；
我们利用叉积来判断每四个点是否由点在其余三个点形成的三角形内，若是，则由定理描述，不是凸包上的点；
把满足要求的点按极坐标大小排列，然后利用画图展示效果
"""
def drawGraph(x,y):
    pl.figure()
    pl.title("ConvexHull-BruteForce")
    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)   #画最后一条封闭线

#蛮力法判断叉积，返回ABxAC的向量中j的系数
"""
a(x1,y1) b(x2,y2)=> axb=x1y2-x2y1
ABxAC=(B[0]-A[0])*(C[1]-A[1])-(C[0]-A[0])*(B[1]-A[1])
"""
def Product(A,B,C):
    return (B[0]-A[0])*(C[1]-A[1])-(C[0]-A[0])*(B[1]-A[1])
#判断M是否在三角形ABC中
"""
利用叉乘判断M是否在三角形ABC中，即
若M在ABC中，则MA X MB,MB X MC,MC X MA
此三组的向量叉乘的结果都是同号的（或都正，或都负），即方向相同的，则说明点M在三角形每条边的同侧，即内部。否则必在外部！
"""
def isInTriangle(A,B,C,M):
    if Product(A,B,M)>=0 and Product(B,C,M)>=0 and Product(C,A,M)>=0:
        return 1
    return 0

#存点的矩阵,每行一个点,列0->x坐标,列1->y坐标,列2->代表极角or访问标记
def matrix(rows,cols):
    cols=3
    mat = [[0 for col in range(cols)]
              for row in range(rows)]
    #0是对这个矩阵赋予的初值
    return mat

#凸包蛮力算法
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
                    # 虽然mat[a/b/c/d]有三个值，但传入到函数中使用，只用了两个，所以不影响，且函数参数的构建也并没有规定参数的类型，不影响；
                    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 #tmp记录y最低的点
            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))
        #print(x,y)
    #print(stack) 列表中的每一个元素是一对x,y坐标；
    return stack


#main
#numpy.random.random(size=None) 生成[0,1)之间的浮点数
max_num = 70#最大点数
x=100*numpy.random.random(max_num)# 生成一百个[0,100)之间的数
y=100*numpy.random.random(max_num)# 生成一百个[0,100)之间的数
#（x,y)组合成了我们需要的100个随机点的坐标
#drawGraph(x,y)  #原始图

CHQ=BruteForce(x,y)
#pl.subplot(132)
drawGraph(x,y)
drawCH(CHQ)     #BruteForce凸包
pl.show()

GrahamScan方法实现凸包

GrahamScan算法:
首先要对点集S进行预处理，选纵坐标最小的点作为P_0，其余点在以P_0为极点、水平方向为极轴的极坐标系下按极角从小到大排序

在极坐标系下按照极角大小排列，然后逆时针方向漫游点集，去除非Convex hull顶点（非左转点）

伪代码：
Graham-Scan(S)：
输入：平面n个点的集合S，经预处理后存储于P[0:m]
输出：S的凸包Q
If m<=1
Then 返回“凸包为空”
将P[0]、P[1]、P[2]依次压入栈Q
For i=3 To n-1 Do
   While Top(Q) ∈∆P[0]P[i]NextToTop(Q) Do
      Pop(Q)
   Push(Q,P[i])

枚举代码实现：

'''
GrahamScan算法:
首先要对点集S进行预处理，选纵坐标最小的点作为P_0，其余点在以P_0为极点、水平方向为极轴的极坐标系下按极角从小到大排序


伪代码：
Graham-Scan(S)：
输入：平面n个点的集合S，经预处理后存储于P[0:m]
输出：S的凸包Q
If m<=1
Then 返回“凸包为空”
将P[0]、P[1]、P[2]依次压入栈Q
For i=3 To n-1 Do
   While Top(Q) ∈∆P[0]P[i]NextToTop(Q) Do
      Pop(Q)
   Push(Q,P[i])
'''
import math
import numpy
import pylab as pl
def drawGraph(x,y):
    pl.figure()
    pl.title("ConvexHull-GrahamScan")
    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])
    #a(x1,y1) b(x2,y2)=> axb=x1y2-x2y1
    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]) #按极角由小到大排序
    #G_C需要按照极角顺序的由小到大来进行遍历求解
    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)
        #在我们人眼识别中是判断是否为左转点 不是就排除，而计算机判断则是利用向量p1p3，p1p2叉积判断 叉积>0=>向量夹角小于90，即是图中反映的左转点；<0 =>向量夹角大于90
    return stack

#main
#numpy.random.random(size=None) 生成[0,1)之间的浮点数
max_num = 30#最大点数
x=100*numpy.random.random(max_num)# 生成一百个[0,100)之间的数
y=100*numpy.random.random(max_num)# 生成一百个[0,100)之间的数
#（x,y)组合成了我们需要的100个随机点的坐标
#drawGraph(x,y)  #原始图

CHQ=GrahamScan(x,y)
#pl.subplot(132)
drawGraph(x,y)
drawCH(CHQ)    #Graham凸包
pl.show()

分治方法实现凸包

Preprocess：找到横坐标最大的点A和最小的点B，标记每个点是否访问，若全部访问则算法停止。(时间复杂度为O(n))；
Divide：作直线AB，把凸包分为SL和SU上下两个子集，对每个部分求得点Pmax,使得SABPmax最大。将三角形三个点和三角形内部的点标记为已访问，删去所有三角形内部及边上的点。(时间复杂度为O(n))；
Conquer: 进一步依据△ABPmax划分成左右两个部分，当作SL和SU，分治 递归、不断重复。(时间复杂度为2T(n/2))。
其中三角形面积计算如下图所示，为S=(x1y2-x1y3+x2y3-x2y1+x3y1-x2y2)。
由master定理解得算法总时间复杂度为O(nlogn)，算法的原理示意图如图3所示。

分治思想图

分治代码实现：

import math
import numpy
import random
import pylab as pl
import matplotlib.pyplot as plt
import convexhull_main

#s生成子图移到main函数
# fig = plt.figure()
# ax = fig.add_subplot(1,1,1)

#坐标生成
all = []
def random_points_generation():
    n=30
    while(len(all)<n):
        x=100*numpy.random.random()
        y=100*numpy.random.random()
        if([x,y] not in all):
                all.append([x,y])
    return all

#计算三角形面积
def triangle_calculate(A,B,C):
        x1=A[0]
        y1=A[1]
        x2=B[0]
        y2=B[1]
        x3=C[0]
        y3=C[1]
        return x1 * y2 + x3 * y1 + x2 * y3 - x3 * y2 - x2 * y1 - x1 * y3

# 将所有的点绘制在坐标图上
def draw(all_points):
    font2 = {'family': 'Times New Roman',
                 'weight': 'normal',
                 'size': 15,
            }
    list_all_x = []
    list_all_y = []
    for item in all_points:
            a, b = item
            list_all_x.append(a)
            list_all_y.append(b)
    ax.set_title("ConvexHull_DivideConquer_Construction",font2)
    ax.set_xlabel("x axis", font2)
    ax.set_ylabel("y axis", font2)
    ax.scatter(list_all_x, list_all_y, c='k', marker='*')

    # pl.figure()
    # pl.title("ConvexHull_DivideConquer_Construction",font2)
    # pl.xlabel("x axis", font2)
    # pl.ylabel("y axis", font2)
    # pl.plot(list_all_x, list_all_y, 'k*')

class Divideconvenhull:
    def __init__(self,left_point, right_point, lists, border_points):
        self.area_max=0
        self.point_max = ()
        self.left_point = left_point
        self.right_point = right_point
        self.lists = lists
        self.border_points = border_points

    # 分为上部分和下部分处理
    def up(self):
        for item in self.lists:
            if item == self.left_point or item == self.right_point:
                continue
            else:
                # abc计算三角形面积 若c在ab的左侧为正
                new_area = triangle_calculate(self.left_point, self.right_point, item)
                if new_area > self.area_max:
                    self.point_max = item
                    self.area_max = new_area

        if self.area_max != 0:
            # 分治处理 上三角形的左边 和上三角形的右边
            self.border_points.append(self.point_max)
            updi=Divideconvenhull(self.left_point, self.point_max, self.lists, self.border_points)
            updi.up()
            updi2=Divideconvenhull(self.point_max,self.right_point,  self.lists, self.border_points)
            updi2.up()

    def down(self):
        for item in self.lists:
            if item == self.left_point or item == self.right_point:
                continue
            else:
                new_area = triangle_calculate(self.left_point, self.right_point, item)
                # 由于ABC，C在AB的右侧，则三角形面积应为负，越小越大；
                if new_area < self.area_max:
                    self.point_max = item
                    self.area_max = new_area

        if self.area_max != 0:
            self.border_points.append(self.point_max)
            downdi = Divideconvenhull(self.left_point, self.point_max, self.lists, self.border_points)
            downdi.down()
            downdi2=Divideconvenhull(self.point_max,self.right_point,  self.lists, self.border_points)
            downdi2.down()

    def convexhulldraw(self,):
        self.border_points.sort()
        first_x, first_y = self.border_points[0]  # 最左边的点
        last_x, last_y = self.border_points[-1]  # 最右边的点
        list_border_up = []  # 上半边界
        for item in self.border_points:
            x, y = item
            if y > max(first_y, last_y):
                list_border_up.append(item)
            if min(first_y, last_y) < y < max(first_y, last_y):
                if triangle_calculate(self.border_points[0], self.border_points[-1], item) > 0:
                    list_border_up.append(item)
                else:
                    continue
        list_border_down = [_ for _ in self.border_points if _ not in list_border_up]  # 下半边界
        list_end = list_border_up + list_border_down[::-1]  # 最终顺时针输出的边界点
        list_end.append(list_end[0])
        for i in range(len(list_end) - 1):
            one_, oneI = list_end[i]
            two_, twoI = list_end[i + 1]
            ax.plot([one_, two_], [oneI, twoI],'c--',linewidth=3)

if __name__ == "__main__":
        fig = plt.figure()
        ax = fig.add_subplot(1, 1, 1)
        all_points =convexhull_main.all_points

        #按横坐标升序排序
        all_points.sort()
        #横坐标最左边和最右边
        left=all_points[0]
        right=all_points[-1]
        border_points=[]#边界点集
        draw(all_points)
        eg=Divideconvenhull(left,right,all_points,border_points)
        eg.up()
        eg.down()
        border_points.append(left)
        border_points.append(right)
        print(border_points)  # 输出边界点
        eg.convexhulldraw()
        plt.show()

运行效果图

下图是为100测试点的效果图