点集凸包算法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
    点赞
  • 46
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
平面点集凸包(Convex Hull)是指在一个平面上,所有点集中点的最小子集,构成一个不包含内部点的凸多边形。这个概念在计算机图形学几何算法和数据分析中非常常见。要编写一个计算凸包的代码,通常会使用诸如 Graham 扫描法、快速 hull 算法(如 Jarvis March 或 Andrew 算法)或 Gift Wrapping 算法等高效的算法。 这里提供一个基于 Gift Wrapping 算法的简单 Python 示例,它假设你有一个二维列表(二维数组)存储了点集: ```python def convex_hull(points): if len(points) < 3: # 凸包至少需要3个点 return points # 将点按 X 坐标排序 points.sort(key=lambda x: x) hull = [points, points] # 初始化凸包,包含第一个点和第二个点 for i in range(2, len(points)): while len(hull) >= 2 and cross_product(hull[-2], hull[-1], points[i]) <= 0: hull.pop() # 如果当前点在凸包内部,则移除最后一个点 hull.append(points[i]) # 因为最后两个点可能共线,所以需要再次检查并添加最后一个点 if hull and cross_product(hull[-2], hull[-1], points) <= 0: hull.pop() return hull def cross_product(p1, p2, p3): return (p2 - p1) * (p3[1] - p1) - (p2 - p1) * (p3 - p1) # 使用示例 points = [[1, 1], [2, 2], [3, 3], [4, 4], [1, 4]] convex_hull_points = convex_hull(points) print("凸包为:", convex_hull_points) ``` 这个代码定义了一个 `convex_hull` 函数,其中包含 `cross_product` 函数用于计算向量的叉积,这在判断三点是否构成凸包方向上很重要。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值