From 6f7805c78dbf652b72e50e9b5ecfe0f62ba0ff5c Mon Sep 17 00:00:00 2001 From: hesuicong Date: Mon, 9 Nov 2026 10:05:17 +0800 Subject: [PATCH] =?UTF-8?q?=E4=B8=8A=E4=BC=A0=E4=BB=A3=E7=A0=81?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .gitignore | 1 - fill2sc/generate_fill_mesh.py | 1177 +++++++++++++++++++++++++++++++++ 2 files changed, 1177 insertions(+), 1 deletion(-) create mode 100644 fill2sc/generate_fill_mesh.py diff --git a/.gitignore b/.gitignore index 9f8ea06..726066d 100644 --- a/.gitignore +++ b/.gitignore @@ -82,6 +82,5 @@ fill2sc/632723.jpg fill2sc/632723.mtl fill2sc/fill_preview.png fill2sc/fill_report.txt -fill2sc/generate_fill_mesh.py .gitignore fill2sc/index.txt diff --git a/fill2sc/generate_fill_mesh.py b/fill2sc/generate_fill_mesh.py new file mode 100644 index 0000000..1fa551d --- /dev/null +++ b/fill2sc/generate_fill_mesh.py @@ -0,0 +1,1177 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +generate_fill_mesh.py +===================== +为 3D 打印模型生成"眼镜 <-> 脸部"之间镂空区域的填充网格。 + +输入（与本脚本同目录）: + 632723.obj 完整模型(封闭网格) + index.txt 眼镜区域的顶点索引(每行一个, 1-based) + +输出(同目录): + fill_mesh.obj 填充网格(独立、封闭、水密) + 632723_filled.obj 原模型 + 填充 合并后的模型(可直接用于切片打印) + fill_preview.png 前视图深度预览(供人工检查) + +原理: + 1. 栅格化眼镜区域三角形 -> 眼镜足迹 mask(补洞=镜片开口, 腐蚀收边) + 2. 对足迹内每个 (x,y) 竖直射线求交: + 镜框管前/后面, 脸面前表面, 网格最后表面 + 3. 填充体 = 高度场实体: + 前表面 h_front: 贴住镜框前面(留 FLUSH 嵌入余量, 与镜框平齐), + 镜片开口区域 = 框前面的稳健二次拟合趋势面 + + 拟合残差的调和延拓(交界处严格连续, 无折痕), + 且不超过附近镜框前面(小半径扩散最大值钳制) + 后表面 h_back : 有真实脸面处 = 脸面 - 嵌入深度(填满镂空), + 无真实脸面处(误标区)不低于 hf - SHELL_MAX(填充只在镜框附近), + 且不低于后脑表面 + 4. 足迹边界提取 -> 重采样 -> 闭合 Laplacian 平滑(外边界平滑) + 5. 区域三角化: 耳切法 + 内部网格点逐点插入(拆分所在三角形, + 拓扑上不可能翻转), 前后表面 + 侧墙组成水密封闭实体 + 6. 输出前做数值校验(水密/定向/无翻转/覆盖/嵌入深度/边界平滑度) + +仅依赖 numpy 与 PIL。 +""" + +import numpy as np +from collections import deque +from PIL import Image, ImageDraw + +# ----------------------------- 参数 ----------------------------- +IN_OBJ = "632723.obj" +IN_IDX = "index.txt" +OUT_FILL = "fill_mesh.obj" +OUT_MERGED = "632723_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 裁剪(只处理正面眼镜部分) +SEP = 0.25 # 区分前后表面的最小 z 间距 +FLUSH = 0.15 # 填充前表面距镜框前面的距离(与镜框平齐, 留嵌入余量) +EMBED_BOT = 0.6 # 填充底面嵌入脸部深度 +SHELL_MAX = 3.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(补偿边界平滑收缩) +BSTEP = 0.18 # 边界重采样步长 +BITER = 25 # 边界 Laplacian 平滑迭代 +BLAM = 0.6 # 边界平滑系数 +TOL = 1e-4 # Jacobi 收敛阈值 + +report = [] +def log(s): + print(s) + report.append(str(s)) + +# ----------------------------- IO ----------------------------- +def parse_obj(path): + vs, fs = [], [] + with open(path) as f: + for line in f: + if line.startswith("v "): + p = line.split() + vs.append((float(p[1]), float(p[2]), float(p[3]))) + elif line.startswith("f "): + p = line.split() + fs.append((int(p[1].split("/")[0]) - 1, + int(p[2].split("/")[0]) - 1, + int(p[3].split("/")[0]) - 1)) + return np.asarray(vs), np.asarray(fs, dtype=np.int64) + +log("== 1. 读入模型 ==") +V, F = parse_obj(IN_OBJ) +idx = np.loadtxt(IN_IDX, dtype=np.int64) - 1 +ism = np.zeros(len(V), bool) +ism[idx] = True +face_g = ism[F].any(1) +log(f" 顶点 {len(V)}, 面 {len(F)}, 眼镜标记顶点 {len(idx)}") + +# 真眼镜物体提取: 标记面的"边连通"主分量. 眼镜与头部网格互不共边, +# 误标皮肤(鼻梁/额头/后脑, 分割噪声)在边连通下是独立小分量 -> 一次性剥离; +# 剥离后的皮肤面自然归入脸面栅格 fz(它们本来就是头部表面). +# 这样足迹/外边框直接来自真眼镜几何, 无需 2D 启发式清洗. +emap = {} +for _fi in np.where(face_g)[0]: + _a, _b, _c = F[_fi] + for _u, _v in ((_a, _b), (_b, _c), (_c, _a)): + _e = (_u, _v) if _u < _v else (_v, _u) + if _e in emap: + emap[_e].append(_fi) + else: + emap[_e] = [_fi] +_seen = set() +_best = [] +for _fi in np.where(face_g)[0]: + if _fi in _seen: + continue + _comp, _stk = [], [_fi] + _seen.add(_fi) + while _stk: + _f0 = _stk.pop() + _comp.append(_f0) + _a, _b, _c = F[_f0] + for _u, _v in ((_a, _b), (_b, _c), (_c, _a)): + _e = (_u, _v) if _u < _v else (_v, _u) + for _f2 in emap[_e]: + if _f2 not in _seen: + _seen.add(_f2) + _stk.append(_f2) + if len(_comp) > len(_best): + _best = _comp +face_g = np.zeros(len(F), bool) +face_g[_best] = True +log(f" 眼镜物体 = 边连通主分量 {len(_best)} 面, " + f"剥离误标皮肤 {int((ism[F].any(1)).sum() - len(_best))} 面") + +# ----------------------------- 栅格 ----------------------------- +W = int(round((X1 - X0) / RES)) + 1 +H = int(round((Y1 - Y0) / RES)) + 1 +log(f"== 2. 栅格化 (窗口 x[{X0},{X1}] y[{Y0},{Y1}], res={RES}, {W}x{H}) ==") + +gz1 = np.full((H, W), -np.inf) # 眼镜最前表面 +gz2 = np.full((H, W), -np.inf) # 眼镜第二表面(镜框管背面) +fz = np.full((H, W), -np.inf) # 非眼镜最前表面(脸) +bz = np.full((H, W), np.inf) # 全网格最后表面(后脑) +mvert = np.zeros((H, W), np.int32) # 眼镜顶点计数 + +yi = np.clip(((V[ism, 1] - Y0) / RES).astype(int), 0, H - 1) +xi = np.clip(((V[ism, 0] - X0) / RES).astype(int), 0, W - 1) +np.add.at(mvert, (yi, xi), 1) + +inwin = ((V[:, 0] > X0 - 0.5) & (V[:, 0] < X1 + 0.5) & + (V[:, 1] > Y0 - 0.5) & (V[:, 1] < Y1 + 0.5)) +fwin = inwin[F].all(1) +log(f" 窗口内三角形 {int(fwin.sum())} (其中眼镜 {(face_g & fwin).sum()})") + +A, B, C = V[F[:, 0]], V[F[:, 1]], V[F[:, 2]] +fi = np.where(fwin)[0] +for t in fi: + ax, ay, az = A[t]; bx, by, bz_ = B[t]; cx, cy, cz = C[t] + g = face_g[t] + i0 = max(int((min(ax, bx, cx) - X0) / RES), 0) + i1 = min(int((max(ax, bx, cx) - X0) / RES), W - 1) + j0 = max(int((min(ay, by, cy) - Y0) / RES), 0) + j1 = min(int((max(ay, by, cy) - Y0) / RES), H - 1) + if i1 < i0 or j1 < j0: + continue + v0x, v0y = cx - ax, cy - ay + v1x, v1y = bx - ax, by - ay + d00 = v0x * v0x + v0y * v0y + d01 = v0x * v1x + v0y * v1y + d11 = v1x * v1x + v1y * v1y + den = d00 * d11 - d01 * d01 + if abs(den) < 1e-14: + continue + for j in range(j0, j1 + 1): + py = Y0 + j * RES + for i in range(i0, i1 + 1): + px = X0 + i * RES + v2x, v2y = px - ax, py - ay + d02 = v0x * v2x + v0y * v2y + d12 = v1x * v2x + v1y * v2y + u = (d11 * d02 - d01 * d12) / den # 权重 C + v = (d00 * d12 - d01 * d02) / den # 权重 B + if u < 0 or v < 0 or u + v > 1: + continue + z = az + u * (cz - az) + v * (bz_ - az) + if g: + if z > gz1[j, i] + 1e-9: + if z - gz1[j, i] > SEP and gz1[j, i] > -np.inf: + gz2[j, i] = max(gz2[j, i], gz1[j, i]) + gz1[j, i] = z + elif z < gz1[j, i] - SEP and z > gz2[j, i]: + gz2[j, i] = z + else: + if z > fz[j, i]: + fz[j, i] = z + if z < bz[j, i]: + bz[j, i] = z + +have_g1 = np.isfinite(gz1) & (gz1 > -1e30) +have_g2 = np.isfinite(gz2) & (gz2 > -1e30) +have_f = np.isfinite(fz) & (fz > -1e30) +have_b = np.isfinite(bz) & (bz < 1e30) +# 眼镜表面场: face_g 已是边连通主分量(真眼镜物体), 无需再按 z 剪杂散. +# 镜圈尾部向耳朵方向包裹折返(最高面 z 从 9.1 一路降到 0.5), +# 这些格必须进入足迹/框带, 否则尾部下镂空填不上(出现缝隙) +g1f = have_g1 +log(f" 有眼镜表面格 {have_g1.sum()}, 管状(双面) {have_g2.sum()}, " + f"有脸面 {have_f.sum()}, 有背面 {have_b.sum()}") + +# ----------------------- 足迹 mask 与填充区域 R ----------------------- +log("== 3. 眼镜足迹 / 填充区域 R ==") +# 足迹 = 真眼镜物体(边连通主分量)覆盖的全部格子(含尾部包裹段). +# 误标杂散已在 face_g 阶段剥离, 无需 z 高度过滤; 旧 gz1>=7 过滤会剪掉 +# 尾部折返带(z 低于 7) -> 尾部镂空填不上 +mask = g1f +_labs = np.zeros(mask.shape, np.int32) +_cur = 0 +from collections import deque as _dq +for _j in range(mask.shape[0]): + for _i in range(mask.shape[1]): + if mask[_j, _i] and _labs[_j, _i] == 0: + _cur += 1 + _q = _dq([(_j, _i)]) + _labs[_j, _i] = _cur + 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 mask[_yy, _xx] and _labs[_yy, _xx] == 0): + _labs[_yy, _xx] = _cur + _q.append((_yy, _xx)) +_sizes = np.bincount(_labs.ravel()) +_sizes[0] = 0 +_main = int(np.argmax(_sizes)) +_nstray = int(mask.sum() - _sizes[_main]) +mask = _labs == _main +log(f" 足迹取主分量: 剔除杂散块 {_nstray} 格 (误标到镜框外的皮肤/后脑)") + +# 中央区剪除(鼻侧多余填充的根源): 镜框外边框在鼻侧 = 镜桥底缘. +# 镜片开口(足迹的大内部洞)之间的中央区里, 开口底缘以下的"眼镜覆盖" +# 是误标鼻面皮肤悬挂半岛 -> 剪掉, 填充严格按镜框外边框 +def _flood_out(msk): + out = np.zeros_like(msk, bool) + q = _dq() + for i in range(msk.shape[1]): + for j in (0, msk.shape[0] - 1): + if not msk[j, i] and not out[j, i]: + out[j, i] = True; q.append((j, i)) + for j in range(msk.shape[0]): + for i in (0, msk.shape[1] - 1): + if not msk[j, i] and not out[j, i]: + out[j, i] = True; q.append((j, i)) + while q: + y, x = q.popleft() + for dj, di in ((1,0),(-1,0),(0,1),(0,-1)): + yy, xx = y + dj, x + di + if 0 <= yy < msk.shape[0] and 0 <= xx < msk.shape[1] and not msk[yy, xx] and not out[yy, xx]: + out[yy, xx] = True; q.append((yy, xx)) + 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) +if len(_big) >= 2: + # 按 x 中心排序, 区分左右开口 + _big.sort(key=len, reverse=True) + h_l, h_r = sorted(_big[:2], key=lambda pts: sum(p[1] for p in pts) / len(pts)) + # 每个开口的 x 范围与底缘 y + def _bbox(pts): + return (min(p[1] for p in pts), max(p[1] for p in pts), min(p[0] for p in pts)) + bl, br = _bbox(h_l), _bbox(h_r) + x_in_l = bl[1] # 左开口右缘(内缘) + x_in_r = 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 + _ncut = int((mask & _cut).sum()) + mask &= ~_cut + log(f" 中央区剪除(两镜片开口之间, 开口底缘 y>{Y0 + y_cut * RES:.1f} 以下): " + f"剔除 {_ncut} 格误标鼻面") +else: + log(" ! 未找到两个镜片开口, 跳过中央区剪除") +# 膨胀 3 格(~0.3) 闭合采样缝隙, 防止镜片开口与外部连通; +# 与腐蚀 2 格(0.2) 配合净外扩 0.1: 边界平滑会压扁凸部(尾部舌形)~0.2, +# 预留外扩让平滑后的边界落在足迹外缘, 镜框下镂空填到边(真外轮廓 +# 比栅格覆盖宽~0.1-0.2, 不会超出镜框) +for _ in range(3): + nm = mask.copy() + nm[1:, :] |= mask[:-1, :]; nm[:-1, :] |= mask[1:, :] + nm[:, 1:] |= mask[:, :-1]; nm[:, :-1] |= mask[:, 1:] + nm[1:, 1:] |= mask[:-1, :-1]; nm[:-1, :-1] |= mask[1:, 1:] + nm[1:, :-1] |= mask[:-1, 1:]; nm[:-1, 1:] |= mask[1:, :-1] + mask = nm + +def flood_outside(mask): + """从窗口边界洪泛, 标记外部; 未标记且非 mask 即为洞""" + outside = np.zeros_like(mask, bool) + dq = deque() + for i in range(mask.shape[1]): + for j in (0, mask.shape[0] - 1): + if not mask[j, i] and not outside[j, i]: + outside[j, i] = True; dq.append((j, i)) + for j in range(mask.shape[0]): + for i in (0, mask.shape[1] - 1): + if not mask[j, i] and not outside[j, i]: + outside[j, i] = True; dq.append((j, i)) + while dq: + j, i = dq.popleft() + for dj, di in ((1,0),(-1,0),(0,1),(0,-1)): + nj, ni = j + dj, i + di + if 0 <= nj < mask.shape[0] and 0 <= ni < mask.shape[1]: + if not mask[nj, ni] and not outside[nj, ni]: + outside[nj, ni] = True; dq.append((nj, ni)) + return ~outside & ~mask + +holes = flood_outside(mask) +log(f" 足迹 {mask.sum()} 格, 补内部洞 {holes.sum()} 格 (镜片开口)") +mask = mask | holes + +# 腐蚀(膨胀0.3 + 腐蚀0.2 => 净外扩0.1): 预补偿边界平滑对凸部的压扁, +# 让最终边界贴足迹外缘, 镜框下镂空(尾部包裹段等)填满到边 +er = mask.copy() +for _ in range(ERODE): + ne = er.copy() + ne[1:, :] &= er[:-1, :]; ne[:-1, :] &= er[1:, :] + ne[:, 1:] &= er[:, :-1]; ne[:, :-1] &= er[:, 1:] + er = ne + +# y 裁剪 + 最大连通块 +jj = np.arange(H) +er = er & (jj[:, None] >= int((RY0 - Y0) / RES)) & (jj[:, None] <= int((RY1 - Y0) / RES)) +R = np.zeros_like(er) +seen = np.zeros_like(er, bool) +best = None +for j0 in range(H): + for i0 in range(W): + if er[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 er[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 +log(f" R: {R.sum()} 格 ({R.sum()*RES*RES:.2f} 面积单位), " + f"触碰窗口边界={touch} (触碰=侧墙切在头内部, 正常)") +if R.sum() < 200: + raise SystemExit("R 过小") + +# ----------------------- 高度场 ----------------------- +log("== 4. 构造填充前/后高度场 ==") + +def smooth_valid(Fld, valid, sigma=SIGMA, it=1): + """归一化高斯卷积(只在 valid 上)""" + rad = int(3 * sigma) + k = np.exp(-0.5 * (np.arange(-rad, rad + 1) / sigma) ** 2) + Fld = Fld.copy() + for _ in range(it): + acc = np.zeros_like(Fld); accw = np.zeros_like(Fld) + val = valid & np.isfinite(Fld) + G = np.where(val, Fld, 0.0); M = val.astype(float) + tmp = np.zeros_like(G); tmpw = np.zeros_like(G) + for d in range(-rad, rad + 1): + if d == 0: + tmp += G * k[rad]; tmpw += M * k[rad] + elif d > 0: + tmp[d:, :] += G[:-d, :] * k[rad + d]; tmpw[d:, :] += M[:-d, :] * k[rad + d] + tmp[:-d, :] += G[d:, :] * k[rad + d]; tmpw[:-d, :] += M[d:, :] * k[rad + d] + acc += tmp; accw += tmpw + tmp = np.zeros_like(G); tmpw = np.zeros_like(G) + for d in range(-rad, rad + 1): + if d == 0: + tmp += acc * k[rad]; tmpw += accw * k[rad] + elif d > 0: + tmp[:, d:] += acc[:, :-d] * k[rad + d]; tmpw[:, d:] += accw[:, :-d] * k[rad + d] + tmp[:, :-d] += acc[:, d:] * k[rad + d]; tmpw[:, :-d] += accw[:, d:] * k[rad + d] + out = np.where(tmpw > 1e-9, tmp / np.maximum(tmpw, 1e-9), Fld) + Fld = np.where(val, out, Fld) + return Fld + +def laplace_fill(Fld, known, iters=6000): + """Dirichlet: known 固定, 其余解调和方程 + 缺失邻居(网格边缘)用反射(补 F 自身), 恒除以 4 保证系数和=1(收敛)""" + F = Fld.copy() + with np.errstate(all='ignore'): + mean_known = np.nanmean(F[known]) if known.any() else 0.0 + if not np.isfinite(mean_known): + mean_known = 0.0 + # 非 known 的格子: 有限值保留, NaN/inf 用已知均值初始化 + F[~known] = np.where(np.isfinite(F[~known]), F[~known], mean_known) + for it in range(iters): + avg = np.zeros_like(F) + for dj, di in ((1,0),(-1,0),(0,1),(0,-1)): + s = np.roll(F, (dj, di), axis=(0, 1)) + ok = np.ones_like(F, bool) + if dj == 1: ok[0, :] = False + if dj == -1: ok[-1, :] = False + if di == 1: ok[:, 0] = False + if di == -1: ok[:, -1] = False + avg += np.where(ok, s, F) # 缺失邻居补 F 自身(反射) + G = np.where(known, F, avg / 4.0) # 恒除 4: 系数和 = 1 + df = np.abs(G - F)[~known] + delta = np.nanmax(df) if df.size and np.isfinite(df).any() else 0.0 + F = G + if delta < TOL * RES: + break + return F + +gx, gy = np.meshgrid(X0 + np.arange(W) * RES, Y0 + np.arange(H) * RES) + +# 补 g1f 栅格化小洞(陡峭角落的覆盖缺口): 闭运算(膨胀2再腐蚀2) +g1c = g1f.copy() +for _ in range(2): + gd = g1c.copy() + for dj, di in ((1,0),(-1,0),(0,1),(0,-1),(1,1),(1,-1),(-1,1),(-1,-1)): + if dj == 1: gd[1:, :] |= g1c[:-1, :] + elif dj == -1: gd[:-1, :] |= g1c[1:, :] + elif di == 1: gd[:, 1:] |= g1c[:, :-1] + elif di == -1: gd[:, :-1] |= g1c[:, 1:] + elif dj == 1 and di == 1: gd[1:, 1:] |= g1c[:-1, :-1] + elif dj == 1 and di == -1: gd[1:, :-1] |= g1c[:-1, 1:] + elif dj == -1 and di == 1: gd[:-1, 1:] |= g1c[1:, :-1] + else: gd[:-1, :-1] |= g1c[1:, 1:] + g1c = gd +for _ in range(2): + ge = g1c.copy() + for dj, di in ((1,0),(-1,0),(0,1),(0,-1),(1,1),(1,-1),(-1,1),(-1,-1)): + if dj == 1: ge[1:, :] &= g1c[:-1, :] + elif dj == -1: ge[:-1, :] &= g1c[1:, :] + elif di == 1: ge[:, 1:] &= g1c[:, :-1] + elif di == -1: ge[:, :-1] &= g1c[:, 1:] + elif dj == 1 and di == 1: ge[1:, 1:] &= g1c[:-1, :-1] + elif dj == 1 and di == -1: ge[1:, :-1] &= g1c[:-1, 1:] + elif dj == -1 and di == 1: ge[:-1, 1:] &= g1c[1:, :-1] + else: ge[:-1, :-1] &= g1c[1:, 1:] + g1c = ge +g1c |= g1f +# 洞格的 gz1 用邻域最大值填充 —— 仅补两类缺口, 不覆盖正常原值: +# 1) g1c 缺口(栅格化角落缝隙, 原值 -inf); +# 2) 孤立针孔: 原值远低于全部邻居中位数(顶面近水平处三角形投影漏采 +# 格心, 射线打到的其实是更低的背面) -> 用邻域最大值顶回真实顶面. +# 尾部包裹段是陡斜坡(有低侧邻居, 单调趋势)不是针孔, 必须保留原值跟随本地框面 +gz1f = np.where(g1f, gz1, -np.inf) +_nb = [] +for dj, di in ((1,0),(-1,0),(0,1),(0,-1),(1,1),(1,-1),(-1,1),(-1,-1)): + s = np.full_like(gz1, np.nan) + if dj == 1: s[1:, :] = gz1[:-1, :] + elif dj == -1: s[:-1, :] = gz1[1:, :] + elif di == 1: s[:, 1:] = gz1[:, :-1] + elif di == -1: s[:, :-1] = gz1[:, 1:] + elif dj == 1 and di == 1: s[1:, 1:] = gz1[:-1, :-1] + elif dj == 1 and di == -1: s[1:, :-1] = gz1[:-1, 1:] + elif dj == -1 and di == 1: s[:-1, 1:] = gz1[1:, :-1] + else: s[:-1, :-1] = gz1[1:, 1:] + _nb.append(s) +_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) +gz1d = gz1f.copy() +for _ in range(3): + gf = gz1d.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(gz1d, -np.inf) + if dj == 1: s[1:, :] = gz1d[:-1, :] + elif dj == -1: s[:-1, :] = gz1d[1:, :] + elif di == 1: s[:, 1:] = gz1d[:, :-1] + elif di == -1: s[:, :-1] = gz1d[:, 1:] + elif dj == 1 and di == 1: s[1:, 1:] = gz1d[:-1, :-1] + elif dj == 1 and di == -1: s[1:, :-1] = gz1d[:-1, 1:] + elif dj == -1 and di == 1: s[:-1, 1:] = gz1d[1:, :-1] + else: s[:-1, :-1] = gz1d[1:, 1:] + gf = np.maximum(gf, s) + gz1d = gf +gz1f = np.where(g1f & ~_pin, gz1, gz1d) +log(f" 顶面针孔修补: {int(_pin.sum())} 格 (孤立低谷, 用邻域最大值补回顶面)") +nfill = int((g1c & ~g1f & R & (gz1f >= 7.0)).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 内真实覆盖格的最小值钳制(保守: 前表面绝不高于本地框面), +# 再轻平滑消除形态学噪声 +_mdil = np.where(g1f, gz1, np.inf) +for _ in range(3): + _mn = _mdil.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(_mdil, np.inf) + if dj == 1: s[1:, :] = _mdil[:-1, :] + elif dj == -1: s[:-1, :] = _mdil[1:, :] + elif di == 1: s[:, 1:] = _mdil[:, :-1] + elif di == -1: s[:, :-1] = _mdil[:, 1:] + elif dj == 1 and di == 1: s[1:, 1:] = _mdil[:-1, :-1] + elif dj == 1 and di == -1: s[1:, :-1] = _mdil[:-1, 1:] + elif dj == -1 and di == 1: s[:-1, 1:] = _mdil[1:, :-1] + 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) +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}") + +# 残差场边界只用拟合内点: 尾部包裹段(偏离前面平面 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为镜框自身几何, 并集后不暴露)") + +hb = np.full((H, W), np.nan) +# 真实脸面: fz >= FACE_MIN. 低于此值的"脸面"是误标导致的深伪面(后脑/口腔), +# 不能作为填充底面(否则镜桥下会出现十几单位深的伪填充柱) +real_f = have_f & (fz >= FACE_MIN) +hb[real_f] = fz[real_f] - EMBED_BOT +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) +log(f" 厚度上限(无真实脸面处, 距前面>{SHELL_MAX})生效格数: {ncap}") + +# 缝隙修补: 平滑可能把底缘底面抬到脸面之前 -> 钳回脸面内(嵌入>=0.25, 并集无缝) +hb = np.where(R & real_f, np.minimum(hb, fz - 0.25), hb) + +bzs = smooth_valid(np.where(have_b, bz, np.nan), have_b) +hb = np.where(R, hb, np.nan) +hb = np.where(have_b, np.maximum(hb, bzs + CLAMP_BACK), hb) +hb = np.minimum(hb, hf - GAP_MIN) + +# 外扩一格填充 NaN(边界采样用) +for _ in range(60): + nanm = ~np.isfinite(hf) + if not nanm.any(): break + for dj, di in ((1,0),(-1,0),(0,1),(0,-1)): + s = np.roll(np.where(nanm, 0, hf), (dj, di), axis=(0, 1)) + ok = ~np.roll(nanm, (dj, di), axis=(0, 1)) & np.isfinite(s) + hf = np.where(nanm & ok, s, hf) +for _ in range(60): + nanm = ~np.isfinite(hb) + if not nanm.any(): break + for dj, di in ((1,0),(-1,0),(0,1),(0,-1)): + 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) + +thickm = np.where(R, hf - hb, np.nan) +log(f" 填充厚度: min {np.nanmin(thickm):.2f}, 中位 {np.nanmedian(thickm):.2f}, " + f"mean {np.nanmean(thickm):.2f}, max {np.nanmax(thickm):.2f}") + +# ----------------------- 边界: 提取/平滑 ----------------------- +log("== 5. 边界提取与平滑 ==") +# 裂缝边(crack edge)链成回路 +edges = {} +def add_edge(a, b): + edges.setdefault(a, []).append(b) +for j in range(H): + for i in range(W): + if not R[j, i]: + continue + if i == 0 or not R[j, i - 1]: + add_edge((i, j), (i, j + 1)) + if i == W - 1 or not R[j, i + 1]: + add_edge((i + 1, j + 1), (i + 1, j)) + if j == 0 or not R[j - 1, i]: + add_edge((i + 1, j), (i, j)) + if j == H - 1 or not R[j + 1, i]: + add_edge((i, j + 1), (i + 1, j + 1)) + +if not edges: + raise SystemExit("无边界") +start = next(iter(edges)) +loop = [start] +cur = edges[start].pop(0) +used = {(start, cur)} +while cur != start: + loop.append(cur) + nxts = [e for e in edges.get(cur, []) if (cur, e) not in used] + if not nxts: + break + nxt = nxts[0] + used.add((cur, nxt)) + cur = nxt +P = np.array([[X0 + i * RES, Y0 + j * RES] for (i, j) in loop], float) +area2 = np.sum(P[:, 0] * np.roll(P[:, 1], -1) - np.roll(P[:, 0], -1) * P[:, 1]) +if area2 < 0: + P = P[::-1] +log(f" 边界回路 {len(P)} 点, 周长 {np.linalg.norm(np.roll(P,-1,0)-P,axis=1).sum():.2f}") + +def resample(P, step): + d = np.linalg.norm(np.roll(P, -1, 0) - P, axis=1) + L = np.concatenate([[0], np.cumsum(d)]) + total = L[-1] + n = max(8, int(round(total / step))) + t = np.linspace(0, total, n, endpoint=False) + x = np.interp(t, L, np.concatenate([P[:, 0], P[:1, 0]])) + y = np.interp(t, L, np.concatenate([P[:, 1], P[:1, 1]])) + return np.stack([x, y], 1) + +def roughness(P): + d = np.roll(P, -1, 0) - P + a = np.arctan2(d[:, 1], d[:, 0]) + da = np.abs((a - np.roll(a, -1) + np.pi) % (2 * np.pi) - np.pi) + return float(da.mean()), float(da.std()) + +r0 = roughness(P) +B = resample(P, BSTEP) +for _ in range(BITER): + B = (1 - BLAM) * B + BLAM * 0.5 * (np.roll(B, 1, 0) + np.roll(B, -1, 0)) +# 面积恢复: 抵消平滑收缩(目标 = R 面积的 97%) +areaR = R.sum() * RES * RES +for _ in range(30): + aB = 0.5 * abs(np.sum(B[:, 0] * np.roll(B[:, 1], -1) - np.roll(B[:, 0], -1) * B[:, 1])) + if aB >= 0.97 * areaR: + break + c = B.mean(0) + s = np.sqrt(0.97 * areaR / max(aB, 1e-9)) + s = min(s, 1.05) + B = c + (B - c) * s +B = resample(B, BSTEP) +r1 = roughness(B) +log(f" 边界粗糙度(平均转折角/标准差): 平滑前 {np.degrees(r0[0]):.2f}°/{np.degrees(r0[1]):.2f}° " + f"-> 平滑后 {np.degrees(r1[0]):.2f}°/{np.degrees(r1[1]):.2f}°") + +def inside_mask(pts, msk): + ii = np.clip(((pts[:, 0] - X0) / RES).astype(int), 0, W - 1) + jj2 = np.clip(((pts[:, 1] - Y0) / RES).astype(int), 0, H - 1) + return msk[jj2, ii] + +bad = ~inside_mask(B, mask) +if bad.any(): + c = B.mean(0) + for k in np.where(bad)[0]: + p = B[k] + for s in np.linspace(0.1, 0.9, 30): + q = c + s * (p - c) + if inside_mask(q[None, :], mask)[0]: + B[k] = q + break + for _ in range(10): + B = (1 - BLAM) * B + BLAM * 0.5 * (np.roll(B, 1, 0) + np.roll(B, -1, 0)) + B = resample(B, BSTEP) +# 边界须压在镜框带上(存在眼镜表面), 统计覆盖率 +bii = np.clip(((B[:, 0] - X0) / RES).astype(int), 0, W - 1) +bjj = np.clip(((B[:, 1] - Y0) / RES).astype(int), 0, H - 1) +nb = np.zeros(len(B), bool) +for dj, di in ((0,0),(1,0),(-1,0),(0,1),(0,-1),(1,1),(1,-1),(-1,1),(-1,-1)): + jj2 = np.clip(bjj + dj, 0, H - 1); ii2 = np.clip(bii + di, 0, W - 1) + nb |= g1f[jj2, ii2] +log(f" 平滑边界 {len(B)} 点, 超出足迹 {int((~inside_mask(B, mask)).sum())}, " + f"压在镜框带上 {nb.mean()*100:.1f}%") + +# ----------------------- 区域三角化: 耳切 + 内部点插入 ----------------------- +log("== 6. 区域三角化 (耳切法 + 内部网格点插入, 无翻转风险) ==") +Nb = len(B) + +def cross3(o, a, b): + return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0]) + +def ear_clip(P): + """P 为 CCW 简单多边形, 返回三角形索引列表""" + idx = list(range(len(P))) + tris = [] + guard = 0 + while len(idx) > 3 and guard < 5000: + guard += 1 + m = len(idx) + clipped = False + best_k, best_score = None, None + for k in range(m): + i0, i1, i2 = idx[(k - 1) % m], idx[k], idx[(k + 1) % m] + a, b, c = P[i0], P[i1], P[i2] + cr = cross3(a, b, c) + if cr <= 1e-10: + continue + inside_any = False + for j in idx: + if j in (i0, i1, i2): + continue + p = P[j] + if (cross3(a, b, p) >= -1e-12 and cross3(b, c, p) >= -1e-12 + and cross3(c, a, p) >= -1e-12): + inside_any = True + break + if not inside_any: + # 质量分: 面积/边长平方 越大越好 + lab = np.hypot(*(np.subtract(b, a))); lac = np.hypot(*(np.subtract(c, a))) + lbc = np.hypot(*(np.subtract(c, b))) + score = cr / max(lab * lac, lab * lbc, lac * lbc, 1e-12) + if best_score is None or score > best_score: + best_k, best_score = k, score + if best_k is None: + # 兜底: 耳切失败(理论不应发生), 取最大正面积顶点 + m = len(idx) + areas = [(cross3(P[idx[(k-1)%m]], P[idx[k]], P[idx[(k+1)%m]]), k) + for k in range(m)] + areas.sort(reverse=True) + best_k = areas[0][1] + log(" ! 耳切兜底触发") + i0, i1, i2 = (idx[(best_k - 1) % len(idx)], idx[best_k], + idx[(best_k + 1) % len(idx)]) + tris.append((i0, i1, i2)) + idx.pop(best_k) + tris.append((idx[0], idx[1], idx[2])) + return tris + +T = np.asarray(ear_clip(B.tolist()), np.int64) +log(f" 耳切三角化: {len(T)} 三角形") + +# 内部点: R 深处格子中心, 间隔 ~0.3 +inner2 = R.copy() +for _ in range(2): + ne = inner2.copy() + 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))) +cand = [] +for j in range(0, H, step): + for i in range(0, W, step): + if inner2[j, i]: + cand.append((X0 + i * RES, Y0 + j * RES)) +cand = np.asarray(cand) +log(f" 内部插入点候选 {len(cand)}") + +# 逐点插入(拆分包含它的三角形 -> 拓扑恒不翻转) +xy = np.vstack([B, cand]) +n_cand = len(cand) +extra = 0 +Ta = [list(map(int, t)) for t in T] +for ci in range(n_cand): + p = xy[Nb + ci] + m = len(Ta) + TA = np.asarray(Ta) + a = xy[TA[:, 0]]; b = xy[TA[:, 1]]; c = xy[TA[:, 2]] + d = ((b[:, 1] - c[:, 1]) * (a[:, 0] - c[:, 0]) + (c[:, 0] - b[:, 0]) * (a[:, 1] - c[:, 1])) + with np.errstate(divide='ignore', invalid='ignore'): + l1 = ((b[:, 1] - c[:, 1]) * (p[0] - c[:, 0]) + (c[:, 0] - b[:, 0]) * (p[1] - c[:, 1])) / d + l2 = ((c[:, 1] - a[:, 1]) * (p[0] - c[:, 0]) + (a[:, 0] - c[:, 0]) * (p[1] - c[:, 1])) / d + l3 = 1 - l1 - l2 + hit = np.where((l1 > 5e-3) & (l2 > 5e-3) & (l3 > 5e-3) + & (l1 < 1 - 5e-3) & (l2 < 1 - 5e-3) & (l3 < 1 - 5e-3))[0] + if len(hit) == 0: + continue + t = int(hit[0]) + i0, i1, i2 = Ta[t] + vnew = Nb + ci + Ta[t] = [i0, i1, vnew] + Ta.append([i1, i2, vnew]) + Ta.append([i2, i0, vnew]) + extra += 1 +T = np.asarray([t for t in Ta if t[0] != t[1] and t[1] != t[2] and t[0] != t[2]], np.int64) +xy = xy[:Nb + n_cand] + +# 内部顶点 XY 轻度 Laplacian 平滑(边界固定), 消除细长三角形; 有翻转则回退 +es = np.unique(np.sort(np.concatenate([T[:, [0, 1]], T[:, [1, 2]], T[:, [2, 0]]]), axis=1), axis=0) +src = np.concatenate([es[:, 0], es[:, 1]]) +dst = np.concatenate([es[:, 1], es[:, 0]]) +deg = np.bincount(src, minlength=len(xy)).astype(float) +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 + +# 边翻转优化(Lawson): 公共边换成严格更短的对角线(总边长严格递减 -> 必收敛), +# 不新增顶点即可消除耳切弦两侧的细长退化三角形(表面粗糙的根源) +def flip_pass(Tl, XY): + em = {} + for fi, t in enumerate(Tl): + a, b, c = t + for u, v in ((a, b), (b, c), (c, a)): + em.setdefault((min(u, v), max(u, v)), []).append(fi) + eset = set(em) + n = 0 + for (u, v), fs in em.items(): + if len(fs) != 2: + continue + f1, f2 = fs + if not (u in Tl[f1] and v in Tl[f1] and u in Tl[f2] and v in Tl[f2]): + continue # 本趟已被改动的面, 下趟再处理 + w1 = next(x for x in Tl[f1] if x != u and x != v) + w2 = next(x for x in Tl[f2] if x != u and x != v) + if w1 == w2 or (min(w1, w2), max(w1, w2)) in eset: + continue + pu, pv = XY[u], XY[v] + s1 = (pv[0]-pu[0])*(XY[w1][1]-pu[1]) - (pv[1]-pu[1])*(XY[w1][0]-pu[0]) + s2 = (pv[0]-pu[0])*(XY[w2][1]-pu[1]) - (pv[1]-pu[1])*(XY[w2][0]-pu[0]) + if s1 * s2 >= 0: + continue # w1,w2 须在 uv 两侧(凸四边形) + q1, q2 = XY[w1], XY[w2] + s3 = (q2[0]-q1[0])*(pu[1]-q1[1]) - (q2[1]-q1[1])*(pu[0]-q1[0]) + s4 = (q2[0]-q1[0])*(pv[1]-q1[1]) - (q2[1]-q1[1])*(pv[0]-q1[0]) + if s3 * s4 >= 0: + 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 拆分 + if s1 > 0: # w1 左, w2 右 + Tl[f1] = [u, w2, w1] + Tl[f2] = [v, w1, w2] + else: # w1 右, w2 左 + Tl[f1] = [u, w1, w2] + Tl[f2] = [v, w2, w1] + eset.add((min(w1, w2), max(w1, w2))) + n += 1 + return n + +Tlist = [list(map(int, t)) for t in T] +nflips = 0 +for _p in range(200): + k = flip_pass(Tlist, xy) + nflips += k + if k == 0: + break +log(f" 边翻转优化: {nflips} 次, 收敛于 {_p+1} 趟") + +# 残余长边轻度细分(翻转后仅少量, 不会级联) +MAXEDGE = 0.5 +nsplit = 0 +for _p in range(6): + emap = {} + for fi, (a, b, c) in enumerate(Tlist): + for u, v in ((a, b), (b, c), (c, a)): + emap.setdefault((min(u, v), max(u, v)), []).append(fi) + longe = [(u, v) for (u, v) in emap if np.hypot(*(xy[v] - xy[u])) > MAXEDGE] + if not longe: + break + mid = {} + newpts = [] + for u, v in longe: + mid[(u, v)] = len(xy) + len(newpts) + newpts.append((xy[u] + xy[v]) / 2.0) + xy = np.vstack([xy, np.asarray(newpts)]) + nsplit += len(newpts) + newT = [] + for (a, b, c) in Tlist: + mab = mid.get((min(a, b), max(a, b))) + mbc = mid.get((min(b, c), max(b, c))) + mca = mid.get((min(c, a), max(c, a))) + if mab is None and mbc is None and mca is None: + newT.append([a, b, c]) + elif mbc is None and mca is None: + newT += [[a, mab, c], [mab, b, c]] + elif mab is None and mca is None: + newT += [[a, b, mbc], [a, mbc, c]] + elif mab is None and mbc is None: + newT += [[a, b, mca], [mca, b, c]] + elif mca is None: + newT += [[a, mab, mbc], [mab, b, mbc], [a, mbc, c]] + elif mab is None: + newT += [[a, b, mbc], [a, mbc, mca], [mbc, c, mca]] + elif mbc is None: + newT += [[a, mab, mca], [mab, b, c], [mab, c, mca]] + else: + newT += [[a, mab, mca], [mab, b, mbc], [mca, mbc, c], [mab, mbc, mca]] + Tlist = newT +T = np.asarray(Tlist, np.int64) +log(f" 残余长边细分(>{MAXEDGE}): 新增 {nsplit} 中点, 总三角形 {len(T)}") + +cr = xy[T[:, 1]] - xy[T[:, 0]] +cs = xy[T[:, 2]] - xy[T[:, 0]] +area = cr[:, 0] * cs[:, 1] - cr[:, 1] * cs[:, 0] +if area.min() <= 0: + raise SystemExit(f"三角形定向异常 (min area {area.min():.2e})") +log(f" 实际插入 {extra} 点, 总三角形 {len(T)}, 最小有向面积 {area.min():.5f} (>0 = 无翻转)") +nV = len(xy) + +# ----------------------- 组装实体 ----------------------- +log("== 7. 组装填充实体 ==") +def bilinear(Grd, pts): + fx = (pts[:, 0] - X0) / RES + fy = (pts[:, 1] - Y0) / RES + i = np.clip(np.floor(fx).astype(int), 0, W - 2) + j = np.clip(np.floor(fy).astype(int), 0, H - 2) + tx = np.clip(fx - i, 0, 1); ty = np.clip(fy - j, 0, 1) + v = (Grd[j, i] * (1 - tx) * (1 - ty) + Grd[j, i + 1] * tx * (1 - ty) + + Grd[j + 1, i] * (1 - tx) * ty + Grd[j + 1, i + 1] * tx * ty) + return v + +zf = bilinear(hf, xy) +zb = bilinear(hb, xy) +zb = np.minimum(zb, zf - 0.15) +zb = np.maximum(zb, zf - 50) + +Pf = np.column_stack([xy, zf]) +Pb = np.column_stack([xy, zb]) +faces = [] +for a, b, c in T: + faces.append((a, b, c)) # 前表面(法向 +z, 向外) +for a, b, c in T: + faces.append((a + nV, c + nV, b + nV)) # 后表面(翻转, 法向 -z, 向外) +# 侧墙: 用细分后的三角化边界环(1-面有向边链成 CCW 环, 与盘边界严格一致) +dired = {} +unded = {} +for a, b, c in T: + for u, v in ((a, b), (b, c), (c, a)): + dired[(u, v)] = dired.get((u, v), 0) + 1 +for (u, v) in dired: + unded[(min(u, v), max(u, v))] = unded.get((min(u, v), max(u, v)), 0) + 1 +bnext = {} +for (u, v) in dired: + if unded[(min(u, v), max(u, v))] == 1: + bnext[u] = v +ring = [next(iter(bnext))] +while bnext[ring[-1]] != ring[0]: + ring.append(bnext[ring[-1]]) +Nr = len(ring) +for i in range(Nr): + u, v = ring[i], ring[(i + 1) % Nr] + faces.append((u, v + nV, v)) # (F_i, B_{i+1}, F_{i+1}) + faces.append((u, u + nV, v + nV)) # (F_i, B_i, B_{i+1}) +Ff = np.asarray(faces, np.int64) +Pall = np.vstack([Pf, Pb]) + +# 定向: 与原模型一致 +def signed_vol(P, Tt): + a = P[Tt[:, 0]]; b = P[Tt[:, 1]]; c = P[Tt[:, 2]] + return float(np.sum(np.einsum('ij,ij->i', a, np.cross(b, c))) / 6.0) +vfill = signed_vol(Pall, Ff) +vorig = signed_vol(V, F) +if vfill * vorig < 0: + Ff = Ff[:, [0, 2, 1]] + vfill = -vfill +# 期望体积 = 三角形 2D 面积 x 平均厚度 (验证定向正确性) +exp_vol = float(np.sum(0.5 * np.abs(area) * ((zf[T[:, 0]] + zf[T[:, 1]] + zf[T[:, 2]] + - zb[T[:, 0]] - zb[T[:, 1]] - zb[T[:, 2]]) / 3.0))) +log(f" 填充实体体积 {vfill:.2f} (期望 {exp_vol:.2f}), 原模型 {vorig:.0f}, 定向一致") +if abs(vfill - exp_vol) > 0.02 * abs(exp_vol) + 1e-6: + raise SystemExit("体积校验失败: 面片定向不一致") + +# 水密性: 无向边计数 + 有向边方向一致性 +e2 = np.sort(np.concatenate([Ff[:, [0, 1]], Ff[:, [1, 2]], Ff[:, [2, 0]]]), axis=1) +ue, cnt = np.unique(e2, axis=0, return_counts=True) +de = np.concatenate([Ff[:, [0, 1]], Ff[:, [1, 2]], Ff[:, [2, 0]]]) +ude, udcnt = np.unique(de, axis=0, return_counts=True) +log(f" 水密性: 非流形边(计数!=2) {int((cnt != 2).sum())} / {len(ue)}, " + f"方向冲突边 {int((udcnt != 1).sum())} (两者都应为 0)") + +# 镜片光滑度: 前表面相邻面二面角(打印外观指标) +nv = np.cross(Pall[Ff[:, 1]] - Pall[Ff[:, 0]], Pall[Ff[:, 2]] - Pall[Ff[:, 0]]) +nv = nv / np.linalg.norm(nv, axis=1, keepdims=True) +fcen = Pall[Ff].mean(1) +emap2 = {} +for fi, (a, b, c) in enumerate(Ff): + for u, v in ((a, b), (b, c), (c, a)): + emap2.setdefault((min(u, v), max(u, v)), []).append(fi) +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: # 左右镜片开口 + lang.append(np.degrees(np.arccos(np.clip(nv[fs[0]] @ nv[fs[1]], -1, 1)))) + lloc.append(m) +lang = np.asarray(lang) +if len(lang): + k = np.argmax(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)") +bmargin = (hb - bzs)[R & have_b] +log(f" 后表面高于最后表面: 最小余量 {np.nanmin(bmargin):.3f} (应>0)") +embed = (fz - hb)[R & real_f] +log(f" 底面嵌入真实脸面: 最小 {np.nanmin(embed):.3f} (应>-0.05, 负=缝隙), " + f"中位 {np.nanmedian(embed):.2f}") +inR = inside_mask(xy, R) +log(f" 圆盘顶点落在 R 内比例 {inR.mean()*100:.1f}%") + +# ----------------------- 输出 ----------------------- +log("== 8. 输出 ==") +with open(OUT_FILL, "w") as f: + f.write("# fill mesh between glasses and face\n") + f.write(f"o fill\n") + for p in Pall: + f.write(f"v {p[0]:.5f} {p[1]:.5f} {p[2]:.5f}\n") + for a, b, c in Ff: + f.write(f"f {a+1} {b+1} {c+1}\n") +log(f" {OUT_FILL}: {len(Pall)} 顶点, {len(Ff)} 面") + +with open(IN_OBJ) as f: + orig = f.read() +with open(OUT_MERGED, "w") as f: + f.write(orig.rstrip("\n") + "\n") + f.write("o fill\n") + off = len(V) + for p in Pall: + f.write(f"v {p[0]:.5f} {p[1]:.5f} {p[2]:.5f}\n") + for a, b, c in Ff: + f.write(f"f {a+1+off} {b+1+off} {c+1+off}\n") +log(f" {OUT_MERGED}: 原模型 + 填充 (顶点偏移 {off})") + +# ----------------------- 预览 ----------------------- +front = np.where(g1f, gz1, np.nan) +mix = np.where(R & np.isfinite(hf), np.where(np.isfinite(front), np.maximum(front, hf), hf), front) +front = mix +panels = [] +span = (4, 12) +def to_img(Grd, span): + G = np.where(np.isfinite(Grd), Grd, span[0]) + t = np.clip((G - span[0]) / (span[1] - span[0]), 0, 1) + return (t * 255).astype(np.uint8) +for name, Grd in (("front-depth", front), ("fill-front", np.where(R, hf, np.nan)), + ("thickness", np.where(R, hf - hb, np.nan))): + im = Image.fromarray(to_img(np.flipud(Grd), span if "thick" not in name else (0, 4))) + im = im.resize((W * 4, H * 4), Image.NEAREST) + d = ImageDraw.Draw(im) + pts = [((b[0] - X0) / RES * 4, (Y1 - b[1]) / RES * 4) for b in B] + d.line(pts + pts[:1], fill=255, width=2) + panels.append((name, im)) +cw = sum(p.width for _, p in panels) + 20 +ch = max(p.height for _, p in panels) +combo = Image.new("L", (cw, ch), 60) +xoff = 0 +for name, p in panels: + combo.paste(p, (xoff, 0)); xoff += p.width + 10 +combo.save(OUT_PREVIEW) +log(f" {OUT_PREVIEW}: 3 面板(前深度/填充前表面/厚度), 红线=平滑边界") + +log("\n== 完成 ==") +with open("fill_report.txt", "w") as f: + f.write("\n".join(report) + "\n") +print("报告已写入 fill_report.txt")