1778 lines
75 KiB
Python
1778 lines
75 KiB
Python
#!/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
|
||
|
||
import subprocess
|
||
import os
|
||
import sys
|
||
import argparse
|
||
|
||
def auto_bevel_in_blender(input_obj_path, output_obj_path, bevel_width=0.5, angle_deg=15):
|
||
blender_script = (
|
||
"import bpy, math\n"
|
||
"bpy.ops.object.select_all(action='SELECT'); bpy.ops.object.delete()\n"
|
||
f"bpy.ops.wm.obj_import(filepath=r'{input_obj_path}')\n"
|
||
"obj = bpy.context.active_object\n"
|
||
"bpy.context.view_layer.objects.active = obj\n"
|
||
"bpy.ops.object.transform_apply(scale=True)\n"
|
||
"bpy.ops.object.mode_set(mode='EDIT')\n"
|
||
"bpy.ops.mesh.select_all(action='DESELECT')\n"
|
||
# 选中非流形边（通常是交界线/边界）\n"
|
||
"bpy.ops.mesh.select_non_manifold(extend=False)\n"
|
||
# 若上面选不出，改用：bpy.ops.mesh.edges_select_sharp(sharpness=math.radians(5))\n"
|
||
f"bpy.ops.mesh.bevel(offset={bevel_width}, segments=4, profile=0.5, clamp_overlap=True)\n"
|
||
"bpy.ops.object.mode_set(mode='OBJECT')\n"
|
||
"bpy.ops.object.shade_smooth()\n"
|
||
f"bpy.ops.wm.obj_export(filepath=r'{output_obj_path}')\n"
|
||
"print('###BEVEL_OK###')\n"
|
||
)
|
||
|
||
script_path = "_temp_bevel.py"
|
||
with open(script_path, 'w') as f:
|
||
f.write(blender_script)
|
||
|
||
try:
|
||
result = subprocess.run(
|
||
['blender', '--background', '--python-exit-code', '1', '--python', script_path],
|
||
capture_output=True, text=True, timeout=120
|
||
)
|
||
# 合并 stdout 和 stderr 一起检查（Blender 喜欢把正常信息输出到 stderr）
|
||
all_output = result.stdout + "\n" + result.stderr
|
||
print("STDOUT:", result.stdout)
|
||
|
||
# 只检查成功标记，忽略 stderr 的 "Writing to" 信息
|
||
if '###BEVEL_OK###' not in all_output:
|
||
print("⚠️ 脚本未执行完")
|
||
print("STDERR:", result.stderr)
|
||
return False
|
||
if result.returncode != 0:
|
||
print("⚠️ Blender 返回码异常:", result.returncode)
|
||
return False
|
||
print("✅ 倒角成功完成")
|
||
return True
|
||
except FileNotFoundError:
|
||
print("❌ 未找到 Blender")
|
||
return False
|
||
except subprocess.TimeoutExpired:
|
||
print("❌ 超时")
|
||
return False
|
||
finally:
|
||
if os.path.exists(script_path):
|
||
os.remove(script_path)
|
||
|
||
class FillMeshGenerator:
|
||
"""眼镜-脸部填充网格生成器"""
|
||
|
||
# ======================================================================
|
||
# 参数配置
|
||
# ======================================================================
|
||
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 = -16.0, 16.0 # 原来 -10.2, 10.2 → 扩大到覆盖 x[-13.61, 13.19]
|
||
Y0, Y1 = 0.0, 62.0
|
||
RY0, RY1 = 1.0, 61.0
|
||
SEP = 0.25 # 区分前后表面的最小 z 间距
|
||
FLUSH = 0.15 # 填充前表面距镜框前面的距离(与镜框平齐, 留嵌入余量)
|
||
EMBED_BOT = 0.8 # 填充底面嵌入脸部深度
|
||
SHELL_MAX = 3.0 # 无真实脸面处的填充厚度上限(相对前表面)
|
||
FACE_MIN = 1.5 # 脸面低于此值视为误标伪面(后脑等), 不作填充底面
|
||
CLAMP_BACK = 0.3 # 距最后表面的安全距离(防穿出后脑)
|
||
GAP_MIN = 0.25 # 填充最小厚度
|
||
SIGMA = 2.0 # 高度场平滑 sigma(格)
|
||
ERODE = 4 # 足迹腐蚀(格) 0.2, 与膨胀0.3配合净外扩0.1(补偿边界平滑收缩)
|
||
BSTEP = 0.18 # 边界重采样步长
|
||
BITER = 25 # 边界 Laplacian 平滑迭代
|
||
BLAM = 0.6 # 边界平滑系数
|
||
TOL = 1e-4 # Jacobi 收敛阈值
|
||
|
||
def __init__(self, in_obj=None, in_idx=None, out_fill=None):
|
||
if in_obj:
|
||
self.IN_OBJ = in_obj
|
||
if in_idx:
|
||
self.IN_IDX = in_idx
|
||
if out_fill:
|
||
self.OUT_FILL = out_fill
|
||
# 根据 OUT_FILL 自动推导 OUT_MERGED 和 OUT_PREVIEW
|
||
base, ext = os.path.splitext(self.OUT_FILL)
|
||
self.OUT_MERGED = base + "_filled" + ext
|
||
self.OUT_PREVIEW = base.replace("mesh", "preview") + ".png"
|
||
if self.OUT_PREVIEW == base + ".png":
|
||
self.OUT_PREVIEW = base + "_preview.png"
|
||
self.report = []
|
||
self.lens_mask = None
|
||
|
||
def log(self, s):
|
||
"""日志记录"""
|
||
print(s)
|
||
self.report.append(str(s))
|
||
|
||
# ======================================================================
|
||
# IO 工具
|
||
# ======================================================================
|
||
def parse_obj(self, 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(self, 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(self, 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(self, 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(self, 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(self, 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(self, 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(self, 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(self, 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(self, 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(self, pts, msk):
|
||
"""判断点是否在 mask 内"""
|
||
ii = np.clip(((pts[:, 0] - self.X0) / self.RES).astype(int), 0, msk.shape[1] - 1)
|
||
jj = np.clip(((pts[:, 1] - self.Y0) / self.RES).astype(int), 0, msk.shape[0] - 1)
|
||
return msk[jj, ii]
|
||
|
||
def bilinear(self, grid, pts):
|
||
"""双线性插值"""
|
||
fx = (pts[:, 0] - self.X0) / self.RES
|
||
fy = (pts[:, 1] - self.Y0) / self.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(self, field, valid, sigma=None, it=1):
|
||
"""归一化高斯卷积(只在 valid 上)"""
|
||
if sigma is None:
|
||
sigma = self.SIGMA
|
||
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(self, field, known, iters=6000):
|
||
"""Dirichlet: known 固定, 其余解调和方程"""
|
||
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 < self.TOL * self.RES:
|
||
break
|
||
return F
|
||
|
||
# ======================================================================
|
||
# 步骤 1: 加载模型并分类眼镜面
|
||
# ======================================================================
|
||
def load_and_classify(self):
|
||
"""加载模型, 通过边连通主分量提取真眼镜物体"""
|
||
self.log("== 1. 读入模型 ==")
|
||
V, F = self.parse_obj(self.IN_OBJ)
|
||
idx = np.loadtxt(self.IN_IDX, dtype=np.int64) - 1
|
||
ism = np.zeros(len(V), bool)
|
||
ism[idx] = True
|
||
face_g = ism[F].any(1)
|
||
self.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
|
||
self.log(f" 眼镜物体 = 边连通主分量 {len(best)} 面, "
|
||
f"剥离误标皮肤 {int((ism[F].any(1)).sum() - len(best))} 面")
|
||
return V, F, face_g
|
||
|
||
# ======================================================================
|
||
# 步骤 2: 栅格化
|
||
# ======================================================================
|
||
def rasterize(self, V, F, face_g):
|
||
"""栅格化三角形到高度场"""
|
||
W = int(round((self.X1 - self.X0) / self.RES)) + 1
|
||
H = int(round((self.Y1 - self.Y0) / self.RES)) + 1
|
||
self.log(f"== 2. 栅格化 (窗口 x[{self.X0},{self.X1}] y[{self.Y0},{self.Y1}], "
|
||
f"res={self.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] > self.X0 - 0.5) & (V[:, 0] < self.X1 + 0.5) &
|
||
(V[:, 1] > self.Y0 - 0.5) & (V[:, 1] < self.Y1 + 0.5))
|
||
fwin = inwin[F].all(1)
|
||
self.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) - self.X0) / self.RES), 0)
|
||
i1 = min(int((max(ax, bx, cx) - self.X0) / self.RES), W - 1)
|
||
j0 = max(int((min(ay, by, cy) - self.Y0) / self.RES), 0)
|
||
j1 = min(int((max(ay, by, cy) - self.Y0) / self.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 = self.Y0 + j * self.RES
|
||
for i in range(i0, i1 + 1):
|
||
px = self.X0 + i * self.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] > self.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] - self.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
|
||
self.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
|
||
|
||
def build_footprint(self, g1f):
|
||
"""构建眼镜足迹 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} 格 (误标到镜框外的皮肤/后脑)")
|
||
|
||
# 中央区剪除(鼻侧多余填充)
|
||
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} 格误标鼻面")
|
||
|
||
# ===== 新增: 创建 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)
|
||
|
||
# 补内部洞 (镜片开口)
|
||
holes = self.flood_outside(mask)
|
||
self.log(f" 足迹 {mask.sum()} 格, 补内部洞 {holes.sum()} 格 (镜片开口)")
|
||
mask = mask | holes
|
||
|
||
# 腐蚀
|
||
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))
|
||
|
||
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} (触碰=侧墙切在头内部, 正常)")
|
||
|
||
# ===== [抽取自原版] 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: 构造填充前/后高度场
|
||
# ======================================================================
|
||
def build_height_fields(self, gz1, gz2, fz, bz, have_g1, have_f, have_b, g1f, mask, R):
|
||
"""构造填充前/后高度场"""
|
||
self.log("== 4. 构造填充前/后高度场 ==")
|
||
H, W = gz1.shape
|
||
gx, gy = np.meshgrid(self.X0 + np.arange(W) * self.RES,
|
||
self.Y0 + np.arange(H) * self.RES)
|
||
|
||
# ---- 修补眼镜表面 ----
|
||
g1c = self.dilate_8(g1f, 2)
|
||
g1c = self.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 = self.dilate_8_inf(gz1f, 3)
|
||
gz1f = np.where(g1f & ~pin, gz1, gz1d)
|
||
self.log(f" 顶面针孔修补: {int(pin.sum())} 格 (孤立低谷, 用邻域最大值补回顶面)")
|
||
nfill = int((g1c & ~g1f & R & (gz1f >= 7.0)).sum())
|
||
self.log(f" 镜框带补栅格化角落缺口: {nfill} 格")
|
||
|
||
# ---- 框前面拟合 ----
|
||
mfit = g1c & R & (gz1f >= 7.0)
|
||
|
||
# 框带定义
|
||
eR = self.erode_4(R, 3)
|
||
outer_band = R & ~eR
|
||
frame_band = (g1c & R) | (outer_band & (gz1f > -1e29))
|
||
|
||
# 陡斜坡外缘钳制
|
||
mdil = np.where(g1f, gz1, np.inf)
|
||
mdil = self.erode_4_inf(mdil, 3)
|
||
fringe = frame_band & ~g1f
|
||
gz1f = np.where(fringe & (mdil < 1e29), np.minimum(gz1f, mdil), gz1f)
|
||
gz1f = np.where(fringe, self.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 = self.pbasis(xp, yp) * wt[:, None]
|
||
coef, *_ = np.linalg.lstsq(Aw, zp * wt, rcond=None)
|
||
rfit = zp - self.pbasis(xp, yp) @ coef
|
||
wt = 1.0 / (1.0 + (np.abs(rfit) / 0.5) ** 2)
|
||
zfit = (self.pbasis(gx.ravel(), gy.ravel()) @ coef).reshape(H, W)
|
||
rf = np.abs((zfit - gz1f)[mfit])
|
||
self.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 = self.smooth_valid(resf, mfit2, sigma=2.5)
|
||
Rext = self.laplace_fill(resf, mfit2)
|
||
Rext = np.where(R, Rext, np.nan)
|
||
|
||
# ---- 前表面 ----
|
||
hf = np.full((H, W), np.nan)
|
||
hf[frame_band] = gz1f[frame_band] - self.FLUSH
|
||
openf = R & ~frame_band
|
||
hfo = np.where(openf, zfit + Rext - self.FLUSH, np.nan)
|
||
|
||
# 开口区安全钳制
|
||
gmax2 = np.where(mfit, gz1f, -np.inf)
|
||
gmax2 = self.dilate_8_inf(gmax2, 3)
|
||
hfo = np.where(openf & (gmax2 > -1e29), np.minimum(hfo, gmax2 - self.FLUSH), hfo)
|
||
hfo = self.smooth_valid(hfo, openf, sigma=2.0, it=2)
|
||
hf = np.where(openf & np.isfinite(hfo), hfo, hf)
|
||
hf = self.smooth_valid(hf, R, sigma=2.0, it=3)
|
||
# 可选：用周围有效值填充孤立NaN（若还有洞）
|
||
hf = np.where(R & ~np.isfinite(hf), np.nanmean(hf[R]), hf)
|
||
|
||
# 墙区双向钳制（不缩回，与镜面齐平）
|
||
edge_in = self.erode_4(R, 4)
|
||
wall_zone = R & ~edge_in
|
||
|
||
# 计算钳制上下界
|
||
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)])
|
||
|
||
# 上界：陡降格防鳍，平缓格不设上限（允许微凸）
|
||
ceil = np.where((gz1f > -1e29) & (gzfm < gmax2 - 1.2), gzfm + 0.30, np.inf)
|
||
# 下界：真实覆盖格不低于本地框面 - 0.15（允许微凹但不缩回）
|
||
bc = R & g1c & (gz1f > -1e29)
|
||
floor = np.where(bc, gz1f - 0.15, -np.inf)
|
||
|
||
# 双向钳制
|
||
hf = np.minimum(np.maximum(hf, floor), ceil)
|
||
|
||
# ===== 镜片弧面凸出 =====
|
||
_nbr = (np.roll(R, 1, 0) | np.roll(R, -1, 0) | np.roll(R, 1, 1) | np.roll(R, -1, 1))
|
||
_ring = ~R & _nbr
|
||
_jr, _ir = np.where(R)
|
||
_js, _is = np.where(_ring)
|
||
_dx = (_ir[:, None] - _is[None, :]) * self.RES
|
||
_dy = (_jr[:, None] - _js[None, :]) * self.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 = self.smooth_valid(mprof, R, sigma=2.5)
|
||
BULGE = 0.45
|
||
hf = np.where(R, hf + BULGE * mprof, np.nan)
|
||
# ===== 结束镜片弧面 =====
|
||
|
||
# 交界连续性检查
|
||
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]), (self.X0 + ja[1] * self.RES, self.Y0 + ja[0] * self.RES)
|
||
self.log(f" 交界处高度跳变最大 {jn:.3f} @ {jn_at} (平面前框带接缝应<0.3; "
|
||
"镜圈圆角处~0.4为镜框自身几何, 并集后不暴露)")
|
||
|
||
# 交界连续性检查
|
||
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]), (self.X0 + ja[1] * self.RES, self.Y0 + ja[0] * self.RES)
|
||
self.log(f" 交界处高度跳变最大 {jn:.3f} @ {jn_at} (平面前框带接缝应<0.3; "
|
||
"镜圈圆角处~0.4为镜框自身几何, 并集后不暴露)")
|
||
|
||
# 在"交界连续性检查"之后，增加：
|
||
nose_zone = R & frame_band & (np.abs(hf - gz1f) < 0.5) # 鼻梁过渡区
|
||
hf = np.where(nose_zone, self.smooth_valid(hf, nose_zone, sigma=3.0, it=4), hf)
|
||
|
||
# ===== 在这里添加双向钳制逻辑 =====
|
||
# 计算 gzfm 和 gmax2（用于防鳍天花板）
|
||
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)])
|
||
|
||
# 下界：允许微凹但不缩回（贴合框面）
|
||
bc = R & g1c & (gz1f > -1e29)
|
||
floor = np.where(bc, gz1f - 0.15, -np.inf)
|
||
# 上界：平缓格不设上限（允许微凸），陡降格防鳍
|
||
ceil = np.where((gz1f > -1e29) & (gzfm < gmax2 - 1.2), gzfm + 0.30, np.inf)
|
||
|
||
# 双向钳制（替代原来的单向缩回）
|
||
hf = np.minimum(np.maximum(hf, floor), ceil)
|
||
# ===== 结束双向钳制 =====
|
||
|
||
# ===== 镜框边界融合（消除棱角） =====
|
||
blend_margin = 0.15 # 越小越窄（0.15~0.3）
|
||
blend_zone = R & frame_band & (np.abs(hf - gz1f) < blend_margin)
|
||
hf = np.where(blend_zone, 0.75 * hf + 0.25 * gz1f, hf) # 混合比例可调
|
||
# ===== 结束边界融合 =====
|
||
|
||
# ==========================================
|
||
# 【新增】镜框内边缘强制填沟（放在这里！）
|
||
# ==========================================
|
||
_inner = R & (
|
||
np.roll(frame_band, 1, 0) | np.roll(frame_band, -1, 0) |
|
||
np.roll(frame_band, 1, 1) | np.roll(frame_band, -1, 1)
|
||
) & ~frame_band
|
||
|
||
if np.any(_inner):
|
||
# 填平深沟：0.5:0.5 可尝试调为 0.6*hf + 0.4*gz1f（更偏面部）或 0.4:0.6（更偏镜框）
|
||
hf = np.where(_inner, 0.5 * hf + 0.5 * gz1f, hf)
|
||
# 强平滑：sigma=2.0, it=4 可微调，越大越平滑但可能损失细节
|
||
hf = np.where(_inner, self.smooth_valid(hf, _inner, sigma=2.0, it=4), hf)
|
||
# ==========================================
|
||
# 【结束新增】
|
||
# ==========================================
|
||
|
||
# ---- 后表面 ----
|
||
hb = np.full((H, W), np.nan)
|
||
real_f = have_f & (fz >= self.FACE_MIN)
|
||
hb[real_f] = fz[real_f] - self.EMBED_BOT
|
||
knownb = R & real_f
|
||
hb = self.laplace_fill(hb, knownb)
|
||
hb = self.smooth_valid(hb, R, sigma=2.0, it=2)
|
||
|
||
cap = hf - self.SHELL_MAX
|
||
ncap = int((R & ~real_f & np.isfinite(hb) & (hb < cap)).sum())
|
||
hb = np.where(R & ~real_f, np.maximum(hb, cap), hb)
|
||
self.log(f" 厚度上限(无真实脸面处, 距前面>{self.SHELL_MAX})生效格数: {ncap}")
|
||
|
||
hb = np.where(R & real_f, np.minimum(hb, fz - 0.25), hb)
|
||
bzs = self.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 + self.CLAMP_BACK), hb)
|
||
hb = np.minimum(hb, hf - self.GAP_MIN)
|
||
hb = self.smooth_valid(hb, R, sigma=2.0, it=2)
|
||
|
||
# ===== 在这里添加后表面钳制 =====
|
||
# 后表面下界（真实脸面 - 嵌入量）
|
||
floor_b = np.where(R & (real_f > -1e29), real_f - self.EMBED_BOT, -np.inf)
|
||
# 后表面上界（防穿出后脑，但允许贴脸）
|
||
ceil_b = np.where(R & (real_f > -1e29), real_f + self.CLAMP_BACK, np.inf)
|
||
hb = np.minimum(np.maximum(hb, floor_b), ceil_b)
|
||
# ===== 结束后表面钳制 =====
|
||
|
||
# ===== 墙区额外压低贴脸 =====
|
||
edge_in = self.erode_4(R, 4)
|
||
wall_zone = R & ~edge_in
|
||
hb = np.where(wall_zone, np.minimum(hb, real_f - self.EMBED_BOT * 0.8), hb)
|
||
# ===== 结束墙区压低 =====
|
||
|
||
# 外扩填充 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)
|
||
self.log(f" 填充厚度: min {np.nanmin(thickm):.2f}, 中位 {np.nanmedian(thickm):.2f}, "
|
||
f"mean {np.nanmean(thickm):.2f}, max {np.nanmax(thickm):.2f}")
|
||
|
||
# # ===== 新增：高度场平滑（三角化前）=====
|
||
# # 只对填充区域 R 内部做平滑，不影响边界
|
||
# from scipy.ndimage import gaussian_filter
|
||
# hf_smooth = hf.copy()
|
||
# hb_smooth = hb.copy()
|
||
|
||
# # 前表面平滑（sigma=1.0 对应 1 格半径，即约 0.1 单位的空间平滑）
|
||
# hf_smooth[R] = gaussian_filter(np.where(R, hf, 0), sigma=1.0, mode='reflect')[R]
|
||
# # 后表面平滑
|
||
# hb_smooth[R] = gaussian_filter(np.where(R, hb, 0), sigma=1.0, mode='reflect')[R]
|
||
|
||
# # 用平滑后的值替换
|
||
# hf = hf_smooth
|
||
# hb = hb_smooth
|
||
# # ===== 结束新增 =====
|
||
|
||
return hf, hb, frame_band, wall_zone, gz1f, bzs, real_f, fz
|
||
|
||
# ======================================================================
|
||
# 步骤 5: 边界提取与平滑
|
||
# ======================================================================
|
||
def extract_and_smooth_boundary(self, R, mask, g1f):
|
||
"""提取并平滑边界"""
|
||
self.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([[self.X0 + i * self.RES, self.Y0 + j * self.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]
|
||
self.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, self.BSTEP)
|
||
for _ in range(self.BITER):
|
||
B = (1 - self.BLAM) * B + self.BLAM * 0.5 * (np.roll(B, 1, 0) + np.roll(B, -1, 0))
|
||
|
||
# 面积恢复
|
||
areaR = R.sum() * self.RES * self.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, self.BSTEP)
|
||
r1 = roughness(B)
|
||
self.log(f" 边界粗糙度(平均转折角/标准差): 平滑前 {np.degrees(r0[0]):.2f}°/{np.degrees(r0[1]):.2f}° "
|
||
f"-> 平滑后 {np.degrees(r1[0]):.2f}°/{np.degrees(r1[1]):.2f}°")
|
||
|
||
# 边界须压在镜框带上
|
||
bad = ~self.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 self.inside_mask(q[None, :], mask)[0]:
|
||
B[k] = q
|
||
break
|
||
for _ in range(10):
|
||
B = (1 - self.BLAM) * B + self.BLAM * 0.5 * (
|
||
np.roll(B, 1, 0) + np.roll(B, -1, 0))
|
||
B = resample(B, self.BSTEP)
|
||
|
||
bii = np.clip(((B[:, 0] - self.X0) / self.RES).astype(int), 0, W - 1)
|
||
bjj = np.clip(((B[:, 1] - self.Y0) / self.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]
|
||
self.log(f" 平滑边界 {len(B)} 点, 超出足迹 {int((~self.inside_mask(B, mask)).sum())}, "
|
||
f"压在镜框带上 {nb.mean()*100:.1f}%")
|
||
|
||
return B, loop
|
||
|
||
# ======================================================================
|
||
# 步骤 6: 三角剖分 (耳切法 + 内部点插入)
|
||
# ======================================================================
|
||
def ear_clip(self, 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 = self.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 (self.cross3(a, b, p) >= -1e-12 and self.cross3(b, c, p) >= -1e-12
|
||
and self.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 = [(self.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]
|
||
self.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(self, B, R, hf, hb):
|
||
"""对填充区域进行三角剖分（含边界平滑，使侧墙圆滑）"""
|
||
self.log("== 6. 区域三角化 (耳切法 + 内部网格点插入, 无翻转风险) ==")
|
||
H, W = R.shape
|
||
Nb = len(B)
|
||
|
||
# ===== 【新增】边界环 Laplacian 平滑（仅平滑 XY，保持闭合） =====
|
||
B_orig = B.copy()
|
||
B_xy = B.copy()
|
||
if Nb > 3:
|
||
smooth_iter = 15
|
||
blend = 0.6
|
||
for _ in range(smooth_iter):
|
||
B_xy_new = 0.5 * B_xy + 0.25 * (np.roll(B_xy, 1, 0) + np.roll(B_xy, -1, 0))
|
||
B_xy = B_xy_new
|
||
# 混合：保留 40% 原始形状（防收缩），60% 平滑结果
|
||
B = (1 - blend) * B_orig + blend * B_xy
|
||
self.log(f" 边界平滑: iter={smooth_iter}, blend={blend}, 点数={Nb}")
|
||
# ================================================================
|
||
|
||
# 耳切三角化（使用平滑后的 B）
|
||
T = np.asarray(self.ear_clip(B.tolist()), np.int64)
|
||
self.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 / self.RES)))
|
||
cand = []
|
||
for j in range(0, H, step):
|
||
for i in range(0, W, step):
|
||
if inner2[j, i]:
|
||
cand.append((self.X0 + i * self.RES, self.Y0 + j * self.RES))
|
||
cand = np.asarray(cand)
|
||
self.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.7 * 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
|
||
self.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)
|
||
self.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})")
|
||
self.log(f" 实际插入 {extra} 点, 总三角形 {len(T)}, 最小有向面积 {area.min():.5f} (>0 = 无翻转)")
|
||
|
||
# 计算 z 坐标
|
||
zf = self.bilinear(self.smooth_valid(hf, R, sigma=1.5, it=2), xy)
|
||
zb = self.bilinear(hb, xy)
|
||
zb = np.minimum(zb, zf - 0.15)
|
||
zb = np.maximum(zb, zf - 50)
|
||
|
||
# ===== 【新增】边界 Z 协同平滑（防侧墙高度突变） =====
|
||
ring_idx = list(range(Nb))
|
||
zf_ring = zf[ring_idx].copy()
|
||
zb_ring = zb[ring_idx].copy()
|
||
for _ in range(8):
|
||
zf_ring = 0.55 * zf_ring + 0.225 * (np.roll(zf_ring, 1) + np.roll(zf_ring, -1))
|
||
zb_ring = 0.55 * zb_ring + 0.225 * (np.roll(zb_ring, 1) + np.roll(zb_ring, -1))
|
||
# 混合：保留 70% 原始 Z（保厚度），30% 平滑（消折角）
|
||
zf[ring_idx] = 0.7 * zf[ring_idx] + 0.3 * zf_ring
|
||
zb[ring_idx] = 0.7 * zb[ring_idx] + 0.3 * zb_ring
|
||
# ==========================================================
|
||
|
||
Pf = np.column_stack([xy, zf])
|
||
Pb = np.column_stack([xy, zb])
|
||
|
||
return Pf, Pb, T, xy
|
||
|
||
import numpy as np
|
||
def assemble_mesh(self, Pf, Pb, T):
|
||
nV = len(Pf)
|
||
faces = []
|
||
|
||
# 1. 前表面（CCW，法线朝前）
|
||
for a, b, c in T:
|
||
faces.append((a, b, c))
|
||
|
||
# 2. 后表面（CW，法线朝后）
|
||
for a, b, c in T:
|
||
faces.append((a + nV, c + nV, b + nV))
|
||
|
||
# 3. 提取边界环
|
||
from collections import Counter, defaultdict
|
||
edge_count = Counter()
|
||
for f in T:
|
||
for a, b in [(f[0], f[1]), (f[1], f[2]), (f[2], f[0])]:
|
||
key = (min(a, b), max(a, b))
|
||
edge_count[key] += 1
|
||
|
||
boundary_edges = [e for e, c in edge_count.items() if c == 1]
|
||
adj = defaultdict(list)
|
||
for a, b in boundary_edges:
|
||
adj[a].append(b)
|
||
adj[b].append(a)
|
||
|
||
ring = []
|
||
visited = set()
|
||
start = next(iter(adj))
|
||
cur = start
|
||
while True:
|
||
ring.append(cur)
|
||
visited.add(cur)
|
||
neighbors = [n for n in adj[cur] if n not in visited]
|
||
if not neighbors:
|
||
break
|
||
cur = neighbors[0]
|
||
|
||
Nr = len(ring)
|
||
|
||
# ========== 4. 计算曲率和外法线 ==========
|
||
curvature = np.zeros(Nr)
|
||
for i in range(Nr):
|
||
prev = ring[(i - 1) % Nr]
|
||
curr = ring[i]
|
||
next_pt = ring[(i + 1) % Nr]
|
||
|
||
v1 = Pf[curr] - Pf[prev]
|
||
v2 = Pf[next_pt] - Pf[curr]
|
||
v1[2] = 0
|
||
v2[2] = 0
|
||
|
||
len1 = np.linalg.norm(v1)
|
||
len2 = np.linalg.norm(v2)
|
||
|
||
if len1 > 1e-6 and len2 > 1e-6:
|
||
v1 /= len1
|
||
v2 /= len2
|
||
cross_z = v1[0] * v2[1] - v1[1] * v2[0]
|
||
dot = np.clip(v1 @ v2, -1, 1)
|
||
angle = np.arccos(dot)
|
||
curvature[i] = angle * (1 if cross_z > 0 else -1)
|
||
|
||
outward_list = []
|
||
for i in range(Nr):
|
||
curr = ring[i]
|
||
next_pt = ring[(i + 1) % Nr]
|
||
prev = ring[(i - 1) % Nr]
|
||
|
||
edge1 = Pf[curr] - Pf[prev]
|
||
edge2 = Pf[next_pt] - Pf[curr]
|
||
edge1[2] = 0
|
||
edge2[2] = 0
|
||
|
||
normal1 = np.array([-edge1[1], edge1[0], 0])
|
||
normal2 = np.array([-edge2[1], edge2[0], 0])
|
||
|
||
len1 = np.linalg.norm(normal1)
|
||
len2 = np.linalg.norm(normal2)
|
||
|
||
if len1 > 1e-6 and len2 > 1e-6:
|
||
outward = (normal1 / len1 + normal2 / len2) / 2
|
||
outward_norm = np.linalg.norm(outward)
|
||
if outward_norm > 1e-6:
|
||
outward /= outward_norm
|
||
else:
|
||
outward = np.array([0, 0, 1])
|
||
else:
|
||
outward = np.array([0, 0, 1])
|
||
|
||
outward_list.append(outward)
|
||
|
||
# ========== 5. 圆角参数 ==========
|
||
base_radius = 0.12
|
||
min_radius = 0.02
|
||
max_radius = 0.15
|
||
segments = 4
|
||
|
||
is_concave = np.zeros(Nr, dtype=bool)
|
||
is_convex = np.zeros(Nr, dtype=bool)
|
||
for i in range(Nr):
|
||
if curvature[i] < -0.05:
|
||
is_concave[i] = True
|
||
elif curvature[i] > 0.05:
|
||
is_convex[i] = True
|
||
|
||
local_radius = np.ones(Nr) * base_radius
|
||
for i in range(Nr):
|
||
if is_concave[i]:
|
||
local_radius[i] = max(min_radius, base_radius * 0.1)
|
||
elif is_convex[i]:
|
||
factor = min(1.3, 1.0 + curvature[i] * 0.5)
|
||
local_radius[i] = min(max_radius, base_radius * factor)
|
||
|
||
# ========== 6. 生成圆角环顶点 ==========
|
||
new_verts = []
|
||
new_base = 2 * nV
|
||
|
||
for i in range(Nr):
|
||
u = ring[i]
|
||
r = local_radius[i]
|
||
for k in range(1, segments + 1):
|
||
t = k / segments
|
||
angle = t * np.pi / 2
|
||
offset_h = r * np.sin(angle)
|
||
offset_v = r * (1 - np.cos(angle))
|
||
pt = Pf[u].copy()
|
||
|
||
pt[:2] -= outward_list[i][:2] * offset_h
|
||
pt[2] -= offset_v
|
||
new_verts.append(pt)
|
||
|
||
Pall = np.vstack([Pf, Pb] + new_verts)
|
||
|
||
def fillet_idx(i, k):
|
||
if k == 0:
|
||
return ring[i]
|
||
return new_base + i * segments + (k - 1)
|
||
|
||
# ========== 7. 圆角环三角形（修正：法线朝外）==========
|
||
for i in range(Nr):
|
||
j = (i + 1) % Nr
|
||
|
||
for k in range(segments):
|
||
a = fillet_idx(i, k)
|
||
b = fillet_idx(i, k + 1)
|
||
c = fillet_idx(j, k)
|
||
d = fillet_idx(j, k + 1)
|
||
|
||
# 关键修正：反转绕序
|
||
# 原来是 (a,c,b) 和 (b,c,d)，现在改为：
|
||
faces.append((a, b, c))
|
||
faces.append((b, d, c))
|
||
|
||
# ========== 8. 直墙（修正：连接圆角底层到对应的后表面位置）==========
|
||
for i in range(Nr):
|
||
j = (i + 1) % Nr
|
||
|
||
a = fillet_idx(i, segments)
|
||
b = fillet_idx(j, segments)
|
||
|
||
# 后表面也做同样的向内偏移
|
||
r_i = local_radius[i]
|
||
r_j = local_radius[j]
|
||
|
||
# 后表面对应点的偏移量（和圆角底层相同）
|
||
offset_h_i = r_i * np.sin(1.0 * np.pi / 2) # segments/segments = 1.0
|
||
offset_v_i = r_i * (1 - np.cos(1.0 * np.pi / 2))
|
||
offset_h_j = r_j * np.sin(1.0 * np.pi / 2)
|
||
offset_v_j = r_j * (1 - np.cos(1.0 * np.pi / 2))
|
||
|
||
a_back_pt = Pb[ring[i]].copy()
|
||
a_back_pt[:2] -= outward_list[i][:2] * offset_h_i
|
||
a_back_pt[2] -= offset_v_i
|
||
|
||
b_back_pt = Pb[ring[j]].copy()
|
||
b_back_pt[:2] -= outward_list[j][:2] * offset_h_j
|
||
b_back_pt[2] -= offset_v_j
|
||
|
||
# 添加后表面偏移后的顶点
|
||
a_back_idx = len(Pall)
|
||
b_back_idx = len(Pall) + 1
|
||
Pall = np.vstack([Pall, a_back_pt, b_back_pt])
|
||
|
||
# 法线朝外
|
||
faces.append((a, b_back_idx, b))
|
||
faces.append((a, a_back_idx, b_back_idx))
|
||
|
||
Ff = np.asarray(faces, np.int64)
|
||
|
||
# ========== 9. 诊断输出 ==========
|
||
front_part = Ff[:len(T)]
|
||
back_part = Ff[len(T):2*len(T)]
|
||
side_part = Ff[2*len(T):]
|
||
|
||
def calc_signed_vol(faces_sub):
|
||
v = 0.0
|
||
for f in faces_sub:
|
||
p0 = Pall[f[0]]
|
||
p1 = Pall[f[1]]
|
||
p2 = Pall[f[2]]
|
||
v += np.dot(p0, np.cross(p1, p2))
|
||
return v / 6.0
|
||
|
||
svol_front = calc_signed_vol(front_part)
|
||
svol_back = calc_signed_vol(back_part)
|
||
svol_side = calc_signed_vol(side_part)
|
||
svol_total = calc_signed_vol(Ff)
|
||
|
||
self.log(f" [诊断] 有向体积: 前={svol_front:.2f} 后={svol_back:.2f} 侧={svol_side:.2f} 总={svol_total:.2f}")
|
||
|
||
# 定向校正
|
||
if svol_total < 0:
|
||
Ff = Ff[:, [0, 2, 1]]
|
||
self.log(f" [诊断] 全局翻转面片方向")
|
||
svol_total = -svol_total
|
||
|
||
# ========== 10. 最终输出 ==========
|
||
vol = abs(svol_total)
|
||
n_side = len(Ff) - 2 * len(T)
|
||
self.log(f" 前表面: {len(T)} 三角形")
|
||
self.log(f" 后表面: {len(T)} 三角形")
|
||
self.log(f" 圆角环+直墙: {n_side} 三角形")
|
||
self.log(f" 总三角形: {len(Ff)}")
|
||
self.log(f" 体积: {vol:.2f}")
|
||
|
||
# 水密性检查
|
||
from collections import defaultdict as dd
|
||
edge_cnt = dd(list)
|
||
for fi, f in enumerate(Ff):
|
||
for j in range(3):
|
||
a, b = min(f[j], f[(j+1)%3]), max(f[j], f[(j+1)%3])
|
||
edge_cnt[(a, b)].append(fi)
|
||
nm = sum(1 for v in edge_cnt.values() if len(v) != 2)
|
||
self.log(f" 水密性: 非流形边 {nm}/{len(edge_cnt)}")
|
||
|
||
n_concave = np.sum(is_concave)
|
||
n_convex = np.sum(is_convex)
|
||
self.log(f" 凹凸分布: 内凹 {n_concave}, 外凸 {n_convex}, 平坦 {Nr - n_concave - n_convex}")
|
||
|
||
return Pall, Ff, vol
|
||
|
||
# ======================================================================
|
||
# 步骤 8: 校验
|
||
# ======================================================================
|
||
def validate_mesh(self, Pall, Ff, Pf, Pb, T, V, F, hf, hb, gz1f, frame_band,
|
||
wall_zone, bzs, have_b, real_f, fz, R, xy):
|
||
"""对生成的网格进行各项校验"""
|
||
self.log("== 8. 校验 ==")
|
||
|
||
# 定向校验
|
||
vfill = self.signed_volume(Pall, Ff)
|
||
vorig = self.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)))
|
||
|
||
self.log(f" 填充实体体积 {vfill:.2f} (期望 {exp_vol:.2f}), 原模型 {vorig:.0f}, 定向一致")
|
||
self.log(f" 体积: {vfill:.2f} (期望 {exp_vol:.2f})")
|
||
|
||
if abs(vfill - exp_vol) > 0.5 * 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)
|
||
self.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)
|
||
self.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]
|
||
self.log(f" 前表面低于本地框面(墙区): 最小余量 {np.nanmin(gmargin):.3f} (应>0)")
|
||
gmargin2 = (gz1f - hf)[R & frame_band]
|
||
self.log(f" 前表面低于本地框面(内部允许贴齐): 最小余量 {np.nanmin(gmargin2):.3f} (应>-0.10)")
|
||
bmargin = (hb - bzs)[R & have_b]
|
||
self.log(f" 后表面高于最后表面: 最小余量 {np.nanmin(bmargin):.3f} (应>0)")
|
||
embed = (fz - hb)[R & real_f]
|
||
self.log(f" 底面嵌入真实脸面: 最小 {np.nanmin(embed):.3f} (应>-0.05, 负=缝隙), "
|
||
f"中位 {np.nanmedian(embed):.2f}")
|
||
inR = self.inside_mask(xy, R)
|
||
self.log(f" 圆盘顶点落在 R 内比例 {inR.mean()*100:.1f}%")
|
||
|
||
return Ff, vfill
|
||
|
||
# ======================================================================
|
||
# 步骤 9: 输出
|
||
# ======================================================================
|
||
def output_mesh(self, Pall, Ff, V, F):
|
||
"""输出填充网格和合并后的模型"""
|
||
self.log("== 9. 输出 ==")
|
||
|
||
# 输出填充网格
|
||
with open(self.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")
|
||
self.log(f" {self.OUT_FILL}: {len(Pall)} 顶点, {len(Ff)} 面")
|
||
|
||
# 输出合并模型
|
||
with open(self.IN_OBJ) as f:
|
||
orig = f.read()
|
||
with open(self.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")
|
||
self.log(f" {self.OUT_MERGED}: 原模型 + 填充 (顶点偏移 {off})")
|
||
|
||
# ======================================================================
|
||
# 步骤 10: 生成预览图
|
||
# ======================================================================
|
||
def generate_preview(self, gz1, hf, hb, R, B):
|
||
"""生成前视图深度预览"""
|
||
self.log("== 10. 生成预览图 ==")
|
||
H, W = self.H, self.W
|
||
|
||
# 合成前表面深度图
|
||
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] - self.X0) / self.RES * 4, (self.Y1 - b[1]) / self.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(self.OUT_PREVIEW)
|
||
self.log(f" {self.OUT_PREVIEW}: 3 面板(前深度/填充前表面/厚度), 红线=平滑边界")
|
||
|
||
# ======================================================================
|
||
# 主流程 run()
|
||
# ======================================================================
|
||
def run(self):
|
||
"""执行完整填充网格生成流程"""
|
||
# 步骤 1: 加载模型
|
||
V, F, face_g = self.load_and_classify()
|
||
|
||
# 步骤 2: 栅格化
|
||
gz1, gz2, fz, bz, have_g1, have_g2, have_f, have_b, g1f, W, H = \
|
||
self.rasterize(V, F, face_g)
|
||
self.W, self.H = W, H
|
||
|
||
# 步骤 3: 构建足迹
|
||
mask, R = self.build_footprint(g1f)
|
||
|
||
# 步骤 4: 构造高度场
|
||
hf, hb, frame_band, wall_zone, gz1f, bzs, real_f, fz = \
|
||
self.build_height_fields(gz1, gz2, fz, bz, have_g1, have_f, have_b,
|
||
g1f, mask, R)
|
||
|
||
# 步骤 5: 边界提取与平滑
|
||
B, loop = self.extract_and_smooth_boundary(R, mask, g1f)
|
||
|
||
# 步骤 6: 三角剖分
|
||
Pf, Pb, T, xy = self.triangulate_fill(B, R, hf, hb)
|
||
|
||
# 步骤 7: 组装实体
|
||
Pall, Ff, nV = self.assemble_mesh(Pf, Pb, T)
|
||
|
||
# 步骤 8: 校验
|
||
Ff, vfill = self.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: 输出
|
||
self.output_mesh(Pall, Ff, V, F)
|
||
|
||
# 步骤 10: 预览图
|
||
self.generate_preview(gz1, hf, hb, R, B)
|
||
|
||
# # 自动调用 Blender 做倒角
|
||
# success = auto_bevel_in_blender(
|
||
# input_obj_path="fill_mesh.obj",
|
||
# output_obj_path="./final_smooth.obj",
|
||
# bevel_width=1.0 # 倒角宽度，可调
|
||
# )
|
||
|
||
# if success:
|
||
# print("✅ 最终平滑模型已生成: final_smooth.obj")
|
||
# else:
|
||
# print("⚠️ 倒角失败，请检查 Blender 安装")
|
||
|
||
self.log("\n== 完成 ==")
|
||
with open("fill_report.txt", "w") as f:
|
||
f.write("\n".join(self.report) + "\n")
|
||
print("报告已写入 fill_report.txt")
|
||
|
||
|
||
# ======================================================================
|
||
# 诊断工具
|
||
# ======================================================================
|
||
def diagnose_indices(obj_path, idx_path):
|
||
"""诊断索引文件是否匹配模型"""
|
||
gen = FillMeshGenerator()
|
||
V, F = gen.parse_obj(obj_path)
|
||
idx = np.loadtxt(idx_path, dtype=np.int64) - 1
|
||
|
||
# 检查索引范围
|
||
print(f"顶点总数: {len(V)}")
|
||
print(f"索引总数: {len(idx)}")
|
||
print(f"索引范围: {idx.min()} ~ {idx.max()}")
|
||
print(f"超出范围的索引: {(idx >= len(V)).sum()}")
|
||
|
||
# 检查眼镜顶点坐标范围
|
||
glasses_V = V[idx]
|
||
print(f"\n眼镜顶点坐标范围:")
|
||
print(f" x: [{glasses_V[:,0].min():.2f}, {glasses_V[:,0].max():.2f}]")
|
||
print(f" y: [{glasses_V[:,1].min():.2f}, {glasses_V[:,1].max():.2f}]")
|
||
print(f" z: [{glasses_V[:,2].min():.2f}, {glasses_V[:,2].max():.2f}]")
|
||
|
||
# 检查是否在栅格化窗口内
|
||
gen2 = FillMeshGenerator()
|
||
in_window = (
|
||
(glasses_V[:,0] > gen2.X0) & (glasses_V[:,0] < gen2.X1) &
|
||
(glasses_V[:,1] > gen2.Y0) & (glasses_V[:,1] < gen2.Y1)
|
||
)
|
||
print(f"\n在栅格化窗口内的顶点: {in_window.sum()} / {len(idx)}")
|
||
print(f" 窗口范围: x[{gen2.X0}, {gen2.X1}], y[{gen2.Y0}, {gen2.Y1}]")
|
||
|
||
# 检查眼镜面数
|
||
ism = np.zeros(len(V), bool)
|
||
ism[idx] = True
|
||
face_g = ism[F].any(1)
|
||
print(f"\n眼镜面数: {face_g.sum()} / {len(F)}")
|
||
|
||
# ======================================================================
|
||
# 入口
|
||
# ======================================================================
|
||
if __name__ == "__main__":
|
||
parser = argparse.ArgumentParser(description="眼镜-脸部填充网格生成器")
|
||
parser.add_argument("--in_obj", type=str, default="632723.obj",
|
||
help="输入 OBJ 模型文件路径 (默认: 632723.obj)")
|
||
parser.add_argument("--in_idx", type=str, default="index.txt",
|
||
help="眼镜顶点索引文件路径 (默认: index.txt)")
|
||
parser.add_argument("--out_fill", type=str, default="fill_mesh.obj",
|
||
help="输出填充网格文件路径 (默认: fill_mesh.obj)")
|
||
parser.add_argument("--diagnose", action="store_true",
|
||
help="仅运行诊断，不生成填充网格")
|
||
|
||
args = parser.parse_args()
|
||
|
||
if args.diagnose:
|
||
diagnose_indices(args.in_obj, args.in_idx)
|
||
else:
|
||
gen = FillMeshGenerator(
|
||
in_obj=args.in_obj,
|
||
in_idx=args.in_idx,
|
||
out_fill=args.out_fill
|
||
)
|
||
gen.run()
|
||
|
||
# python generate_fill_mesh_hsc.py --in_obj 632723.obj --in_idx index.txt --out_fill fill_mesh.obj |