Python 多边形和圆集合之间的形状差异无法正常工作

Question

我正在尝试测量一组圆圈内的空位。 空位是表面积大于目标面积的空白区域 我的目标算法应该像这样工作：

1. 我首先用坐标 x,y 定义圆。
2. 然后定义包围圆簇的凸包。
3. 我统一这些圆，然后从凸包中减去它们。
4. 我识别出表面积大于目标区域的空缺，并丢弃其余的（将它们填充为背景颜色。

为了帮助完成此练习，我每次都会绘制结果。

这是我写的测试代码：

from shapely.geometry import Polygon, Point
from shapely.ops import unary_union, polygonize
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)

    # Creating a hexagonal grid of points within the hexagon
    row = 0
    while (row * vertical_distance) < (2 * x_max):
        if row % 2 == 0:
            x_offset = diameter
        else:
            x_offset = 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)
    ]
    hexagon = Polygon(vertices)
    return hexagon
def fill_hexagon_with_circles(hexagon, points):
    # Initialize a plot
    fig, ax = plt.subplots()
    
    # Set the background color of the plot to black
    ax.set_facecolor('black')

    # List to store circle geometries
    circles = []
    occupied_area_center = []

    # Plot circles within the hexagon with black color
    i = 0
    for x, y in points:

        if i!=int(len(points)/2):
            point = Point(x, y)
            circle = point.buffer(0.5)  # Circle with diameter 1
            if circle.centroid not in occupied_area_center:
                if hexagon.contains(circle):
                    occupied_area_center.append(circle.centroid)
                    circles.append(circle)
                    x, y = circle.exterior.xy
                    ax.fill(x, y, alpha=1, fc='yellow', edgecolor='black')
        i+=1
    occupied_area = unary_union(circles)
    
    # Create a polygon that surrounds the circles by creating a convex hull
    surrounding_polygon = unary_union(circles).convex_hull

    # Calculate the surface area of a bead
    bead_area = math.pi * (0.5) ** 2

    # Calculate the target area
    max_packing_fraction = math.pi / (2 * math.sqrt(3))
    target_area = bead_area / max_packing_fraction
    vacant_space = surrounding_polygon.difference(occupied_area)
    vacant_areas = list(polygonize(vacant_space.boundary))
    valid_vacancies = [poly for poly in vacant_areas if poly.area > target_area]
    total_vacant_area = sum(poly.area for poly in valid_vacancies)
    for i, vacancy in enumerate(valid_vacancies):
        x, y = vacancy.exterior.xy
        ax.fill(x, y, alpha=0.5, fc='red', label=f'Vacancy {i + 1}')
    plt.show()


hex_diameter = 12
circle_diameter = 1.0

# Calculate the side length of the hexagon
side_length = hex_diameter * np.sqrt(3) / 2

# Generate a grid of points within the hexagon
points = generate_hexagonal_grid(side_length, circle_diameter)

# Create hexagon and fill it with circles
hexagon = draw_flat_topped_hexagon(hex_diameter)
fill_hexagon_with_circles(hexagon, points)

结果如下所示：
我期待更多这样的事情：

一条重要信息：为了生成所需的结果，我修改了这一行：

circle = point.buffer(0.5)  # Circle with diameter 1 

到

circle = point.buffer(0.501)  # Circle with diameter 1

这些都是假位置，实际位置会更远，所以这个技巧不起作用。

任何帮助将不胜感激。

Answer 1

主要问题是许多/一些相邻的空置区域保持连接，因此它们的区域被合并......

为了避免这种情况，我会寻找一种结构化的方法来将“空白区域”细分为逻辑部分。

我尝试使用基于圆的中点的德劳内三角形，这样就有可能达到你预期的结果：

import shapely
import shapely.plotting as sh_plot
from shapely.geometry import MultiLineString, MultiPoint, Polygon, Point
from shapely.ops import unary_union, polygonize
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)

    # Creating a hexagonal grid of points within the hexagon
    row = 0
    while (row * vertical_distance) < (2 * x_max):
        if row % 2 == 0:
            x_offset = diameter
        else:
            x_offset = 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),
    ]
    hexagon = Polygon(vertices)
    return hexagon


def fill_hexagon_with_circles(hexagon, points):
    # Initialize a plot
    fig, ax = plt.subplots()

    # Set the background color of the plot to black
    ax.set_facecolor("black")

    # List to store circle geometries
    circles = []
    occupied_area_center = []

    # Plot circles within the hexagon with black color
    i = 0
    for x, y in points:
        if i != int(len(points) / 2):
            point = Point(x, y)
            circle = point.buffer(0.5)  # Circle with diameter 1
            if circle.centroid not in occupied_area_center:
                if hexagon.contains(circle):
                    occupied_area_center.append(circle.centroid)
                    circles.append(circle)
                    x, y = circle.exterior.xy
                    ax.fill(x, y, alpha=1, fc="yellow", edgecolor="black")
        i += 1
    occupied_area = unary_union(circles)

    # Create a polygon that surrounds the circles by creating a convex hull
    surrounding_polygon = unary_union(circles).convex_hull

    # Calculate the surface area of a bead
    bead_area = math.pi * (0.5) ** 2

    # Calculate the target area
    max_packing_fraction = math.pi / (2 * math.sqrt(3))
    target_area = bead_area / max_packing_fraction
    vacant_space = surrounding_polygon.difference(occupied_area)

    # Subdivide the vacant space using delaunay triangles. The target area needs to
    # be divided by the number of triangles
    delaunays = shapely.get_parts(shapely.delaunay_triangles(MultiPoint(points)))
    vacant_areas = shapely.intersection(vacant_space, delaunays)
    # vacant_areas = list(polygonize(vacant_space.boundary))
    target_area /= 6

    # Plot points and delaunay edges to illustrate what is happening
    # delaunay_edges = shapely.delaunay_triangles(MultiPoint(points), only_edges=True)
    # sh_plot.plot_line(MultiLineString(delaunay_edges), color="green")
    # sh_plot.plot_points(MultiPoint(points), color="orange")

    # Only retain sufficiently large areas
    valid_vacancies = [poly for poly in vacant_areas if poly.area > target_area]
    total_vacant_area = sum(poly.area for poly in valid_vacancies)
    for i, vacancy in enumerate(valid_vacancies):
        x, y = vacancy.exterior.xy
        ax.fill(x, y, alpha=0.5, fc="red", label=f"Vacancy {i + 1}")
    plt.show()


hex_diameter = 12
circle_diameter = 1.0

# Calculate the side length of the hexagon
side_length = hex_diameter * np.sqrt(3) / 2

# Generate a grid of points within the hexagon
points = generate_hexagonal_grid(side_length, circle_diameter)

# Create hexagon and fill it with circles
hexagon = draw_flat_topped_hexagon(hex_diameter)
fill_hexagon_with_circles(hexagon, points)

为了便于说明，橙色点是点，绿色线是用于细分空白区域的德劳奈三角形。正如你所看到的，较大的空白区域被细分为 6 部分，这就是我将 target_area 除以 6 的原因。