谁离你最近？

pip install matplotlib

import math
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon


def clip_halfplane(polygon, a, b, c):
    """
    将多边形裁剪到 ax + by <= c 的半平面内。

    polygon: [(x, y), ...]
    返回裁剪后的多边形顶点。
    """
    if not polygon:
        return []

    def f(p):
        return a * p[0] + b * p[1] - c

    def inside(p):
        return f(p) <= 1e-9

    result = []

    for i in range(len(polygon)):
        p = polygon[i]
        q = polygon[(i + 1) % len(polygon)]

        fp, fq = f(p), f(q)
        pin, qin = inside(p), inside(q)

        # 一条边穿过半平面边界
        if pin != qin:
            t = fp / (fp - fq)

            intersection = (
                p[0] + t * (q[0] - p[0]),
                p[1] + t * (q[1] - p[1])
            )

            result.append(intersection)

        # 边的终点在允许区域中
        if qin:
            result.append(q)

    return result


def voronoi_cell(site, sites, bounds):
    """
    计算某个站点的裁剪后沃罗诺伊单元。

    bounds: (xmin, ymin, xmax, ymax)
    """

    xmin, ymin, xmax, ymax = bounds

    polygon = [
        (xmin, ymin),
        (xmax, ymin),
        (xmax, ymax),
        (xmin, ymax)
    ]

    sx, sy = site

    for tx, ty in sites:
        if (sx, sy) == (tx, ty):
            continue

        # 距离平方比较：
        # (x-sx)^2+(y-sy)^2
        # <= (x-tx)^2+(y-ty)^2

        a = 2 * (tx - sx)
        b = 2 * (ty - sy)

        c = (
            tx * tx + ty * ty
            - sx * sx - sy * sy
        )

        polygon = clip_halfplane(
            polygon, a, b, c
        )

        if not polygon:
            break

    return polygon


def polygon_area(polygon):
    if len(polygon) < 3:
        return 0.0

    total = 0.0

    for i, p in enumerate(polygon):
        q = polygon[(i + 1) % len(polygon)]

        total += (
            p[0] * q[1]
            - q[0] * p[1]
        )

    return abs(total) / 2


# 四家虚构医院
sites = {
    "A": (20, 24),
    "B": (74, 25),
    "C": (29, 78),
    "D": (77, 72),
}

bounds = (0, 0, 100, 100)

all_positions = list(sites.values())

fig, ax = plt.subplots(figsize=(7, 7))

total_area = 0

for name, site in sites.items():
    region = voronoi_cell(
        site,
        all_positions,
        bounds
    )

    area = polygon_area(region)
    total_area += area

    ax.add_patch(
        Polygon(
            region,
            closed=True,
            alpha=0.3,
            edgecolor="black"
        )
    )

    ax.scatter(*site)
    ax.annotate(
        name,
        (site[0] + 2, site[1] + 2)
    )

    print(f"医院 {name}：区域面积 {area:.2f}")

ax.set_xlim(0, 100)
ax.set_ylim(0, 100)
ax.set_aspect("equal")

ax.set_title("Voronoi Diagram")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.grid(alpha=0.15)

plt.show()

print("全部区域面积之和:", total_area)

assert math.isclose(
    total_area,
    10000,
    abs_tol=1e-6
)

import numpy as np
from scipy.spatial import Voronoi

points = np.array([
    [20, 24],
    [74, 25],
    [29, 78],
    [77, 72],
    [49, 49],
])

vor = Voronoi(points)

print("沃罗诺伊顶点：")
print(vor.vertices)

print("相邻站点对：")
print(vor.ridge_points)

print("区域顶点索引：")
print(vor.regions)