什么是凸包？

凸包定义

点集p的凸包是指一个最小凸多边形（内角均小于180°），满足p中的点或者在多边形边上或者在其内

下图中的红色线段表示的多边形就是点集p={p0,p1,p2,p3,…………,p12}的凸包

通俗理解

• 一组平面上的点，求一个包含所有点的最小的凸多边形
• 这可以形象地想成这样：在地上放置一些不可移动的木桩(代表点集中的点)，用一根绳子把他们尽量紧地圈起来，这就是凸包

凸包有什么特点？

• 整个凸包都在任意一条边的一侧
• 凸包任意两点的中点都在凸包内
• 凸包内的任意点集的加权平均（凸组合）都在凸包内

凸包有什么用途？

• 从点集中抽象出一个唯一确定的凸多边形，即一组点集的凸包是唯一的，凹包并不唯一
• 用尽量少的点来描述一个点集的边界
• 使点集有序
• 对复杂多边形进行化简
• 为其它算法做预处理

绘制凸包的常见算法？

• Jarvis march（包裹法）

• Graham Scan（扫描法）

• Divide and conquer（分治法）

• Divide and conquer（分治法）

这里只讲解最简单常用的Jarvis march（包裹法）

Jarvis march（包裹法）

这种算法利用了凸包的一个特性：整个凸包都在任意一条边的一侧

基本思想如下：

• 选取必定在凸包上的一个点作为起点p0，可以是最低点，最高点，最左侧点，最右侧点其中一个，以最低点为例，将其放入凸包点集convexPoints=[p0]
• 从点集剩余点中选取一点p1，使得点集中的所有点都在p0p1连线的同侧，方法如下：

点集中任意一点与p0的连线与x轴的夹角为alfa，选择使得夹角alfa最小的点为p1，此时点集中的所有点都在p0p1连线的同侧,将p1放入凸包点集中convexPoints=[p0,p1],并从点集p中删除点p1

• 以点集p1为起点startPoint，p0为上一点previousPoint，在点集p的剩余点中寻找下一个凸包点p2=nextPoint，使得从startPoint出发到previousPoint 和nextPoint的夹角最大，此时该点为下一凸包点convexPoints=[p0,p1,p2]，从点集p中删除点nextPoint

• 以p2为新起点startPoint，p1为新的上一点previousPoints,重复上一步骤，在p剩余点集中寻找下一点p3=nextPoint，直到寻找的nextPoint在convexPoints中已经出现，此时凸包点已经首位闭合

算法流程图

补充

如何计算与x轴的逆时针夹角alfa

• 第一象限

alfa=arctan(dy/dx)

• 第二象限

alfa=arctan(dy/dx)+pi

• 第三象限

alfa=arctan(dy/dx)+pi

• 第四象限

alfa=arctan(dy/dx)+2pi

python代码

defgetAlfa(startPoint, endPoint):
    dx= endPoint[0]- startPoint[0]
    dy= endPoint[1]- startPoint[1]if dx==0:if dy>0:return np.pi/2elif dy<0:return3* np.pi/2else:return np.pi
    alfa= np.arctan(dy/ dx)if dx<0:
        alfa+= np.pielif dy<0:
        alfa+= np.pi*2return alfa

如何计算startPoint与previousPoint，nextPoint连线的夹角dalfa

• 记startPoint与previousPoint连线与x轴的夹角为alfa1
• 记startPoint与nextPoint连线与x轴的夹角为alfa2

dlfa=abs（alfa2-alfa1）

如果dalfa大于pi，由于凸多边形的内角小于180°，因此此时还需要用dalfa=2pi-dalfa

完整代码

基于上述算法，用gda库计算点SHP文件的凸包

