diff --git a/generate_fill_mesh_hsc.py b/generate_fill_mesh_hsc.py deleted file mode 100644 index 7957495..0000000 --- a/generate_fill_mesh_hsc.py +++ /dev/null @@ -1,1249 +0,0 @@ -#!/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 前视图深度预览(供人工检查) - -仅依赖 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): - """解析 OBJ 文件, 返回顶点数组和面数组""" - 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) - - -# ====================================================================== -# 形态学操作 -# ====================================================================== -def dilate_8(mask, iterations=1): - """二值图像膨胀 (8-邻域)""" - result = mask.copy() - for _ in range(iterations): - nm = result.copy() - nm[1:, :] |= result[:-1, :]; nm[:-1, :] |= result[1:, :] - nm[:, 1:] |= result[:, :-1]; nm[:, :-1] |= result[:, 1:] - nm[1:, 1:] |= result[:-1, :-1]; nm[:-1, :-1] |= result[1:, 1:] - nm[1:, :-1] |= result[:-1, 1:]; nm[:-1, 1:] |= result[1:, :-1] - result = nm - return result - - -def erode_4(mask, iterations=1): - """二值图像腐蚀 (4-邻域)""" - result = mask.copy() - for _ in range(iterations): - ne = result.copy() - ne[1:, :] &= result[:-1, :]; ne[:-1, :] &= result[1:, :] - ne[:, 1:] &= result[:, :-1]; ne[:, :-1] &= result[:, 1:] - result = ne - return result - - -def dilate_8_inf(field, iterations=1): - """浮点数场膨胀 (8-邻域取最大值)""" - result = field.copy() - for _ in range(iterations): - gf = result.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(result, -np.inf) - if dj == 1: s[1:, :] = result[:-1, :] - elif dj == -1: s[:-1, :] = result[1:, :] - elif di == 1: s[:, 1:] = result[:, :-1] - elif di == -1: s[:, :-1] = result[:, 1:] - elif dj == 1 and di == 1: s[1:, 1:] = result[:-1, :-1] - elif dj == 1 and di == -1: s[1:, :-1] = result[:-1, 1:] - elif dj == -1 and di == 1: s[:-1, 1:] = result[1:, :-1] - else: s[:-1, :-1] = result[1:, 1:] - gf = np.maximum(gf, s) - result = gf - return result - - -def erode_4_inf(field, iterations=1): - """浮点数场腐蚀 (4-邻域取最小值)""" - result = field.copy() - for _ in range(iterations): - mn = result.copy() - for dj, di in ((1,0),(-1,0),(0,1),(0,-1)): - s = np.full_like(result, np.inf) - if dj == 1: s[1:, :] = result[:-1, :] - elif dj == -1: s[:-1, :] = result[1:, :] - elif di == 1: s[:, 1:] = result[:, :-1] - elif di == -1: s[:, :-1] = result[:, 1:] - mn = np.minimum(mn, s) - result = mn - return result - - -def flood_outside(mask): - """从窗口边界洪泛, 标记外部; 未标记且非 mask 即为洞""" - outside = np.zeros_like(mask, bool) - dq = deque() - h, w = mask.shape - for i in range(w): - for j in (0, h - 1): - if not mask[j, i] and not outside[j, i]: - outside[j, i] = True; dq.append((j, i)) - for j in range(h): - for i in (0, w - 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 < h and 0 <= ni < w: - if not mask[nj, ni] and not outside[nj, ni]: - outside[nj, ni] = True; dq.append((nj, ni)) - return ~outside & ~mask - - -def largest_component_8(mask): - """取 mask 的最大连通分量 (8-邻域)""" - labels = np.zeros(mask.shape, np.int32) - cur = 0 - for j in range(mask.shape[0]): - for i in range(mask.shape[1]): - if mask[j, i] and labels[j, i] == 0: - cur += 1 - q = deque([(j, i)]) - labels[j, i] = cur - 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 mask[ny, nx] and labels[ny, nx] == 0): - labels[ny, nx] = cur - q.append((ny, nx)) - sizes = np.bincount(labels.ravel()) - sizes[0] = 0 - main_label = int(np.argmax(sizes)) - stray = int(mask.sum() - sizes[main_label]) - return labels == main_label, stray - - -# ====================================================================== -# 数值工具 -# ====================================================================== -def pbasis(x, y): - """二次多项式基函数: [1, x, y, x^2, xy, y^2]""" - return np.stack([np.ones_like(x), x, y, x * x, x * y, y * y], 1) - - -def cross3(o, a, b): - """2D 叉积 (o->a) x (o->b)""" - return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0]) - - -def signed_volume(verts, faces): - """计算有符号体积""" - a = verts[faces[:, 0]]; b = verts[faces[:, 1]]; c = verts[faces[:, 2]] - return float(np.sum(np.einsum('ij,ij->i', a, np.cross(b, c))) / 6.0) - - -def inside_mask(pts, msk): - """判断点是否在 mask 内""" - ii = np.clip(((pts[:, 0] - X0) / RES).astype(int), 0, msk.shape[1] - 1) - jj = np.clip(((pts[:, 1] - Y0) / RES).astype(int), 0, msk.shape[0] - 1) - return msk[jj, ii] - - -def bilinear(grid, pts): - """双线性插值""" - fx = (pts[:, 0] - X0) / RES - fy = (pts[:, 1] - Y0) / RES - i = np.clip(np.floor(fx).astype(int), 0, grid.shape[1] - 2) - j = np.clip(np.floor(fy).astype(int), 0, grid.shape[0] - 2) - tx = np.clip(fx - i, 0, 1); ty = np.clip(fy - j, 0, 1) - return (grid[j, i] * (1 - tx) * (1 - ty) + - grid[j, i + 1] * tx * (1 - ty) + - grid[j + 1, i] * (1 - tx) * ty + - grid[j + 1, i + 1] * tx * ty) - - -# ====================================================================== -# 平滑与插值 -# ====================================================================== -def smooth_valid(field, valid, sigma=SIGMA, it=1): - """归一化高斯卷积(只在 valid 上)""" - rad = int(3 * sigma) - k = np.exp(-0.5 * (np.arange(-rad, rad + 1) / sigma) ** 2) - field = field.copy() - for _ in range(it): - acc = np.zeros_like(field); accw = np.zeros_like(field) - val = valid & np.isfinite(field) - G = np.where(val, field, 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), field) - field = np.where(val, out, field) - return field - - -def laplace_fill(field, known, iters=6000): - """Dirichlet: known 固定, 其余解调和方程 - 缺失邻居(网格边缘)用反射(补 F 自身), 恒除以 4 保证系数和=1(收敛)""" - F = field.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 - 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) - G = np.where(known, F, avg / 4.0) - 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 - - -# ====================================================================== -# 步骤 1: 加载模型并分类眼镜面 -# ====================================================================== -def load_and_classify(): - """加载模型, 通过边连通主分量提取真眼镜物体""" - 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)}") - - # 边连通主分量提取 - 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) - emap.setdefault(e, []).append(fi) - - seen = set() - best = [] - for fi in np.where(face_g)[0]: - if fi in seen: - continue - comp, stack = [], [fi] - seen.add(fi) - while stack: - f0 = stack.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) - stack.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))} 面") - return V, F, face_g - - -# ====================================================================== -# 步骤 2: 栅格化 -# ====================================================================== -def rasterize(V, F, face_g): - """栅格化三角形到高度场""" - 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) # 全网格最后表面(后脑) - - 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]] - for t in np.where(fwin)[0]: - 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 - v = (d00 * d12 - d01 * d02) / den - 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) - g1f = have_g1 - log(f" 有眼镜表面格 {have_g1.sum()}, 管状(双面) {have_g2.sum()}, " - f"有脸面 {have_f.sum()}, 有背面 {have_b.sum()}") - - return gz1, gz2, fz, bz, have_g1, have_g2, have_f, have_b, g1f, W, H - - -# ====================================================================== -# 步骤 3: 构建足迹 mask 与填充区域 R -# ====================================================================== -def build_footprint(g1f): - """构建眼镜足迹 mask 和填充区域 R""" - log("== 3. 眼镜足迹 / 填充区域 R ==") - mask = g1f.copy() - - # 取主分量 - mask, n_stray = largest_component_8(mask) - log(f" 足迹取主分量: 剔除杂散块 {n_stray} 格 (误标到镜框外的皮肤/后脑)") - - # 中央区剪除(鼻侧多余填充) - holes = ~mask & ~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 - log(f" 中央区剪除(两镜片开口之间, 开口底缘 y>{Y0 + y_cut * RES:.1f} 以下): " - f"剔除 {n_cut} 格误标鼻面") - else: - log(" ! 未找到两个镜片开口, 跳过中央区剪除") - - # 膨胀 3 格 - mask = dilate_8(mask, 3) - - # 补内部洞 (镜片开口) - holes = flood_outside(mask) - log(f" 足迹 {mask.sum()} 格, 补内部洞 {holes.sum()} 格 (镜片开口)") - mask = mask | holes - - # 腐蚀 - eroded = erode_4(mask, ERODE) - - # y 裁剪 + 最大连通块 - H, W = eroded.shape - jj = np.arange(H) - eroded = eroded & (jj[:, None] >= int((RY0 - Y0) / RES)) & \ - (jj[:, None] <= int((RY1 - Y0) / 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) - - 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 过小") - - return mask, R - - -# ====================================================================== -# 步骤 4: 构造填充前/后高度场 -# ====================================================================== -def build_height_fields(gz1, gz2, fz, bz, have_g1, have_f, have_b, g1f, mask, R): - """构造填充前/后高度场""" - log("== 4. 构造填充前/后高度场 ==") - H, W = gz1.shape - gx, gy = np.meshgrid(X0 + np.arange(W) * RES, Y0 + np.arange(H) * RES) - - # ---- 修补眼镜表面 ---- - g1c = dilate_8(g1f, 2) - g1c = erode_4(g1c, 2) - g1c |= g1f - - # 修补针孔 - 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 = dilate_8_inf(gz1f, 3) - gz1f = np.where(g1f & ~pin, gz1, gz1d) - log(f" 顶面针孔修补: {int(pin.sum())} 格 (孤立低谷, 用邻域最大值补回顶面)") - nfill = int((g1c & ~g1f & R & (gz1f >= 7.0)).sum()) - log(f" 镜框带补栅格化角落缺口: {nfill} 格") - - # ---- 框前面拟合 ---- - mfit = g1c & R & (gz1f >= 7.0) - - # 框带定义 - eR = erode_4(R, 3) - outer_band = R & ~eR - frame_band = (g1c & R) | (outer_band & (gz1f > -1e29)) - - # 陡斜坡外缘钳制 - mdil = np.where(g1f, gz1, np.inf) - mdil = erode_4_inf(mdil, 3) - fringe = frame_band & ~g1f - gz1f = np.where(fringe & (mdil < 1e29), np.minimum(gz1f, mdil), gz1f) - 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}") - - # 残差场调和延拓 - mfit2 = mfit & (np.abs(gz1f - zfit) <= 1.0) - 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) - - # ---- 前表面 ---- - 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) - - # 开口区安全钳制 - gmax2 = np.where(mfit, gz1f, -np.inf) - gmax2 = dilate_8_inf(gmax2, 3) - 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) - - # 墙区安全钳制 - edge_in = erode_4(R, 4) - wall_zone = R & ~edge_in - 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) - - # 交界连续性检查 - 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) - 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) - - 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}") - - 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}") - - return hf, hb, frame_band, wall_zone, gz1f, bzs, real_f, fz - - -# ====================================================================== -# 步骤 5: 边界提取与平滑 -# ====================================================================== -def extract_and_smooth_boundary(R, mask, g1f): - """提取并平滑边界""" - log("== 5. 边界提取与平滑 ==") - H, W = R.shape - - # 裂缝边链成回路 - 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)) - - # 面积恢复 - 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}°") - - # 边界须压在镜框带上 - 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}%") - - return B, loop - - -# ====================================================================== -# 步骤 6: 三角剖分 (耳切法 + 内部点插入) -# ====================================================================== -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) - 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 - - -def triangulate_fill(B, R, hf, hb): - """对填充区域进行三角剖分""" - log("== 6. 区域三角化 (耳切法 + 内部网格点插入, 无翻转风险) ==") - H, W = R.shape - Nb = len(B) - - # 耳切三角化 - T = np.asarray(ear_clip(B.tolist()), np.int64) - log(f" 耳切三角化: {len(T)} 三角形") - - # 内部点候选 - 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] - - # 内部顶点轻度 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 - - # 边翻转优化 - 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 - 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 - if s1 > 0: - Tl[f1] = [u, w2, w1] - Tl[f2] = [v, w1, w2] - else: - 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 = 无翻转)") - - # 计算 z 坐标 - 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]) - - return Pf, Pb, T, xy - -# ====================================================================== -# 步骤 7: 组装实体 (前/后表面 + 侧墙) -# ====================================================================== -def assemble_mesh(Pf, Pb, T): - """组装填充实体: 前表面 + 后表面 + 侧墙""" - log("== 7. 组装填充实体 ==") - nV = len(Pf) - faces = [] - - # 前表面 (法向 +z, 向外) - 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)) - - # 侧墙: 提取边界环 - 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]) - - return Pall, Ff, nV - - -# ====================================================================== -# 步骤 8: 校验 -# ====================================================================== -def validate_mesh(Pall, Ff, Pf, Pb, T, V, F, hf, hb, gz1f, frame_band, - wall_zone, bzs, have_b, real_f, fz, R, xy): - """对生成的网格进行各项校验""" - log("== 8. 校验 ==") - - # 定向校验 - vfill = signed_volume(Pall, Ff) - vorig = signed_volume(V, F) - if vfill * vorig < 0: - Ff = Ff[:, [0, 2, 1]] - vfill = -vfill - - # 期望体积 - area = np.cross(Pf[T[:, 1]] - Pf[T[:, 0]], Pf[T[:, 2]] - Pf[T[:, 0]])[:, 2] - exp_vol = float(np.sum(0.5 * np.abs(area) * - ((Pf[T[:, 0], 2] + Pf[T[:, 1], 2] + Pf[T[:, 2], 2] - - Pb[T[:, 0], 2] - Pb[T[:, 1], 2] - Pb[T[:, 2], 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}°, " - f"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}%") - - return Ff, vfill - - -# ====================================================================== -# 步骤 9: 输出 -# ====================================================================== -def output_mesh(Pall, Ff, V, F): - """输出填充网格和合并后的模型""" - log("== 9. 输出 ==") - - # 输出填充网格 - with open(OUT_FILL, "w") as f: - f.write("# fill mesh between glasses and face\n") - f.write("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})") - - -# ====================================================================== -# 步骤 10: 生成预览图 -# ====================================================================== -def generate_preview(gz1, hf, hb, R, B): - """生成前视图深度预览""" - log("== 10. 生成预览图 ==") - - # 合成前表面深度图 - front = np.where(np.isfinite(gz1) & (gz1 > -1e30), 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 面板(前深度/填充前表面/厚度), 红线=平滑边界") - - -# ====================================================================== -# 主流程 -# ====================================================================== -if __name__ == "__main__": - # 步骤 1: 加载模型 - V, F, face_g = load_and_classify() - - # 步骤 2: 栅格化 - gz1, gz2, fz, bz, have_g1, have_g2, have_f, have_b, g1f, W, H = \ - rasterize(V, F, face_g) - - # 步骤 3: 构建足迹 - mask, R = build_footprint(g1f) - - # 步骤 4: 构造高度场 - hf, hb, frame_band, wall_zone, gz1f, bzs, real_f, fz = \ - build_height_fields(gz1, gz2, fz, bz, have_g1, have_f, have_b, - g1f, mask, R) - - # 步骤 5: 边界提取与平滑 - B, loop = extract_and_smooth_boundary(R, mask, g1f) - - # 步骤 6: 三角剖分 - Pf, Pb, T, xy = triangulate_fill(B, R, hf, hb) - - # 步骤 7: 组装实体 - Pall, Ff, nV = assemble_mesh(Pf, Pb, T) - - # 步骤 8: 校验 - Ff, vfill = validate_mesh(Pall, Ff, Pf, Pb, T, V, F, hf, hb, - gz1f, frame_band, wall_zone, bzs, - have_b, real_f, fz, R, xy) - - # 步骤 9: 输出 - output_mesh(Pall, Ff, V, F) - - # 步骤 10: 预览图 - generate_preview(gz1, hf, hb, R, B) - - log("\n== 完成 ==") - with open("fill_report.txt", "w") as f: - f.write("\n".join(report) + "\n") - print("报告已写入 fill_report.txt") \ No newline at end of file