点集凸包算法python实现

什么是凸包?

凸包定义

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

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

image-20220422102723122

通俗理解

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

image-20220422103122651

凸包有什么特点?

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

凸包有什么用途?

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

绘制凸包的常见算法?

  • Jarvis march(包裹法)

  • Graham Scan(扫描法)

  • Divide and conquer(分治法)

  • Divide and conquer(分治法)

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

Jarvis march(包裹法)

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

基本思想如下:

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

1

点集中任意一点与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

2

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

算法流程图

未命名文件 (2)

补充

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

  • 第一象限

象限角1

alfa=arctan(dy/dx)

  • 第二象限

象限角2

alfa=arctan(dy/dx)+pi

  • 第三象限

象限角3

alfa=arctan(dy/dx)+pi

  • 第四象限

象限角4

alfa=arctan(dy/dx)+2pi

python代码

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

如何计算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 osgeo import ogr, gdalconst, osr
import numpy as np
import pandas as pd
import os

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:
        if n == 1:
            startPoint = convexPoints[0]
            minAlfa = 0
            index = 0
            for i in range(xys.shape[0]):
                endPoint = xys[i, :]
                alfa = getAlfa(startPoint, endPoint)
                if i == 0:
                    minAlfa = alfa
                    index = 0
                else:
                    if alfa < minAlfa:
                        minAlfa = alfa
                        index = i
            convexPoints.append(xys[index, :])
            xys = np.delete(xys, index, axis=0)
            n += 1
        else:
            startPoint = convexPoints[-1]
            previousPoint = convexPoints[-2]
            alfa0 = getAlfa(startPoint, previousPoint)
            maxAlfa = 0
            index = 0
            for i in range(xys.shape[0]):
                nextPoint = xys[i, :]
                alfa1 = getAlfa(startPoint, nextPoint)
                dalfa = (alfa0 - alfa1) if alfa0 > alfa1 else (alfa1 - alfa0)
                if dalfa > np.pi:
                    dalfa = np.pi * 2 - dalfa
                if i == 0:
                    maxAlfa = dalfa
                    index = 0
                else:
                    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 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!')

结算结果,如下图所示:

image-20220422155423391

  • 4
    点赞
  • 45
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值