179 lines
6.6 KiB
Python
179 lines
6.6 KiB
Python
"""问题一独立参考解：大框 Sutherland–Hodgman 半平面裁剪 + 精确直径/MEC。
|
||
|
||
与生产实现 ``mcm26b.geometry.localization.solve`` 的独立性：
|
||
|
||
- 生产实现枚举"半平面边界直线两两交点"再筛可行点；
|
||
- 本参考解从一个大方框出发，逐条半平面做 Sutherland–Hodgman 裁剪，
|
||
是另一条构造路径，且这里的半平面系数、直径、最小包围圆都在本文件内
|
||
独立编码（不复用生产实现的 ``halfplanes`` / ``diameter`` / ``mec``）。
|
||
|
||
只要区域有界，大框（±1e7 m）裁剪得到的顶点即精确多边形；顶点触到大框
|
||
边界则判无界。最小包围圆用"枚举所有 1/2/3 点支撑圆"的暴力算法，直径用
|
||
顶点对暴力枚举，二者都是 O(m²)/O(m³) 的朴素精确算法，m 为多边形顶点数。
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
import math
|
||
from time import perf_counter
|
||
|
||
from project.common import Method
|
||
|
||
METHOD = Method(
|
||
key="reference",
|
||
name="大框半平面裁剪 + 精确直径/MEC",
|
||
kind="reference",
|
||
note="独立于生产顶点枚举的精确参考解",
|
||
)
|
||
|
||
BOX_M = 1.0e7
|
||
_PARALLEL_EPS = 1e-12
|
||
|
||
|
||
def _planes(case: dict) -> list[tuple[float, float, float]]:
|
||
"""独立重写 a*x+b*y<=c 的半平面构造（与题面角域定义一致）。"""
|
||
planes = []
|
||
for row in case["measurements"]:
|
||
half = math.radians(row.get("half_angle_deg", 1.0))
|
||
theta = math.radians(row["bearing_deg"] % 360.0)
|
||
sx, sy = row["x"], row["y"]
|
||
for angle, sign in ((theta - half, 1.0), (theta + half, -1.0)):
|
||
a = sign * math.sin(angle)
|
||
b = -sign * math.cos(angle)
|
||
planes.append((a, b, a * sx + b * sy))
|
||
# 前向约束：不允许把角域取成反向无限双楔
|
||
a, b = -math.cos(theta), -math.sin(theta)
|
||
planes.append((a, b, a * sx + b * sy))
|
||
return planes
|
||
|
||
|
||
def _clip(polygon, a, b, c, tol=1e-7):
|
||
"""Sutherland–Hodgman：保留 a*x+b*y<=c 的一侧。"""
|
||
if not polygon:
|
||
return []
|
||
output = []
|
||
count = len(polygon)
|
||
for index in range(count):
|
||
p, q = polygon[index], polygon[(index + 1) % count]
|
||
fp = a * p[0] + b * p[1] - c
|
||
fq = a * q[0] + b * q[1] - c
|
||
p_in, q_in = fp <= tol, fq <= tol
|
||
if p_in:
|
||
output.append(p)
|
||
if p_in != q_in:
|
||
ratio = fp / (fp - fq)
|
||
output.append((p[0] + ratio * (q[0] - p[0]),
|
||
p[1] + ratio * (q[1] - p[1])))
|
||
deduped = []
|
||
for point in output:
|
||
if not deduped or math.dist(point, deduped[-1]) > 1e-9:
|
||
deduped.append(point)
|
||
if len(deduped) > 1 and math.dist(deduped[0], deduped[-1]) <= 1e-9:
|
||
deduped.pop()
|
||
return deduped
|
||
|
||
|
||
def _brute_force_diameter(points):
|
||
best = 0.0
|
||
pair = (points[0], points[0])
|
||
for i, pi in enumerate(points):
|
||
for pj in points[i + 1:]:
|
||
value = math.dist(pi, pj)
|
||
if value > best:
|
||
best, pair = value, (pi, pj)
|
||
return best, pair
|
||
|
||
|
||
def _circumcircle(a, b, c):
|
||
ux, uy = b[0] - a[0], b[1] - a[1]
|
||
vx, vy = c[0] - a[0], c[1] - a[1]
|
||
det = 2.0 * (ux * vy - uy * vx)
|
||
if abs(det) <= _PARALLEL_EPS * max(
|
||
math.hypot(ux, uy) * math.hypot(vx, vy), 1e-300):
|
||
return None
|
||
uu, vv = ux * ux + uy * uy, vx * vx + vy * vy
|
||
center = (a[0] + (uu * vy - vv * uy) / det,
|
||
a[1] + (ux * vv - vx * uu) / det)
|
||
return center, math.dist(center, a)
|
||
|
||
|
||
def _brute_force_mec(points):
|
||
"""枚举 1/2/3 点支撑圆的精确最小包围圆（朴素但无迭代）。"""
|
||
xs = [p[0] for p in points]
|
||
ys = [p[1] for p in points]
|
||
best_c = ((min(xs) + max(xs)) / 2.0, (min(ys) + max(ys)) / 2.0)
|
||
best_r = max(math.dist(best_c, p) for p in points)
|
||
candidates: list[tuple[tuple[float, float], float]] = []
|
||
for i, p in enumerate(points):
|
||
for q in points[i + 1:]:
|
||
candidates.append((((p[0] + q[0]) / 2.0, (p[1] + q[1]) / 2.0),
|
||
math.dist(p, q) / 2.0))
|
||
for i, p in enumerate(points):
|
||
for j in range(i + 1, len(points)):
|
||
for k in range(j + 1, len(points)):
|
||
circle = _circumcircle(p, points[j], points[k])
|
||
if circle is not None:
|
||
candidates.append(circle)
|
||
for center, radius in candidates:
|
||
if radius >= best_r:
|
||
continue
|
||
if all(math.dist(center, p) <= radius + 1e-7 for p in points):
|
||
best_c, best_r = center, radius
|
||
return best_c, best_r
|
||
|
||
|
||
def solve_case(case: dict) -> dict:
|
||
start = perf_counter()
|
||
planes = _planes(case)
|
||
polygon = [(-BOX_M, -BOX_M), (BOX_M, -BOX_M), (BOX_M, BOX_M), (-BOX_M, BOX_M)]
|
||
for a, b, c in planes:
|
||
polygon = _clip(polygon, a, b, c)
|
||
if not polygon:
|
||
return _finish(_degenerate("empty"), start)
|
||
boundary_tol = 1e-6 * BOX_M
|
||
if any(abs(point[0]) >= BOX_M - boundary_tol
|
||
or abs(point[1]) >= BOX_M - boundary_tol for point in polygon):
|
||
return _finish(_degenerate("unbounded"), start)
|
||
|
||
if len(polygon) == 1:
|
||
summary = {"status": "point", "diameter_m": 0.0, "covers": True,
|
||
"mec_radius_m": 0.0, "max_excess_m": 0.0,
|
||
"vertex_count": 1}
|
||
else:
|
||
value, (a, b) = _brute_force_diameter(polygon)
|
||
if len(polygon) == 2:
|
||
summary = {"status": "segment", "diameter_m": value, "covers": True,
|
||
"mec_radius_m": value / 2.0, "max_excess_m": 0.0,
|
||
"vertex_count": 2}
|
||
else:
|
||
_mec_center, radius = _brute_force_mec(polygon)
|
||
diameter_center = ((a[0] + b[0]) / 2.0, (a[1] + b[1]) / 2.0)
|
||
max_distance = max(math.dist(diameter_center, p) for p in polygon)
|
||
summary = {
|
||
"status": "polygon",
|
||
"diameter_m": value,
|
||
"covers": bool(radius <= value / 2.0 + 1e-7),
|
||
"mec_radius_m": radius,
|
||
"max_excess_m": round(max_distance - value / 2.0, 9),
|
||
"vertex_count": len(polygon),
|
||
}
|
||
return _finish(summary, start)
|
||
|
||
|
||
def _degenerate(status: str) -> dict:
|
||
return {"status": status, "diameter_m": None, "covers": None,
|
||
"mec_radius_m": None, "max_excess_m": None, "vertex_count": 0}
|
||
|
||
|
||
def _finish(summary: dict, start: float) -> dict:
|
||
result = {
|
||
"method": METHOD.key,
|
||
"algorithm": METHOD.name,
|
||
"resolution_m": 0.0,
|
||
"diameter_infinite": summary["status"] == "unbounded",
|
||
"outside_vertex_count": None,
|
||
"runtime_ms": round((perf_counter() - start) * 1000.0, 4),
|
||
}
|
||
result.update(summary)
|
||
return result
|