python实现MPR（最小扰动半径）计算

# -*- coding: utf-8 -*-
# @Time    : 2024/6/13 下午4:08
# @Author  :Fivem
# @File    : MPR.py
# @Software: PyCharm
# @last modified:2024/6/13 下午4:08
import glob
import os

import geopandas as gpd
import numpy as np
from scipy.spatial import Delaunay

from watermark.vector_process.get_coor import get_coor_nested, get_coor_array
from watermark.vector_process.select_file import select_folder


# 计算矢量地图坐标点的最小扰动区域（MPR）的半径
# s1:获取矢量地图的所有坐标点的数组
# s2:对所有坐标点生成Delaunay三角网
# s3:对所有三角网生成内切圆，并计算最小的半径

def compute_inradius(vertices):
    # 计算三角形的边长
    a = np.linalg.norm(vertices[1] - vertices[0])
    b = np.linalg.norm(vertices[2] - vertices[1])
    c = np.linalg.norm(vertices[0] - vertices[2])

    # 计算半周长
    s = (a + b + c) / 2

    # 计算面积（使用海伦公式）
    area = np.sqrt(s * (s - a) * (s - b) * (s - c))

    # 计算内切圆半径
    inradius = area / s

    return inradius


def compute_incenter(vertices):
    # 计算三角形的边长
    a = np.linalg.norm(vertices[1] - vertices[0])
    b = np.linalg.norm(vertices[2] - vertices[1])
    c = np.linalg.norm(vertices[0] - vertices[2])

    # 计算内切圆中心
    P = a + b + c
    incenter = (a * vertices[2] + b * vertices[0] + c * vertices[1]) / P

    return incenter


def MPR(shpfile_path):
    """
    计算MPR 最小扰动区域
    :param shpfile_path:输入数据的路径
    :return:最小半径
    """

    # -------------------------设置临时目录--------------------------------
    os.environ['TMPDIR'] = 'C:\\Temp'  # 将临时目录设置为 /tmp

    # -------------------------数据读取--------------------------------
    # 读取数据
    shpfile = gpd.read_file(shpfile_path)
    # 获取数据的嵌套坐标组
    coor_nested, shp_type = get_coor_nested(shpfile)
    # 获取数据的坐标数组
    coor_array = get_coor_array(coor_nested, shp_type)
    # 去重
    coor_array = np.unique(coor_array, axis=0)
    coor_array = coor_array.T

    # 使用 SciPy 生成 Delaunay 三角网
    delaunay_tri = Delaunay(coor_array)

    # -------------------------计算内切圆并找到最小半径--------------------------------
    min_inradius = float('inf')
    circles = []

    for simplex in delaunay_tri.simplices:
        vertices = coor_array[simplex]
        inradius = compute_inradius(vertices)
        incenter = compute_incenter(vertices)

        circles.append((incenter, inradius))

        if inradius < min_inradius: min_inradius = inradius

    print(f"最小的内切圆半径: {min_inradius}")

    # # 绘图
    # plt.figure()
    # plt.triplot(coor_array[:, 0], coor_array[:, 1], delaunay_tri.simplices)
    # plt.plot(coor_array[:, 0], coor_array[:, 1], 'o')
    #
    # # 绘制内切圆
    # # for center, radius in circles:
    # #     circle = plt.Circle(center, radius, color='r', fill=False)
    # #     plt.gca().add_patch(circle)
    #
    # plt.xlabel('X coordinate')
    # plt.ylabel('Y coordinate')
    # plt.title('Delaunay Triangulation')
    # plt.gca().set_aspect('equal', adjustable='box')
    # plt.show()
    return min_inradius


if __name__ == "__main__":
    # -------------------------Data process------------------------------
    # 读取数据
    folder_path = select_folder()
    # 获取文件夹所有后缀为shp的文件路径 保存为数组
    shapefiles = glob.glob(os.path.join(folder_path, '*.shp'))
    for shpfile_path in shapefiles:
        print('------------------------------------\n')
        print(f'当前处理的数据是{os.path.basename(shpfile_path)}')
        MPR(shpfile_path)
        print('\n-----------------\n')

    # -------------------------drawing--------------------------------
    # # 创建图形和坐标轴
    # fig, ax = plt.subplots(figsize=(10, 10))
    # # 绘图
    # gpd.read_file(original_shpfile_path).plot(ax=ax, color='red', linewidth=3)
    # processed_shpfile.plot(ax=ax, color='blue', linewidth=1.5)
    # # 添加图例
    # ax.legend(['Original', 'transformed'])
    # # 显示图形
    # plt.show()