新算法

This commit is contained in:
hesuicong committed 2026-11-12 11:27:57 +08:00
1 parent abb81f6bfd
commit 2434c4140d
1 file changed
+482
+482
View File
@@ -0,0 +1,482 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
make_shell.py — 为每个眼镜实例的顶点连通区域生成"最小壳" (minimal tight shell)。
输入支持两种方式:
1. data.json (推荐, 支持多眼镜实例, 每个实例单独处理并各导出一个壳):
{"input": "xxx.obj",
"instances": [
{"id": 0, "type": "ring-pair", "score": 0.91, "count": 1234,
"indices": [12, 13, 14, ...]},
{"id": 1, "type": "ring-pair", "score": 0.87, "count": 1180,
"indices": [5012, 5013, ...]}
]}
"input" 相对 data.json 所在目录解析 (解析不到时按同名 obj 兜底); 也可用 --obj 显式指定 (优先级更高)。
各实例独立计算, 壳最后合并导出为单个 obj (默认 shell.obj)。
2. index.txt (兼容旧的单实例流程): 顶点索引, 每行一个, 输出单个 shell.obj。
不是 AABB, 也不是带 padding 的膨胀包裹, 而是 alpha shape (3D 凹包):
- 壳面由区域最外侧的原始采样点构成, 紧贴外侧轮廓 (含凹陷, 如眼镜框);
- alpha 自动取"使壳保持单一连通的最小值" (二分搜索), 即最小壳;
- 桥接只发生在顶点漏采的空洞处, 桥的粗细由最小 alpha 自然决定;
- 最后做可选的 Taubin 平滑 (不收缩), 兼顾"最小"与"光滑"。
流程 (每个实例独立执行一遍):
1. 读 obj (全局一次) 与实例顶点索引 (默认 0-based), 取区域顶点 +
"三顶点均在区域内"的三角形重心/边中点加密 -> 采样点;
2. 3D Delaunay 四面体化, 计算每个四面体的外接球半径;
3. 二分搜索最小 alpha: alpha 复合体 (外接球半径 <= alpha 的四面体之并)
满足 (a) 连通 (b) 覆盖全部采样点;
4. 若最小 alpha 下边界非流形/不水密, alpha 向上微调至水密;
5. 抽复合体的边界三角形 (外向绕序), 得到紧贴的壳;
6. 可选 Taubin 平滑 + 数值校验 (点偏移, 体积, 水密性);
7. 全部实例的壳合并 (同一坐标系直接拼接) 导出单个 obj。
用法:
python make_shell.py # 读同目录 data.json, 各实例合并导出 shell.obj
python make_shell.py --data other.json # 指定 data.json
python make_shell.py --index index.txt # 兼容旧单索引, 导出 shell.obj
python make_shell.py --alpha 3.0 # 手动指定 alpha (更大=更圆润)
python make_shell.py --smooth 0 # 关闭平滑, 严格贴点
"""
import argparse
import glob
import json
import os
import numpy as np
from scipy.spatial import Delaunay, cKDTree
from scipy.sparse import coo_matrix
from scipy.sparse.csgraph import connected_components
from scipy import ndimage
from skimage import measure
import trimesh
from trimesh.smoothing import filter_taubin
# ---------------------------------------------------------------- IO
def load_obj(path):
verts, faces = [], []
with open(path) as f:
for line in f:
if line.startswith("v "):
verts.append([float(x) for x in line.split()[1:4]])
elif line.startswith("f "):
faces.append([int(t.split("/")[0]) - 1 for t in line.split()[1:4]])
return np.asarray(verts, np.float64), np.asarray(faces, np.int64)
def load_data(args):
"""解析输入, 返回 (instances, input_obj)。
instances: [(实例名, 顶点索引 ndarray), ...]
input_obj: data.json 声明的 obj 路径 (未解析相对路径, 可能为 None)
优先级: --index 单实例 txt > --data 多实例 json。
"""
if args.index:
idx = np.loadtxt(args.index, dtype=np.int64).ravel()
if args.one_based:
idx = idx - 1
return [("index", idx)], None
if not os.path.exists(args.data):
raise SystemExit(f"找不到 {args.data}, 请用 --data 指定 data.json, "
f"或用 --index 指定单实例索引文件")
with open(args.data) as f:
data = json.load(f)
insts = []
for inst in data.get("instances", []):
name = str(inst.get("id", len(insts)))
idx = np.asarray(inst.get("indices", []), dtype=np.int64).ravel()
if args.one_based:
idx = idx - 1
cnt = inst.get("count")
if cnt is not None and cnt != len(idx):
print(f"[warn] 实例 {name}: count={cnt} 与 indices 数量 {len(idx)} 不一致")
if len(idx) == 0:
print(f"[warn] 实例 {name}: 索引为空, 跳过")
continue
insts.append((name, idx))
if not insts:
raise SystemExit(f"{args.data} 中没有可用实例 (instances 为空)")
return insts, data.get("input")
def resolve_obj(args, input_obj):
"""确定 obj 路径: --obj 显式指定 > data.json 的 input 字段 > data 目录内唯一 obj。
input 字段可能是相对其他根目录的路径 (如 "datasets/xxx/xxx.obj"), 本地解析不到时
按 basename 在 data 目录内找同名 obj 兜底。
"""
if args.obj:
return args.obj
base = os.path.dirname(os.path.abspath(args.data))
if input_obj:
p = input_obj if os.path.isabs(input_obj) else os.path.join(base, input_obj)
if os.path.exists(p):
return p
name = os.path.basename(input_obj) # 按文件名在 data 目录内找同名 obj
cands = sorted(glob.glob(os.path.join(base, name)))
if len(cands) == 1:
print(f"[warn] data.json 的 input '{input_obj}' 本地不存在, "
f"使用同名 obj: {cands[0]}")
return cands[0]
print(f"[warn] data.json 的 input '{input_obj}' 本地不存在, 尝试目录内查找 obj")
cands = sorted(glob.glob(os.path.join(base, "*.obj")))
if len(cands) == 1:
return cands[0]
if not cands:
raise SystemExit(f"未找到 obj: 请用 --obj 指定, 或在 {base} 放置 .obj 文件")
raise SystemExit(f"{base} 下有多个 obj, 请用 --obj 指定其一: "
f"{', '.join(os.path.basename(c) for c in cands)}")
def densify_points(verts, faces, region_mask):
"""区域顶点 + 区域内部三角形的重心/边中点采样点。"""
pts = [verts[region_mask]]
inner = faces[np.all(region_mask[faces], axis=1)]
if len(inner):
tris = verts[inner] # (F,3,3)
for a, b, c in ((1/3, 1/3, 1/3),
(0.5, 0.5, 0.0), (0.5, 0.0, 0.5), (0.0, 0.5, 0.5)):
pts.append(a * tris[:, 0] + b * tris[:, 1] + c * tris[:, 2])
pts = np.concatenate(pts, axis=0)
return np.unique(pts, axis=0), inner # 去重, 避免退化四面体
def largest_cluster(pts, r):
"""半径 r 欧氏聚类, 返回最大簇的布尔掩码 (剔除分割噪声)。"""
tree = cKDTree(pts)
pairs = np.array(list(tree.query_pairs(r=r)), dtype=np.int64).reshape(-1, 2)
n = len(pts)
if len(pairs):
m = coo_matrix((np.ones(len(pairs)), (pairs[:, 0], pairs[:, 1])), shape=(n, n))
_, lab = connected_components(m, directed=False)
sizes = np.bincount(lab)
return lab == int(np.argmax(sizes)), int(sizes.max())
return np.ones(n, bool), n
# ------------------------------------------------------ Delaunay 预处理
class AlphaStructure:
"""Delaunay 四面体化 + 面邻接/外接球半径, 供 alpha 二分搜索复用。"""
def __init__(self, pts):
self.pts = pts
tri = Delaunay(pts)
self.tri = tri # find_simplex 用于凸包内部栅格化
tets = tri.simplices.copy()
a = pts[tets[:, 0]]
m = np.stack([pts[tets[:, 1]] - a, pts[tets[:, 2]] - a, pts[tets[:, 3]] - a], axis=1)
det = np.linalg.det(m)
scale = float(np.abs(det).max()) ** (1 / 3) + 1e-30
good = np.abs(det) > 1e-10 * scale # 剔除退化四面体
tets, det, a = tets[good], det[good], a[good]
tets[det < 0] = tets[det < 0][:, [0, 1, 3, 2]] # 统一正体积绕序
self.tets = tets
self.T = len(tets)
# 外接球半径: 解 2(b-a).x = |b|^2-|a|^2 线性方程组
e = [np.einsum("ij,ij->i", pts[tets[:, k]], pts[tets[:, k]]) for k in range(4)]
A = 2.0 * np.stack([pts[tets[:, 1]] - a, pts[tets[:, 2]] - a, pts[tets[:, 3]] - a], axis=1)
rhs = np.stack([e[1] - e[0], e[2] - e[0], e[3] - e[0]], axis=1)
center = np.linalg.solve(A, rhs)
self.circum = np.linalg.norm(center - a, axis=1)
# 每个四面体的 4 个面 (tet 优先排列: 线性位置 t*4+f), 外向绕序 + 面共享关系
faces = np.concatenate([tets[:, [1, 2, 3]], tets[:, [0, 3, 2]],
tets[:, [0, 1, 3]], tets[:, [0, 2, 1]]],
axis=1).reshape(self.T * 4, 3)
fcent = pts[faces].mean(axis=1)
tcent = np.repeat(pts[tets].mean(axis=1), 4, axis=0)
n = np.cross(pts[faces[:, 1]] - pts[faces[:, 0]], pts[faces[:, 2]] - pts[faces[:, 0]])
flip = np.einsum("ij,ij->i", n, fcent - tcent) < 0
faces[flip] = faces[flip][:, [0, 2, 1]] # 法向朝外
self.faces = faces
keys = np.sort(self.faces, axis=1)
_, inv = np.unique(keys, axis=0, return_inverse=True)
self.face_id = np.asarray(inv).reshape(-1)
cnt = np.bincount(self.face_id)
shared = np.flatnonzero(cnt == 2)
# 每个共享面 -> 一对四面体 (面在线性数组中的两个出现位置)
order = np.argsort(self.face_id, kind="stable")
inv_sorted = self.face_id[order]
starts = np.searchsorted(inv_sorted, shared, side="left")
self.pairs = np.stack([order[starts], order[starts + 1]], axis=1) // 4
self.n_pts = len(pts)
# ---- 给定 alpha 的复合体状态 --------------------------------------
def status(self, alpha):
keep = self.circum <= alpha
used = np.zeros(self.n_pts, bool)
used[self.tets[keep].ravel()] = True
ok_used = bool(used.all())
pe = keep[self.pairs[:, 0]] & keep[self.pairs[:, 1]]
g = coo_matrix((np.ones(int(pe.sum()), bool),
(self.pairs[pe, 0], self.pairs[pe, 1])), shape=(self.T, self.T))
nc, lab = connected_components(g, directed=False)
sizes = np.bincount(lab, weights=keep.astype(np.int64))
ok_conn = nc - int((~keep).sum()) == 1 and int(sizes.max()) == int(keep.sum())
return keep, ok_used, ok_conn
def boundary_faces(self, keep):
"""复合体的边界三角形 (恰属于一个保留四面体的面, 外向绕序)。"""
mask = np.repeat(keep, 4) # faces 按 [tet 的 4 面] 排列
fid = self.face_id[mask]
cnt = np.bincount(fid) # 该面在保留四面体中的出现次数
return self.faces[mask][cnt[fid] == 1]
# ------------------------------------------------------ 眼镜盒式外壳
def build_case(pts, tri, res, sigma, faces_target):
"""凸包(最外侧点连成的最小凸外壳) -> 内部栅格化 -> 高斯平滑 -> MC。
壳贴着最外侧点(误差 ~ 1 voxel), 边角以 sigma*voxel 尺度圆润, 均匀光滑, 水密。
"""
hull = trimesh.Trimesh(pts).convex_hull
lo, hi = pts.min(axis=0), pts.max(axis=0)
voxel = (hi - lo).max() / res
pad = 6 * sigma * voxel + 0.5
lo2 = lo - pad
dims = np.maximum(np.ceil((hi + pad - lo2) / voxel).astype(int) + 1, 2)
# 体素中心是否在凸包内 (Delaunay find_simplex, C 实现)
grids = np.meshgrid(*[lo2[i] + (np.arange(dims[i]) + 0.5) * voxel for i in range(3)], indexing="ij")
centers = np.stack([g.ravel() for g in grids], axis=1)
occ = (tri.find_simplex(centers) >= 0).reshape(dims) # 体素中心在凸包内
field = ndimage.gaussian_filter(occ.astype(np.float32), sigma=sigma)
v, f, _, _ = measure.marching_cubes(field, level=0.5)
v = (v + 0.5) * voxel + lo2 # 网格索引 -> 世界坐标
m = trimesh.Trimesh(vertices=v, faces=f[:, ::-1], process=True)
m = max(m.split(only_watertight=False), key=lambda x: x.volume)
filter_taubin(m, lamb=0.5, nu=0.5, iterations=5)
if faces_target and len(m.faces) > faces_target: # pymeshlab 简化 + 轻平滑
import pymeshlab as pml
ms = pml.MeshSet()
ms.add_mesh(pml.Mesh(vertex_matrix=m.vertices, face_matrix=m.faces.astype(np.int64)))
ms.meshing_decimation_quadric_edge_collapse(targetfacenum=int(faces_target),
preservenormal=True)
out = ms.current_mesh()
m = trimesh.Trimesh(vertices=np.array(out.vertex_matrix()),
faces=np.array(out.face_matrix()), process=False)
filter_taubin(m, lamb=0.5, nu=0.5, iterations=3)
return m, hull
# ------------------------------------------------------ 单实例流水线
def run_instance(verts, faces, idx, args):
"""处理一个眼镜实例 (一组顶点索引) 的完整流程, 返回该实例的壳 mesh (由调用方合并导出)。"""
assert idx.min() >= 0 and idx.max() < len(verts), \
"index 超出顶点范围, 请确认索引含义 (0/1-based, --one-based 可切换)"
region_mask = np.zeros(len(verts), bool)
region_mask[idx] = True # Blender 0-based 顶点索引
# ---- 噪声过滤: 只保留最大空间簇 (index 中混有 SAM3 多视角假阳性) ----
if args.filter == "main":
ridx = np.flatnonzero(region_mask)
mkeep, nkeep = largest_cluster(verts[ridx], args.cluster_r)
region_mask[ridx[~mkeep]] = False
print(f"[filter] 聚类(r={args.cluster_r}) 保留最大簇 {nkeep}/{len(ridx)} 顶点, "
f"剔除 {len(ridx) - nkeep} 个噪声点")
pts, inner = densify_points(verts, faces, region_mask)
print(f"[region] 顶点 {int(region_mask.sum())}, 内部面 {len(inner)}, "
f"采样点 {len(pts)}, bbox {np.round(pts.max(0) - pts.min(0), 3)}")
print("[delaunay] 四面体化中 ...")
st = AlphaStructure(pts)
print(f"[delaunay] 四面体 {st.T}, 采样密度决定的最小可行 alpha ~ "
f"{np.percentile(st.circum, 50):.3f}")
if args.mode == "case":
# ---- 眼镜盒模式: 最外侧点连成凸外壳 + 栅格化高斯隐式面 + MC -----
mesh, hull = build_case(pts, st.tri, args.res, args.sigma, args.case_faces)
kd = cKDTree(mesh.vertices)
dd, _ = kd.query(pts, k=1)
print(f"[case] 凸包体积 {hull.volume:.1f} -> 壳体积 {mesh.volume:.1f} "
f"({mesh.volume / hull.volume * 100:.1f}%)")
print(f"[shell] 顶点 {len(mesh.vertices)}, 面 {len(mesh.faces)}, "
f"水密 {mesh.is_watertight}, 欧拉数 {mesh.euler_number}, "
f"bbox {np.round(mesh.extents, 3)} (区域 bbox {np.round(pts.max(0) - pts.min(0), 3)})")
print(f"[fit] 采样点到壳最近距离: mean={dd.mean():.4f}, "
f"p99={np.percentile(dd, 99):.4f}, max={dd.max():.4f}")
return mesh
# ---- 最小 alpha 二分搜索 ------------------------------------------
if args.alpha is not None:
alpha = args.alpha
keep, ok_used, ok_conn = st.status(alpha)
print(f"[alpha] 手动 alpha={alpha:.4f}, 连通={ok_conn}, 覆盖={ok_used}")
else:
cand = np.unique(st.circum)
lo, hi = 0, len(cand) - 1
keep, ok_used, ok_conn = st.status(cand[hi])
assert ok_conn and ok_used, "最大 alpha 仍不连通, 数据异常"
while lo < hi:
mid = (lo + hi) // 2
k, u, c = st.status(cand[mid])
if u and c:
hi, keep, ok_used, ok_conn = mid, k, u, c
else:
lo = mid + 1
alpha = float(cand[lo])
print(f"[alpha] 最小连通 alpha = {alpha:.4f}")
# ---- 边界抽取 + 掐点拓扑修复 (仅局部补四面体, 不牺牲最小性) ---------
def try_extract(keep):
"""抽边界; 若有非 2 重边(掐点), 把该边周围的四面体扇补进复合体再试。"""
keep = keep.copy()
for _ in range(12):
bf = st.boundary_faces(keep)
m = trimesh.Trimesh(vertices=pts, faces=bf, process=False)
uniq, cnt = np.unique(m.edges_sorted, axis=0, return_counts=True)
bad = uniq[cnt != 2]
if len(bad) == 0:
return keep, m, m.is_watertight
added = 0
for a, b in bad:
fan = (st.tets == a).any(axis=1) & (st.tets == b).any(axis=1) & ~keep
keep |= fan
added += int(fan.sum())
if added == 0: # 无可补的四面体, 修复失败
return keep, m, False
return keep, m, m.is_watertight
keep, mesh, ok = try_extract(keep)
if not ok:
for _ in range(6): # 兜底: 逐步放大 alpha 重试
alpha = alpha * 1.1
k1, _, _ = st.status(alpha)
keep, mesh, ok = try_extract(k1)
print(f"[alpha] 局部修复不足, alpha 上调 10% -> {alpha:.4f}")
if ok:
break
if not mesh.is_watertight:
print("[warn] 仍未水密, 直接导出当前结果")
# ---- 外壳模式: 搜最小 alpha 使边界为单一 genus-0 闭合外壳 ----------
if args.mode == "outer" and mesh.euler_number != 2:
def pred(a):
"""alpha=a 时边界(修复掐点后)是否已是单一 genus-0 外壳。"""
k, _, _ = st.status(a)
_, m, _ = try_extract(k)
parts = m.split(only_watertight=False)
outer = [p for p in parts if p.volume > 0]
if len(outer) != 1:
return False, None
m = outer[0]
if m.is_watertight and m.euler_number == 2:
return True, m
return False, None
a_lo, a_hi, found = alpha, None, None # a_lo: 当前非外壳
a = alpha
for _ in range(20): # 指数上扫找可行 alpha
a = a * 1.3
okf, mf = pred(a)
if okf:
a_hi, found = a, mf
break
if a_hi is None:
print("[warn] 未找到 genus-0 外壳, 保持紧凑壳")
else:
for _ in range(8): # 二分细化到最小外壳 alpha
mid = (a_lo + a_hi) / 2
okf, mf = pred(mid)
if okf:
a_hi, found = mid, mf
else:
a_lo = mid
alpha, mesh = a_hi, found
print(f"[alpha] 外壳模式: 最小 genus-0 外壳 alpha = {alpha:.4f}")
vol0 = mesh.volume
# 去掉内腔壳 (负体积分量), 只留外壳
parts = mesh.split(only_watertight=False)
outer = [p for p in parts if p.volume > 0]
if len(outer) != len(parts):
print(f"[shell] 去除 {len(parts) - len(outer)} 个内腔面片")
mesh = max(outer, key=lambda m: m.volume) if len(outer) > 1 else outer[0]
# ---- Taubin 平滑 (不收缩) ------------------------------------------
if args.smooth > 0:
before = mesh.vertices.copy()
filter_taubin(mesh, lamb=0.5, nu=0.5, iterations=args.smooth)
drift = np.linalg.norm(mesh.vertices - before, axis=1)
print(f"[smooth] Taubin x{args.smooth}, 顶点漂移 mean={drift.mean():.4f} max={drift.max():.4f}")
# ---- 数值校验 ------------------------------------------------------
prox = trimesh.proximity.ProximityQuery(mesh)
d = np.abs(prox.signed_distance(pts))
print(f"[shell] 顶点 {len(mesh.vertices)}, 面 {len(mesh.faces)}, "
f"水密 {mesh.is_watertight}, 欧拉数 {mesh.euler_number}")
print(f"[shell] 体积 {mesh.volume:.1f} (平滑前 {vol0:.1f}), "
f"bbox {np.round(mesh.extents, 3)}")
print(f"[fit] 采样点到壳距离 |d|: mean={d.mean():.4f}, "
f"p99={np.percentile(d, 99):.4f}, max={d.max():.4f}")
return mesh
# ------------------------------------------------------ 主流程
def main():
ap = argparse.ArgumentParser(description="生成每个眼镜实例的最小壳 (alpha shape)")
here = os.path.dirname(os.path.abspath(__file__))
ap.add_argument("--obj", default=None,
help="obj 路径; 缺省依次取 data.json 的 input 字段 / data 目录内唯一 obj")
ap.add_argument("--data", default=os.path.join(here, "data.json"),
help="多实例 data.json (instances[].indices), 每实例单独处理")
ap.add_argument("--index", default=None,
help="兼容旧的单实例 index.txt 顶点索引文件; 指定时优先于 --data")
ap.add_argument("--out", default=os.path.join(here, "shell.obj"),
help="输出路径, 所有实例的壳合并导出为这一个文件 (默认 shell.obj)")
ap.add_argument("--mode", choices=["case", "outer", "tight"], default="case",
help="case=眼镜盒式光滑凸外壳(最外侧连成一体), outer=跟随外形的 genus-0 外壳, "
"tight=贴表面的紧身壳 (默认 case)")
ap.add_argument("--res", type=int, default=128, help="case 模式体素分辨率 (默认 128)")
ap.add_argument("--sigma", type=float, default=2.5, help="case 模式高斯平滑 sigma, 体素单位 (默认 2.5)")
ap.add_argument("--case-faces", type=int, default=25000, help="case 模式简化目标面数, 0=不简化 (默认 25000)")
ap.add_argument("--alpha", type=float, default=None,
help="outer/tight 模式手动指定 alpha; 缺省自动搜索")
ap.add_argument("--smooth", type=int, default=10,
help="Taubin 平滑迭代数, 0=关闭取严格最小壳 (默认 10)")
ap.add_argument("--filter", choices=["main", "all"], default="main",
help="main=只保留最大空间簇(去 SAM3 噪声), all=用全部 index 点 (默认 main)")
ap.add_argument("--cluster-r", type=float, default=2.5,
help="噪声过滤的聚类半径, 世界单位 (默认 2.5)")
ap.add_argument("--one-based", action="store_true",
help="索引为 1-based (默认按 Blender 导出的 0-based 读取)")
args = ap.parse_args()
instances, input_obj = load_data(args)
obj_path = resolve_obj(args, input_obj)
verts, faces = load_obj(obj_path)
print(f"[load] {obj_path}: 顶点 {len(verts)}, 面 {len(faces)}, "
f"待处理实例 {len(instances)} 个")
meshes, n_ok = [], 0
for name, idx in instances:
print(f"\n===== 实例 {name}: {len(idx)} 顶点 =====")
try:
meshes.append(run_instance(verts, faces, idx, args))
n_ok += 1
except Exception as e:
print(f"[error] 实例 {name} 处理失败, 跳过: {e}")
if not meshes:
raise SystemExit("[error] 所有实例处理失败, 未导出任何壳")
# ---- 合并所有实例的壳 (同一坐标系, 直接拼接), 导出单个 obj ----
combined = meshes[0] if len(meshes) == 1 else trimesh.util.concatenate(meshes)
combined.export(args.out)
print(f"\n[summary] 完成 {n_ok}/{len(instances)} 个实例, 合并导出: {args.out} "
f"(顶点 {len(combined.vertices)}, 面 {len(combined.faces)}, "
f"连通分量 {len(combined.split())})")
if __name__ == "__main__":
main()