1249 lines
53 KiB
Python
1249 lines
53 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
|
|
|
|
|
|
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 = -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 收敛阈值
|
|
|
|
def __init__(self, in_obj=None, in_idx=None):
|
|
if in_obj:
|
|
self.IN_OBJ = in_obj
|
|
if in_idx:
|
|
self.IN_IDX = in_idx
|
|
self.report = []
|
|
|
|
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
|
|
|
|
# ======================================================================
|
|
# 步骤 3: 构建足迹 mask 与填充区域 R
|
|
# ======================================================================
|
|
def build_footprint(self, g1f):
|
|
"""构建眼镜足迹 mask 和填充区域 R"""
|
|
self.log("== 3. 眼镜足迹 / 填充区域 R ==")
|
|
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} 格误标鼻面")
|
|
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} (触碰=侧墙切在头内部, 正常)")
|
|
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=1.5)
|
|
|
|
# 墙区安全钳制
|
|
edge_in = self.erode_4(R, 4)
|
|
wall_zone = R & ~edge_in
|
|
hf = np.where(frame_band & wall_zone, np.minimum(hf, gz1f - self.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]), (self.X0 + ja[1] * self.RES, self.Y0 + ja[0] * self.RES)
|
|
self.log(f" 交界处高度跳变最大 {jn:.3f} @ {jn_at} (平面前框带接缝应<0.3; "
|
|
"镜圈圆角处~0.4为镜框自身几何, 并集后不暴露)")
|
|
|
|
# ---- 后表面 ----
|
|
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)
|
|
|
|
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)
|
|
|
|
# 外扩填充 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}")
|
|
|
|
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)
|
|
|
|
# 耳切三角化
|
|
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.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
|
|
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(hf, xy)
|
|
zb = self.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(self, Pf, Pb, T):
|
|
"""组装填充实体: 前表面 + 后表面 + 侧墙"""
|
|
self.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(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}, 定向一致")
|
|
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)
|
|
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)
|
|
|
|
self.log("\n== 完成 ==")
|
|
with open("fill_report.txt", "w") as f:
|
|
f.write("\n".join(self.report) + "\n")
|
|
print("报告已写入 fill_report.txt")
|
|
|
|
|
|
# ======================================================================
|
|
# 入口
|
|
# ======================================================================
|
|
if __name__ == "__main__":
|
|
gen = FillMeshGenerator()
|
|
gen.run() |