from osgeoimport ogr, gdalconst, osrimport numpyas npimport pandasas pdimport osdefgetAlfa(startPoint, endPoint):
    dx= endPoint[0]- startPoint[0]
    dy= endPoint[1]- startPoint[1]if dx==0:if dy>0:return np.pi/2elif dy<0:return3* np.pi/2else:return np.pi
    alfa= np.arctan(dy/ dx)if dx<0:
        alfa+= np.pielif dy<0:
        alfa+= np.pi*2return alfadefconvexHull(xys: np.array):
    convexPoints=[]
    ymin= np.min(xys, axis=0)[1]
    index= np.where(xys[:,1]== ymin)[0][0]
    convexPoints.append(xys[index,:])
    n=1whileTrue:if n==1:
            startPoint= convexPoints[0]
            minAlfa=0
            index=0for iinrange(xys.shape[0]):
                endPoint= xys[i,:]
                alfa= getAlfa(startPoint, endPoint)if i==0:
                    minAlfa= alfa
                    index=0else:if alfa< minAlfa:
                        minAlfa= alfa
                        index= i
            convexPoints.append(xys[index,:])
            xys= np.delete(xys, index, axis=0)
            n+=1else:
            startPoint= convexPoints[-1]
            previousPoint= convexPoints[-2]
            alfa0= getAlfa(startPoint, previousPoint)
            maxAlfa=0
            index=0for iinrange(xys.shape[0]):
                nextPoint= xys[i,:]
                alfa1= getAlfa(startPoint, nextPoint)
                dalfa=(alfa0- alfa1)if alfa0> alfa1else(alfa1- alfa0)if dalfa> np.pi:
                    dalfa= np.pi*2- dalfaif i==0:
                    maxAlfa= dalfa
                    index=0else:if dalfa> maxAlfa:
                        maxAlfa= dalfa
                        index= i
            convexPoints.append(xys[index,:])
            xys= np.delete(xys, index, axis=0)
            n+=1# 判断首位是否重合
        firstPoint= convexPoints[0]
        lastPoint= convexPoints[-1]if firstPoint[0]== lastPoint[0]and firstPoint[1]== lastPoint[1]:break
    wktPolygon=""for iinrange(n):
        x= convexPoints[i][0]
        y= convexPoints[i][1]
        wktPolygon='{} {},'.format(x, y)+ wktPolygon
    wktPolygon= wktPolygon[0:-1]
    wktPolygon="POLYGON(({}))".format(wktPolygon)return wktPolygonif __name__=="__main__":
    ogr.RegisterAll()
    ds= ogr.Open("./point.shp", gdalconst.GA_ReadOnly)
    oLay= ogr.DataSource.GetLayer(ds,0)
    ogr.Layer.ResetReading(oLay)
    xys=[]whileTrue:
        oFea= ogr.Layer.GetNextFeature(oLay)if oFea==None:break
        oPoi= ogr.Feature.GetGeometryRef(oFea)
        x= ogr.Geometry.GetX(oPoi)
        y= ogr.Geometry.GetY(oPoi)
        xys.append([x, y])
    xys= np.array(xys)
    wktPolygon= convexHull(xys)
    driver= ogr.GetDriverByName("ESRI Shapefile")
    convexds= ogr.Driver.CreateDataSource(driver,"convexHull.shp")
    srs= ogr.Layer.GetSpatialRef(oLay)
    convexLay= ogr.DataSource.CreateLayer(convexds,"convexhull", srs, ogr.wkbPolygon)
    labelField= ogr.FieldDefn("label", ogr.OFTInteger)
    ogr.Layer.CreateField(convexLay, labelField)
    convexFea= ogr.Feature(ogr.Layer.GetLayerDefn(convexLay))
    ogr.Feature.SetField(convexFea,"label",1)
    convexPolygon=ogr.CreateGeometryFromWkt(wktPolygon)
    ogr.Feature.SetGeometry(convexFea,convexPolygon)
    ogr.Layer.CreateFeature(convexLay,convexFea)
    ogr.DataSource.Destroy(convexds)
    ogr.DataSource.Destroy(ds)print('ok!')

结算结果，如下图所示：