#!/usr/bin/env python
# -*- coding: utf-8 -*-
import re
import math
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from collections import defaultdict
# 设置中文字体支持
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题
plt.rcParams['figure.figsize'] = (12, 8)
plt.rcParams['figure.dpi'] = 100
def edge_heading(p, q):
"""计算两点 p->q 的方向角(弧度,范围 -pi~pi)"""
return math.atan2(q[1]-p[1], q[0]-p[0])
def parse_points_from_cell(val):
"""
解析单元格中的{x,y,z}格式数据
val: 单元格内容
返回: 点列表 [(x, y, z), ...]
"""
if pd.isna(val):
return []
s = str(val).strip()
# 如是 '{x,y,z}',去掉大括号
if s.startswith('{') and s.endswith('}'):
s = s[1:-1]
# 分割字符串
parts = re.split(r'[,\s]+', s)
numbers = []
for part in parts:
try:
numbers.append(float(part))
except ValueError:
continue
points = []
# 每三个数字组成一个点
for i in range(0, len(numbers), 3):
if i + 2 < len(numbers):
x = numbers[i]
y = numbers[i + 1]
z = numbers[i + 2] if i + 2 < len(numbers) else 0.0
points.append((x, y, z))
return points
def compute_turning_points_angle(points, angle_threshold=30):
"""
基于相邻两段的方向角变化,判断转折点
points: 路径上的点列表,形如 [(x0,y0,z0), (x1,y1,z1), ...]
angle_threshold: 阈值,单位为度
返回: 转折点列表和转折点数量
"""
if len(points) < 3:
return [], 0
# 只用 2D 坐标
pts2d = [(p[0], p[1]) for p in points]
heads = [edge_heading(pts2d[i-1], pts2d[i]) for i in range(1, len(pts2d))]
turning_points = []
angles = []
th_rad = math.radians(angle_threshold)
for i in range(1, len(heads)):
a = (heads[i] - heads[i-1] + math.pi) % (2*math.pi) - math.pi
if abs(a) > th_rad:
turning_points.append(pts2d[i])
angles.append(abs(a))
# 去重并返回
return turning_points, len(turning_points)
def compute_turning_points_sliding(points, window=1, threshold=30):
"""
基于滑动窗口的转折点检测
points: 路径上的点列表
window: 左右各取的点数
threshold: 角度阈值,单位为度
返回: 转折点列表和转折点数量
"""
if len(points) < 2*window + 1:
return [], 0
pts2d = [(p[0], p[1]) for p in points]
turning = []
th = math.radians(threshold)
for i in range(window, len(pts2d) - window):
p_left = pts2d[i - window]
p_mid = pts2d[i]
p_right = pts2d[i + window]
v1 = (p_mid[0] - p_left[0], p_mid[1] - p_left[1])
v2 = (p_right[0] - p_mid[0], p_right[1] - p_mid[1])
n1 = math.hypot(v1[0], v1[1])
n2 = math.hypot(v2[0], v2[1])
if n1 < 1e-10 or n2 < 1e-10:
continue
dot = v1[0]*v2[0] + v1[1]*v2[1]
cosang = max(-1.0, min(1.0, dot/(n1*n2)))
ang = math.acos(cosang)
if ang > th:
turning.append(pts2d[i])
return turning, len(turning)
def curvature_three_points(p0, p1, p2):
"""
用三点近似曲率 κ = 4Δ/(a*b*c)
p0, p1, p2: 三个点 (x,y)
返回 κ(曲率,单位 1/长度),若存在零长度则返回0
Δ:三角形面积
a = |p1 - p0|, b = |p2 - p1|, c = |p2 - p0|
"""
a = math.hypot(p1[0]-p0[0], p1[1]-p0[1])
b = math.hypot(p2[0]-p1[0], p2[1]-p1[1])
c = math.hypot(p2[0]-p0[0], p2[1]-p0[1])
if a < 1e-12 or b < 1e-12 or c < 1e-12:
return 0.0
# Δ = |(p1-p0) x (p2-p0)| / 2
cross = (p1[0]-p0[0])*(p2[1]-p0[1]) - (p1[1]-p0[1])*(p2[0]-p0[0])
delta = abs(cross) / 2.0
kappa = 4.0 * delta / (a*b*c)
return kappa
def compute_turning_points_curvature(points, window=1, kappa_thresh=None, percentile=None):
"""
基于曲率的转折点检测
points: 路径点列表,((x,y,z))
window: 使用的滑动窗口宽度 w(左右各取 w 点),实际算的是 p[i-w], p[i], p[i+w]
kappa_thresh: 直接给定的阈值(单位 1/长度),若提供则使用它
percentile: 使用百分位阈值(如 95),若提供则取该分位数作为阈值
返回: 转折点列表和转折点数量
"""
if len(points) < 2*window + 1:
return [], 0
coords2d = [(p[0], p[1]) for p in points]
kappas = []
candidates = []
for i in range(window, len(coords2d) - window):
p0 = coords2d[i - window]
p1 = coords2d[i]
p2 = coords2d[i + window]
k = curvature_three_points(p0, p1, p2)
kappas.append(k)
candidates.append((i, k))
thresh = 0.0
if kappa_thresh is not None:
thresh = kappa_thresh
elif percentile is not None and len(kappas) > 0:
thresh = float(np.percentile(kappas, percentile))
else:
thresh = 0.0
turning = []
for idx, k in candidates:
if k > thresh:
turning.append(coords2d[idx])
return turning, len(turning)
def line_intersection(p1, p2, q1, q2, eps=1e-9):
"""
计算两条线段 p1-p2 与 q1-q2 的交点,若有返回点 (x,y),否则返回 None
使用参数方程方法
"""
x1, y1 = p1
x2, y2 = p2
x3, y3 = q1
x4, y4 = q2
denom = (y4 - y3) * (x2 - x1) - (x4 - x3) * (y2 - y1)
if abs(denom) < eps:
return None # 平行或共线
ua = ((x4 - x3) * (y1 - y3) - (y4 - y3) * (x1 - x3)) / denom
ub = ((x2 - x1) * (y1 - y3) - (y2 - y1) * (x1 - x3)) / denom
if 0 <= ua <= 1 and 0 <= ub <= 1:
x = x1 + ua * (x2 - x1)
y = y1 + ub * (y2 - y1)
return (x, y)
return None
def dedup_points(points, decimals=4):
"""对点集合去重,按给定小数位进行圆整比较"""
seen = set()
uniq = []
for p in points:
key = (round(p[0], decimals), round(p[1], decimals))
if key not in seen:
seen.add(key)
uniq.append((p[0], p[1]))
return uniq
def compute_segments(road_paths):
"""
将 road_paths 转换为线段集合
road_paths: list of paths,路径为 [(x,y,z), ...]
返回:segments = [(p1, p2, path_idx, seg_idx), ...]
"""
segments = []
for pid, path in enumerate(road_paths):
for i in range(len(path) - 1):
p1 = (float(path[i][0]), float(path[i][1]))
p2 = (float(path[i+1][0]), float(path[i+1][1]))
segments.append((p1, p2, pid, i))
return segments
def calculate_intersections(road_paths, eps=1e-9):
"""
计算所有不同线段之间的交点,避免相邻同一路径的线段
road_paths: list of paths,路径为 [(x,y,z), ...]
返回: intersections (list[(x,y)]), count
"""
intersections = []
segments = compute_segments(road_paths)
n = len(segments)
for i in range(n):
p1, p2, pid1, idx1 = segments[i]
for j in range(i+1, n):
q1, q2, pid2, idx2 = segments[j]
# 同一路径且相邻的线段不计为交点
if pid1 == pid2 and abs(idx1 - idx2) == 1:
continue
pt = line_intersection(p1, p2, q1, q2, eps=eps)
if pt is not None:
intersections.append(pt)
# 去重
intersections = dedup_points(intersections, decimals=4)
return intersections, len(intersections)
def visualize_garden_2d(excel_path, sheet_name=0, output_path=None,
method='angle', window=1, threshold=30, percentile=None,
kappa_thresh=None, road_col=0,
entrance=(131100.0, 40400.0), exit_point=(-5100.0, 21800.0)):
"""
可视化花园布局的二维图,计算道路的转折点和交叉点
excel_path: Excel文件路径
sheet_name: 工作表名称或索引
output_path: 输出图像路径(可选)
method: 转折点检测方法:'angle'、'sliding'、'curvature'
window: 滑动窗口宽度(对 sliding 方法有效)
threshold: 角度阈值(对 angle 或 sliding 有效,单位度)
percentile: 当 method 为 'curvature' 时,使用 κ 的 percentile 作为阈值
kappa_thresh: 当 method 为 'curvature' 时,直接给定 κ 的阈值(单位 1/长度)
road_col: 道路数据所在列索引(Excel 第几列存放道路数据)
entrance/exit_point: 入口、出口坐标
"""
try:
print(f"正在读取Excel文件: {excel_path}")
df = pd.read_excel(excel_path, sheet_name=sheet_name, header=None)
print(f"成功读取数据,数据形状: {df.shape}")
print(f"前几行数据:\n{df.head()}")
# 读取道路数据(第一列作为道路数据示例,若有多列可扩展)
road_paths = []
road_col_idx = int(road_col)
for cell in df.iloc[:, road_col_idx]:
pts = parse_points_from_cell(cell)
if pts:
road_paths.append(pts)
if not road_paths:
print("未在指定列解析到有效的道路数据,请检查格式。")
return
# 绘图初始化
fig, ax = plt.subplots(figsize=(12, 8))
# 绘制道路点与路径
road_points_all = []
for path in road_paths:
road_points_all.extend([(p[0], p[1]) for p in path])
xs = [p[0] for p in path]
ys = [p[1] for p in path]
if len(path) >= 2:
ax.plot(xs, ys, color='gray', linewidth=2, alpha=0.7)
if road_points_all:
xs, ys = zip(*road_points_all)
ax.scatter(xs, ys, s=5, c='gray', alpha=0.6, label='道路点')
# 入口/出口
ax.scatter(entrance[0], entrance[1], s=200, c='red', marker='>',
label='入口', alpha=0.9, edgecolors='black')
ax.annotate('入口', xy=entrance, xytext=(10, 10),
textcoords='offset points', fontsize=12, color='red',
bbox=dict(boxstyle="round,pad=0.3", fc="white", ec="red", alpha=0.7))
ax.scatter(exit_point[0], exit_point[1], s=200, c='blue', marker='<',
label='出口', alpha=0.9, edgecolors='black')
ax.annotate('出口', xy=exit_point, xytext=(10, 10),
textcoords='offset points', fontsize=12, color='blue',
bbox=dict(boxstyle="round,pad=0.3", fc="white", ec="blue", alpha=0.7))
# 计算转折点与交叉点
turning_points_all = []
turning_count = 0
bifurcation_points = []
bifurcation_count = 0
if method == 'angle':
for path in road_paths:
tps, tcount = compute_turning_points_angle(path, angle_threshold=threshold)
turning_points_all.extend(tps)
turning_count += tcount
elif method == 'sliding':
for path in road_paths:
tps, tcount = compute_turning_points_sliding(path, window=window, threshold=threshold)
turning_points_all.extend(tps)
turning_count += tcount
elif method == 'curvature':
for path in road_paths:
tps, tcount = compute_turning_points_curvature(
path, window=window, kappa_thresh=kappa_thresh, percentile=percentile
)
turning_points_all.extend(tps)
turning_count += tcount
else:
raise ValueError("method must be one of 'angle', 'sliding', 'curvature'")
# 去重
turning_points_all = dedup_points(turning_points_all, decimals=4)
# 绘制转折点
if turning_points_all:
txs, tys = zip(*turning_points_all)
ax.scatter(txs, tys, c='red', marker='x', s=80, label='转折点')
# 计算并绘制交叉点
intersections, inter_count = calculate_intersections(road_paths)
if intersections:
ix, iy = zip(*intersections)
ax.scatter(ix, iy, c='purple', marker='+', s=60, label='交叉点')
# 显示统计信息
ax.text(0.02, 0.92, f'转折点数量: {turning_count}', transform=ax.transAxes,
fontsize=12, color='red', bbox=dict(boxstyle="round,pad=0.3", fc="white", ec="red", alpha=0.8))
ax.text(0.02, 0.86, f'交叉点数量: {inter_count}', transform=ax.transAxes,
fontsize=12, color='purple', bbox=dict(boxstyle="round,pad=0.3", fc="white", ec="purple", alpha=0.8))
ax.set_xlabel('X坐标')
ax.set_ylabel('Y坐标')
ax.set_title('花园布局二维可视化(含转折点与交叉点)')
ax.legend(loc='best')
ax.grid(True, linestyle='--', alpha=0.7)
ax.axis('equal')
# 保存或显示图像
if output_path:
plt.savefig(output_path, dpi=300, bbox_inches='tight')
print(f"\n图像已保存至: {output_path}")
plt.show()
print("\n=== 统计结果 ===")
print(f"道路转折点数量: {turning_count}")
print(f"道路交叉点数量: {inter_count}")
except Exception as e:
print(f"发生错误: {str(e)}")
import traceback
traceback.print_exc()
print("请检查文件路径和数据格式是否正确")
if __name__ == "__main__":
excel_file = r"D:\数模01\数模01.xlsx" # 替换为您的实际文件路径
visualize_garden_2d(excel_file, output_path="garden_turn_points_curvature_sliding.png",
method='curvature', window=2, percentile=95)
最新发布