diff --git a/fill2sc/generate_fill_mesh.py b/fill2sc/generate_fill_mesh.py index 1fa551d..61d6bd9 100644 --- a/fill2sc/generate_fill_mesh.py +++ b/fill2sc/generate_fill_mesh.py @@ -34,34 +34,47 @@ generate_fill_mesh.py 仅依赖 numpy 与 PIL。 """ +import sys import numpy as np from collections import deque from PIL import Image, ImageDraw # ----------------------------- 参数 ----------------------------- -IN_OBJ = "632723.obj" -IN_IDX = "index.txt" +# 用法: python3 generate_fill_mesh.py [模型.obj] [眼镜顶点索引.txt] +# (缺省用下方默认值; 输出 fill_mesh.obj / <模型名>_filled.obj / +# fill_preview.png / fill_report.txt 写在当前目录) +IN_OBJ = sys.argv[1] if len(sys.argv) > 1 else "632723.obj" +IN_IDX = sys.argv[2] if len(sys.argv) > 2 else "index.txt" OUT_FILL = "fill_mesh.obj" -OUT_MERGED = "632723_filled.obj" +OUT_MERGED = (IN_OBJ[:-4] if IN_OBJ.endswith(".obj") else IN_OBJ) + "_filled.obj" OUT_PREVIEW = "fill_preview.png" RES = 0.1 # 栅格分辨率 -X0, X1 = -10.2, 10.2 -Y0, Y1 = 9.2, 19.8 -RY0, RY1 = 9.7, 19.3 # 填充区域 y 裁剪(只处理正面眼镜部分) +# 栅格窗口不写死(旧值 x[-10.2,10.2] y[9.2,19.8] 只适配 632723 的头部位置, +# 换模型会切掉镜圈 -> 镜片开口连通窗口边界 -> 检测不到开口 -> 镜片镂空不填): +# 由眼镜主分量 XY bbox 自动外扩推导, 见 face_g 提取之后 +WMARG = 2.0 # 窗口边距(眼镜 bbox 外扩, 容纳足迹膨胀/边界平滑/背景) +RY_PAD = 0.5 # 填充区域 y 裁剪内缩(避开窗口边缘栅格化截断) SEP = 0.25 # 区分前后表面的最小 z 间距 -FLUSH = 0.15 # 填充前表面距镜框前面的距离(与镜框平齐, 留嵌入余量) +FLUSH = 0.05 # 前表面低于本地框面的余量(≈平齐; 局部平滑抬高=轻微凸出) +BULGE = 0.45 # 镜片弧面凸出幅值(抛物面, 中心最高, 边缘归零平滑衔接) EMBED_BOT = 0.6 # 填充底面嵌入脸部深度 -SHELL_MAX = 3.0 # 无真实脸面处的填充厚度上限(相对前表面) +SHELL_MAX = 7.0 # 无真实脸面处厚度上限(透明支撑体, 尾部厚属预期) FACE_MIN = 2.0 # 脸面低于此值视为误标伪面(后脑等), 不作填充底面 CLAMP_BACK = 0.4 # 距最后表面的安全距离(防穿出后脑) GAP_MIN = 0.25 # 填充最小厚度 SIGMA = 2.0 # 高度场平滑 sigma(格) -ERODE = 2 # 足迹腐蚀(格) 0.2, 与膨胀0.3配合净外扩0.1(补偿边界平滑收缩) +ERODE = 3 # 足迹腐蚀(格) 0.3, 与膨胀0.3配合净外扩0(XY严格不超镜框轮廓) BSTEP = 0.18 # 边界重采样步长 BITER = 25 # 边界 Laplacian 平滑迭代 BLAM = 0.6 # 边界平滑系数 TOL = 1e-4 # Jacobi 收敛阈值 +GEO_MIN_COMP = 500 # 标记边连通主分量面数低于此值 -> 判定 index 无效, 切几何模式 +FUSE_BLOB = 25 # 贴合脸面的标记格最大 8 连通团块格数 >= 此值 -> 疑似融合 +FUSE_SKIN_F = 2500 # 且主分量中贴合皮肤面数 >= 此值(或占比>=10%) -> 判定真融合, 切几何模式 + # (眼镜搭鼻上的轻接触 ~800 面/4-5% 无害: 该处填充前面=脸面高度, + # 标记模式天然兼容; 只有分割把大块皮肤涂进标记才需要几何兜底) +GRES = 0.2 # 几何预扫描分辨率(仅用于定位眼镜, 主栅格仍为 RES) report = [] def log(s): @@ -85,11 +98,32 @@ def parse_obj(path): log("== 1. 读入模型 ==") V, F = parse_obj(IN_OBJ) -idx = np.loadtxt(IN_IDX, dtype=np.int64) - 1 +_idx_raw = np.loadtxt(IN_IDX, dtype=np.int64) +# index 基址自适应: 不同批次数据的 index 有 0 基(直接顶点编号)与 1 基 +# (OBJ 行号)两种约定. 判据 = "3 顶点全标记面"数: 正确基址下标记落在 +# 眼镜表面的连续顶点块上, 全标记面多; 错位基址把标记挪到邻近无关顶点 +# (顶点编号与空间不相关时彻底散开, 如 20261102), 全标记面骤减. +def _full_marked(_ids): + _ids = _ids[(_ids >= 0) & (_ids < len(V))] + if len(_ids) == 0: + return 0 + _m = np.zeros(len(V), bool) + _m[_ids] = True + return int((_m[F].sum(1) == 3).sum()) +_n1 = _full_marked(_idx_raw - 1) +_n0 = _full_marked(_idx_raw) +if _n1 > _n0: + idx = _idx_raw - 1 + log(f" 顶点 {len(V)}, 面 {len(F)}, 眼镜标记顶点 {len(_idx_raw)} " + f"(判 1 基: 全标记面 {_n1} > {_n0})") +else: + idx = _idx_raw + log(f" 顶点 {len(V)}, 面 {len(F)}, 眼镜标记顶点 {len(_idx_raw)} " + f"(判 0 基: 全标记面 {_n0} >= {_n1})") +idx = idx[(idx >= 0) & (idx < len(V))] ism = np.zeros(len(V), bool) ism[idx] = True face_g = ism[F].any(1) -log(f" 顶点 {len(V)}, 面 {len(F)}, 眼镜标记顶点 {len(idx)}") # 真眼镜物体提取: 标记面的"边连通"主分量. 眼镜与头部网格互不共边, # 误标皮肤(鼻梁/额头/后脑, 分割噪声)在边连通下是独立小分量 -> 一次性剥离; @@ -128,6 +162,265 @@ face_g[_best] = True log(f" 眼镜物体 = 边连通主分量 {len(_best)} 面, " f"剥离误标皮肤 {int((ism[F].any(1)).sum() - len(_best))} 面") +# ------------- 粗栅格管对扫描(几何兜底 + 融合检测共用) ------------- +# 眼镜 = 悬浮在脸面前的"管状"结构. 沿 z 射线在 (x,y) 粗格内出现一对 +# 相邻 hit(镜框前后两面, 间距 0.1~2.0), 且这对身后还有表面(落后>=0.15, +# 即脸) -> 管格; 全模型最前方的管格 8 连通块 = 镜圈/镜桥正面, 再沿 +# 管带 8 连通生长带上镜腿/尾部折返段(它们自身也是管对, 身后是头部). +# 该扫描无条件执行, 服务两种用途: +# a) 几何兜底: index 无效(标记主分量 < GEO_MIN_COMP)或标记与脸融合 +# 时, 眼镜面片 = 质心落在眼镜区域且 z 落在邻域(3x3 管格)管对 +# 极值窗口 [min(zb)-0.1, max(zf)+0.15] 内: 检测间隙 0.15 > 窗口 +# 余量 0.1, 结构上保证管身后的皮肤/脸面进不了窗口 -> 无接触泄漏 +# (无需洪泛, 眼镜与头部共边也不受影响); 镜桥 pinch 格邻域无管格 +# -> 窗口空, 不选面(由调和膜插值). 后续管线不变. +# b) 融合检测: 有些数据(如 20261102)的眼镜与脸部部分连在一起, 分割 +# 把接触区皮肤整块标了进来 -> 主分量 = 眼镜+融合皮肤, 直接用会把 +# 填充泄到脸上. 皮肤面就贴在脸上: 主分量覆盖格的 5x5 内非主分量 +# 面最高 z(脸面参考, 排除比覆盖面低 4.0+ 的后脑面)与覆盖面最高 z +# 落差 <=0.15 = 贴合格; 纯模型只有眉梁/镜桥贴触细条(几个格), +# 融合模型的皮肤是大块贴合团 -> 最大 8 连通贴合格团块 >= FUSE_BLOB +# -> 判定融合, 切几何模式. +# 1) 全模型粗栅格, 每格收集全部 z hit +# 全身模型(躯干/手臂/手)会在身体上产生大量伪管对, 经闭运算污染 +# 眼镜区域 -> 有有效标记主分量时把扫描限定在标记 bbox 外扩 3.0 +# (覆盖标记漏标的框上缘 ~1.2 与镜尾折返段; 标记无效仍全模型扫描) +_gx0 = float(np.floor(V[:, 0].min() / GRES)) * GRES +_gx1 = float(np.ceil(V[:, 0].max() / GRES)) * GRES +_gy0 = float(np.floor(V[:, 1].min() / GRES)) * GRES +_gy1 = float(np.ceil(V[:, 1].max() / GRES)) * GRES +if len(_best) >= GEO_MIN_COMP: + _bv = V[np.unique(F[_best])] + _gx0 = max(_gx0, float(np.floor((_bv[:, 0].min() - 3.0) / GRES)) * GRES) + _gx1 = min(_gx1, float(np.ceil((_bv[:, 0].max() + 3.0) / GRES)) * GRES) + _gy0 = max(_gy0, float(np.floor((_bv[:, 1].min() - 3.0) / GRES)) * GRES) + _gy1 = min(_gy1, float(np.ceil((_bv[:, 1].max() + 3.0) / GRES)) * GRES) + log(f" 管格扫描限定于标记 bbox 外扩 3.0: " + f"x[{_gx0:.1f},{_gx1:.1f}] y[{_gy0:.1f},{_gy1:.1f}]") +_gw = int(round((_gx1 - _gx0) / GRES)) + 1 +_gh = int(round((_gy1 - _gy0) / GRES)) + 1 +_ghit = [[[] for _ in range(_gw)] for _ in range(_gh)] +for _t in range(len(F)): + _a, _b, _c = F[_t] + _ax, _ay, _az = V[_a] + _bx, _by, _bz2 = V[_b] + _cx, _cy, _cz2 = V[_c] + _det = (_bx - _ax) * (_cy - _ay) - (_by - _ay) * (_cx - _ax) + if _det == 0: + continue + _i0 = int((min(_ax, _bx, _cx) - _gx0) / GRES) + _i1 = int((max(_ax, _bx, _cx) - _gx0) / GRES) + _j0 = int((min(_ay, _by, _cy) - _gy0) / GRES) + _j1 = int((max(_ay, _by, _cy) - _gy0) / GRES) + if _i1 < 0 or _j1 < 0 or _i0 >= _gw or _j0 >= _gh: + continue # 限定扫描范围外的面 + _i0 = max(_i0, 0) + _j0 = max(_j0, 0) + _i1 = min(_i1, _gw - 1) + _j1 = min(_j1, _gh - 1) + for _j in range(_j0, _j1 + 1): + _py = _gy0 + (_j + 0.5) * GRES + for _i in range(_i0, _i1 + 1): + _px = _gx0 + (_i + 0.5) * GRES + _u = ((_px - _ax) * (_cy - _ay) + - (_py - _ay) * (_cx - _ax)) / _det + if _u < -1e-9: + continue + _v = ((_bx - _ax) * (_py - _ay) + - (_by - _ay) * (_px - _ax)) / _det + if _v < -1e-9 or _u + _v > 1 + 1e-9: + continue + _ghit[_j][_i].append( + _az + _u * (_bz2 - _az) + _v * (_cz2 - _az)) +# 2) 管对检测: 从最前一对向后找第一个满足 +# (对间距 0.1~2.0 且身后还有表面落后>=0.15) 的悬浮对 +_tube = np.zeros((_gh, _gw), bool) +_gzf = np.full((_gh, _gw), -np.inf) # 管前面 z +_gzb = np.full((_gh, _gw), -np.inf) # 管背面 z +for _j in range(_gh): + for _i in range(_gw): + _zs = sorted(_ghit[_j][_i]) + if len(_zs) < 3: + continue + for _k in range(len(_zs) - 2, 0, -1): + _dz = _zs[_k + 1] - _zs[_k] + if 0.1 <= _dz <= 2.0 and _zs[_k - 1] <= _zs[_k] - 0.15: + _tube[_j, _i] = True + _gzf[_j, _i] = _zs[_k + 1] + _gzb[_j, _i] = _zs[_k] + break +# 3) 管带闭运算连通: 镜桥压在鼻上无悬浮间隙 -> 无管格, 两镜圈管带 +# 相距 ~1.4; 尾部折返段与镜圈相距 ~0.5-1.0. 圆盘闭运算(半径 1.0, +# 填 <=2.0 的 pinch, 不外扩外轮廓)把镜圈/镜桥带/尾部并成一个 +# 连通域; 头发/手等杂散管 >=3.2 不会并入. 取含全局最前管格的 +# 连通域 = 眼镜区域(管格 + 镜桥 pinch 填充) +_offs = [(_dj, _di) for _dj in range(-5, 6) for _di in range(-5, 6) + if _dj * _dj + _di * _di <= 25] +_td = np.zeros_like(_tube) +for _dj, _di in _offs: + _td |= np.roll(_tube, (_dj, _di), (0, 1)) +_tc = np.ones_like(_tube) +for _dj, _di in _offs: + _tc &= np.roll(_td, (_dj, _di), (0, 1)) +# 有效标记时: 种子取标记足迹(外扩 2.0)内最前的管格(全身模型上头发/ +# 耳朵伪管可能比镜框更前), 且区域生长完成后与标记足迹外扩 2.0 相交, +# 排除远处伪管(头发距标记中位 2.8+)并入 -> 窗口选面污染 +_mfootd = None +if len(_best) >= GEO_MIN_COMP: + _mv = V[np.unique(F[_best])] + _mf = np.zeros((_gh, _gw), bool) + _mi = np.clip(((_mv[:, 0] - _gx0) / GRES).astype(int), 0, _gw - 1) + _mj = np.clip(((_mv[:, 1] - _gy0) / GRES).astype(int), 0, _gh - 1) + _mf[_mj, _mi] = True + _mfootd = _mf.copy() + for _ in range(10): # 外扩 2.0 单位 + _g = _mfootd.copy() + for _dj, _di in ((1, 0), (-1, 0), (0, 1), (0, -1)): + _g |= np.roll(_mfootd, (_dj, _di), (0, 1)) + _mfootd = _g + _tt = _tube & _mfootd +else: + _tt = _tube +if _tt.any(): + _jj, _ii = np.unravel_index( + np.argmax(np.where(_tt, _gzf, -np.inf)), _gzf.shape) +else: + _jj, _ii = np.unravel_index( + np.argmax(np.where(_tube, _gzf, -np.inf)), _gzf.shape) +_region = np.zeros_like(_tube) +_region[_jj, _ii] = True +_stk = [(_jj, _ii)] +while _stk: + _y, _x = _stk.pop() + for _dy in (-1, 0, 1): + for _dx in (-1, 0, 1): + _ny, _nx = _y + _dy, _x + _dx + if (0 <= _ny < _gh and 0 <= _nx < _gw + and _tc[_ny, _nx] and not _region[_ny, _nx]): + _region[_ny, _nx] = True + _stk.append((_ny, _nx)) +if _mfootd is not None: + _region &= _mfootd +_ry, _rx = np.where(_region) +log(f" 管格 {int(_tube.sum())}, 眼镜区域(闭运算) {int(_region.sum())} 格 " + f"x[{_gx0 + _rx.min() * GRES:.1f},{_gx0 + _rx.max() * GRES:.1f}] " + f"y[{_gy0 + _ry.min() * GRES:.1f},{_gy0 + _ry.max() * GRES:.1f}], " + f"管前面 z 顶 {_gzf[_jj, _ii]:.2f}") +# ---- 融合检测: 标记主分量是否混入贴合脸面的皮肤 ---- +# 皮肤面就贴在脸上 -> 与周围脸面(非主分量面)同高(落差 |dz|<=0.15); +# 眼镜面悬浮在脸前(正落差), 镜尾折返面低于头侧表面(大负落差) -> +# 都不是贴合. 贴合格聚成大 8 连通团块 = 融合皮肤(纯模型只有眉梁/ +# 镜桥/轻度皮肤 fringe 的贴触细条, 十几个格). +_fused = False +if len(_best) >= GEO_MIN_COMP and _tube.any(): + _fcs = V[F[_best]].mean(1) + _fci = np.clip(((_fcs[:, 0] - _gx0) / GRES).astype(int), 0, _gw - 1) + _fcj = np.clip(((_fcs[:, 1] - _gy0) / GRES).astype(int), 0, _gh - 1) + _fcc = np.zeros((_gh, _gw), bool) + _fcc[_fcj, _fci] = True + _fcz = np.full((_gh, _gw), -np.inf) + np.maximum.at(_fcz, (_fcj, _fci), _fcs[:, 2]) + # 脸面参考: 非主分量面质心的每格最高 z, 逐格排除比覆盖面低 4.0+ + # 的后脑面, 7x7 邻域传播(皮肤块深处的参考取自块外脸面的延续) + _allc = V[F].mean(1) + _aci = np.clip(((_allc[:, 0] - _gx0) / GRES).astype(int), 0, _gw - 1) + _acj = np.clip(((_allc[:, 1] - _gy0) / GRES).astype(int), 0, _gh - 1) + _nc = ~face_g + _nz = np.full((_gh, _gw), -np.inf) + np.maximum.at(_nz, (_acj[_nc], _aci[_nc]), _allc[_nc, 2]) + _nv = np.where(_nz > _fcz - 4.0, _nz, -np.inf) + _nzp = _nv.copy() + for _dj in range(-3, 4): + for _di in range(-3, 4): + _nzp = np.maximum(_nzp, np.roll(_nv, (_dj, _di), (0, 1))) + _nzq = np.where(_nzp > -np.inf, _nzp, -1e9) + _att = _fcc & (np.abs(_fcz - _nzq) <= 0.15) + _skinf = int(_att[_fcj, _fci].sum()) + _avis = set(map(tuple, np.argwhere(_att))) + _blobmax = 0 + while _avis: + _seed = _avis.pop() + _stk2, _sz = [_seed], 0 + while _stk2: + _cj2, _ci2 = _stk2.pop() + _sz += 1 + for _dy in (-1, 0, 1): + for _dx in (-1, 0, 1): + _nb = (_cj2 + _dy, _ci2 + _dx) + if _nb in _avis: + _avis.remove(_nb) + _stk2.append(_nb) + _blobmax = max(_blobmax, _sz) + # 融合 = 大块贴合团 且 皮肤面量大(绝对或占比). 只有团块时可能只是 + # 鼻托轻接触(如 617892: 团 26/皮肤 755 面 4.5%), 标记模式即可; + # 真融合(20261102: 团 38/皮肤 5000+ 面 17.8%)才切几何模式. + _fused = (_blobmax >= FUSE_BLOB + and (_skinf >= FUSE_SKIN_F or _skinf >= 0.10 * len(_best))) + log(f" 融合检测: 主分量覆盖 {int(_fcc.sum())} 格, 贴合格 " + f"{int(_att.sum())}, 最大贴合团块 {_blobmax} 格 (阈值 {FUSE_BLOB}), " + f"贴合皮肤面 {_skinf}/{len(_best)}" + f" -> {'融合, 切几何模式' if _fused else '标记纯净(轻接触无害), 走标记模式'}") + +# ---- 路由: index 无效 / 标记与脸融合 -> 几何模式; 否则标记模式 ---- +if len(_best) < GEO_MIN_COMP: + log(f"! 标记主分量仅 {len(_best)} 面 (<{GEO_MIN_COMP}), index 无效" + f" -> 切换几何检测模式") +elif _fused: + log("! 标记与脸融合(大块贴合皮肤团) -> 切换几何检测模式") +if len(_best) < GEO_MIN_COMP or _fused: + if not _tube.any(): + raise SystemExit("几何模式: 未找到任何悬浮管对, 无法定位眼镜") + # 4) 眼镜面片 = 质心落在眼镜区域内, 且 z 落在邻域(3x3 管格)管对 + # 极值窗口 [min(zb)-0.1, max(zf)+0.15] 内. + # 单格 ±0.1 窗口会漏选大量管面(陡斜管壁上质心与格中心 z 错位>0.1) + # -> face_g 覆盖稀疏 -> 环角/尾部出现大片调和膜口袋, 前表面 + # 斜降到头后穿出. 邻域极值跟随管壁倾斜; 皮肤至少低于每个管背 + # 0.15(检测间隙), 也低于邻域 min(zb), -0.1 余量仍保证不吸皮肤. + # 镜桥 pinch 格邻域无管格 -> 窗口空, 不选面(由调和膜插值) + _gzb3 = np.full((_gh, _gw), np.inf) + _gzf3 = np.full((_gh, _gw), -np.inf) + for _dj in (-1, 0, 1): + for _di in (-1, 0, 1): + _gzb3 = np.minimum( + _gzb3, np.roll(np.where(_tube, _gzb, np.inf), (_dj, _di), (0, 1))) + _gzf3 = np.maximum( + _gzf3, np.roll(np.where(_tube, _gzf, -np.inf), (_dj, _di), (0, 1))) + _gzlo = np.where(_region, _gzb3 - 0.1, np.inf) + _gzhi = np.where(_region, _gzf3 + 0.15, -np.inf) + _gfc = V[F].mean(1) + _gci = np.clip(((_gfc[:, 0] - _gx0) / GRES).astype(int), 0, _gw - 1) + _gcj = np.clip(((_gfc[:, 1] - _gy0) / GRES).astype(int), 0, _gh - 1) + _gcz = _gfc[:, 2] + _gsel = (_gcz >= _gzlo[_gcj, _gci]) & (_gcz <= _gzhi[_gcj, _gci]) + _gfaces = np.where(_gsel)[0] + if not 300 <= len(_gfaces) <= 60000: + raise SystemExit( + f"几何模式: 管对窗口选出 {len(_gfaces)} 面, 超出合理范围" + f" [300, 60000], 定位失败") + face_g = np.zeros(len(F), bool) + face_g[_gfaces] = True + _gvz = V[np.unique(F[_gfaces])] + log(f" 几何模式选定眼镜面片: {len(_gfaces)} 面, " + f"bbox x[{_gvz[:,0].min():.1f},{_gvz[:,0].max():.1f}] " + f"y[{_gvz[:,1].min():.1f},{_gvz[:,1].max():.1f}] " + f"z[{_gvz[:,2].min():.1f},{_gvz[:,2].max():.1f}]") + _GEO_MODE = True + _GEO_FOOT = _region +else: + _GEO_MODE = False + +# 栅格窗口自动推导: 眼镜主分量顶点 XY bbox 外扩 WMARG(对齐 RES), +# 任何模型/任何头部位置通用 +_gv = V[np.unique(F[face_g])] +X0 = float(np.floor((_gv[:, 0].min() - WMARG) / RES)) * RES +X1 = float(np.ceil((_gv[:, 0].max() + WMARG) / RES)) * RES +Y0 = float(np.floor((_gv[:, 1].min() - WMARG) / RES)) * RES +Y1 = float(np.ceil((_gv[:, 1].max() + WMARG) / RES)) * RES +RY0, RY1 = Y0 + RY_PAD, Y1 - RY_PAD # 填充区域 y 裁剪(只处理正面眼镜部分) +log(f" 窗口自动推导: x[{X0:.1f},{X1:.1f}] y[{Y0:.1f},{Y1:.1f}] " + f"(眼镜 bbox x[{_gv[:,0].min():.1f},{_gv[:,0].max():.1f}] " + f"y[{_gv[:,1].min():.1f},{_gv[:,1].max():.1f}] + {WMARG} 边距)") + # ----------------------------- 栅格 ----------------------------- W = int(round((X1 - X0) / RES)) + 1 H = int(round((Y1 - Y0) / RES)) + 1 @@ -208,7 +501,20 @@ log("== 3. 眼镜足迹 / 填充区域 R ==") # 足迹 = 真眼镜物体(边连通主分量)覆盖的全部格子(含尾部包裹段). # 误标杂散已在 face_g 阶段剥离, 无需 z 高度过滤; 旧 gz1>=7 过滤会剪掉 # 尾部折返带(z 低于 7) -> 尾部镂空填不上 -mask = g1f +if _GEO_MODE: + # 几何模式: 足迹 = 眼镜区域(管带闭运算, 已含镜桥 pinch 填充)上采样. + # 桥位格无眼镜面片 -> 无 gz1, 与镜片开口一样由调和膜插值前表面 + mask = np.zeros((H, W), bool) + for _j, _i in zip(*np.where(_GEO_FOOT)): + _mj = int(np.floor(((_gy0 + _j * GRES) - Y0) / RES + 1e-9)) + _mi = int(np.floor(((_gx0 + _i * GRES) - X0) / RES + 1e-9)) + for _dj in (0, 1): + for _di in (0, 1): + _yy, _xx = _mj + _dj, _mi + _di + if 0 <= _yy < H and 0 <= _xx < W: + mask[_yy, _xx] = True +else: + mask = g1f _labs = np.zeros(mask.shape, np.int32) _cur = 0 from collections import deque as _dq @@ -256,25 +562,51 @@ def _flood_out(msk): return out _holes = ~mask & ~_flood_out(mask) -_hl = np.zeros(mask.shape, np.int32) -_big = [] -for _j in range(mask.shape[0]): - for _i in range(mask.shape[1]): - if _holes[_j, _i] and _hl[_j, _i] == 0: - _hl[_j, _i] = 1 - _q = _dq([(_j, _i)]) - _pts = [(_j, _i)] - while _q: - _y, _x = _q.popleft() - for _dj, _di in ((1,0),(-1,0),(0,1),(0,-1),(1,1),(1,-1),(-1,1),(-1,-1)): - _yy, _xx = _y + _dj, _x + _di - if (0 <= _yy < mask.shape[0] and 0 <= _xx < mask.shape[1] - and _holes[_yy, _xx] and _hl[_yy, _xx] == 0): - _hl[_yy, _xx] = 1 - _q.append((_yy, _xx)) - _pts.append((_yy, _xx)) - if len(_pts) >= 500: - _big.append(_pts) + +def _label_comps(msk): + """8 连通分量标注, 返回所有分量的格点列表(按扫描顺序)""" + _lab = np.zeros(msk.shape, np.int32) + _comps = [] + for _j in range(msk.shape[0]): + for _i in range(msk.shape[1]): + if msk[_j, _i] and _lab[_j, _i] == 0: + _lab[_j, _i] = 1 + _q = _dq([(_j, _i)]) + _pts = [(_j, _i)] + while _q: + _y, _x = _q.popleft() + for _dj, _di in ((1,0),(-1,0),(0,1),(0,-1),(1,1),(1,-1),(-1,1),(-1,-1)): + _yy, _xx = _y + _dj, _x + _di + if (0 <= _yy < msk.shape[0] and 0 <= _xx < msk.shape[1] + and msk[_yy, _xx] and _lab[_yy, _xx] == 0): + _lab[_yy, _xx] = 1 + _q.append((_yy, _xx)) + _pts.append((_yy, _xx)) + _comps.append(_pts) + return _comps + +_comps = _label_comps(_holes) +_hsz = [len(_p) for _p in _comps] +_big = [_p for _p in _comps if len(_p) >= 250] +# 镜片开口掩码(供后面镜片光滑度校验区域使用, 不再写死坐标): +# 取两大开口洞并腐蚀 2 格(0.2) -> 只留镜片膜主体, 自动排除镜圈/ +# 镜腿折叠处的陡降与狭缝袋(洞 bbox 会经狭缝袋延伸到镜圈外侧) +lens_boxes = [] +lens_mask = np.zeros(mask.shape, bool) +for _pts in _big[:2]: + _js = [p[0] for p in _pts]; _is = [p[1] for p in _pts] + lens_boxes.append((X0 + min(_is) * RES, X0 + max(_is) * RES, + Y0 + min(_js) * RES, Y0 + max(_js) * RES)) + for p in _pts: + lens_mask[p[0], p[1]] = True +for _ in range(2): + _nm = lens_mask.copy() + _nm[1:, :] &= lens_mask[:-1, :]; _nm[:-1, :] &= lens_mask[1:, :] + _nm[:, 1:] &= lens_mask[:, :-1]; _nm[:, :-1] &= lens_mask[:, 1:] + lens_mask = _nm +if lens_boxes: + log(" 镜片开口框: " + "; ".join( + "x[%.1f,%.1f] y[%.1f,%.1f]" % b for b in lens_boxes)) if len(_big) >= 2: # 按 x 中心排序, 区分左右开口 _big.sort(key=len, reverse=True) @@ -294,11 +626,11 @@ if len(_big) >= 2: log(f" 中央区剪除(两镜片开口之间, 开口底缘 y>{Y0 + y_cut * RES:.1f} 以下): " f"剔除 {_ncut} 格误标鼻面") else: - log(" ! 未找到两个镜片开口, 跳过中央区剪除") + log(" ! 未找到两个镜片开口, 跳过中央区剪除" + + (f"; 洞分量尺寸: {sorted(_hsz, reverse=True)[:8]}" if _hsz else "; 无内部洞(开口与外部连通)")) # 膨胀 3 格(~0.3) 闭合采样缝隙, 防止镜片开口与外部连通; -# 与腐蚀 2 格(0.2) 配合净外扩 0.1: 边界平滑会压扁凸部(尾部舌形)~0.2, -# 预留外扩让平滑后的边界落在足迹外缘, 镜框下镂空填到边(真外轮廓 -# 比栅格覆盖宽~0.1-0.2, 不会超出镜框) +# 与腐蚀 3 格(0.3) 配合净外扩 0: 边界平滑后落在栅格足迹内, +# XY 严格不超镜框轮廓(真外轮廓比栅格覆盖宽~0.1-0.2) for _ in range(3): nm = mask.copy() nm[1:, :] |= mask[:-1, :]; nm[:-1, :] |= mask[1:, :] @@ -328,12 +660,66 @@ def flood_outside(mask): outside[nj, ni] = True; dq.append((nj, ni)) return ~outside & ~mask -holes = flood_outside(mask) +# 几何模式: 镜圈内上角管对检测失败的楔形缺口把镜片开口与外部连通, +# flood 找不到开口 -> 开口镂空不填. 对已膨胀掩码再做 0.5 闭运算(仅用于 +# 洞检测, 掩码本身不动, 外轮廓逻辑不变): 封闭 <=1.0 的缺口, 两镜片开口 +# 成为内部洞 -> 正常补上(开口边界仅角部圆化 <=0.5). marks 模式掩码为 +# 密集覆盖, 无此问题 +if _GEO_MODE: + _clm = mask.copy() + for _ in range(5): + _nm = _clm.copy() + _nm[1:, :] |= _clm[:-1, :]; _nm[:-1, :] |= _clm[1:, :] + _nm[:, 1:] |= _clm[:, :-1]; _nm[:, :-1] |= _clm[:, 1:] + _nm[1:, 1:] |= _clm[:-1, :-1]; _nm[:-1, :-1] |= _clm[1:, 1:] + _nm[1:, :-1] |= _clm[:-1, 1:]; _nm[:-1, 1:] |= _clm[1:, :-1] + _clm = _nm + for _ in range(5): + _nm = _clm.copy() + _nm[1:, :] &= _clm[:-1, :]; _nm[:-1, :] &= _clm[1:, :] + _nm[:, 1:] &= _clm[:, :-1]; _nm[:, :-1] &= _clm[:, 1:] + _nm[1:, 1:] &= _clm[:-1, :-1]; _nm[:-1, :-1] &= _clm[1:, 1:] + _nm[1:, :-1] &= _clm[:-1, 1:]; _nm[:-1, 1:] &= _clm[1:, :-1] + _clm = _nm + holes = flood_outside(_clm) +else: + holes = flood_outside(mask) log(f" 足迹 {mask.sum()} 格, 补内部洞 {holes.sum()} 格 (镜片开口)") mask = mask | holes +# 镜片开口检测兜底(所有模式): 原始掩码上开口可能经细通道(镜圈栅格 +# 缝隙/真实镜片面片打断)与外部连通 -> 洪泛把它当外部, 开口漏检 -> +# 该区无 lens_mask 豁免, 被穿皮钳压到眼面(填充凹陷). 在膨胀掩码 +# (缝隙已闭合)的洞上重检测补齐 +if len(lens_boxes) < 2: + _comps2 = _label_comps(holes) + _big2 = sorted((_p for _p in _comps2 if len(_p) >= 250), key=len, reverse=True) + if _big2: + for _pts in _big2[:2]: + # 跳过已检出的开口(与现有 lens_mask 有重叠的洞) + _ov = 0 + for p in _pts: + if lens_mask[p[0], p[1]]: + _ov += 1 + if _ov > len(_pts) // 3: + continue + _js = [p[0] for p in _pts]; _is = [p[1] for p in _pts] + lens_boxes.append((X0 + min(_is) * RES, X0 + max(_is) * RES, + Y0 + min(_js) * RES, Y0 + max(_js) * RES)) + for p in _pts: + lens_mask[p[0], p[1]] = True + for _ in range(2): + _nm = lens_mask.copy() + _nm[1:, :] &= lens_mask[:-1, :]; _nm[:-1, :] &= lens_mask[1:, :] + _nm[:, 1:] &= lens_mask[:, :-1]; _nm[:, :-1] &= lens_mask[:, 1:] + lens_mask = _nm + log(" 镜片开口兜底(膨胀后重检测): " + + ("; ".join("x[%.1f,%.1f] y[%.1f,%.1f]" % b for b in lens_boxes) + + "; 洞分量尺寸: " + str(sorted((len(_p) for _p in _comps2), reverse=True)[:8]) + if lens_boxes else + "未检出; 洞分量尺寸: " + str(sorted((len(_p) for _p in _comps2), reverse=True)[:8]))) -# 腐蚀(膨胀0.3 + 腐蚀0.2 => 净外扩0.1): 预补偿边界平滑对凸部的压扁, -# 让最终边界贴足迹外缘, 镜框下镂空(尾部包裹段等)填满到边 +# 腐蚀(膨胀0.3 + 腐蚀0.3 => 净外扩0, 闭运算): 边界平滑收缩后 +# 仍落在栅格足迹内, XY 不超镜框轮廓 er = mask.copy() for _ in range(ERODE): ne = er.copy() @@ -374,6 +760,45 @@ log(f" R: {R.sum()} 格 ({R.sum()*RES*RES:.2f} 面积单位), " if R.sum() < 200: raise SystemExit("R 过小") +# ---- R 裁剪到眼镜真实轮廓(填充严格不超镜框范围) ---- +# 上面的区域构造会外扩: 几何模式粗栅格闭运算外扩 ~1 单位; 标记模式 +# 开口泄漏通道(镜圈栅格缝隙 0.4-0.6)让开口绕过洞检测直接并入. +# 轮廓 = dilate(g1f, 0.6) 填内部洞后整体腐蚀 0.6(闭运算, 作用在填 +# 洞后的实心区域上: 薄框带存活、开口为内部无边界、外边界光滑), +# 再并入镜片开口区(脚本自身开口检测的 lens_mask, 应对框环上角 +# 选面失败的楔形缺口>0.6 时洪泛泄漏). 闭运算是扩展运算, 轮廓必含 +# 全部覆盖格(镜腿等细结构不裁), XY 超出眼镜覆盖 <= 0.2. +_sil = g1f.copy() +for _ in range(6): + _sm = _sil.copy() + _sm[1:, :] |= _sil[:-1, :]; _sm[:-1, :] |= _sil[1:, :] + _sm[:, 1:] |= _sil[:, :-1]; _sm[:, :-1] |= _sil[:, 1:] + _sm[1:, 1:] |= _sil[:-1, :-1]; _sm[:-1, :-1] |= _sil[1:, 1:] + _sm[1:, :-1] |= _sil[:-1, 1:]; _sm[:-1, 1:] |= _sil[1:, :-1] + _sil = _sm +_sil = _sil | flood_outside(_sil) # 填内部洞(镜片开口) -> 实心轮廓 +for _ in range(6): + _sm = _sil.copy() + _sm[1:, :] &= _sil[:-1, :]; _sm[:-1, :] &= _sil[1:, :] + _sm[:, 1:] &= _sil[:, :-1]; _sm[:, :-1] &= _sil[:, 1:] + _sm[1:, 1:] &= _sil[:-1, :-1]; _sm[:-1, :-1] &= _sil[1:, 1:] + _sm[1:, :-1] &= _sil[:-1, 1:]; _sm[:-1, 1:] &= _sil[1:, :-1] + _sil = _sm +_lexr = lens_mask.copy() +for _ in range(2): + _sm = _lexr.copy() + _sm[1:, :] |= _lexr[:-1, :]; _sm[:-1, :] |= _lexr[1:, :] + _sm[:, 1:] |= _lexr[:, :-1]; _sm[:, :-1] |= _lexr[:, 1:] + _sm[1:, 1:] |= _lexr[:-1, :-1]; _sm[:-1, :-1] |= _lexr[1:, 1:] + _sm[1:, :-1] |= _lexr[:-1, 1:]; _sm[:-1, 1:] |= _lexr[1:, :-1] + _lexr = _sm +_sil = _sil | _lexr +_nclip = int((R & ~_sil).sum()) +R &= _sil +log(f" R 裁剪到眼镜轮廓: 剔除轮廓外 {_nclip} 格, R={R.sum()} 格") +if R.sum() < 200: + raise SystemExit("R 过小(轮廓裁剪后)") + # ----------------------- 高度场 ----------------------- log("== 4. 构造填充前/后高度场 ==") @@ -484,7 +909,11 @@ _nb = np.stack(_nb) _valid_nb = _nb > -1e29 with np.errstate(all='ignore'): _nbmed = np.nanmedian(np.where(_valid_nb, _nb, np.nan), axis=0) -_pin = g1f & np.isfinite(_nbmed) & (gz1 < _nbmed - 0.4) + _nbmin = np.nanmin(np.where(_valid_nb, _nb, np.inf), axis=0) +# 针孔 = 低于"所有"覆盖邻格最小值 0.4+ 的孤点. 旧判据(低于邻域中位数 +# 0.4)会误杀镜腿尾端的真实陡降格(周围皆高, 如 632667 尾尖 6.3 被补成 +# 11.0 -> 填充前端浮起 4+ 成鳍); 真陡降总有相近低的邻格(连成下降路径) +_pin = g1f & np.isfinite(_nbmed) & (gz1 < _nbmin - 0.4) gz1d = gz1f.copy() for _ in range(3): gf = gz1d.copy() @@ -502,39 +931,20 @@ for _ in range(3): gz1d = gf gz1f = np.where(g1f & ~_pin, gz1, gz1d) log(f" 顶面针孔修补: {int(_pin.sum())} 格 (孤立低谷, 用邻域最大值补回顶面)") -nfill = int((g1c & ~g1f & R & (gz1f >= 7.0)).sum()) +# 镜框带补栅格化角落缺口(阈值相对镜框管前面中位高度, 适配任意模型; +# 原写死 z>=7.0 只适配 632723) +_ztube = gz1[g1f & (gz2 > -1e29)] +_zcut = float(np.median(_ztube)) - 2.0 if len(_ztube) else -np.inf +nfill = int((g1c & ~g1f & R & (gz1f >= _zcut)).sum()) log(f" 镜框带补栅格化角落缺口: {nfill} 格") -# ---- 镜片开口前表面: 框前面曲面 = 稳健二次拟合 + 残差调和延拓 ---- -# zfit = 框前面二次趋势面(IRLS 降权误标/离群格); Rext = 拟合残差的调和延拓 -# -> 交界处自动严格连续 (zfit + (gz1 - zfit) = gz1), 内部沿框前面趋势(与框平齐) -def pbasis(x, y): - return np.stack([np.ones_like(x), x, y, x * x, x * y, y * y], 1) - -# mfit = 框前面(镜片平面方向)拟合范围: 前面区 z>=7 (框前面实际范围 8~10). -# 尾部包裹段(z 更低的陡斜坡)不参与镜片曲面拟合, 但参与框带钉扎(见 frame_band) -mfit = g1c & R & (gz1f >= 7.0) -# frame_band = 完整框带(含尾部折返段): 前面钉扎到本地最高面 - FLUSH. -# 净外扩+0.1 后, R 最外 ~0.3 带内存在"无栅格覆盖但 gz1f(最大填充)有值" -# 的外缘格(镜框轮廓边). 它们必须归框带: 若按"无覆盖"误入镜片开口, -# 前面拟合面在尾部折返斜坡处会高出包裹面 1~2 (凸出镜框). -_eR = R.copy() -for _ in range(3): - _e = _eR.copy() - for dj, di in ((1,0),(-1,0),(0,1),(0,-1),(1,1),(1,-1),(-1,1),(-1,-1)): - if dj == 1: _e[1:, :] &= _eR[:-1, :] - elif dj == -1: _e[:-1, :] &= _eR[1:, :] - elif di == 1: _e[:, 1:] &= _eR[:, :-1] - elif di == -1: _e[:, :-1] &= _eR[:, 1:] - elif dj == 1 and di == 1: _e[1:, 1:] &= _eR[:-1, :-1] - elif dj == 1 and di == -1: _e[1:, :-1] &= _eR[:-1, 1:] - elif dj == -1 and di == 1: _e[:-1, 1:] &= _eR[1:, :-1] - else: _e[:-1, :-1] &= _eR[1:, 1:] - _eR = _e -outer_band = R & ~_eR -frame_band = (g1c & R) | (outer_band & (gz1f > -1e29)) -# 陡斜坡外缘: 最大填充取邻域最高, 在折返坡会高出本地框面 ~1. -# 用 0.3 内真实覆盖格的最小值钳制(保守: 前表面绝不高于本地框面), +# ---- 前表面(v6): 跟随镜框前面 + 开口调和膜 + 整体平滑 ---- +# 有栅格覆盖处(框前面/鼻梁/镜尾包裹斜坡) 前面 = 本地框面(平齐-FLUSH); +# 无覆盖处(镜片开口/中央缝隙/外缘) = 覆盖值的调和延拓膜. 调和解内部无 +# 局部极值 => 镜片面天然无凹陷; 镜尾跟随包裹段下折(侧视无鳍状多余凸出), +# 局部低谷(鼻梁等)被整体平滑略微抬高 = 允许的轻微凸出(不突兀) +# 角落缺口格(g1c & ~g1f)的 gz1f 为形态学最大延拓, 在折返坡会高出本地框面 ~1. +# 用 0.3 内真实覆盖格的最小值钳制(让覆盖值跟随本地框面), # 再轻平滑消除形态学噪声 _mdil = np.where(g1f, gz1, np.inf) for _ in range(3): @@ -551,96 +961,195 @@ for _ in range(3): else: s[:-1, :-1] = _mdil[1:, 1:] _mn = np.minimum(_mn, s) _mdil = _mn -_fringe = frame_band & ~g1f -gz1f = np.where(_fringe & (_mdil < 1e29), np.minimum(gz1f, _mdil), gz1f) +_fringe = g1c & ~g1f +# 最小值钳制扩展到所有"无覆盖但 0.3 内有覆盖"的格(不止角落缺口格): +# 折叠处(镜圈-镜腿仅隔 0.1~0.3, 如 632667)闭合掩码 g1c 有洞, 未钳制的 +# 格保留形态学最大延拓(取自镜圈 10.9), 与相邻钳制格(取自镜腿 7.9)交替 +# -> 种子/天花板基准场本身振荡 -> 前表面出现 1+ 深的凹坑 +gz1f = np.where((~g1f) & (_mdil < 1e29), np.minimum(gz1f, _mdil), gz1f) with np.errstate(all='ignore'): gz1f = np.where(_fringe, smooth_valid(gz1f, _fringe, sigma=2.0, it=1), gz1f) -xp, yp, zp = gx[mfit], gy[mfit], gz1f[mfit] -wt = np.ones(len(zp)) -for _ in range(3): - Aw = pbasis(xp, yp) * wt[:, None] - coef, *_ = np.linalg.lstsq(Aw, zp * wt, rcond=None) - rfit = zp - pbasis(xp, yp) @ coef - wt = 1.0 / (1.0 + (np.abs(rfit) / 0.5) ** 2) -zfit = (pbasis(gx.ravel(), gy.ravel()) @ coef).reshape(H, W) -rf = np.abs((zfit - gz1f)[mfit]) -log(f" 框前面拟合(稳健): 残差 mean {rf.mean():.3f}, p90 {np.percentile(rf, 90):.3f}, " - f"max {rf.max():.3f}") +if _GEO_MODE: + # 陡降段基准软化: 几何模式选面在折返框陡壁上的逐格最大 z 是台阶状, + # 而种子/走廊钳制(天花板 gz1f+0.3 / 下界 gz1f-0.15)每轮平滑后都把 + # 前表面拉回该台阶基准 -> 前表面台阶化(片二面角 p90 ~3 倍于标记模式, + # 632723/632667 无此问题: 稠密标记的框面本来就平滑). 对陡降格 + # (1.0 邻域内框面高出本地 >=1.2, 与走廊天花板同判据)的基准 gz1f 做 + # 掩码高斯平滑: 高斯保线性 -> 光滑斜坡原样保留, 只消台阶抖动; 再用 + # 3x3 原始最大值封顶, 基准不超出近旁真实框前面 -> 防鳍性质保留 + # (填充不伸到框面前方之外, 与原 +0.3 天花板余量同级). 平缓格/镜片 + # 膜/标记模式不受影响. + _gzv = gz1f > -1e29 + _gzm = np.where(_gzv, gz1f, np.nan) + _gmax0 = np.where(_gzv, gz1f, -np.inf) + _gmx = _gmax0.copy() + for _ in range(10): + _gmx = np.maximum.reduce([_gmx, np.roll(_gmx, 1, 0), np.roll(_gmx, -1, 0), + np.roll(_gmx, 1, 1), np.roll(_gmx, -1, 1)]) + _steepb = _gzv & (_gmax0 < _gmx - 1.2) + if _steepb.any(): + _gzs = smooth_valid(_gzm, _gzv, sigma=2.0, it=1) + _m3 = np.full((H, W), -np.inf) + for _dj in (-1, 0, 1): + for _di in (-1, 0, 1): + _m3 = np.maximum(_m3, np.roll(_gmax0, (_dj, _di), (0, 1))) + _dlt = np.where(_steepb & np.isfinite(_gzs), _gzs - gz1f, np.nan) + gz1f = np.where(_steepb & np.isfinite(_gzs), + np.minimum(_gzs, _m3), gz1f) + log(f" 陡降段基准软化(几何模式): {int(_steepb.sum())} 格, " + f"基准改动量中位 {np.nanmedian(np.abs(_dlt)):.2f} " + f"max {np.nanmax(np.abs(_dlt)):.2f}") +# 前表面: 边界条件 = 覆盖格(含钳制后的角落缺口格)的框前面, 轻平滑去栅格 +# 噪声(σ 取小, 防止把镜尾陡崖格拉向邻值), 调和延拓到无覆盖区 +# (镜片开口/缝隙/外缘) —— 调和解内部无局部极值 => 无凹陷 +seed = np.full((H, W), np.nan) +bc = R & g1c & (gz1f > -1e29) +seed[bc] = gz1f[bc] +seed = smooth_valid(seed, bc, sigma=1.0) +hf = laplace_fill(seed, bc) +hf = np.where(R, hf - FLUSH, np.nan) +# 镜片弧面(轻微凸出): 以镜框前方两侧(镜圈前面)为缘, 向开口内部连成小 +# 弧度抛物面: 凸出量 = BULGE * (1-(1-d/D)^2), d = 离填充区边界的精确 +# 欧氏距离(城市块距离的水平集是菱形, 沿对角线有折角, 会让弧面起折痕), +# D = 全区最大深度. 镜片中心约 +0.35(整体平滑后), 边缘归零 -> 与框带 +# 平滑衔接(不突兀); 盖住镜框前方属于允许行为. 其后走廊钳制仍生效 +# (镜尾陡降格防鳍不受影响) +_nbr = (np.roll(R, 1, 0) | np.roll(R, -1, 0) | np.roll(R, 1, 1) | np.roll(R, -1, 1)) +_ring = ~R & _nbr # R 边界外圈格(距离源) +_jr, _ir = np.where(R) +_js, _is = np.where(_ring) +_dx = (_ir[:, None] - _is[None, :]) * RES +_dy = (_jr[:, None] - _js[None, :]) * RES +_dist = np.full((H, W), np.inf) +_dist[_jr, _ir] = np.sqrt(_dx * _dx + _dy * _dy).min(1) +dmax = float(np.max(_dist[R])) +mprof = np.where(R, 1.0 - (1.0 - np.minimum(_dist / max(dmax, 1e-9), 1.0)) ** 2, + np.nan) +mprof = smooth_valid(mprof, R, sigma=2.5) +hf = np.where(R, hf + BULGE * mprof, np.nan) +# 整体平滑: 柔化框带/镜片膜交界折角(打印外观光滑) +hf = smooth_valid(hf, R, sigma=2.0, it=2) +# 走廊钳制: 平滑/延拓会把陡崖格与缺口格拉向邻值, 造成镜尾鳍状凸出(向上) +# 或管底缝隙(向下, 如管顶 6.0 邻格 1.8 的悬崖). 钳制范围与量: +# 上界: 仅陡降格(1.0 邻域内框面高出本地 >=1.2: 镜尾尖端/对角缺口/陡崖) +# 钳到 gz1f+0.3 防鳍; 平缓格不设天花板 -- 弧面允许盖住镜框前方(用户 +# 许可), 单格最大值代理会压出尖点凹坑. 邻域取 1.0: 尾尖孤立低点(四周 +# 皆低)会被整体平滑从远处(镜圈在 ~1.0 外)拉高 ~3, 窄邻域看不见 -> 悬浮 +# 鳍(632667 尾尖曾 fill 高出框面 3.0/皮肤 2.5) +# 下界(仅真实覆盖格) = gz1f - 0.15: ① 真实镜框管深(前后表面差)>=0.27, +# 故该下界恒在管背之上 -> 与框体交叠无缝(旧值 -0.5 按粗框管半径 0.67 +# 调, 细框模型管深 0.3 时填充顶会低于管背 -> 内部空隙); ② 防止整体 +# 平滑把覆盖带前表面拉离本地框面超过 0.15(陡崖邻域的拖拽, 632667 的 +# 镜圈-镜腿折叠处曾拉低 0.5 造成交界折角) +# 钳后轻平滑柔化折痕, 再钳(两轮收敛) +# 镜片邻近带(镜片开口+0.6): 防鳍天花板豁免区. 开口边缘无覆盖格的延拓 +# gz1f 取自镜圈内侧面(比框前沿低 ~1), 与 1.0 内框前沿差 >1.2 会被误判 +# "陡降格"钳到 gzfm+0.3, 在穹顶边缘压出 0.4-0.7 深的沟(617892 右镜下缘 +# ceil=9.36 恰等于沟底). 穹顶/开口边缘允许升到框面, 不属鳍. +_lex = lens_mask.copy() +for _ in range(6): + _lg = _lex.copy() + for _dj, _di in ((1, 0), (-1, 0), (0, 1), (0, -1), (1, 1), (1, -1), (-1, 1), (-1, -1)): + _lg |= np.roll(_lex, (_dj, _di), (0, 1)) + _lex = _lg +gzfm = np.where(gz1f > -1e29, gz1f, -np.inf) +gmax2 = gzfm.copy() +for _ in range(10): + gmax2 = np.maximum.reduce([gmax2, np.roll(gmax2, 1, 0), np.roll(gmax2, -1, 0), + np.roll(gmax2, 1, 1), np.roll(gmax2, -1, 1)]) +# 防鳍豁免镜片邻近带: 见上方 _lex 注释(镜片区误钳出沟) +ceil = np.where((gz1f > -1e29) & (gzfm < gmax2 - 1.2) & ~_lex, gzfm + 0.30, np.inf) +floor = np.where(bc, gz1f - 0.15, -np.inf) +if _GEO_MODE: + # 贴脸锚连续权重版: 走廊/覆盖格的真实前界 = 该列最前表面 = max(镜框延 + # 拓, 脸面). 脸高于框参考 gap>0.35 的格进入贴脸影响带, 地板 = 贴脸目标 + # T - (0.15 + 0.5*(1-w)), w = clip((gap-0.35)/1.65) 随 gap 线性收紧 + # [T-0.65, T-0.15] —— 连续无阈值悬崖(阈值版在边缘一格差 2.5, 83°二 + # 面角); 填充贴到脸面(镂空无气袋), 过渡带留 0.65 圆角余量. gap>=2(深 + # 臂覆盖格 gz1f 可为负深度, 旧版前表面潜入峡谷 80°; 走廊袋)时天花板 + # 直接覆盖反鳍天花板(负深度会削回峡谷); 0.35 -1e29, gz1f, -np.inf) + _hasf = R & (gz1f > -1e29) & have_f & (fz > _gzf + 0.35) + _pcl = _hasf & (fz > _gzf + 2.0) + if _hasf.any(): + _fzs = smooth_valid(np.where(have_f, fz, np.nan), have_f, sigma=1.5, it=1) + _T = np.where(_hasf, np.minimum(_fzs, fz) - 0.05, np.nan) + _gap = np.where(_hasf, fz - _gzf, 0.0) + _w = np.clip((_gap - 0.35) / 1.65, 0.0, 1.0) + _fl = _T - 0.15 - 0.5 * (1.0 - _w) + floor = np.maximum(floor, np.where(_hasf, np.nan_to_num(_fl, nan=-np.inf), -np.inf)) + ceil = np.where(_pcl & np.isfinite(_T), _T, ceil) + ceil = np.maximum(ceil, np.where(_hasf & np.isfinite(_T), _T, -np.inf)) + # 贴脸带 0.5 邻域内无脸格(镜尾走廊)天花板放开到贴脸目标, 铺斜坡; + # 有脸面处以 min(平滑脸面, 脸面)-0.05 封顶(不穿皮, 不出脸) + _bd = np.where(_pcl & np.isfinite(_T), _T, -np.inf) + for _ in range(5): + _bd = np.maximum.reduce([_bd, np.roll(_bd, 1, 0), np.roll(_bd, -1, 0), + np.roll(_bd, 1, 1), np.roll(_bd, -1, 1)]) + _capr = np.where(have_f & np.isfinite(_fzs), np.minimum(_fzs, fz) - 0.05, np.inf) + _raise = np.where(_bd > -1e29, np.minimum(_bd, _capr), -np.inf) + ceil = np.maximum(ceil, np.where(R & ~g1f & ~_hasf, _raise, -np.inf)) + log(f" 贴脸锚(连续权重): 影响带 {int(_hasf.sum())} 格, " + f"其中贴皮带(gap>=2) {int(_pcl.sum())} 格") -# 残差场边界只用拟合内点: 尾部包裹段(偏离前面平面 1 以上)不参与, -# 其低残差经平滑/调和延拓渗入镜片开口会拉低镜片面 -mfit2 = mfit & (np.abs(gz1f - zfit) <= 1.0) - -# 残差场: 镜框带上 = gz1 - zfit; 先平滑边界条件(消除格尺度噪声 -> 延拓解光滑), -# 再调和延拓到开口内部 -resf = np.full((H, W), np.nan) -resf[mfit2] = gz1f[mfit2] - zfit[mfit2] -resf = smooth_valid(resf, mfit2, sigma=2.5) -Rext = laplace_fill(resf, mfit2) -Rext = np.where(R, Rext, np.nan) - -# 前表面: 完整框带(含尾部)贴住本地最高面(FLUSH 嵌入余量, 尾部跟随包裹斜坡); -# 开口 = zfit + Rext - FLUSH -hf = np.full((H, W), np.nan) -hf[frame_band] = gz1f[frame_band] - FLUSH -openf = R & ~frame_band -hfo = np.where(openf, zfit + Rext - FLUSH, np.nan) - -# 开口区安全钳制(不得超过 3 格半径内框前沿), 在平滑前施加 -> 折角被平滑柔化 -gmax2 = np.where(mfit, gz1f, -np.inf) -for _ in range(3): - gm = gmax2.copy() - for dj, di in ((1,0),(-1,0),(0,1),(0,-1),(1,1),(1,-1),(-1,1),(-1,-1)): - s = np.full_like(gmax2, -np.inf) - if dj == 1: s[1:, :] = gmax2[:-1, :] - elif dj == -1: s[:-1, :] = gmax2[1:, :] - elif di == 1: s[:, 1:] = gmax2[:, :-1] - elif di == -1: s[:, :-1] = gmax2[:, 1:] - elif dj == 1 and di == 1: s[1:, 1:] = gmax2[:-1, :-1] - elif dj == 1 and di == -1: s[1:, :-1] = gmax2[:-1, 1:] - elif dj == -1 and di == 1: s[:-1, 1:] = gmax2[1:, :-1] - else: s[:-1, :-1] = gmax2[1:, 1:] - gm = np.maximum(gm, s) - gmax2 = gm -hfo = np.where(openf & (gmax2 > -1e29), np.minimum(hfo, gmax2 - FLUSH), hfo) -# 开口区整体平滑(镜片表面打印光滑) -hfo = smooth_valid(hfo, openf, sigma=2.0, it=2) -hf = np.where(openf & np.isfinite(hfo), hfo, hf) - -# 整体轻度平滑柔化交界折角 -hf = smooth_valid(hf, R, sigma=1.5) - -# 墙区(R 外缘 4 格): 侧墙顶所在, 需严格嵌入框实体 -edge_in = R.copy() -for _ in range(4): - ne = edge_in.copy() - ne[1:, :] &= edge_in[:-1, :]; ne[:-1, :] &= edge_in[1:, :] - ne[:, 1:] &= edge_in[:, :-1]; ne[:, :-1] &= edge_in[:, 1:] - edge_in = ne -wall_zone = R & ~edge_in - -# 安全钳制(平滑后仅作用于框带; 开口的钳制已在平滑前施加并被柔化): -# 框带墙区 -> 本地框面 - FLUSH(严格嵌入, 尾部跟随包裹斜坡) -# 框带内部 -> 至多贴到框面 + 0.10(交界平齐, 并集无缝) -hf = np.where(frame_band & wall_zone, np.minimum(hf, gz1f - FLUSH), hf) -hf = np.where(frame_band & ~wall_zone, np.minimum(hf, gz1f + 0.10), hf) - -# 交界连续性: 镜片开口格与相邻"平面前框带"格(与镜片面近共面, |gz1f-zfit|<=0.3) -# 的前表面高度差. 陡壁/包裹段格(镜圈管内壁等)与镜片面的台阶是镜框自身几何, -# 并集后不暴露到外表面, 不计入此项 -jn, jn_at = 0.0, None -bflat = mfit2 & (gz1f > zfit - 0.3) -for dj, di in ((1,0),(-1,0),(0,1),(0,-1)): - a = np.roll(bflat, (dj, di), axis=(0, 1)) & openf & np.isfinite(hf) - if a.any(): - b = np.roll(a, (-dj, -di), axis=(0, 1)) - dd = np.abs(np.where(a, hf, 0.0)[a] - np.where(np.isfinite(hf), hf, 0.0)[b]) - k = int(np.nanargmax(dd)) - if dd[k] > jn: - ja = np.argwhere(a)[k] - jn, jn_at = float(dd[k]), (X0 + ja[1] * RES, Y0 + ja[0] * RES) -log(f" 交界处高度跳变最大 {jn:.3f} @ {jn_at} (平面前框带接缝应<0.3; " - "镜圈圆角处~0.4为镜框自身几何, 并集后不暴露)") +# ---- 穿皮防护(全模式): 无镜框覆盖且非镜片穹顶的真实脸面格, 前面不高出脸面 ---- +# 走廊/外缘格的膜面跟随镜框平面(防鳍天花板来自延拓框面), 而该处皮肤 +# 退到框平面下方 1-3 -> 填充高出皮肤, 正面/侧面可见"穿透脸皮" +# (20261102 外上角最显著 +2.8). 钳到 T=min(平滑脸,脸)-0.05; 有镜框 +# 直接覆盖(g1c)或镜片穹顶区(含 0.4 边缘环, 穹顶高于其后眼面属镜片 +# 设计)豁免. 距覆盖带 4 格内线性松弛(每格 0.45), 与钳制-平滑循环 +# 一起铺成连续斜坡, 避免框缘一步崖(v9 教训: 硬钳=103° 二面角). +_fzs2 = smooth_valid(np.where(have_f, fz, np.nan), have_f, sigma=1.5, it=1) +_rf = have_f & (fz >= FACE_MIN) +_T2 = np.where(_rf, np.minimum(np.nan_to_num(_fzs2, nan=np.inf), fz) - 0.05, + np.inf) +# _lex(镜片开口+0.6 膨胀带)已在防鳍天花板前计算(上方), 此处直接引用 +# 豁免 = 镜框直接覆盖(g1c) 或 镜片穹顶区(含 0.4 边缘环; 穹顶高于其后 +# 眼面属镜片设计). 注意不能用延拓 gz1f 判"框在上方": 延拓只说明框在 +# 0.2-0.3 外侧, 该列正上方无框, 填充高出照样可见. +_exempt = g1c | _lex +_dcv = np.full((H, W), 99) +_dcv[g1c] = 0 +_dcur = g1c.copy() +for _s in range(1, 6): + _dg = _dcur.copy() + for _dj, _di in ((1, 0), (-1, 0), (0, 1), (0, -1)): + _dg |= np.roll(_dcur, (_dj, _di), (0, 1)) + _nw = _dg & ~_dcur + _dcv[_nw] = _s + _dcur = _dg +# 局部框面前沿参考(6 格延拓, 供落差计算): 斜坡从框面线性降到脸面 +_gref = gzfm.copy() +for _ in range(6): + _gg = _gref.copy() + for _dj, _di in ((1, 0), (-1, 0), (0, 1), (0, -1)): + _gg = np.maximum(_gg, np.roll(_gref, (_dj, _di), (0, 1))) + _gref = _gg +_drop = np.where(np.isfinite(_T2) & np.isfinite(_gref), + np.clip(_gref + 0.15 - _T2, 0.0, None), 0.0) +_slack = _drop * np.maximum(0, 3 - np.minimum(_dcv, 3)) / 3.0 +# 穿皮钳(R 内)停用: R 已裁剪到眼镜轮廓+0.2, 轮廓内填充前面高出脸面 = +# 填满镜框-脸面之间镂空的支撑体本体(设计意图, 镜片穹顶/镜桥同理). +# 旧版把无框覆盖但有脸面的格(镜桥 pinch/镜下袋)钳到脸面-0.05 -> 前表面 +# 凹陷(617892 鼻梁塌 2.8, 20261102 塌 0.6-0.9, 用户明确要求镜片区无凹陷). +# 原防穿皮目标(填充伸出镜框轮廓外高出脸面)已由 R 轮廓裁剪根治; +# 防鳍天花板(陡降格 gzfm+0.3)与贴脸锚(脸面在框延拓前方的融合带)仍生效. +_skincap = np.zeros_like(R) +log(" 穿皮防护(内部): 已停用 (R 限轮廓内; 镜桥/镜下袋/镜片穹顶统一填到框面)") +hf = np.minimum(np.maximum(hf, floor), ceil) +for _ in range(2): + hf = smooth_valid(hf, R, sigma=1.0, it=1) + hf = np.where(R, np.minimum(np.maximum(hf, floor), ceil), np.nan) +# 几何模式袋状斜坡追加收敛: 主循环 2 轮 σ1.0 不足以把 2+ 单位悬崖完全 +# 铺开, 追加 2 轮(钳制同前, 袋状贴皮带/防穿皮帽持续生效) +if _GEO_MODE and _pcl.any(): + for _ in range(2): + hf = smooth_valid(hf, R, sigma=1.0, it=1) + hf = np.where(R, np.minimum(np.maximum(hf, floor), ceil), np.nan) +log(f" 前表面(框面跟随+调和膜)范围: z [{np.nanmin(hf):.2f}, {np.nanmax(hf):.2f}]") hb = np.full((H, W), np.nan) # 真实脸面: fz >= FACE_MIN. 低于此值的"脸面"是误标导致的深伪面(后脑/口腔), @@ -651,8 +1160,8 @@ knownb = R & real_f hb = laplace_fill(hb, knownb) hb = smooth_valid(hb, R) -# 无真实脸面处(镜桥下/误标区): 填充厚度不超过 hf - SHELL_MAX, -# 保证填充只在镜框附近的 z 范围内; 有真实脸面处仍嵌到脸面(填满镂空) +# 无真实脸面处(误标区): 厚度上限仅作安全钳(透明支撑体尾部较厚属预期); +# 有真实脸面处仍嵌到脸面(填满镂空) cap = hf - SHELL_MAX ncap = int((R & ~real_f & np.isfinite(hb) & (hb < cap)).sum()) hb = np.where(R & ~real_f, np.maximum(hb, cap), hb) @@ -681,6 +1190,24 @@ for _ in range(60): s = np.roll(np.where(nanm, 0, hb), (dj, di), axis=(0, 1)) ok = ~np.roll(nanm, (dj, di), axis=(0, 1)) & np.isfinite(s) hb = np.where(nanm & ok, s, hb) +# 外扩带坡度延拓: 逐格复制外扩是平檐, 与陡坡(包裹带/袋状斜坡)交界成折角 +# (前表面外缘片二面角 p90 ~35°). 对 R 外的有限格做掩码平滑, 坡面平滑延续 +# 到边界环(仅影响边界采样, 不动 R 内字段); 并保持外扩带厚度 >= 0.1 防 +# 前后面交叠. +_fring = np.isfinite(hf) & ~R +if _fring.any(): + with np.errstate(all='ignore'): + _hfs = smooth_valid(hf, np.isfinite(hf), sigma=1.0, it=2) + _hbs = smooth_valid(hb, np.isfinite(hb), sigma=1.0, it=2) + hb = np.where(_fring, _hbs, hb) + hf = np.where(_fring, np.maximum(_hfs, hb + 0.10), hf) + # 穿皮防护延伸到外扩带: fringe 会被 σ 平滑从内部覆盖带拉高, 而 + # 边界环顶点采样的是 fringe 值 -> 边界环仍高出皮肤(外上角 +2.3, + # 内部钳制对此无效). 用同一 T+距框松弛钳外扩带, 边界环收口到皮肤. + # 外扩带穿皮钳停用: 它把边界环压到皮肤, 前表面边缘与底面在皮肤处 + # 捏成刀刃(脸面距框面 1-3 时侧视尖角). 停用后边界环 = 坡度延拓 + # (≈本地框面), 侧墙从框面平滑落到脸面 = 实心支撑体外形. + log(" 穿皮防护(外扩带): 已停用 (边界环=坡度延拓, 侧墙平滑落脸)") thickm = np.where(R, hf - hb, np.nan) log(f" 填充厚度: min {np.nanmin(thickm):.2f}, 中位 {np.nanmedian(thickm):.2f}, " @@ -852,7 +1379,9 @@ for _ in range(2): ne[1:, :] &= inner2[:-1, :]; ne[:-1, :] &= inner2[1:, :] ne[:, 1:] &= inner2[:, :-1]; ne[:, :-1] &= inner2[:, 1:] inner2 = ne -step = max(1, int(round(0.3 / RES))) +# 网格步长 0.2(原 0.3): 弧面鞍形曲率处交替对角线扭折 ~ 间距^2, 加密后 +# 前表面二面角明显下降(镜片光滑) +step = max(1, int(round(0.2 / RES))) cand = [] for j in range(0, H, step): for i in range(0, W, step): @@ -899,19 +1428,55 @@ def tri_areas(xy, T): cr = xy[T[:, 1]] - xy[T[:, 0]] cs = xy[T[:, 2]] - xy[T[:, 0]] return cr[:, 0] * cs[:, 1] - cr[:, 1] * cs[:, 0] -xy_best = xy.copy() -for _ in range(8): - sx = np.bincount(src, weights=xy[dst, 0], minlength=len(xy)) - sy = np.bincount(src, weights=xy[dst, 1], minlength=len(xy)) - new = np.stack([sx, sy], 1) / np.maximum(deg, 1)[:, None] - xy2 = xy.copy() - xy2[Nb:] = 0.5 * xy[Nb:] + 0.5 * new[Nb:] - if tri_areas(xy2, T).min() > 1e-6: - xy = xy2 - xy_best = xy2.copy() - else: - break -xy = xy_best +if _GEO_MODE: + # 几何模式逐点 Laplacian 平滑: 全局整轮接受/回退在楔形缺口(细长三 + # 角形)处第一轮就出现反演 -> 完全不平滑 -> 网格噪声大(镜片二面角 + # 中位数翻倍). 逐点接受: 只移动使所有关联三角形保持正面积的顶点, + # 其余留原地, 细长三角形逐步被规整化 + _vf = [[] for _ in range(len(xy))] + for _fi, (_a, _b, _c) in enumerate(T): + _vf[_a].append(_fi); _vf[_b].append(_fi); _vf[_c].append(_fi) + for _r in range(8): + _moved = 0 + for _v in range(Nb, len(xy)): + _nb = set() + for _fi in _vf[_v]: + for _u in T[_fi]: + if _u != _v: + _nb.add(int(_u)) + if not _nb: + continue + _tgt = xy[list(_nb)].mean(0) + _newp = 0.5 * xy[_v] + 0.5 * _tgt + _ok = True + for _fi in _vf[_v]: + _t = T[_fi] + _p0 = _newp if _t[0] == _v else xy[_t[0]] + _p1 = _newp if _t[1] == _v else xy[_t[1]] + _p2 = _newp if _t[2] == _v else xy[_t[2]] + _cr = _p1 - _p0; _cs = _p2 - _p0 + if _cr[0] * _cs[1] - _cr[1] * _cs[0] <= 1e-9: + _ok = False + break + if _ok: + xy[_v] = _newp + _moved += 1 + if _moved == 0: + break +else: + xy_best = xy.copy() + for _ in range(8): + sx = np.bincount(src, weights=xy[dst, 0], minlength=len(xy)) + sy = np.bincount(src, weights=xy[dst, 1], minlength=len(xy)) + new = np.stack([sx, sy], 1) / np.maximum(deg, 1)[:, None] + xy2 = xy.copy() + xy2[Nb:] = 0.5 * xy[Nb:] + 0.5 * new[Nb:] + if tri_areas(xy2, T).min() > 1e-6: + xy = xy2 + xy_best = xy2.copy() + else: + break + xy = xy_best # 边翻转优化(Lawson): 公共边换成严格更短的对角线(总边长严格递减 -> 必收敛), # 不新增顶点即可消除耳切弦两侧的细长退化三角形(表面粗糙的根源) @@ -943,6 +1508,11 @@ def flip_pass(Tl, XY): s4 = (q2[0]-q1[0])*(pv[1]-q1[1]) - (q2[1]-q1[1])*(pv[0]-q1[0]) if s3 * s4 >= 0: continue + # 共线格点容差: 顶点恰好落在对角线上时(|s|~1e-15, 浮点噪声使严格 + # 符号判断失效), 翻转会生成零面积/翻反面三角形(定向校验失败); + # 要求四个点-线有向面积都显著大于零 + if min(abs(s1), abs(s2), abs(s3), abs(s4)) < 1e-9: + continue if np.hypot(*(XY[w2] - XY[w1])) >= np.hypot(*(XY[v] - XY[u])): continue # 新对角线须更短 # 凸四边形 u->a->v->b (a 在 u->v 右侧, b 在左侧) 沿对角线 a-b 拆分 @@ -978,9 +1548,20 @@ for _p in range(6): break mid = {} newpts = [] + _pos = {(round(float(p[0]), 6), round(float(p[1]), 6)): k + for k, p in enumerate(xy)} for u, v in longe: - mid[(u, v)] = len(xy) + len(newpts) - newpts.append((xy[u] + xy[v]) / 2.0) + m = (xy[u] + xy[v]) / 2.0 + _key = (round(float(m[0]), 6), round(float(m[1]), 6)) + _hit = _pos.get(_key) + if _hit is not None: + # 中点与已有顶点重合(边界/插入点落在对角格点上): 复用旧点, + # 否则产生零长度边 -> 零面积三角形(定向校验失败) + mid[(u, v)] = _hit + else: + mid[(u, v)] = len(xy) + len(newpts) + _pos[_key] = mid[(u, v)] + newpts.append(m) xy = np.vstack([xy, np.asarray(newpts)]) nsplit += len(newpts) newT = [] @@ -1033,6 +1614,62 @@ zb = bilinear(hb, xy) zb = np.minimum(zb, zf - 0.15) zb = np.maximum(zb, zf - 50) +# 前表面对角线优化: 高度场本身光滑, 但网格四边形对角线方向交替会在鞍形 +# 曲率处产生小扭折(前表面二面角 ~6°). 对"前表面-前表面"内部边尝试翻转 +# 对角线, 二面角显著减小则翻(边界边不动 => 侧墙环不受影响; 面片绕向按 +# XY 正面积统一, 边界有向环自动保持一致) +def _tri_nrm(P3, t): + n = np.cross(P3[t[1]] - P3[t[0]], P3[t[2]] - P3[t[0]]) + return n / max(np.linalg.norm(n), 1e-12) + +def _pair_ang(P3, t0, t1): + n0, n1 = _tri_nrm(P3, t0), _tri_nrm(P3, t1) + return np.degrees(np.arccos(np.clip(n0 @ n1, -1, 1))) + +P3f = np.column_stack([xy, zf]) +for _round in range(4): + em = {} + for fi, (a, b, c) in enumerate(T): + for u, v in ((a, b), (b, c), (c, a)): + em.setdefault((min(u, v), max(u, v)), []).append(fi) + nfl = 0 + for (u, v), fsl in list(em.items()): + if len(fsl) != 2: + continue + f0, f1 = fsl + # em 在轮内会过期: 翻转会把四边形边界边在两面间重新分配, 后续条目 + # 引用的面可能已不含该边. 要求两面各恰含 u,v 一次(第三点唯一), 否则 + # 视为过期条目跳过(下一轮重建 em 后再处理) + _p = [t for t in T[f0] if t != u and t != v] + _q = [t for t in T[f1] if t != u and t != v] + if len(_p) != 1 or len(_q) != 1 or _p[0] == _q[0]: + continue + p, q = _p[0], _q[0] + # 凸四边形检验: 两条对角线在 XY 严格相交才可翻 + cr = ((xy[v, 0] - xy[u, 0]) * (xy[p, 1] - xy[u, 1]) - + (xy[v, 1] - xy[u, 1]) * (xy[p, 0] - xy[u, 0])) + cr2 = ((xy[v, 0] - xy[u, 0]) * (xy[q, 1] - xy[u, 1]) - + (xy[v, 1] - xy[u, 1]) * (xy[q, 0] - xy[u, 0])) + cr3 = ((xy[q, 0] - xy[p, 0]) * (xy[u, 1] - xy[p, 1]) - + (xy[q, 1] - xy[p, 1]) * (xy[u, 0] - xy[p, 0])) + cr4 = ((xy[q, 0] - xy[p, 0]) * (xy[v, 1] - xy[p, 1]) - + (xy[q, 1] - xy[p, 1]) * (xy[v, 0] - xy[p, 0])) + if not (cr * cr2 < 0 and cr3 * cr4 < 0): + continue + t1 = (p, q, u) if ((xy[q, 0] - xy[p, 0]) * (xy[u, 1] - xy[p, 1]) - + (xy[q, 1] - xy[p, 1]) * (xy[u, 0] - xy[p, 0])) > 0 else (u, q, p) + t2 = (p, q, v) if ((xy[q, 0] - xy[p, 0]) * (xy[v, 1] - xy[p, 1]) - + (xy[q, 1] - xy[p, 1]) * (xy[v, 0] - xy[p, 0])) > 0 else (v, q, p) + a_before = _pair_ang(P3f, tuple(T[f0]), tuple(T[f1])) + a_after = _pair_ang(P3f, t1, t2) + # 翻转后二面角改善 >0.3° 才翻(防振荡) + if a_after < a_before - 0.3: + T[f0], T[f1] = t1, t2 + nfl += 1 + if nfl == 0: + break +log(f" 前表面对角线优化: 翻转 {nfl} 边(减小鞍形曲率处折角)") + Pf = np.column_stack([xy, zf]) Pb = np.column_stack([xy, zb]) faces = [] @@ -1099,7 +1736,10 @@ lang, lloc = [], [] for (u, v), fs in emap2.items(): if len(fs) == 2 and nv[fs[0], 2] > 0.5 and nv[fs[1], 2] > 0.5: m = fcen[fs[0]] - if 2.2 < abs(m[0]) < 6.3 and 12.5 < m[1] < 16.2: # 左右镜片开口 + # 左右镜片膜主体(开口洞腐蚀 0.2 后的掩码, 自动适配任意模型) + _gi = int((m[0] - X0) / RES); _gj = int((m[1] - Y0) / RES) + if (0 <= _gj < lens_mask.shape[0] and 0 <= _gi < lens_mask.shape[1] + and lens_mask[_gj, _gi]): lang.append(np.degrees(np.arccos(np.clip(nv[fs[0]] @ nv[fs[1]], -1, 1)))) lloc.append(m) lang = np.asarray(lang) @@ -1108,11 +1748,17 @@ if len(lang): log(f" 镜片光滑度(相邻面二面角): 中位 {np.median(lang):.2f}°, p90 {np.percentile(lang, 90):.2f}°, " f"max {lang[k]:.2f}° @ ({lloc[k][0]:.1f},{lloc[k][1]:.1f}) (中位应<2°)") -# 嵌入校验 -gmargin = (gz1f - hf)[R & frame_band & wall_zone] -log(f" 前表面低于本地框面(墙区): 最小余量 {np.nanmin(gmargin):.3f} (应>0)") -gmargin2 = (gz1f - hf)[R & frame_band] -log(f" 前表面低于本地框面(内部允许贴齐): 最小余量 {np.nanmin(gmargin2):.3f} (应>-0.10)") +# 前表面校验: 陡降格(镜尾尖端/缺口)走廊 +0.3 内防鳍; 平缓格允许弧面 +# 盖住框前方(轻微凸出); 低于框面属正常(框体覆盖). 贴皮带(前表面在本 +# 地脸面 0.70 内: 贴皮锚/连续权重过渡带, 几何模式)以皮为外边界, 凸出 +# 非鳍, 不计入鳍检查(过渡带允许 0.65 圆角下沉) +_skinc = R & have_f & (fz >= hf - 0.70) +_mchk = R & g1c & (gz1f > -1e29) & ~_skinc +_steep = _mchk & (gz1f < gmax2 - 1.2) +gf = (hf - gz1f)[_mchk] +gs = (hf - gz1f)[_steep] +log(f" 前表面相对本地框面: 陡降格最大凸出 {np.nanmax(gs):.3f} (应<0.35), " + f"平缓格最大 {np.nanmax(gf):.3f} (弧面盖框, 允许), 最大低于 {np.nanmin(gf):.3f}") bmargin = (hb - bzs)[R & have_b] log(f" 后表面高于最后表面: 最小余量 {np.nanmin(bmargin):.3f} (应>0)") embed = (fz - hb)[R & real_f] diff --git a/fill2sc/generate_fill_mesh_hsc.py b/fill2sc/generate_fill_mesh_hsc.py index ee6a722..f7913d7 100644 --- a/fill2sc/generate_fill_mesh_hsc.py +++ b/fill2sc/generate_fill_mesh_hsc.py @@ -58,6 +58,7 @@ class FillMeshGenerator: if in_idx: self.IN_IDX = in_idx self.report = [] + self.lens_mask = None def log(self, s): """日志记录""" @@ -401,110 +402,153 @@ class FillMeshGenerator: return gz1, gz2, fz, bz, have_g1, have_g2, have_f, have_b, g1f, W, H - # ====================================================================== - # 步骤 3: 构建足迹 mask 与填充区域 R - # ====================================================================== def build_footprint(self, g1f): - """构建眼镜足迹 mask 和填充区域 R""" - self.log("== 3. 眼镜足迹 / 填充区域 R ==") - mask = g1f.copy() + """构建眼镜足迹 mask 和填充区域 R""" + self.log("== 3. 眼镜足迹 / 填充区域 R ==") + self.lens_mask = None # <-- 新增: 初始化 + mask = g1f.copy() - # 取主分量 - mask, n_stray = self.largest_component_8(mask) - self.log(f" 足迹取主分量: 剔除杂散块 {n_stray} 格 (误标到镜框外的皮肤/后脑)") + # 取主分量 + mask, n_stray = self.largest_component_8(mask) + self.log(f" 足迹取主分量: 剔除杂散块 {n_stray} 格 (误标到镜框外的皮肤/后脑)") - # 中央区剪除(鼻侧多余填充) - holes = ~mask & ~self.flood_outside(mask) - hl = np.zeros(mask.shape, np.int32) - big_holes = [] - for j in range(mask.shape[0]): - for i in range(mask.shape[1]): - if holes[j, i] and hl[j, i] == 0: - hl[j, i] = 1 - q = deque([(j, i)]) - pts = [(j, i)] - while q: - y, x = q.popleft() - for dy, dx in ((1,0),(-1,0),(0,1),(0,-1), - (1,1),(1,-1),(-1,1),(-1,-1)): - ny, nx = y + dy, x + dx - if (0 <= ny < mask.shape[0] and 0 <= nx < mask.shape[1] - and holes[ny, nx] and hl[ny, nx] == 0): - hl[ny, nx] = 1 - q.append((ny, nx)) - pts.append((ny, nx)) - if len(pts) >= 500: - big_holes.append(pts) + # 中央区剪除(鼻侧多余填充) + holes = ~mask & ~self.flood_outside(mask) + hl = np.zeros(mask.shape, np.int32) + big_holes = [] + for j in range(mask.shape[0]): + for i in range(mask.shape[1]): + if holes[j, i] and hl[j, i] == 0: + hl[j, i] = 1 + q = deque([(j, i)]) + pts = [(j, i)] + while q: + y, x = q.popleft() + for dy, dx in ((1,0),(-1,0),(0,1),(0,-1), + (1,1),(1,-1),(-1,1),(-1,-1)): + ny, nx = y + dy, x + dx + if (0 <= ny < mask.shape[0] and 0 <= nx < mask.shape[1] + and holes[ny, nx] and hl[ny, nx] == 0): + hl[ny, nx] = 1 + q.append((ny, nx)) + pts.append((ny, nx)) + if len(pts) >= 500: + big_holes.append(pts) - if len(big_holes) >= 2: - big_holes.sort(key=len, reverse=True) - h_l, h_r = sorted(big_holes[:2], - key=lambda pts: sum(p[1] for p in pts) / len(pts)) - bl = (min(p[1] for p in h_l), max(p[1] for p in h_l), min(p[0] for p in h_l)) - br = (min(p[1] for p in h_r), max(p[1] for p in h_r), min(p[0] for p in h_r)) - x_in_l, x_in_r = bl[1], br[0] - y_cut = min(bl[2], br[2]) - if x_in_l < x_in_r: - cut = np.zeros_like(mask, bool) - cut[:y_cut, x_in_l + 1:x_in_r] = True - n_cut = int((mask & cut).sum()) - mask &= ~cut - self.log(f" 中央区剪除(两镜片开口之间, 开口底缘 y>{self.Y0 + y_cut * self.RES:.1f} 以下): " - f"剔除 {n_cut} 格误标鼻面") - else: - self.log(" ! 未找到两个镜片开口, 跳过中央区剪除") + if len(big_holes) >= 2: + big_holes.sort(key=len, reverse=True) + h_l, h_r = sorted(big_holes[:2], + key=lambda pts: sum(p[1] for p in pts) / len(pts)) + bl = (min(p[1] for p in h_l), max(p[1] for p in h_l), min(p[0] for p in h_l)) + br = (min(p[1] for p in h_r), max(p[1] for p in h_r), min(p[0] for p in h_r)) + x_in_l, x_in_r = bl[1], br[0] + y_cut = min(bl[2], br[2]) + if x_in_l < x_in_r: + cut = np.zeros_like(mask, bool) + cut[:y_cut, x_in_l + 1:x_in_r] = True + n_cut = int((mask & cut).sum()) + mask &= ~cut + self.log(f" 中央区剪除(两镜片开口之间, 开口底缘 y>{self.Y0 + y_cut * self.RES:.1f} 以下): " + f"剔除 {n_cut} 格误标鼻面") + + # ===== 新增: 创建 lens_mask ===== + self.lens_mask = np.zeros(mask.shape, bool) + for pts in big_holes[:2]: + for p in pts: + self.lens_mask[p[0], p[1]] = True + for _ in range(2): # 腐蚀 2 格 + _nm = self.lens_mask.copy() + _nm[1:, :] &= self.lens_mask[:-1, :]; _nm[:-1, :] &= self.lens_mask[1:, :] + _nm[:, 1:] &= self.lens_mask[:, :-1]; _nm[:, :-1] &= self.lens_mask[:, 1:] + self.lens_mask = _nm + # ================================ + else: + self.log(" ! 未找到两个镜片开口, 跳过中央区剪除") - # 膨胀 3 格 - mask = self.dilate_8(mask, 3) + # 膨胀 3 格 + mask = self.dilate_8(mask, 3) - # 补内部洞 (镜片开口) - holes = self.flood_outside(mask) - self.log(f" 足迹 {mask.sum()} 格, 补内部洞 {holes.sum()} 格 (镜片开口)") - mask = mask | holes + # 补内部洞 (镜片开口) + holes = self.flood_outside(mask) + self.log(f" 足迹 {mask.sum()} 格, 补内部洞 {holes.sum()} 格 (镜片开口)") + mask = mask | holes - # 腐蚀 - eroded = self.erode_4(mask, self.ERODE) + # 腐蚀 + eroded = self.erode_4(mask, self.ERODE) - # y 裁剪 + 最大连通块 - H, W = eroded.shape - jj = np.arange(H) - eroded = eroded & (jj[:, None] >= int((self.RY0 - self.Y0) / self.RES)) & \ - (jj[:, None] <= int((self.RY1 - self.Y0) / self.RES)) + # y 裁剪 + 最大连通块 + H, W = eroded.shape + jj = np.arange(H) + eroded = eroded & (jj[:, None] >= int((self.RY0 - self.Y0) / self.RES)) & \ + (jj[:, None] <= int((self.RY1 - self.Y0) / self.RES)) - R = np.zeros_like(eroded) - seen = np.zeros_like(eroded, bool) - best = None - for j0 in range(H): - for i0 in range(W): - if eroded[j0, i0] and not seen[j0, i0]: - comp = [] - dq = deque([(j0, i0)]) - seen[j0, i0] = True - touch = False - while dq: - j, i = dq.popleft() - comp.append((j, i)) - if i in (0, W - 1) or j in (0, H - 1): - touch = True - for dj, di in ((1,0),(-1,0),(0,1),(0,-1)): - nj, ni = j + dj, i + di - if 0 <= nj < H and 0 <= ni < W and eroded[nj, ni] and not seen[nj, ni]: - seen[nj, ni] = True - dq.append((nj, ni)) - if best is None or len(comp) > len(best[0]): - best = (comp, touch) + R = np.zeros_like(eroded) + seen = np.zeros_like(eroded, bool) + best = None + for j0 in range(H): + for i0 in range(W): + if eroded[j0, i0] and not seen[j0, i0]: + comp = [] + dq = deque([(j0, i0)]) + seen[j0, i0] = True + touch = False + while dq: + j, i = dq.popleft() + comp.append((j, i)) + if i in (0, W - 1) or j in (0, H - 1): + touch = True + for dj, di in ((1,0),(-1,0),(0,1),(0,-1)): + nj, ni = j + dj, i + di + if 0 <= nj < H and 0 <= ni < W and eroded[nj, ni] and not seen[nj, ni]: + seen[nj, ni] = True + dq.append((nj, ni)) + if best is None or len(comp) > len(best[0]): + best = (comp, touch) - if best is None: - raise SystemExit("R 为空: 足迹检测失败") - comp, touch = best - for j, i in comp: - R[j, i] = True - self.log(f" R: {R.sum()} 格 ({R.sum()*self.RES*self.RES:.2f} 面积单位), " - f"触碰窗口边界={touch} (触碰=侧墙切在头内部, 正常)") - if R.sum() < 200: - raise SystemExit("R 过小") + if best is None: + raise SystemExit("R 为空: 足迹检测失败") + comp, touch = best + for j, i in comp: + R[j, i] = True + self.log(f" R: {R.sum()} 格 ({R.sum()*self.RES*self.RES:.2f} 面积单位), " + f"触碰窗口边界={touch} (触碰=侧墙切在头内部, 正常)") - return mask, R + # ===== [抽取自原版] R 裁剪到眼镜真实轮廓 ===== + _sil = g1f.copy() + for _ in range(6): + _sm = _sil.copy() + _sm[1:, :] |= _sil[:-1, :]; _sm[:-1, :] |= _sil[1:, :] + _sm[:, 1:] |= _sil[:, :-1]; _sm[:, :-1] |= _sil[:, 1:] + _sm[1:, 1:] |= _sil[:-1, :-1]; _sm[:-1, :-1] |= _sil[1:, 1:] + _sm[1:, :-1] |= _sil[:-1, 1:]; _sm[:-1, 1:] |= _sil[1:, :-1] + _sil = _sm + _sil = _sil | self.flood_outside(_sil) + for _ in range(6): + _sm = _sil.copy() + _sm[1:, :] &= _sil[:-1, :]; _sm[:-1, :] &= _sil[1:, :] + _sm[:, 1:] &= _sil[:, :-1]; _sm[:, :-1] &= _sil[:, 1:] + _sm[1:, 1:] &= _sil[:-1, :-1]; _sm[:-1, :-1] &= _sil[1:, 1:] + _sm[1:, :-1] &= _sil[:-1, 1:]; _sm[:-1, 1:] &= _sil[1:, :-1] + _sil = _sm + if self.lens_mask is not None: + _lexr = self.lens_mask.copy() + for _ in range(2): + _sm = _lexr.copy() + _sm[1:, :] |= _lexr[:-1, :]; _sm[:-1, :] |= _lexr[1:, :] + _sm[:, 1:] |= _lexr[:, :-1]; _sm[:, :-1] |= _lexr[:, 1:] + _sm[1:, 1:] |= _lexr[:-1, :-1]; _sm[:-1, :-1] |= _lexr[1:, 1:] + _sm[1:, :-1] |= _lexr[:-1, 1:]; _sm[:-1, 1:] |= _lexr[1:, :-1] + _lexr = _sm + _sil = _sil | _lexr + _nclip = int((R & ~_sil).sum()) + R &= _sil + self.log(f" R 裁剪到眼镜轮廓: 剔除轮廓外 {_nclip} 格, R={R.sum()} 格") + # =========================================== + + if R.sum() < 200: + raise SystemExit("R 过小") + + return mask, R # ====================================================================== # 步骤 4: 构造填充前/后高度场 @@ -644,6 +688,7 @@ class FillMeshGenerator: hb = np.minimum(hb, hf - self.GAP_MIN) hb = self.smooth_valid(hb, R, sigma=2.0, it=2) + # 外扩填充 NaN for _ in range(60): nanm = ~np.isfinite(hf) diff --git a/fill2sc/vacuum_wrap.py b/fill2sc/vacuum_wrap.py new file mode 100644 index 0000000..45f90b0 --- /dev/null +++ b/fill2sc/vacuum_wrap.py @@ -0,0 +1,1046 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +vacuum_wrap.py —— 眼镜“抽真空包装” mesh 生成脚本 +================================================ +参考 youshui_result.obj 中新增的 o Transparent 部分(一个完全包裹眼镜外轮廓、 +并填充镜片镂空处的水密闭合 mesh), 本脚本从 youshui.obj + index.txt 出发重新生成它。 + +输入/输出(自动探测, 也可用 --obj/--index/--out 显式指定): + 模型 本目录最大的非输出 *.obj(排除 *result*/*wrap*/*_fill*) + 索引 模型同目录的 index.txt(眼镜顶点索引, 每行一个, 0-based) + 输出 模型同目录 <模型名>_wrap_result.obj = 原模型 + 追加 o Transparent + 处理其它数据集时, 把本脚本拷到该数据目录下直接运行即可。 + +思路(“抽真空包装”): + 1. 由 index.txt 提取眼镜子网格(三个顶点都在索引里的面); + 2. 判定正面视线轴后, 只取“朝前的面”(法线与视线轴余弦 > 阈值)投影 + 栅格化, 形态学闭运算补缝, 取含镜片镂空的最大连通组件, + 内部封闭镂空自动填满 → 得到镜框正面“外轮廓”多边形 + (镜片镂空被盖住, 鼻梁处的凹口保留; 向后弯折的镜腿被排除, + 包装不超出眼镜正面外框); + 3. 在外轮廓内铺一层网格点, 作为“包装薄膜”的采样点; + 4. 每个网格点从正前方贴到眼镜表面(取该处最靠前的表面高度); + 镜片镂空处没有表面 → 用周边镜框高度拟合二次趋势面外推 + (跟随镜框倾斜/缓弯, 平齐不下陷), 残差在孔内做调和衔接, + 平滑无台阶; footprint 按镜架整体投影剪裁, 不超出镜架; + 5. 薄膜沿视线后退 thickness 得到背面膜(最深不超过镜架最后缘), 周边侧壁缝合 + → 一个水密(watertight)闭合的包装 mesh; + 6. 棋盘式交替对角线剖分; 薄膜边缘沿边界曲线做一维拉普拉斯收边, + 磨掉 footprint 栅格化的锯齿轮廓; + 7. 原样拷贝原 obj(仅修正 mtllib/usemtl 指向本目录的 youshui.mtl), + 把新 mesh 追加在文件末尾(o Transparent / usemtl Transparent)。 + +只依赖 numpy; 所有读写都限制在本脚本所在目录。 + +与参考结果 youshui_result.obj 的对比: 结构一致(原模型 + o Transparent 水密闭合 +包装 mesh); 本脚本只包装镜框正面核心(含镜片镂空填充), 不包裹向后弯折的镜腿, +因此与参考一样不超出眼镜正面外框、z 向深度一致; 参考在镜腿后另有一小片 +手工处理的收拢薄膜, 此处不复现。 + +用法: + python3 vacuum_wrap.py + python3 vacuum_wrap.py --spacing 0.06 --thickness 3.4 --cell 0.05 --front 0.25 +""" + +import os +import glob +import argparse +import numpy as np + +SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) + + +def resolve(p): + """路径一律相对本脚本目录解析, 不碰其它目录。""" + return p if os.path.isabs(p) else os.path.join(SCRIPT_DIR, p) + + +# ----------------------------------------------------------------- OBJ 读取 +def load_obj(path): + """读取顶点与三角面(多于三个顶点的面做三角扇拆分)。""" + verts, faces = [], [] + with open(path, 'r', errors='ignore') as f: + for line in f: + if line.startswith('v '): + p = line.split() + verts.append((float(p[1]), float(p[2]), float(p[3]))) + elif line.startswith('f '): + ids = [int(x.split('/')[0]) - 1 for x in line.split()[1:]] + for k in range(1, len(ids) - 1): + faces.append((ids[0], ids[k], ids[k + 1])) + return np.asarray(verts, np.float64), np.asarray(faces, np.int64) + + +def load_index(path, nverts): + idx = [] + with open(path) as f: + for line in f: + line = line.strip() + if line: + idx.append(int(line)) + idx = np.asarray(idx, np.int64) + if len(idx) == 0: + raise RuntimeError('index.txt 为空') + if idx.min() < 0 or idx.max() >= nverts: + raise RuntimeError('index.txt 中索引越界(0..%d)' % (nverts - 1)) + return idx + + +# ------------------------------------------------------------ 正面朝向判定 +def pick_front_axis(verts, sub, wind, coarse): + """判定“抽真空”的视线轴: 沿该轴投影时, 眼镜区域内部会出现封闭的 + 镜片镂空(孔洞), 孔洞面积最大的轴即包装轴(如从正面看眼镜)。 + wind: 整mesh体积符号(+1=法线朝外), 用于确定正面在轴的正/负侧。 + 返回 (k, sgn, u轴, v轴), (u, v, w=sgn*第k轴) 构成右手系。""" + p0, p1, p2 = verts[sub[:, 0]], verts[sub[:, 1]], verts[sub[:, 2]] + n = np.cross(p1 - p0, p2 - p0) + sgn_area = wind * n.sum(axis=0) # 各轴带符号投影面积 + best_hole, best_k = -1.0, 2 + rc = max(2, int(round(0.6 / coarse))) # 较大闭运算: 桥接处覆盖稀疏时也能封住镂空 + for k in range(3): + u, v = (k + 1) % 3, (k + 2) % 3 + tri_uv = verts[sub][:, :, [u, v]] + tri_w = verts[sub][:, :, k] + cover, _, _ = rasterize(tri_uv, tri_w, coarse) + foot = fill_enclosed(binary_closing(cover, rc)) + hole = float((foot & ~cover).sum()) * coarse * coarse + print(' 投影轴 %s: 表面覆盖 %.1f, 封闭镂空 %.1f' + % ('xyz'[k], cover.sum() * coarse * coarse, hole)) + if hole > best_hole: + best_hole, best_k = hole, k + k = best_k + sgn = 1.0 if sgn_area[k] >= 0 else -1.0 + return k, sgn, (k + 1) % 3, (k + 2) % 3 + + +# ------------------------------------------------------------------ 栅格化 +def rasterize(tri_uv, tri_w, h, keep=None, collect=False): + """把子网格三角形投影到 (u,v) 平面栅格化。 + keep: 可选的三角形布尔掩码(如只保留朝前的面)。 + collect: 额外收集 每格覆盖的三角形编号表(用于后续精确求高)。 + 返回 cover(是否被覆盖) 与 topz(该格内最靠前的表面高度 w); + collect=True 时再附带 cell_tris。""" + umin, vmin = tri_uv[..., 0].min(), tri_uv[..., 1].min() + umax, vmax = tri_uv[..., 0].max(), tri_uv[..., 1].max() + # 留出边距(>= 闭运算半径+2), 保证膨胀/洪泛不会越过网格边界 + mg = 8 + umin -= mg * h; vmin -= mg * h + nu = int(np.ceil((umax - umin) / h)) + 1 + mg + nv = int(np.ceil((vmax - vmin) / h)) + 1 + mg + if keep is not None: # 网格 bbox 一律取自全集, 保证可对齐 + tri_uv = tri_uv[keep] + tri_w = tri_w[keep] + nu -= nu % 2; nv -= nv % 2 # 偶数尺寸, 红黑迭代对齐 + cover = np.zeros((nu, nv), bool) + topz = np.full((nu, nv), np.nan) + cell_tris = {} if collect else None + ntri = tri_uv.shape[0] + for t in range(ntri): + a = tri_uv[t, 0]; b = tri_uv[t, 1]; c = tri_uv[t, 2] + wa, wb, wc = tri_w[t] + i0 = max(0, int(np.floor((min(a[0], b[0], c[0]) - umin) / h)) - 1) + i1 = min(nu - 1, int(np.floor((max(a[0], b[0], c[0]) - umin) / h)) + 1) + j0 = max(0, int(np.floor((min(a[1], b[1], c[1]) - vmin) / h)) - 1) + j1 = min(nv - 1, int(np.floor((max(a[1], b[1], c[1]) - vmin) / h)) + 1) + if i1 < i0 or j1 < j0: + continue + cx = umin + (np.arange(i0, i1 + 1) + 0.5) * h + cy = vmin + (np.arange(j0, j1 + 1) + 0.5) * h + CX, CY = np.meshgrid(cx, cy, indexing='ij') + v0 = c - a + v1 = b - a + d00 = v0 @ v0; d01 = v0 @ v1; d11 = v1 @ v1 + den = d00 * d11 - d01 * d01 + if abs(den) < 1e-14: # 投影退化(近似垂直面), 跳过 + continue + px, py = CX - a[0], CY - a[1] + d20 = px * v0[0] + py * v0[1] + d21 = px * v1[0] + py * v1[1] + bu = (d11 * d20 - d01 * d21) / den + bv = (d00 * d21 - d01 * d20) / den + ins = (bu >= -1e-9) & (bv >= -1e-9) & (bu + bv <= 1 + 1e-9) + if not ins.any(): + continue + z = wa + bu * (wb - wa) + bv * (wc - wa) + sub = topz[i0:i1 + 1, j0:j1 + 1] + topz[i0:i1 + 1, j0:j1 + 1] = np.where(ins, np.fmax(sub, z), sub) + cover[i0:i1 + 1, j0:j1 + 1] |= ins + if collect: + ii, jj = np.nonzero(ins) + for ci_, cj_ in zip(ii, jj): + cell_tris.setdefault((int(ci_) + i0, int(cj_) + j0), []).append(t) + if collect: + return cover, topz, (umin, vmin), cell_tris + return cover, topz, (umin, vmin) + + +# ------------------------------------------------------- 精确薄膜高度 +def exact_film_heights(pts, tri_uv, tri_w, cell_tris, u0, v0, h): + """每个薄膜网格点沿视线轴与子网格求精确交点, 取最靠前的高度。 + 相比栅格最近格值, 没有分段常数的阶梯量化, 薄膜贴合更光滑。""" + n = len(pts) + Z = np.full(n, np.nan) + for p in range(n): + x, y = pts[p, 0], pts[p, 1] + i = int((x - u0) / h) + j = int((y - v0) / h) + ids = cell_tris.get((i, j)) + if not ids: + continue + best = -np.inf + for t in ids: + a = tri_uv[t, 0]; b = tri_uv[t, 1]; c = tri_uv[t, 2] + e0 = c - a + e1 = b - a + d00 = e0 @ e0; d01 = e0 @ e1; d11 = e1 @ e1 + den = d00 * d11 - d01 * d01 + if abs(den) < 1e-14: + continue + d20 = (x - a[0]) * e0[0] + (y - a[1]) * e0[1] + d21 = (x - a[0]) * e1[0] + (y - a[1]) * e1[1] + bu = (d11 * d20 - d01 * d21) / den + bv = (d00 * d21 - d01 * d20) / den + if bu >= -1e-9 and bv >= -1e-9 and bu + bv <= 1 + 1e-9: + wa, wb, wc = tri_w[t] + z = wa + bu * (wb - wa) + bv * (wc - wa) + if z > best: + best = z + if best > -np.inf: + Z[p] = best + return Z + + +# ---------------------------------------------------------------- 形态学 +def dilate(m): + r = m.copy() + r[1:, :] |= m[:-1, :] + r[:-1, :] |= m[1:, :] + r[:, 1:] |= m[:, :-1] + r[:, :-1] |= m[:, 1:] + return r + + +def erode(m): + r = m.copy() + r[1:, :] &= m[:-1, :] + r[:-1, :] &= m[1:, :] + r[:, 1:] &= m[:, :-1] + r[:, :-1] &= m[:, 1:] + return r + + +def binary_closing(mask, r): + if r <= 0: + return mask.copy() + m = np.pad(mask, r, constant_values=False) # pad 足够大, 膨胀不会越界 + for _ in range(r): + m = dilate(m) + for _ in range(r): + m = erode(m) + return m[r:-r, r:-r] + + +def fill_enclosed(mask): + """从外部洪泛, 补上内部封闭的镂空(镜片孔洞), 开向边界的凹口(鼻梁)保留。""" + outside = np.zeros_like(mask) + outside[0, :] = ~mask[0, :] + outside[-1, :] = ~mask[-1, :] + outside[:, 0] = ~mask[:, 0] + outside[:, -1] = ~mask[:, -1] + while True: + new = dilate(outside) & ~mask + if np.array_equal(new, outside): + break + outside = new + return ~outside + + +def largest_filled_component(mask): + """先补内部封闭孔洞, 再取面积最大的连通组件(含镜片镂空的镜框核心), + 丢掉镜腿等散碎覆盖 → 包装不会越过镜框正面外框。""" + filled = fill_enclosed(mask) + BIG = filled.size + 1 + lab = np.where(filled, np.arange(filled.size).reshape(filled.shape), BIG) + while True: # 4邻域最小标签传播 → 连通组件标记 + new = lab.copy() + for (s1, s2) in (((slice(1, None), slice(None)), (slice(None, -1), slice(None))), + ((slice(None, -1), slice(None)), (slice(1, None), slice(None))), + ((slice(None), slice(1, None)), (slice(None), slice(None, -1))), + ((slice(None), slice(None, -1)), (slice(None), slice(1, None)))): + cand = np.full_like(lab, BIG) + cand[s1] = lab[s2] + cand[~filled] = BIG # 标签只在组件内传播, 不漏进空格 + np.minimum(new, cand, out=new) + if np.array_equal(new, lab): + break + lab = new + flat = lab.ravel() + labels = flat[flat < BIG] + if len(labels) == 0: + raise RuntimeError('没有可用连通组件') + best_lab, best_area = None, -1 + for L in np.unique(labels): + a = int((flat == L).sum()) + if a > best_area: + best_area, best_lab = a, L + return lab == best_lab + + +# ----------------------------------------------------------- 外轮廓追踪 +def trace_all_cycles(mask): + """沿 cell 边缘追踪 mask 的所有边界环(内部在行进方向左侧)。""" + nu, nv = mask.shape + edges = {} + + def add(a, b): + edges.setdefault(a, []).append(b) + + I, J = np.nonzero(mask) + for i, j in zip(I, J): + if i == 0 or not mask[i - 1, j]: # 左侧外部 → 向 -v 走 + add((i, j + 1), (i, j)) + if i == nu - 1 or not mask[i + 1, j]: # 右侧外部 → 向 +v 走 + add((i + 1, j), (i + 1, j + 1)) + if j == 0 or not mask[i, j - 1]: # 下侧外部 → 向 +u 走 + add((i, j), (i + 1, j)) + if j == nv - 1 or not mask[i, j + 1]: # 上侧外部 → 向 -u 走 + add((i + 1, j + 1), (i, j + 1)) + + cycles, visited = [], set() + for start in edges: + if start in visited: + continue + cyc = [start]; visited.add(start); cur = start + while True: + nxts = [e for e in edges.get(cur, ()) if e not in visited] + if not nxts: + break + nxt = nxts[0] + cyc.append(nxt); visited.add(nxt); cur = nxt + cycles.append(cyc) + return cycles + + +def signed_area(cyc): + a = 0.0 + for k in range(len(cyc)): + x1, y1 = cyc[k]; x2, y2 = cyc[(k + 1) % len(cyc)] + a += x1 * y2 - x2 * y1 + return a / 2.0 + + +def trace_outer_cycle(mask): + """返回 (最大面积的边界环, 显著环个数)。 + 环的方向统一为逆时针(内部在左侧)。""" + cycles = trace_all_cycles(mask) + if not cycles: + raise RuntimeError('外轮廓追踪失败: mask 无边界') + cycles.sort(key=lambda c: abs(signed_area(c)), reverse=True) + main = cycles[0] + if signed_area(main) < 0: + main = main[::-1] + ref = abs(signed_area(main)) + n_sig = 1 + sum(1 for c in cycles[1:] + if abs(signed_area(c)) > 0.02 * ref) + return main, n_sig + + +# -------------------------------------------------------------- 点在多边形内 +def points_in_poly(pts, poly): + """even-odd 规则, pts:(N,2), poly:(K,2) 隐式闭合。""" + x, y = pts[:, 0], pts[:, 1] + inside = np.zeros(len(pts), bool) + x1, y1 = poly[:, 0], poly[:, 1] + x2, y2 = np.roll(x1, -1), np.roll(y1, -1) + with np.errstate(divide='ignore', invalid='ignore'): + for k in range(len(poly)): + cond = (y1[k] > y) != (y2[k] > y) + xint = (x2[k] - x1[k]) * (y - y1[k]) / (y2[k] - y1[k]) + x1[k] + inside ^= cond & (x < xint) + return inside + + +# --------------------------------------------- 镂空处高度场(与镜框平齐) +def _neighbor_avg4(Z, ok_mask): + """4 邻域均值: 只统计 ok_mask 为真且值非 nan 的格子。 + 返回 (均值, 计数)。""" + acc = np.zeros(Z.shape) + cnt = np.zeros(Z.shape) + nanZ = np.isnan(Z) + for sl_src, sl_dst in ((slice(None, -1), slice(1, None)), + (slice(1, None), slice(None, -1))): + have = ok_mask[sl_src, :] & ~nanZ[sl_src, :] + acc[sl_dst, :] += np.where(have, Z[sl_src, :], 0.0) + cnt[sl_dst, :] += have + have = ok_mask[:, sl_src] & ~nanZ[:, sl_src] + acc[:, sl_dst] += np.where(have, Z[:, sl_src], 0.0) + cnt[:, sl_dst] += have + return np.where(cnt > 0, acc / np.maximum(cnt, 1), np.nan), cnt + + +def _fit_surface(di, dj, z): + """最小二乘拟合 z ≈ a + b·di + c·dj (平面趋势)。 + + 采样点只分布在镂空周围的环带上: 二次趋势面对这种采样是病态 + 外推(环带上常数项与二次项近似共线, 中心值可飘出 ±2), 平面 + 趋势则稳定; 曲率起伏留给残差调和项(边界值插值, 良定)。 + 点太少/退化时回退到常数。返回 (系数, 项数)。""" + n = len(z) + di = di.astype(np.float64) + dj = dj.astype(np.float64) + if n >= 7: + A = np.stack([np.ones(n), di, dj], 1) + G = A.T @ A + G[np.diag_indices_from(G)] += 1e-6 * max(1.0, G.diagonal().max()) + try: + return np.linalg.solve(G, A.T @ z), 3 + except np.linalg.LinAlgError: + pass + return np.array([float(np.mean(z))]), 1 + + +def _eval_surface(coef, nterm, di, dj): + di = np.asarray(di, np.float64) + dj = np.asarray(dj, np.float64) + if nterm == 3: + return coef[0] + coef[1] * di + coef[2] * dj + return np.full(di.shape, coef[0], float) + + +def trend_fill(Z, valid, select=None, blend_passes=500, omega=0.8): + """镜片镂空处高度场填充: 全局平面趋势 + 残差调和, 平齐且平滑。 + + 薄膜要"平齐"(跟随镜框正面的倾斜)也要平滑(无台阶无尖峰): + 1) 用正面表面(select)上的贴合点拟合平面趋势 T(镜框正面倾角)。 + 采样点铺满整片 footprint, 平面拟合稳定(二次面在只环绕 + 镜片镂空的环带采样下是病态外推, 中心值可飘出 ±2); + 2) 残差 R = 贴合高度 - T 在**所有**贴合点锚定(含外缘陡边/ + 镜托/镜腿根处的深锚点——薄膜要包住它们), 镂空处解 + Laplace 方程: 调和解无内部极值、不会过冲, 天然平滑; + 镜片中央 ≈ 趋势面(平齐), 孔缘精确衔接, 深锚点周围形成 + 光滑下凹而不是尖 crater。 + blend_passes 为残差调和迭代次数; 迭代用欠松弛(ω≤1)—— + 超松弛会放大棋盘模式留下二阶差分尖峰。""" + Z = Z.copy() + unknown = valid & np.isnan(Z) + if not unknown.any(): + return Z + known = valid & ~np.isnan(Z) + # 1) 平面趋势: 只用正面表面高度, 稳健剔除离群锚点 + K = known if select is None else (known & select) + if K.sum() >= 7: + ki, kj = np.nonzero(K) + for _ in range(2): + coef, nterm = _fit_surface(ki - ki.mean(), kj - kj.mean(), + Z[ki, kj]) + if len(ki) < 16: + break + res = Z[ki, kj] - _eval_surface( + coef, nterm, ki - ki.mean(), kj - kj.mean()) + keep = np.abs(res) <= max(0.12, 4.0 * np.median(np.abs(res))) + if keep.all() or keep.sum() < 8: + break + ki, kj = ki[keep], kj[keep] + ci, cj = ki.mean(), kj.mean() + elif K.sum() > 0: + coef, nterm = np.array([float(np.mean(Z[K]))]), 1 + ci = cj = 0.0 + else: + coef, nterm = np.array([float(np.nanmean(Z[known]))]), 1 + ci = cj = 0.0 + ii, jj = np.indices(Z.shape) + T = _eval_surface(coef, nterm, ii - ci, jj - cj) + # 2) 残差调和填充: 全体镂空处一起解, 锚定所有贴合点 + R = np.where(known, Z - T, 0.0) + for _ in range(blend_passes): + avg, cnt = _neighbor_avg4(R, valid) + upd = unknown & (cnt > 0) + if not upd.any(): + break + R[upd] += omega * (avg[upd] - R[upd]) + Z[unknown] = T[unknown] + R[unknown] + return Z + + +# ------------------------------------------------------------ 贴合带平滑 +def smooth_adherence(Z, known, valid, thr, spacing, + alpha=0.6, passes=60, near=0.3, near_w=0.4): + """贴合带受限平滑(消除镜框陡坡处高度场的棱面)。 + + 高度场 z(u,v) 表示坡度为 s 的表面时, 其曲率会被放大 ~s^2 倍, + 镜框侧弯坡度大处逐格贴合会出现 30-60 度的相邻面法线跳变。 + 本函数把与 8 邻域均值偏差超过 thr 的点向邻域均值加权松弛: + - 贴合点(known) 权重 1: 磨掉曲率突变棱面, 薄膜在陡坡处"绷而不贴"; + 坡度平缓处/均匀斜面的偏差本来就小于 thr, 不会被改动; + - 绷紧区离贴合边界 near 以内的点权重 near_w: 轻推, 软化接缝褶皱; + - 更远处不参与(保持真空 Laplace 解)。 + thr 的量纲是高度: 相邻面法线跳变约 atan(4*thr/spacing)。""" + band = known.copy() + for _ in range(max(1, int(round(near / spacing)))): + band = dilate(band) + w = np.where(known, 1.0, np.where(band & valid, near_w, 0.0)) + Z = Z.copy() + ni, nj = Z.shape + for _ in range(passes): + Zp = np.pad(Z, 1, constant_values=np.nan) + Vp = np.pad(valid, 1, constant_values=False) + acc = np.zeros(Z.shape) + cnt = np.zeros(Z.shape) + for di in (-1, 0, 1): + for dj in (-1, 0, 1): + if di == 0 and dj == 0: + continue + sl = (slice(1 + di, 1 + di + ni), slice(1 + dj, 1 + dj + nj)) + nb = np.where(Vp[sl], Zp[sl], np.nan) + ok = ~np.isnan(nb) + acc += np.where(ok, nb, 0.0) + cnt += ok + avg = np.where(cnt > 0, acc / np.maximum(cnt, 1), np.nan) + res = np.abs(Z - avg) + upd = (w > 0) & (cnt > 0) & (res > thr) + if not upd.any(): + break + Z[upd] += w[upd] * alpha * (avg[upd] - Z[upd]) + return Z, w + + +# ------------------------------------------------------------ 网格清理 +def drop_small_components(faces, tags=None): + """只保留面数最多的连通组件(正常情况下包装 mesh 只有唯一一层壳, + 栅格噪声/异常轮廓可能产生碎块或多壳, 一律丢弃)。 + tags: 面标签(0=顶面 1=底面 2=侧壁), 随面一同过滤; 返回 faces + 或 (faces, tags)。""" + from collections import defaultdict + parent = list(range(len(faces))) + + def find(x): + while parent[x] != x: + parent[x] = parent[parent[x]] + x = parent[x] + return x + + e2f = defaultdict(list) + for fi, f in enumerate(faces): + for x, y in ((f[0], f[1]), (f[1], f[2]), (f[2], f[0])): + e2f[(min(x, y), max(x, y))].append(fi) + for fl in e2f.values(): + for fi in fl[1:]: + parent[find(fi)] = find(fl[0]) + groups = defaultdict(int) + for fi in range(len(faces)): + groups[find(fi)] += 1 + if len(groups) <= 1: + return (faces, tags) if tags is not None else faces + best = max(groups.values()) + if best < len(faces) * 0.9: + print(' 警告: 网格有 %d 个组件, 仅保留最大者(%d/%d 面)' + % (len(groups), best, len(faces))) + keep = [fi for fi in range(len(faces)) if groups[find(fi)] == best] + if len(keep) != len(faces): + keep = np.asarray(keep) + faces = faces[keep] + if tags is not None: + tags = tags[keep] + return (faces, tags) if tags is not None else faces + + +def split_pinch_vertices(verts, faces): + """拆分“夹点”顶点(一个顶点被多片不相连的面共享, 网格对角接触造成), + 拆分后网格为严格的顶点流形。返回 (verts, faces)。""" + from collections import defaultdict + vfaces = defaultdict(list) + for fi, f in enumerate(faces): + for v in f: + vfaces[int(v)].append(fi) + # v -> [ (本组面id列表), ... ] + splits = [] + for v, fl in vfaces.items(): + if len(fl) < 2: + continue + edge_faces = defaultdict(list) + for fi in fl: + f = faces[fi] + for x, y in ((f[0], f[1]), (f[1], f[2]), (f[2], f[0])): + u = y if x == v else (x if y == v else None) + if u is not None: + edge_faces[int(u)].append(fi) + parent = {fi: fi for fi in fl} + + def find(x): + while parent[x] != x: + parent[x] = parent[parent[x]] + x = parent[x] + return x + + for fl2 in edge_faces.values(): + for fi in fl2[1:]: + parent[find(fi)] = find(fl2[0]) + groups = defaultdict(list) + for fi in fl: + groups[find(fi)].append(fi) + if len(groups) > 1: + splits.append((v, [g for g in groups.values()])) + if not splits: + return verts, faces + faces = faces.copy() + extra_verts = [] + for v, groups in splits: + for g in groups[1:]: + nv = len(verts) + len(extra_verts) + extra_verts.append(verts[v]) + for fi in g: + f = faces[fi] + faces[fi] = [nv if x == v else x for x in f] + if extra_verts: + verts = np.vstack([verts, np.asarray(extra_verts)]) + return verts, faces + + +# ------------------------------------------------------------ 网格松弛 +def smooth_closed_curve(P, passes=3, lam=0.5, mu=-0.53): + """闭合折线 Taubin(λ/μ) 低通平滑: 消除栅格楼梯而不收缩轮廓。 + 交替"向邻点均值收缩(λ)"与"反向外推(μ)", 大尺度形状(鼻梁凹口、 + 镜框圆角)保留, 只磨掉 cell 量级的锯齿。用于外轮廓多边形。""" + P = np.array(P, np.float64, copy=True) + for _ in range(passes): + for a in (lam, mu): + avg = 0.5 * (np.roll(P, 1, axis=0) + np.roll(P, -1, axis=0)) + P = P + a * (avg - P) + return P + + +def smooth_rim_curve(verts, faces, tags, passes=20, lam=0.5): + """薄膜周边收边: 沿顶膜边界曲线做 1D 拉普拉斯平滑。 + + footprint 是网格点掩码, 薄膜边缘天生带 spacing 级楼梯。本函数: + - 取顶面(tags==0)的边界边 → 每个边缘顶点恰有两个曲线邻居, + 沿这条闭合曲线反复做拉普拉斯松弛, 把楼梯磨成圆角; + 曲线点数近千、点距 ~spacing, 整体收缩量可忽略; + - 底膜边缘顶点(顶/底膜是同一网格点的上下两个副本, 组装时 + (u,v) 相同)通过侧壁找到配对, 加上与顶膜边缘相同的位移, + 侧壁整体随缘平移, 不会被剪切; + - 薄膜内部完全不被动(贴合不受影响)。 + + 注意: 不能改用对整个闭合 mesh 的拉普拉斯松弛来做收边/平滑—— + 侧壁是 spacing x thickness 的细长三角形(长宽比可达 70:1), + 伞形算子会把膜面边缘顶点拽向侧壁另一端的底膜顶点, 反而在 + 周边产生大幅 z 向锯齿(膜面顶点也会被拉伸错位)。 + tags: build_wrap_mesh 的面标签(0=顶面 1=底面 2=侧壁)。""" + from collections import Counter, defaultdict + topF = faces[tags == 0] + botF = faces[tags == 1] + wallF = faces[tags == 2] + in_bot = np.zeros(len(verts), bool) + in_bot[botF.ravel()] = True + # 顶膜边界无向边(在顶面内出现一次的边) → 曲线邻接表 + cnt = Counter() + for f in topF: + for x, y in ((f[0], f[1]), (f[1], f[2]), (f[2], f[0])): + cnt[(min(x, y), max(x, y))] += 1 + nbr = {} + for (a, b) in (e for e, c in cnt.items() if c == 1): + nbr.setdefault(int(a), []).append(int(b)) + nbr.setdefault(int(b), []).append(int(a)) + # 顶/底膜边缘配对: 与顶膜边缘顶点共享侧壁且 (u,v) 相同的底膜点 + v2wall = defaultdict(list) + for fi, f in enumerate(wallF): + for x in f: + v2wall[int(x)].append(fi) + keys = sorted(nbr.keys()) + pair = {} + for v in keys: + pv = verts[v] + for fi in v2wall.get(v, ()): + for x in wallF[fi]: + x = int(x) + if in_bot[x] and abs(verts[x, 0] - pv[0]) < 1e-9 \ + and abs(verts[x, 1] - pv[1]) < 1e-9: + pair[v] = x + # 沿曲线拉普拉斯松弛(Jacobi) + pos = {v: verts[v].copy() for v in keys} + for _ in range(passes): + new = {} + for v, P in pos.items(): + a, b = nbr[v] + new[v] = P + lam * (0.5 * (pos[a] + pos[b]) - P) + pos = new + # 应用: 顶膜边缘取新位置, 底膜配对点加相同位移(侧壁随缘平移) + for v in keys: + delta = pos[v] - verts[v] + verts[v] = pos[v] + b = pair.get(v) + if b is not None: + verts[b] += delta + return verts + + +# ------------------------------------------------------------ 包装 mesh 组装 +def build_wrap_mesh(U, V, Z, inside, thickness, w_floor=None): + """顶面膜 + 底面膜 + 侧壁 → 闭合网格。 + 顶点/面均为 (u,v,w) 坐标, 返回 (verts, faces, tags), + tags: 每个面的标签(0=顶面 1=底面 2=侧壁), 供收边平滑 + 定位薄膜边缘(周边锯齿所在的侧壁)。 + w_floor: 底面膜最深不超过该值(默认=镜架最后缘再退 0.1), + 防止薄膜沿侧影线下潜时背面膜探出镜架范围。""" + ni, nj = Z.shape + idx = -np.ones((ni, nj), np.int64) + ii, jj = np.nonzero(inside) + n = len(ii) + idx[ii, jj] = np.arange(n) + zt = Z[ii, jj] + zb = zt - thickness + if w_floor is not None: + zb = np.maximum(zb, w_floor) + verts = np.concatenate([ + np.stack([U[ii, jj], V[ii, jj], zt], axis=1), + np.stack([U[ii, jj], V[ii, jj], zb], axis=1), + ]) + # 单元格四角: a=(i,j) b=(i+1,j) c=(i+1,j+1) d=(i,j+1) + a = idx[:-1, :-1]; b = idx[1:, :-1]; c = idx[1:, 1:]; d = idx[:-1, 1:] + ok = (a >= 0) & (b >= 0) & (c >= 0) & (d >= 0) + ci, cj = np.nonzero(ok) + A, B, C, D = a[ci, cj], b[ci, cj], c[ci, cj], d[ci, cj] + # 棋盘式交替对角线: 相邻单元的剖分方向相反, 避免整片同向的机械条纹 + m0 = ((ci + cj) % 2 == 0) + topF = np.concatenate([ + np.stack([A[m0], B[m0], C[m0]], 1), np.stack([A[m0], C[m0], D[m0]], 1), + np.stack([A[~m0], B[~m0], D[~m0]], 1), np.stack([B[~m0], C[~m0], D[~m0]], 1)]) + botF = np.concatenate([ + np.stack([A[m0], C[m0], B[m0]], 1), np.stack([A[m0], D[m0], C[m0]], 1), + np.stack([A[~m0], D[~m0], B[~m0]], 1), np.stack([B[~m0], D[~m0], C[~m0]], 1)]) + n + + # 顶面的边界有向边(内部在左侧) → 缝合到背面 + # 先做无向计数找边界(内部边在相邻两面中方向相反、各出现一次), + # 再取其在面中出现的方向用于侧壁绕向 + from collections import Counter + und = Counter() + directed = [] + for f in topF: + for x, y in ((f[0], f[1]), (f[1], f[2]), (f[2], f[0])): + key = (x, y) if x < y else (y, x) + und[key] += 1 + directed.append((x, y)) + bndset = {e for e, c in und.items() if c == 1} + rim_pairs = [(p, q) for (p, q) in directed + if ((p, q) if p < q else (q, p)) in bndset] + rim = [[q, p, p + n] for (p, q) in rim_pairs] + rim += [[q, p + n, q + n] for (p, q) in rim_pairs] + faces = np.concatenate([topF, botF, np.asarray(rim, np.int64)]) + tags = np.concatenate([np.zeros(len(topF), np.int8), + np.ones(len(botF), np.int8), + np.full(len(rim), 2, np.int8)]) + + # 去掉未被任何面引用的孤立顶点 + faces, tags = drop_small_components(faces, tags) + verts, faces = split_pinch_vertices(verts, faces) + used = np.zeros(len(verts), bool) + used[faces.ravel()] = True + remap = -np.ones(len(verts), np.int64) + remap[used] = np.arange(used.sum()) + return verts[used], remap[faces], tags + + +def check_closed(faces, nverts): + """水密性与拓扑自检: 每条有向边恰好出现一次且其反向也存在 + (方向一致的水密闭合), 欧拉示性数 V - E + F = 2。""" + d = {} + for f in faces: + for x, y in ((f[0], f[1]), (f[1], f[2]), (f[2], f[0])): + d[(x, y)] = d.get((x, y), 0) + 1 + bad = [e for e, c in d.items() if c != 1 or (e[1], e[0]) not in d] + n_e = 3 * len(faces) // 2 + chi = nverts - n_e + len(faces) + return len(bad) == 0, chi + + +# ------------------------------------------------------------------- 主流程 +def default_obj_path(): + """自动探测输入模型: 取脚本所在目录下最大的 *.obj, + 排除历史输出(*_result/wrap_result/*_fill* 等)。""" + cands = [] + for p in glob.glob(os.path.join(SCRIPT_DIR, '*.obj')): + b = os.path.basename(p).lower() + if 'result' in b or 'wrap' in b or '_fill' in b: + continue + cands.append((os.path.getsize(p), p)) + if not cands: + return None + return max(cands)[1] + + +def main(): + ap = argparse.ArgumentParser(description='眼镜“抽真空包装”mesh 生成') + ap.add_argument('--obj', default=None, + help='原始模型 obj(默认: 自动探测本目录最大的非输出 obj)') + ap.add_argument('--index', default=None, + help='眼镜顶点索引(默认: 模型同目录的 index.txt)') + ap.add_argument('--out', default=None, + help='输出路径(默认: 模型同目录 <模型名>_wrap_result.obj)') + ap.add_argument('--cell', type=float, default=0.05, + help='外轮廓栅格分辨率(模型单位, 默认0.05)') + ap.add_argument('--spacing', type=float, default=0.06, + help='包装薄膜网格点间距(默认0.06)') + ap.add_argument('--thickness', type=float, default=3.4, + help='背面膜相对薄膜的后退厚度(默认3.4)') + ap.add_argument('--front', type=float, default=0.25, + help='朝前面阈值(法线与视线轴余弦, 默认0.25); ' + '调小会包含更斜的面乃至镜腿, 调大只贴更正对的面') + ap.add_argument('--closing', type=int, default=4, + help='形态学闭运算半径(格数, 默认4)') + ap.add_argument('--hole-blend', type=int, default=500, + help='镜片镂空处残差调和的迭代次数(默认500): ' + '趋势面跟随镜框平齐, 欠松弛迭代负责孔缘平滑衔接; ' + '越大越平滑') + ap.add_argument('--rim-smooth', type=int, default=20, + help='周边收边平滑次数(0=关闭, 默认20): 沿薄膜边缘曲线 ' + '拉普拉斯平滑, 磨掉轮廓栅格锯齿') + ap.add_argument('--poly-smooth', type=int, default=3, + help='外轮廓多边形 Taubin 平滑次数(0=关闭, 默认3), ' + '从源头去掉 cell 级楼梯') + ap.add_argument('--band-smooth', type=float, default=0.035, + help='贴合带平滑阈值(×spacing, 0=关闭; 磨掉陡坡棱面)') + ap.add_argument('--mtl', default=None, + help='覆盖 mtllib 行的材料文件名(默认自动找本目录 *.mtl)') + args = ap.parse_args() + + if args.obj is None: + args.obj = default_obj_path() + if args.obj is None: + raise RuntimeError('未指定 --obj, 且本目录没有可自动探测的模型文件') + obj_path = resolve(args.obj) + obj_dir = os.path.dirname(obj_path) + if args.index is None: + args.index = os.path.join(obj_dir, 'index.txt') + if args.out is None: + stem = os.path.splitext(os.path.basename(obj_path))[0] + args.out = os.path.join(obj_dir, stem + '_wrap_result.obj') + idx_path, out_path = resolve(args.index), resolve(args.out) + + print('[1/6] 读取 %s ...' % os.path.basename(obj_path)) + verts, faces = load_obj(obj_path) + print(' %d 顶点, %d 三角面' % (len(verts), len(faces))) + gidx = load_index(idx_path, len(verts)) + print(' index.txt: %d 个眼镜顶点' % len(gidx)) + + # ---- 眼镜子网格 & 正面轴向 + in_g = np.zeros(len(verts), bool) + in_g[gidx] = True + sel = in_g[faces].all(axis=1) + sub = faces[sel] + if len(sub) < 10: + raise RuntimeError('眼镜子网格面太少(%d), 请检查 index.txt' % len(sub)) + print(' 眼镜子网格: %d 面' % len(sub)) + # 整 mesh 体积符号 → 面法线是否朝外 + p0, p1, p2 = verts[faces[:, 0]], verts[faces[:, 1]], verts[faces[:, 2]] + vol = np.einsum('ij,ij->i', p0, np.cross(p1, p2)).sum() / 6.0 + wind = 1.0 if vol >= 0 else -1.0 + coarse = max(0.15, args.cell * 1.5) + print(' 判定包装视线轴(封闭镂空最大者):') + k, sgn, ax_u, ax_v = pick_front_axis(verts, sub, wind, coarse) + print(' 正面视线轴 = %s, u=轴%d(%s), v=轴%d(%s)' + % ('xyz'[k] + ('+' if sgn > 0 else '-'), + ax_u, 'xyz'[ax_u], ax_v, 'xyz'[ax_v])) + + tri_uv = verts[sub][:, :, [ax_u, ax_v]] + tri_w = sgn * verts[sub][:, :, k] + + # ---- 朝前面掩码(法线与视线轴余弦) → 只包装镜框正面, 排除向后弯折的镜腿 + sp0, sp1, sp2 = verts[sub[:, 0]], verts[sub[:, 1]], verts[sub[:, 2]] + nrm = np.cross(sp1 - sp0, sp2 - sp0) + ln = np.linalg.norm(nrm, axis=1) + front = (sgn * nrm[:, k] / np.where(ln > 1e-12, ln, 1.0)) > args.front + print(' 朝前面: %d/%d (front>%.2f)' % (front.sum(), len(sub), args.front)) + + # ---- 栅格化 → 外轮廓(闭运算补缝 → 只取含镜片镂空的镜框核心组件, + # 内部封闭孔洞填满; 半径自适应有限加大, 不会封死鼻梁凹口) + print('[2/6] 栅格化朝前表面 (cell=%.2f) ...' % args.cell) + cover, topz, (u0, v0) = rasterize(tri_uv, tri_w, args.cell, front) + cover_all, topz_all, (ua, va), cell_tris = rasterize( + tri_uv, tri_w, args.cell, collect=True) + assert (ua, va) == (u0, v0) # 两次栅格化共用同一网格 + print(' 朝前表面覆盖 %d/%d 格' % (cover.sum(), cover.size)) + r = args.closing + cap = max(args.closing + 4, int(round(1.0 / args.cell))) + cover_area = float(cover.sum()) * args.cell ** 2 + while True: + # 镜片镂空改以"全表面轮缘环"的封闭孔判定: + # 轮缘正面覆盖稀疏/有缺口时, 正面环闭不上会漏填整块镜片 + # (632587 右镜片曾缺失 1/3); 全表面环更完整, 而 开向边界的 + # 鼻梁凹口不属于封闭孔, 不会被误填 + ring = binary_closing(cover_all, r) + holes = fill_enclosed(ring) & ~ring + foot = largest_filled_component(binary_closing(cover, r) | holes) + # 上/下缘的掠射面不在正面覆盖内但属于镜架: + # 正面 footprint 附近 0.2 内的轮缘格也纳入, 包装才铺满镜框外缘 + near = foot.copy() + for _ in range(4): + near = dilate(near) + foot = largest_filled_component(foot | (ring & near)) + cycle, n_sig = trace_outer_cycle(foot) + cyc_area = abs(signed_area(cycle)) * args.cell ** 2 + foot_area = float(foot.sum()) * args.cell ** 2 + # 判据一: 外轮廓单环完整(追踪不破); + # 判据二: 镜片镂空被填上(footprint 达到朝前覆盖的 2 倍以上, + # 覆盖稀疏/有缺口的数据需要更大的闭运算半径) + if (cyc_area >= 0.9 * foot_area and foot_area >= 2.0 * cover_area) \ + or r >= cap: + break + print(' 闭运算半径 %d → 外轮廓覆盖 footprint 的 %.0f%%, ' + '镂空填充率 %.2f, 继续加大...' + % (r, 100 * cyc_area / foot_area, foot_area / cover_area)) + r += 2 + if foot_area < 1.5 * cover_area: + print(' 注意: 镜片镂空可能未被填满(填充率 %.2f), ' + '可尝试加大 --closing' % (foot_area / cover_area)) + print(' 镜框正面 footprint: %d 格(面积 %.1f), 闭运算半径 %d' + % (foot.sum(), foot_area, r)) + poly = np.array([[u0 + i * args.cell, v0 + j * args.cell] for (i, j) in cycle], + dtype=np.float64) + print(' 外轮廓多边形: %d 个角点' % len(poly)) + if args.poly_smooth > 0: + poly = smooth_closed_curve(poly, passes=args.poly_smooth) + print(' 外轮廓 Taubin 平滑 %d 次(消除 cell 级栅格楼梯, 不收缩轮廓)' + % args.poly_smooth) + + # ---- 薄膜网格点 + print('[3/6] 铺薄膜网格 (spacing=%.2f) ...' % args.spacing) + s = args.spacing + gu0 = poly[:, 0].min() - s; gu1 = poly[:, 0].max() + s + gv0 = poly[:, 1].min() - s; gv1 = poly[:, 1].max() + s + ni = int(np.ceil((gu1 - gu0) / s)) + 1 + nj = int(np.ceil((gv1 - gv0) / s)) + 1 + GU, GV = np.meshgrid(gu0 + np.arange(ni) * s, + gv0 + np.arange(nj) * s, indexing='ij') + pts = np.stack([GU.ravel(), GV.ravel()], axis=1) + ins = points_in_poly(pts, poly).reshape(ni, nj) + print(' 轮廓内网格点: %d' % ins.sum()) + if ins.sum() < 4: + raise RuntimeError('轮廓内网格点太少, 请增大 --spacing 或减小 --cell') + + # ---- 按镜架整体投影剪裁: 闭运算/轮缘带会把 footprint 向外胀出 + # 0.1~0.2 的薄裙边(镜架上缘/下缘可见"超出"), 以镜架外轮廓 + # (全表面投影 + 补采样缝 + 填充内部镂空)为准裁掉 + frame_sil = fill_enclosed(binary_closing(cover_all, 2)) + ic = np.clip(np.round((GU - u0) / args.cell).astype(int), + 0, frame_sil.shape[0] - 1) + jc = np.clip(np.round((GV - v0) / args.cell).astype(int), + 0, frame_sil.shape[1] - 1) + n_clip = int((ins & ~frame_sil[ic, jc]).sum()) + if n_clip: + ins &= frame_sil[ic, jc] + print(' 按镜架投影裁掉外轮廓之外的裙边: %d 个网格点' % n_clip) + # 趋势拟合的锚点只取正面表面(cover 膨胀 3 格)上的贴合高度: + # 外缘陡边/掠射面处薄膜本该顺侧影下潜, 高度不能拿来外推镜片平面 + cov_d = cover.copy() + for _ in range(3): + cov_d = dilate(cov_d) + front_sel = cov_d[ic, jc] + + # ---- 顶面高度: 贴表面 + 镂空处与镜框平齐 + print('[4/6] 计算薄膜高度(贴眼镜表面, 镜片处与镜框平齐) ...') + Z = exact_film_heights(pts, tri_uv, tri_w, cell_tris, + u0, v0, args.cell).reshape(ni, nj) + known = ~np.isnan(Z) & ins + n_hole = int((ins & ~known).sum()) + print(' %d 点直接贴到表面, %d 点位于镂空处(趋势平齐填充)' + % (known.sum(), n_hole)) + Z = trend_fill(Z, ins, select=front_sel, blend_passes=args.hole_blend) + + # ---- 贴合带平滑: 镜框侧弯坡度大处, 高度场曲率被坡度放大成棱面 + # (相邻面法线跳变 30-60 度), 受限松弛磨掉棱面并软化接缝; + # 浅坡处与均匀斜面本来就光滑, 不会被改动 + if args.band_smooth > 0: + thr = args.band_smooth * args.spacing + Zs, w = smooth_adherence(Z, known, ins, thr, args.spacing) + n_up = int((np.abs(Zs - Z) > 1e-6).sum()) + if n_up: + # 松弛后的点(含软化接缝带)作为新边界, 镂空处重新平齐外推 + Z = trend_fill(np.where(w > 0, Zs, np.nan), ins, + select=front_sel, blend_passes=args.hole_blend) + print(' 贴合带平滑: %d 个点被松弛(阈值 %.4f), ' + '镂空处重新平齐' % (n_up, thr)) + + # ---- 组装闭合 mesh + print('[5/6] 组装水密包装 mesh (thickness=%.2f) ...' % args.thickness) + # 底面膜深度钳制: 不超过镜架最后缘(留 0.1 薄膜厚度), + # 防止薄膜沿侧影线下潜时背面膜探出镜架范围 + w_floor = float((sgn * verts[sub][:, :, k]).min()) - 0.1 + wverts, wfaces, ftags = build_wrap_mesh(GU, GV, Z, ins, args.thickness, w_floor) + if args.rim_smooth > 0: + # 周边收边: 顶膜边缘沿边界曲线 1D 拉普拉斯, 磨掉栅格楼梯, + # 底膜边缘同步位移(侧壁随缘平移), 薄膜内部不动。 + # 注意: 不做整 mesh 的拉普拉斯松弛——侧壁细长三角形会把 + # 伞形算子变成"跨 3~4 个单位距离的强混合", 边缘反而出锯齿。 + print(' 周边收边: 边缘曲线拉普拉斯平滑 %d 次 ...' % args.rim_smooth) + wverts = smooth_rim_curve(wverts, wfaces, ftags, passes=args.rim_smooth) + closed, chi = check_closed(wfaces, len(wverts)) + print(' %d 顶点, %d 面; 水密=%s, 欧拉示性数=%d' + % (len(wverts), len(wfaces), closed, chi)) + if not closed: + raise RuntimeError('生成的网格不是水密闭合网格, 请检查参数') + + # 变换回世界坐标: (u, v, w) -> 原坐标 + W = sgn * wverts[:, 2] + verts_world = np.zeros_like(wverts) + verts_world[:, ax_u] = wverts[:, 0] + verts_world[:, ax_v] = wverts[:, 1] + verts_world[:, k] = W + + # 顶点法线(面积加权) + p0, p1, p2 = verts_world[wfaces[:, 0]], verts_world[wfaces[:, 1]], verts_world[wfaces[:, 2]] + fn = np.cross(p1 - p0, p2 - p0) + vn = np.zeros_like(verts_world) + for c in range(3): + np.add.at(vn, wfaces[:, c], fn) + ln = np.linalg.norm(vn, axis=1, keepdims=True) + vn = vn / np.where(ln > 1e-12, ln, 1.0) + + # ---- 写出 obj + print('[6/6] 写出 %s ...' % os.path.basename(out_path)) + mtl = args.mtl + if mtl is None: # 先找模型同目录的 mtl, 再退回脚本目录 + found = sorted(glob.glob(os.path.join(obj_dir, '*.mtl'))) or \ + sorted(glob.glob(os.path.join(SCRIPT_DIR, '*.mtl'))) + mtl = os.path.basename(found[0]) if found else None + n_v = n_vn = 0 + with open(out_path, 'w') as out, open(obj_path, 'r', errors='ignore') as src: + for line in src: # 原模型逐行原样拷贝 + if line.startswith('v '): + n_v += 1 + elif line.startswith('vn '): + n_vn += 1 + elif line.startswith('mtllib ') and mtl: + line = 'mtllib %s\n' % mtl # 指向本目录实际存在的 mtl + elif line.startswith('usemtl ') and line.strip().endswith('.001'): + line = 'usemtl %s\n' % line.split()[1][:-4] # Blender 重命名残留 + out.write(line) + # 追加包装 mesh + out.write('o Transparent\n') + out.write('s 0\n') + out.write('usemtl Transparent\n') + for p in verts_world: + out.write('v %.6f %.6f %.6f\n' % (p[0], p[1], p[2])) + base_v, base_n = n_v, n_vn + for nn in vn: + out.write('vn %.6f %.6f %.6f\n' % (nn[0], nn[1], nn[2])) + for f in wfaces: + out.write('f %d//%d %d//%d %d//%d\n' % ( + f[0] + 1 + base_v, f[0] + 1 + base_n, + f[1] + 1 + base_v, f[1] + 1 + base_n, + f[2] + 1 + base_v, f[2] + 1 + base_n)) + + lo, hi = W.min(), W.max() + if sgn < 0: + lo, hi = -hi, -lo + print('完成! 输出: %s' % out_path) + print(' 新增 Transparent mesh: %d 顶点 / %d 面 (参考结果为 9999/19994)' + % (len(wverts), len(wfaces))) + print(' 沿视线轴 %s 的包装范围: [%.2f, %.2f]' % ('xyz'[k], lo, hi)) + print(' 包装 xy 外接盒: x[%.2f, %.2f], y[%.2f, %.2f]' + % (verts_world[:, ax_u].min(), verts_world[:, ax_u].max(), + verts_world[:, ax_v].min(), verts_world[:, ax_v].max())) + + +if __name__ == '__main__': + main()