Shapely中多边形与圆形集合求差异常问题求助

我需要测量圆形集群中面积大于目标值的空区域（vacancy），预期算法流程：

• 定义带x、y坐标的圆形
• 生成包围圆形集群的凸包
• 合并所有圆形后从凸包中减去
• 识别面积大于目标值的vacancy并丢弃其余区域

但使用Shapely的difference方法得到的结果不符合预期，仅临时将圆形buffer值从0.5改为0.501时能得到接近预期的效果，但该方法不适用于实际场景（实际圆形位置间距更大）。

核心问题是Shapely几何运算的浮点数精度限制：

• 当圆形紧密排列、边界相切时，Shapely无法精准识别相切处的拓扑边界，导致difference运算无法正确分割出内部空区域
• 修改buffer值为0.501是通过微小膨胀让圆形从相切变为相交，间接规避了精度问题，但属于临时方案，无法适配间距更大的场景

针对精度问题，采用两种可靠的拓扑修复+正确空区域提取方式：

1. 使用buffer(0)修复几何对象：对合并后的圆形区域和凸包执行buffer(0)，可自动修复因精度导致的拓扑错误（如相切边界的模糊问题）
2. 直接提取空区域几何部件：不再通过polygonize处理边界，而是直接遍历difference返回的Polygon/MultiPolygon对象，避免额外误差

from shapely.geometry import Polygon, Point
from shapely.ops import unary_union
import matplotlib.pyplot as plt
import math
import numpy as np
import matplotlib as mpl
mpl.use("qt5agg")

def generate_hexagonal_grid(side_length, diameter):
    points = []
    radius = diameter / 2
    vertical_distance = np.sqrt(3) * radius
    x_max = side_length * np.cos(np.pi / 6)

    row = 0
    while (row * vertical_distance) < (2 * x_max):
        x_offset = diameter if row % 2 == 0 else radius
        x = -x_max + x_offset
        while x < x_max:
            y = -x_max + row * vertical_distance
            if np.abs(y) <= x_max:
                points.append((x, y))
            x += diameter
        row += 1
    return points

def draw_flat_topped_hexagon(apothem):
    s = 2 * apothem * math.tan(math.pi / 6)
    vertices = [
        (-apothem * math.tan(math.pi / 6), apothem),
        (apothem * math.tan(math.pi / 6), apothem),
        (2 * apothem * math.tan(math.pi / 6), 0),
        (apothem * math.tan(math.pi / 6), -apothem),
        (-apothem * math.tan(math.pi / 6), -apothem),
        (-2 * apothem * math.tan(math.pi / 6), 0)
    ]
    return Polygon(vertices)

def fill_hexagon_with_circles(hexagon, points):
    fig, ax = plt.subplots()
    ax.set_facecolor('black')

    circles = []
    occupied_area_center = []

    for idx, (x, y) in enumerate(points):
        if idx != int(len(points)/2):
            point = Point(x, y)
            circle = point.buffer(0.5)
            if circle.centroid not in occupied_area_center and hexagon.contains(circle):
                occupied_area_center.append(circle.centroid)
                circles.append(circle)
                cx, cy = circle.exterior.xy
                ax.fill(cx, cy, alpha=1, fc='yellow', edgecolor='black')

    # 合并圆形并修复拓扑错误
    occupied_area = unary_union(circles).buffer(0)
    # 生成凸包并修复拓扑错误
    surrounding_polygon = occupied_area.convex_hull.buffer(0)

    bead_area = math.pi * (0.5) ** 2
    max_packing_fraction = math.pi / (2 * math.sqrt(3))
    target_area = bead_area / max_packing_fraction

    # 直接提取空区域，区分单个Polygon和MultiPolygon
    vacant_space = surrounding_polygon.difference(occupied_area)
    valid_vacancies = []
    if vacant_space.geom_type == 'MultiPolygon':
        for poly in vacant_space.geoms:
            if poly.area > target_area:
                valid_vacancies.append(poly)
    elif vacant_space.geom_type == 'Polygon':
        if vacant_space.area > target_area:
            valid_vacancies.append(vacant_space)

    total_vacant_area = sum(poly.area for poly in valid_vacancies)
    for i, vacancy in enumerate(valid_vacancies):
        vx, vy = vacancy.exterior.xy
        ax.fill(vx, vy, alpha=0.5, fc='red', label=f'Vacancy {i + 1}')
    plt.legend()
    plt.show()

hex_diameter = 12
circle_diameter = 1.0
side_length = hex_diameter * np.sqrt(3) / 2
points = generate_hexagonal_grid(side_length, circle_diameter)
hexagon = draw_flat_topped_hexagon(hex_diameter)
fill_hexagon_with_circles(hexagon, points)

• 对合并后的圆形区域和凸包执行buffer(0)，修复相切边界的拓扑精度问题
• 直接处理difference返回的几何对象，替代polygonize处理边界的方式，避免额外误差
• 优化循环逻辑，提升代码可读性

内容的提问来源于stack exchange，提问作者mugenop