点集凸包算法python实现（二）

算法逻辑

在点集凸包算法python实现这篇博客中介绍了一种凸包算法，这种算法中凸包点搜索的过程较为麻烦，主要是因为计算点集连线与X轴的夹角需要考虑到四个不同象限，在这里通过计算向量夹角的方式，对凸包点集的搜索过程进行了优化

算法逻辑如下图所示：

python代码实现

from osgeo import ogr, gdalconst, osr
import numpy as np


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


def convexHull(xys: np.array):
    convexPoints = []
    ymin = np.min(xys, axis=0)[1]
    index = np.where(xys[:, 1] == ymin)[0][0]
    convexPoints.append(xys[index, :])
    n = 1
    while True:
        startPoint=convexPoints[-1]
        if n==1:
            previousPoint=np.array([startPoint[0]+100,startPoint[1]])
        else:
            previousPoint=convexPoints[-2]
        index=0
        minCosAngle=0
        v1=previousPoint-startPoint
        for i in range(xys.shape[0]):
            nextPoint=xys[i,:]
            v2=nextPoint-startPoint
            cosAngle=np.dot(v1,v2)/(np.sqrt(np.sum(v1**2))*np.sqrt(np.sum(v2**2)))
            if i==0:
                index=i
                minCosAngle=cosAngle
            else:
                if cosAngle<minCosAngle:
                    minCosAngle=cosAngle
                    index=i
        convexPoints.append(xys[index,:])
        n+=1
        xys=np.delete(xys,index,axis=0)
        firstPoint=convexPoints[0]
        endPoint=convexPoints[-1]
        if firstPoint[0]==endPoint[0] and firstPoint[1]==endPoint[1]:
            break
    wktPolygon = ""
    for i in range(n):
        x = convexPoints[i][0]
        y = convexPoints[i][1]
        wktPolygon = '{} {},'.format(x, y) + wktPolygon
    wktPolygon = wktPolygon[0:-1]
    wktPolygon = "POLYGON(({}))".format(wktPolygon)
    return wktPolygon


if __name__ == "__main__":
    ogr.RegisterAll()
    ds = ogr.Open("./point.shp", gdalconst.GA_ReadOnly)
    oLay = ogr.DataSource.GetLayer(ds, 0)
    ogr.Layer.ResetReading(oLay)
    xys = []
    while True:
        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!')

代码运行结果