12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394 |
"""演示8: 面 — 长方体最小环基识别 6 个表面 + 平面拟合校验。
已知一个长方体 3D 线框(8 顶点、12 条边), 用图论方法:
1. 构建顶点-边图;
2. 求最小环基(Minimum Cycle Basis), 期望得到 6 个四边形环;
3. 对每个环做空间平面拟合, 验证共面性(判定为平面而非曲面)。
"""
import os
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from figure_style import figures_dir, configure
configure()
OUT = figures_dir()
os.makedirs(OUT, exist_ok=True)
# ---- 长方体 8 顶点(尺寸 2x1x1.5) ----
w, h, d = 2.0, 1.5, 1.0
verts3d = np.array([
[0, 0, 0], [w, 0, 0], [w, h, 0], [0, h, 0], # 底面
[0, 0, d], [w, 0, d], [w, h, d], [0, h, d], # 顶面
])
# ---- 12 条棱(边) ----
edges = [
(0, 1), (1, 2), (2, 3), (3, 0), # 底环
(4, 5), (5, 6), (6, 7), (7, 4), # 顶环
(0, 4), (1, 5), (2, 6), (3, 7), # 垂直棱
]
G = nx.Graph()
G.add_edges_from(edges)
# ---- 找所有四边形表面: 枚举简单 4 环并去重 ----
# 说明: minimum_cycle_basis 返回的是环空间的一组基(维数 E-V+1=5), 而长方体的
# 6 个表面是平面嵌入的"面"; 这里枚举所有简单四边形环, 恰好得到全部 6 个表面。
DG = nx.DiGraph()
for u, v in edges:
DG.add_edge(u, v); DG.add_edge(v, u)
quad_cycles = {frozenset(c) for c in nx.simple_cycles(DG) if len(c) == 4}
cycles = [list(c) for c in quad_cycles]
print(f"四边形环个数: {len(cycles)} (期望 6 = 长方体的 6 个面)")
# ---- 平面拟合校验: 环顶点是否共面 ----
def plane_fit(pts):
"""最小二乘拟合平面, 返回 (法向量, 残差)。"""
centroid = pts.mean(axis=0)
_, s, Vt = np.linalg.svd(pts - centroid)
normal = Vt[-1]
res = np.linalg.norm((pts - centroid) @ normal)
return normal, res
faces = []
for cyc in cycles:
pts = verts3d[list(cyc)]
n, res = plane_fit(pts)
faces.append((list(cyc), n, res))
print(f" 环 {len(faces)}: 顶点 {len(cyc)} 个, 平面拟合残差 {res:.4f}"
f" -> {'平面' if res < 1e-6 else '非平面(曲面?)'}")
# ---- 可视化 ----
fig = plt.figure(figsize=(13, 6))
ax1 = fig.add_subplot(1, 2, 1, projection="3d")
for u, v in edges:
ax1.plot(*zip(verts3d[u], verts3d[v]), "k-", lw=1.5)
ax1.scatter(*verts3d.T, c="tab:blue", s=40)
ax1.set_title("长方体 3D 线框(8 顶点 / 12 条边)")
ax2 = fig.add_subplot(1, 2, 2, projection="3d")
cmap = plt.cm.tab10
for i, (cyc, n, res) in enumerate(faces):
poly = verts3d[cyc]
tri = ax2.add_collection3d(Poly3DCollection([poly], alpha=0.35,
facecolor=cmap(i % 10),
edgecolor="k", lw=1))
for u, v in edges:
ax2.plot(*zip(verts3d[u], verts3d[v]), "k-", lw=1.2)
ax2.set_title(f"四边形环识别: {len(faces)} 个表面(平面拟合校验通过)")
for ax in (ax1, ax2):
ax.set_xlabel("X"); ax.set_ylabel("Y"); ax.set_zlabel("Z")
ax.set_xlim(-0.5, 2.5); ax.set_ylim(-0.5, 2.0); ax.set_zlim(-0.5, 1.6)
fig.suptitle("三维识别·第三步: 面 — 最小环基 + 平面拟合识别长方体表面", fontsize=14)
fig.tight_layout(rect=(0, 0, 1, 0.95))
out_path = os.path.join(OUT, "fig_3d_box_faces.png")
fig.savefig(out_path, dpi=150)
print(f"-> {out_path}")
|
