diff --git a/fill2sc/make_shell.py b/fill2sc/make_shell.py new file mode 100644 index 0000000..f67b138 --- /dev/null +++ b/fill2sc/make_shell.py @@ -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()