openMVS/fill2sc/vacuum_wrap.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
vacuum_wrap.py —— 眼镜“抽真空包装” mesh 生成脚本
================================================
参考 youshui_result.obj 中新增的 o Transparent 部分(一个完全包裹眼镜外轮廓、
并填充镜片镂空处的水密闭合 mesh), 本脚本从 youshui.obj + index.txt 出发重新生成它。

输入/输出(自动探测, 也可用 --obj/--index/--out 显式指定):
    模型  本目录最大的非输出 *.obj(排除 *result*/*wrap*/*_fill*)
    索引  模型同目录的 index.txt(眼镜顶点索引, 每行一个, 0-based)
    输出  模型同目录 <模型名>_wrap_result.obj = 原模型 + 追加 o Transparent
    处理其它数据集时, 把本脚本拷到该数据目录下直接运行即可。

思路(“抽真空包装”):
    1. 由 index.txt 提取眼镜子网格(三个顶点都在索引里的面);
    2. 判定正面视线轴后, 只取“朝前的面”(法线与视线轴余弦 > 阈值)投影
       栅格化, 形态学闭运算补缝, 取含镜片镂空的最大连通组件,
       内部封闭镂空自动填满 → 得到镜框正面“外轮廓”多边形
       (镜片镂空被盖住, 鼻梁处的凹口保留; 向后弯折的镜腿被排除,
       包装不超出眼镜正面外框);
    3. 在外轮廓内铺一层网格点, 作为“包装薄膜”的采样点;
    4. 每个网格点从正前方贴到眼镜表面(取该处最靠前的表面高度);
       镜片镂空处没有表面 → 用周边镜框高度拟合二次趋势面外推
       (跟随镜框倾斜/缓弯, 平齐不下陷), 残差在孔内做调和衔接,
       平滑无台阶; footprint 按镜架整体投影剪裁, 不超出镜架;
    5. 薄膜沿视线后退 thickness 得到背面膜(最深不超过镜架最后缘), 周边侧壁缝合
       → 一个水密(watertight)闭合的包装 mesh;
    6. 棋盘式交替对角线剖分; 薄膜边缘沿边界曲线做一维拉普拉斯收边,
       磨掉 footprint 栅格化的锯齿轮廓;
    7. 原样拷贝原 obj(仅修正 mtllib/usemtl 指向本目录的 youshui.mtl),
       把新 mesh 追加在文件末尾(o Transparent / usemtl Transparent)。

只依赖 numpy; 所有读写都限制在本脚本所在目录。

与参考结果 youshui_result.obj 的对比: 结构一致(原模型 + o Transparent 水密闭合
包装 mesh); 本脚本只包装镜框正面核心(含镜片镂空填充), 不包裹向后弯折的镜腿,
因此与参考一样不超出眼镜正面外框、z 向深度一致; 参考在镜腿后另有一小片
手工处理的收拢薄膜, 此处不复现。

用法:
    python3 vacuum_wrap.py
    python3 vacuum_wrap.py --spacing 0.06 --thickness 3.4 --cell 0.05 --front 0.25
"""

import os
import glob
import argparse
import numpy as np

SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))


def resolve(p):
    """路径一律相对本脚本目录解析, 不碰其它目录。"""
    return p if os.path.isabs(p) else os.path.join(SCRIPT_DIR, p)


# ----------------------------------------------------------------- OBJ 读取
def load_obj(path):
    """读取顶点与三角面(多于三个顶点的面做三角扇拆分)。"""
    verts, faces = [], []
    with open(path, 'r', errors='ignore') as f:
        for line in f:
            if line.startswith('v '):
                p = line.split()
                verts.append((float(p[1]), float(p[2]), float(p[3])))
            elif line.startswith('f '):
                ids = [int(x.split('/')[0]) - 1 for x in line.split()[1:]]
                for k in range(1, len(ids) - 1):
                    faces.append((ids[0], ids[k], ids[k + 1]))
    return np.asarray(verts, np.float64), np.asarray(faces, np.int64)


def load_index(path, nverts):
    idx = []
    with open(path) as f:
        for line in f:
            line = line.strip()
            if line:
                idx.append(int(line))
    idx = np.asarray(idx, np.int64)
    if len(idx) == 0:
        raise RuntimeError('index.txt 为空')
    if idx.min() < 0 or idx.max() >= nverts:
        raise RuntimeError('index.txt 中索引越界(0..%d)' % (nverts - 1))
    return idx


# ------------------------------------------------------------ 正面朝向判定
def pick_front_axis(verts, sub, wind, coarse):
    """判定“抽真空”的视线轴: 沿该轴投影时, 眼镜区域内部会出现封闭的
    镜片镂空(孔洞), 孔洞面积最大的轴即包装轴(如从正面看眼镜)。
    wind: 整mesh体积符号(+1=法线朝外), 用于确定正面在轴的正/负侧。
    返回 (k, sgn, u轴, v轴), (u, v, w=sgn*第k轴) 构成右手系。"""
    p0, p1, p2 = verts[sub[:, 0]], verts[sub[:, 1]], verts[sub[:, 2]]
    n = np.cross(p1 - p0, p2 - p0)
    sgn_area = wind * n.sum(axis=0)          # 各轴带符号投影面积
    best_hole, best_k = -1.0, 2
    rc = max(2, int(round(0.6 / coarse)))    # 较大闭运算: 桥接处覆盖稀疏时也能封住镂空
    for k in range(3):
        u, v = (k + 1) % 3, (k + 2) % 3
        tri_uv = verts[sub][:, :, [u, v]]
        tri_w = verts[sub][:, :, k]
        cover, _, _ = rasterize(tri_uv, tri_w, coarse)
        foot = fill_enclosed(binary_closing(cover, rc))
        hole = float((foot & ~cover).sum()) * coarse * coarse
        print('    投影轴 %s: 表面覆盖 %.1f, 封闭镂空 %.1f'
              % ('xyz'[k], cover.sum() * coarse * coarse, hole))
        if hole > best_hole:
            best_hole, best_k = hole, k
    k = best_k
    sgn = 1.0 if sgn_area[k] >= 0 else -1.0
    return k, sgn, (k + 1) % 3, (k + 2) % 3


# ------------------------------------------------------------------ 栅格化
def rasterize(tri_uv, tri_w, h, keep=None, collect=False):
    """把子网格三角形投影到 (u,v) 平面栅格化。
    keep: 可选的三角形布尔掩码(如只保留朝前的面)。
    collect: 额外收集 每格覆盖的三角形编号表(用于后续精确求高)。
    返回 cover(是否被覆盖) 与 topz(该格内最靠前的表面高度 w);
    collect=True 时再附带 cell_tris。"""
    umin, vmin = tri_uv[..., 0].min(), tri_uv[..., 1].min()
    umax, vmax = tri_uv[..., 0].max(), tri_uv[..., 1].max()
    # 留出边距(>= 闭运算半径+2), 保证膨胀/洪泛不会越过网格边界
    mg = 8
    umin -= mg * h; vmin -= mg * h
    nu = int(np.ceil((umax - umin) / h)) + 1 + mg
    nv = int(np.ceil((vmax - vmin) / h)) + 1 + mg
    if keep is not None:                       # 网格 bbox 一律取自全集, 保证可对齐
        tri_uv = tri_uv[keep]
        tri_w = tri_w[keep]
    nu -= nu % 2; nv -= nv % 2          # 偶数尺寸, 红黑迭代对齐
    cover = np.zeros((nu, nv), bool)
    topz = np.full((nu, nv), np.nan)
    cell_tris = {} if collect else None
    ntri = tri_uv.shape[0]
    for t in range(ntri):
        a = tri_uv[t, 0]; b = tri_uv[t, 1]; c = tri_uv[t, 2]
        wa, wb, wc = tri_w[t]
        i0 = max(0, int(np.floor((min(a[0], b[0], c[0]) - umin) / h)) - 1)
        i1 = min(nu - 1, int(np.floor((max(a[0], b[0], c[0]) - umin) / h)) + 1)
        j0 = max(0, int(np.floor((min(a[1], b[1], c[1]) - vmin) / h)) - 1)
        j1 = min(nv - 1, int(np.floor((max(a[1], b[1], c[1]) - vmin) / h)) + 1)
        if i1 < i0 or j1 < j0:
            continue
        cx = umin + (np.arange(i0, i1 + 1) + 0.5) * h
        cy = vmin + (np.arange(j0, j1 + 1) + 0.5) * h
        CX, CY = np.meshgrid(cx, cy, indexing='ij')
        v0 = c - a
        v1 = b - a
        d00 = v0 @ v0; d01 = v0 @ v1; d11 = v1 @ v1
        den = d00 * d11 - d01 * d01
        if abs(den) < 1e-14:            # 投影退化(近似垂直面), 跳过
            continue
        px, py = CX - a[0], CY - a[1]
        d20 = px * v0[0] + py * v0[1]
        d21 = px * v1[0] + py * v1[1]
        bu = (d11 * d20 - d01 * d21) / den
        bv = (d00 * d21 - d01 * d20) / den
        ins = (bu >= -1e-9) & (bv >= -1e-9) & (bu + bv <= 1 + 1e-9)
        if not ins.any():
            continue
        z = wa + bu * (wb - wa) + bv * (wc - wa)
        sub = topz[i0:i1 + 1, j0:j1 + 1]
        topz[i0:i1 + 1, j0:j1 + 1] = np.where(ins, np.fmax(sub, z), sub)
        cover[i0:i1 + 1, j0:j1 + 1] |= ins
        if collect:
            ii, jj = np.nonzero(ins)
            for ci_, cj_ in zip(ii, jj):
                cell_tris.setdefault((int(ci_) + i0, int(cj_) + j0), []).append(t)
    if collect:
        return cover, topz, (umin, vmin), cell_tris
    return cover, topz, (umin, vmin)


# ------------------------------------------------------- 精确薄膜高度
def exact_film_heights(pts, tri_uv, tri_w, cell_tris, u0, v0, h):
    """每个薄膜网格点沿视线轴与子网格求精确交点, 取最靠前的高度。
    相比栅格最近格值, 没有分段常数的阶梯量化, 薄膜贴合更光滑。"""
    n = len(pts)
    Z = np.full(n, np.nan)
    for p in range(n):
        x, y = pts[p, 0], pts[p, 1]
        i = int((x - u0) / h)
        j = int((y - v0) / h)
        ids = cell_tris.get((i, j))
        if not ids:
            continue
        best = -np.inf
        for t in ids:
            a = tri_uv[t, 0]; b = tri_uv[t, 1]; c = tri_uv[t, 2]
            e0 = c - a
            e1 = b - a
            d00 = e0 @ e0; d01 = e0 @ e1; d11 = e1 @ e1
            den = d00 * d11 - d01 * d01
            if abs(den) < 1e-14:
                continue
            d20 = (x - a[0]) * e0[0] + (y - a[1]) * e0[1]
            d21 = (x - a[0]) * e1[0] + (y - a[1]) * e1[1]
            bu = (d11 * d20 - d01 * d21) / den
            bv = (d00 * d21 - d01 * d20) / den
            if bu >= -1e-9 and bv >= -1e-9 and bu + bv <= 1 + 1e-9:
                wa, wb, wc = tri_w[t]
                z = wa + bu * (wb - wa) + bv * (wc - wa)
                if z > best:
                    best = z
        if best > -np.inf:
            Z[p] = best
    return Z


# ---------------------------------------------------------------- 形态学
def dilate(m):
    r = m.copy()
    r[1:, :] |= m[:-1, :]
    r[:-1, :] |= m[1:, :]
    r[:, 1:] |= m[:, :-1]
    r[:, :-1] |= m[:, 1:]
    return r


def erode(m):
    r = m.copy()
    r[1:, :] &= m[:-1, :]
    r[:-1, :] &= m[1:, :]
    r[:, 1:] &= m[:, :-1]
    r[:, :-1] &= m[:, 1:]
    return r


def binary_closing(mask, r):
    if r <= 0:
        return mask.copy()
    m = np.pad(mask, r, constant_values=False)   # pad 足够大, 膨胀不会越界
    for _ in range(r):
        m = dilate(m)
    for _ in range(r):
        m = erode(m)
    return m[r:-r, r:-r]


def fill_enclosed(mask):
    """从外部洪泛, 补上内部封闭的镂空(镜片孔洞), 开向边界的凹口(鼻梁)保留。"""
    outside = np.zeros_like(mask)
    outside[0, :] = ~mask[0, :]
    outside[-1, :] = ~mask[-1, :]
    outside[:, 0] = ~mask[:, 0]
    outside[:, -1] = ~mask[:, -1]
    while True:
        new = dilate(outside) & ~mask
        if np.array_equal(new, outside):
            break
        outside = new
    return ~outside


def largest_filled_component(mask):
    """先补内部封闭孔洞, 再取面积最大的连通组件(含镜片镂空的镜框核心),
    丢掉镜腿等散碎覆盖 → 包装不会越过镜框正面外框。"""
    filled = fill_enclosed(mask)
    BIG = filled.size + 1
    lab = np.where(filled, np.arange(filled.size).reshape(filled.shape), BIG)
    while True:                                  # 4邻域最小标签传播 → 连通组件标记
        new = lab.copy()
        for (s1, s2) in (((slice(1, None), slice(None)), (slice(None, -1), slice(None))),
                         ((slice(None, -1), slice(None)), (slice(1, None), slice(None))),
                         ((slice(None), slice(1, None)), (slice(None), slice(None, -1))),
                         ((slice(None), slice(None, -1)), (slice(None), slice(1, None)))):
            cand = np.full_like(lab, BIG)
            cand[s1] = lab[s2]
            cand[~filled] = BIG                 # 标签只在组件内传播, 不漏进空格
            np.minimum(new, cand, out=new)
        if np.array_equal(new, lab):
            break
        lab = new
    flat = lab.ravel()
    labels = flat[flat < BIG]
    if len(labels) == 0:
        raise RuntimeError('没有可用连通组件')
    best_lab, best_area = None, -1
    for L in np.unique(labels):
        a = int((flat == L).sum())
        if a > best_area:
            best_area, best_lab = a, L
    return lab == best_lab


# ----------------------------------------------------------- 外轮廓追踪
def trace_all_cycles(mask):
    """沿 cell 边缘追踪 mask 的所有边界环(内部在行进方向左侧)。"""
    nu, nv = mask.shape
    edges = {}

    def add(a, b):
        edges.setdefault(a, []).append(b)

    I, J = np.nonzero(mask)
    for i, j in zip(I, J):
        if i == 0 or not mask[i - 1, j]:          # 左侧外部 → 向 -v 走
            add((i, j + 1), (i, j))
        if i == nu - 1 or not mask[i + 1, j]:     # 右侧外部 → 向 +v 走
            add((i + 1, j), (i + 1, j + 1))
        if j == 0 or not mask[i, j - 1]:          # 下侧外部 → 向 +u 走
            add((i, j), (i + 1, j))
        if j == nv - 1 or not mask[i, j + 1]:     # 上侧外部 → 向 -u 走
            add((i + 1, j + 1), (i, j + 1))

    cycles, visited = [], set()
    for start in edges:
        if start in visited:
            continue
        cyc = [start]; visited.add(start); cur = start
        while True:
            nxts = [e for e in edges.get(cur, ()) if e not in visited]
            if not nxts:
                break
            nxt = nxts[0]
            cyc.append(nxt); visited.add(nxt); cur = nxt
        cycles.append(cyc)
    return cycles


def signed_area(cyc):
    a = 0.0
    for k in range(len(cyc)):
        x1, y1 = cyc[k]; x2, y2 = cyc[(k + 1) % len(cyc)]
        a += x1 * y2 - x2 * y1
    return a / 2.0


def trace_outer_cycle(mask):
    """返回 (最大面积的边界环, 显著环个数)。
    环的方向统一为逆时针(内部在左侧)。"""
    cycles = trace_all_cycles(mask)
    if not cycles:
        raise RuntimeError('外轮廓追踪失败: mask 无边界')
    cycles.sort(key=lambda c: abs(signed_area(c)), reverse=True)
    main = cycles[0]
    if signed_area(main) < 0:
        main = main[::-1]
    ref = abs(signed_area(main))
    n_sig = 1 + sum(1 for c in cycles[1:]
                    if abs(signed_area(c)) > 0.02 * ref)
    return main, n_sig


# -------------------------------------------------------------- 点在多边形内
def points_in_poly(pts, poly):
    """even-odd 规则, pts:(N,2), poly:(K,2) 隐式闭合。"""
    x, y = pts[:, 0], pts[:, 1]
    inside = np.zeros(len(pts), bool)
    x1, y1 = poly[:, 0], poly[:, 1]
    x2, y2 = np.roll(x1, -1), np.roll(y1, -1)
    with np.errstate(divide='ignore', invalid='ignore'):
        for k in range(len(poly)):
            cond = (y1[k] > y) != (y2[k] > y)
            xint = (x2[k] - x1[k]) * (y - y1[k]) / (y2[k] - y1[k]) + x1[k]
            inside ^= cond & (x < xint)
    return inside


# --------------------------------------------- 镂空处高度场(与镜框平齐)
def _neighbor_avg4(Z, ok_mask):
    """4 邻域均值: 只统计 ok_mask 为真且值非 nan 的格子。
    返回 (均值, 计数)。"""
    acc = np.zeros(Z.shape)
    cnt = np.zeros(Z.shape)
    nanZ = np.isnan(Z)
    for sl_src, sl_dst in ((slice(None, -1), slice(1, None)),
                           (slice(1, None), slice(None, -1))):
        have = ok_mask[sl_src, :] & ~nanZ[sl_src, :]
        acc[sl_dst, :] += np.where(have, Z[sl_src, :], 0.0)
        cnt[sl_dst, :] += have
        have = ok_mask[:, sl_src] & ~nanZ[:, sl_src]
        acc[:, sl_dst] += np.where(have, Z[:, sl_src], 0.0)
        cnt[:, sl_dst] += have
    return np.where(cnt > 0, acc / np.maximum(cnt, 1), np.nan), cnt


def _fit_surface(di, dj, z):
    """最小二乘拟合 z ≈ a + b·di + c·dj (平面趋势)。

    采样点只分布在镂空周围的环带上: 二次趋势面对这种采样是病态
    外推(环带上常数项与二次项近似共线, 中心值可飘出 ±2), 平面
    趋势则稳定; 曲率起伏留给残差调和项(边界值插值, 良定)。
    点太少/退化时回退到常数。返回 (系数, 项数)。"""
    n = len(z)
    di = di.astype(np.float64)
    dj = dj.astype(np.float64)
    if n >= 7:
        A = np.stack([np.ones(n), di, dj], 1)
        G = A.T @ A
        G[np.diag_indices_from(G)] += 1e-6 * max(1.0, G.diagonal().max())
        try:
            return np.linalg.solve(G, A.T @ z), 3
        except np.linalg.LinAlgError:
            pass
    return np.array([float(np.mean(z))]), 1


def _eval_surface(coef, nterm, di, dj):
    di = np.asarray(di, np.float64)
    dj = np.asarray(dj, np.float64)
    if nterm == 3:
        return coef[0] + coef[1] * di + coef[2] * dj
    return np.full(di.shape, coef[0], float)


def trend_fill(Z, valid, select=None, blend_passes=500, omega=0.8):
    """镜片镂空处高度场填充: 全局平面趋势 + 残差调和, 平齐且平滑。

    薄膜要"平齐"(跟随镜框正面的倾斜)也要平滑(无台阶无尖峰):
    1) 用正面表面(select)上的贴合点拟合平面趋势 T(镜框正面倾角)。
       采样点铺满整片 footprint, 平面拟合稳定(二次面在只环绕
       镜片镂空的环带采样下是病态外推, 中心值可飘出 ±2);
    2) 残差 R = 贴合高度 - T 在**所有**贴合点锚定(含外缘陡边/
       镜托/镜腿根处的深锚点——薄膜要包住它们), 镂空处解
       Laplace 方程: 调和解无内部极值、不会过冲, 天然平滑;
       镜片中央 ≈ 趋势面(平齐), 孔缘精确衔接, 深锚点周围形成
       光滑下凹而不是尖 crater。
    blend_passes 为残差调和迭代次数; 迭代用欠松弛(ω≤1)——
    超松弛会放大棋盘模式留下二阶差分尖峰。"""
    Z = Z.copy()
    unknown = valid & np.isnan(Z)
    if not unknown.any():
        return Z
    known = valid & ~np.isnan(Z)
    # 1) 平面趋势: 只用正面表面高度, 稳健剔除离群锚点
    K = known if select is None else (known & select)
    if K.sum() >= 7:
        ki, kj = np.nonzero(K)
        for _ in range(2):
            coef, nterm = _fit_surface(ki - ki.mean(), kj - kj.mean(),
                                       Z[ki, kj])
            if len(ki) < 16:
                break
            res = Z[ki, kj] - _eval_surface(
                coef, nterm, ki - ki.mean(), kj - kj.mean())
            keep = np.abs(res) <= max(0.12, 4.0 * np.median(np.abs(res)))
            if keep.all() or keep.sum() < 8:
                break
            ki, kj = ki[keep], kj[keep]
        ci, cj = ki.mean(), kj.mean()
    elif K.sum() > 0:
        coef, nterm = np.array([float(np.mean(Z[K]))]), 1
        ci = cj = 0.0
    else:
        coef, nterm = np.array([float(np.nanmean(Z[known]))]), 1
        ci = cj = 0.0
    ii, jj = np.indices(Z.shape)
    T = _eval_surface(coef, nterm, ii - ci, jj - cj)
    # 2) 残差调和填充: 全体镂空处一起解, 锚定所有贴合点
    R = np.where(known, Z - T, 0.0)
    for _ in range(blend_passes):
        avg, cnt = _neighbor_avg4(R, valid)
        upd = unknown & (cnt > 0)
        if not upd.any():
            break
        R[upd] += omega * (avg[upd] - R[upd])
    Z[unknown] = T[unknown] + R[unknown]
    return Z


# ------------------------------------------------------------ 贴合带平滑
def smooth_adherence(Z, known, valid, thr, spacing,
                     alpha=0.6, passes=60, near=0.3, near_w=0.4):
    """贴合带受限平滑(消除镜框陡坡处高度场的棱面)。

    高度场 z(u,v) 表示坡度为 s 的表面时, 其曲率会被放大 ~s^2 倍,
    镜框侧弯坡度大处逐格贴合会出现 30-60 度的相邻面法线跳变。
    本函数把与 8 邻域均值偏差超过 thr 的点向邻域均值加权松弛:
    - 贴合点(known) 权重 1: 磨掉曲率突变棱面, 薄膜在陡坡处"绷而不贴";
      坡度平缓处/均匀斜面的偏差本来就小于 thr, 不会被改动;
    - 绷紧区离贴合边界 near 以内的点权重 near_w: 轻推, 软化接缝褶皱;
    - 更远处不参与(保持真空 Laplace 解)。
    thr 的量纲是高度: 相邻面法线跳变约 atan(4*thr/spacing)。"""
    band = known.copy()
    for _ in range(max(1, int(round(near / spacing)))):
        band = dilate(band)
    w = np.where(known, 1.0, np.where(band & valid, near_w, 0.0))
    Z = Z.copy()
    ni, nj = Z.shape
    for _ in range(passes):
        Zp = np.pad(Z, 1, constant_values=np.nan)
        Vp = np.pad(valid, 1, constant_values=False)
        acc = np.zeros(Z.shape)
        cnt = np.zeros(Z.shape)
        for di in (-1, 0, 1):
            for dj in (-1, 0, 1):
                if di == 0 and dj == 0:
                    continue
                sl = (slice(1 + di, 1 + di + ni), slice(1 + dj, 1 + dj + nj))
                nb = np.where(Vp[sl], Zp[sl], np.nan)
                ok = ~np.isnan(nb)
                acc += np.where(ok, nb, 0.0)
                cnt += ok
        avg = np.where(cnt > 0, acc / np.maximum(cnt, 1), np.nan)
        res = np.abs(Z - avg)
        upd = (w > 0) & (cnt > 0) & (res > thr)
        if not upd.any():
            break
        Z[upd] += w[upd] * alpha * (avg[upd] - Z[upd])
    return Z, w


# ------------------------------------------------------------ 网格清理
def drop_small_components(faces, tags=None):
    """只保留面数最多的连通组件(正常情况下包装 mesh 只有唯一一层壳,
    栅格噪声/异常轮廓可能产生碎块或多壳, 一律丢弃)。
    tags: 面标签(0=顶面 1=底面 2=侧壁), 随面一同过滤; 返回 faces
    或 (faces, tags)。"""
    from collections import defaultdict
    parent = list(range(len(faces)))

    def find(x):
        while parent[x] != x:
            parent[x] = parent[parent[x]]
            x = parent[x]
        return x

    e2f = defaultdict(list)
    for fi, f in enumerate(faces):
        for x, y in ((f[0], f[1]), (f[1], f[2]), (f[2], f[0])):
            e2f[(min(x, y), max(x, y))].append(fi)
    for fl in e2f.values():
        for fi in fl[1:]:
            parent[find(fi)] = find(fl[0])
    groups = defaultdict(int)
    for fi in range(len(faces)):
        groups[find(fi)] += 1
    if len(groups) <= 1:
        return (faces, tags) if tags is not None else faces
    best = max(groups.values())
    if best < len(faces) * 0.9:
        print('    警告: 网格有 %d 个组件, 仅保留最大者(%d/%d 面)'
              % (len(groups), best, len(faces)))
    keep = [fi for fi in range(len(faces)) if groups[find(fi)] == best]
    if len(keep) != len(faces):
        keep = np.asarray(keep)
        faces = faces[keep]
        if tags is not None:
            tags = tags[keep]
    return (faces, tags) if tags is not None else faces


def split_pinch_vertices(verts, faces):
    """拆分“夹点”顶点(一个顶点被多片不相连的面共享, 网格对角接触造成),
    拆分后网格为严格的顶点流形。返回 (verts, faces)。"""
    from collections import defaultdict
    vfaces = defaultdict(list)
    for fi, f in enumerate(faces):
        for v in f:
            vfaces[int(v)].append(fi)
    # v -> [ (本组面id列表), ... ]
    splits = []
    for v, fl in vfaces.items():
        if len(fl) < 2:
            continue
        edge_faces = defaultdict(list)
        for fi in fl:
            f = faces[fi]
            for x, y in ((f[0], f[1]), (f[1], f[2]), (f[2], f[0])):
                u = y if x == v else (x if y == v else None)
                if u is not None:
                    edge_faces[int(u)].append(fi)
        parent = {fi: fi for fi in fl}

        def find(x):
            while parent[x] != x:
                parent[x] = parent[parent[x]]
                x = parent[x]
            return x

        for fl2 in edge_faces.values():
            for fi in fl2[1:]:
                parent[find(fi)] = find(fl2[0])
        groups = defaultdict(list)
        for fi in fl:
            groups[find(fi)].append(fi)
        if len(groups) > 1:
            splits.append((v, [g for g in groups.values()]))
    if not splits:
        return verts, faces
    faces = faces.copy()
    extra_verts = []
    for v, groups in splits:
        for g in groups[1:]:
            nv = len(verts) + len(extra_verts)
            extra_verts.append(verts[v])
            for fi in g:
                f = faces[fi]
                faces[fi] = [nv if x == v else x for x in f]
    if extra_verts:
        verts = np.vstack([verts, np.asarray(extra_verts)])
    return verts, faces


# ------------------------------------------------------------ 网格松弛
def smooth_closed_curve(P, passes=3, lam=0.5, mu=-0.53):
    """闭合折线 Taubin(λ/μ) 低通平滑: 消除栅格楼梯而不收缩轮廓。
    交替"向邻点均值收缩(λ)"与"反向外推(μ)", 大尺度形状(鼻梁凹口、
    镜框圆角)保留, 只磨掉 cell 量级的锯齿。用于外轮廓多边形。"""
    P = np.array(P, np.float64, copy=True)
    for _ in range(passes):
        for a in (lam, mu):
            avg = 0.5 * (np.roll(P, 1, axis=0) + np.roll(P, -1, axis=0))
            P = P + a * (avg - P)
    return P


def smooth_rim_curve(verts, faces, tags, passes=20, lam=0.5):
    """薄膜周边收边: 沿顶膜边界曲线做 1D 拉普拉斯平滑。

    footprint 是网格点掩码, 薄膜边缘天生带 spacing 级楼梯。本函数:
    - 取顶面(tags==0)的边界边 → 每个边缘顶点恰有两个曲线邻居,
      沿这条闭合曲线反复做拉普拉斯松弛, 把楼梯磨成圆角;
      曲线点数近千、点距 ~spacing, 整体收缩量可忽略;
    - 底膜边缘顶点(顶/底膜是同一网格点的上下两个副本, 组装时
      (u,v) 相同)通过侧壁找到配对, 加上与顶膜边缘相同的位移,
      侧壁整体随缘平移, 不会被剪切;
    - 薄膜内部完全不被动(贴合不受影响)。

    注意: 不能改用对整个闭合 mesh 的拉普拉斯松弛来做收边/平滑——
    侧壁是 spacing x thickness 的细长三角形(长宽比可达 70:1),
    伞形算子会把膜面边缘顶点拽向侧壁另一端的底膜顶点, 反而在
    周边产生大幅 z 向锯齿(膜面顶点也会被拉伸错位)。
    tags: build_wrap_mesh 的面标签(0=顶面 1=底面 2=侧壁)。"""
    from collections import Counter, defaultdict
    topF = faces[tags == 0]
    botF = faces[tags == 1]
    wallF = faces[tags == 2]
    in_bot = np.zeros(len(verts), bool)
    in_bot[botF.ravel()] = True
    # 顶膜边界无向边(在顶面内出现一次的边) → 曲线邻接表
    cnt = Counter()
    for f in topF:
        for x, y in ((f[0], f[1]), (f[1], f[2]), (f[2], f[0])):
            cnt[(min(x, y), max(x, y))] += 1
    nbr = {}
    for (a, b) in (e for e, c in cnt.items() if c == 1):
        nbr.setdefault(int(a), []).append(int(b))
        nbr.setdefault(int(b), []).append(int(a))
    # 顶/底膜边缘配对: 与顶膜边缘顶点共享侧壁且 (u,v) 相同的底膜点
    v2wall = defaultdict(list)
    for fi, f in enumerate(wallF):
        for x in f:
            v2wall[int(x)].append(fi)
    keys = sorted(nbr.keys())
    pair = {}
    for v in keys:
        pv = verts[v]
        for fi in v2wall.get(v, ()):
            for x in wallF[fi]:
                x = int(x)
                if in_bot[x] and abs(verts[x, 0] - pv[0]) < 1e-9 \
                        and abs(verts[x, 1] - pv[1]) < 1e-9:
                    pair[v] = x
    # 沿曲线拉普拉斯松弛(Jacobi)
    pos = {v: verts[v].copy() for v in keys}
    for _ in range(passes):
        new = {}
        for v, P in pos.items():
            a, b = nbr[v]
            new[v] = P + lam * (0.5 * (pos[a] + pos[b]) - P)
        pos = new
    # 应用: 顶膜边缘取新位置, 底膜配对点加相同位移(侧壁随缘平移)
    for v in keys:
        delta = pos[v] - verts[v]
        verts[v] = pos[v]
        b = pair.get(v)
        if b is not None:
            verts[b] += delta
    return verts


# ------------------------------------------------------------ 包装 mesh 组装
def build_wrap_mesh(U, V, Z, inside, thickness, w_floor=None):
    """顶面膜 + 底面膜 + 侧壁 → 闭合网格。
    顶点/面均为 (u,v,w) 坐标, 返回 (verts, faces, tags),
    tags: 每个面的标签(0=顶面 1=底面 2=侧壁), 供收边平滑
    定位薄膜边缘(周边锯齿所在的侧壁)。
    w_floor: 底面膜最深不超过该值(默认=镜架最后缘再退 0.1),
    防止薄膜沿侧影线下潜时背面膜探出镜架范围。"""
    ni, nj = Z.shape
    idx = -np.ones((ni, nj), np.int64)
    ii, jj = np.nonzero(inside)
    n = len(ii)
    idx[ii, jj] = np.arange(n)
    zt = Z[ii, jj]
    zb = zt - thickness
    if w_floor is not None:
        zb = np.maximum(zb, w_floor)
    verts = np.concatenate([
        np.stack([U[ii, jj], V[ii, jj], zt], axis=1),
        np.stack([U[ii, jj], V[ii, jj], zb], axis=1),
    ])
    # 单元格四角: a=(i,j) b=(i+1,j) c=(i+1,j+1) d=(i,j+1)
    a = idx[:-1, :-1]; b = idx[1:, :-1]; c = idx[1:, 1:]; d = idx[:-1, 1:]
    ok = (a >= 0) & (b >= 0) & (c >= 0) & (d >= 0)
    ci, cj = np.nonzero(ok)
    A, B, C, D = a[ci, cj], b[ci, cj], c[ci, cj], d[ci, cj]
    # 棋盘式交替对角线: 相邻单元的剖分方向相反, 避免整片同向的机械条纹
    m0 = ((ci + cj) % 2 == 0)
    topF = np.concatenate([
        np.stack([A[m0], B[m0], C[m0]], 1), np.stack([A[m0], C[m0], D[m0]], 1),
        np.stack([A[~m0], B[~m0], D[~m0]], 1), np.stack([B[~m0], C[~m0], D[~m0]], 1)])
    botF = np.concatenate([
        np.stack([A[m0], C[m0], B[m0]], 1), np.stack([A[m0], D[m0], C[m0]], 1),
        np.stack([A[~m0], D[~m0], B[~m0]], 1), np.stack([B[~m0], D[~m0], C[~m0]], 1)]) + n

    # 顶面的边界有向边(内部在左侧) → 缝合到背面
    # 先做无向计数找边界(内部边在相邻两面中方向相反、各出现一次),
    # 再取其在面中出现的方向用于侧壁绕向
    from collections import Counter
    und = Counter()
    directed = []
    for f in topF:
        for x, y in ((f[0], f[1]), (f[1], f[2]), (f[2], f[0])):
            key = (x, y) if x < y else (y, x)
            und[key] += 1
            directed.append((x, y))
    bndset = {e for e, c in und.items() if c == 1}
    rim_pairs = [(p, q) for (p, q) in directed
                 if ((p, q) if p < q else (q, p)) in bndset]
    rim = [[q, p, p + n] for (p, q) in rim_pairs]
    rim += [[q, p + n, q + n] for (p, q) in rim_pairs]
    faces = np.concatenate([topF, botF, np.asarray(rim, np.int64)])
    tags = np.concatenate([np.zeros(len(topF), np.int8),
                           np.ones(len(botF), np.int8),
                           np.full(len(rim), 2, np.int8)])

    # 去掉未被任何面引用的孤立顶点
    faces, tags = drop_small_components(faces, tags)
    verts, faces = split_pinch_vertices(verts, faces)
    used = np.zeros(len(verts), bool)
    used[faces.ravel()] = True
    remap = -np.ones(len(verts), np.int64)
    remap[used] = np.arange(used.sum())
    return verts[used], remap[faces], tags


def check_closed(faces, nverts):
    """水密性与拓扑自检: 每条有向边恰好出现一次且其反向也存在
    (方向一致的水密闭合), 欧拉示性数 V - E + F = 2。"""
    d = {}
    for f in faces:
        for x, y in ((f[0], f[1]), (f[1], f[2]), (f[2], f[0])):
            d[(x, y)] = d.get((x, y), 0) + 1
    bad = [e for e, c in d.items() if c != 1 or (e[1], e[0]) not in d]
    n_e = 3 * len(faces) // 2
    chi = nverts - n_e + len(faces)
    return len(bad) == 0, chi


# ------------------------------------------------------------------- 主流程
def default_obj_path():
    """自动探测输入模型: 取脚本所在目录下最大的 *.obj,
    排除历史输出(*_result/wrap_result/*_fill* 等)。"""
    cands = []
    for p in glob.glob(os.path.join(SCRIPT_DIR, '*.obj')):
        b = os.path.basename(p).lower()
        if 'result' in b or 'wrap' in b or '_fill' in b:
            continue
        cands.append((os.path.getsize(p), p))
    if not cands:
        return None
    return max(cands)[1]


def main():
    ap = argparse.ArgumentParser(description='眼镜“抽真空包装”mesh 生成')
    ap.add_argument('--obj', default=None,
                    help='原始模型 obj(默认: 自动探测本目录最大的非输出 obj)')
    ap.add_argument('--index', default=None,
                    help='眼镜顶点索引(默认: 模型同目录的 index.txt)')
    ap.add_argument('--out', default=None,
                    help='输出路径(默认: 模型同目录 <模型名>_wrap_result.obj)')
    ap.add_argument('--cell', type=float, default=0.05,
                    help='外轮廓栅格分辨率(模型单位, 默认0.05)')
    ap.add_argument('--spacing', type=float, default=0.06,
                    help='包装薄膜网格点间距(默认0.06)')
    ap.add_argument('--thickness', type=float, default=3.4,
                    help='背面膜相对薄膜的后退厚度(默认3.4)')
    ap.add_argument('--front', type=float, default=0.25,
                    help='朝前面阈值(法线与视线轴余弦, 默认0.25); '
                         '调小会包含更斜的面乃至镜腿, 调大只贴更正对的面')
    ap.add_argument('--closing', type=int, default=4,
                    help='形态学闭运算半径(格数, 默认4)')
    ap.add_argument('--hole-blend', type=int, default=500,
                    help='镜片镂空处残差调和的迭代次数(默认500): '
                         '趋势面跟随镜框平齐, 欠松弛迭代负责孔缘平滑衔接; '
                         '越大越平滑')
    ap.add_argument('--rim-smooth', type=int, default=20,
                    help='周边收边平滑次数(0=关闭, 默认20): 沿薄膜边缘曲线 '
                         '拉普拉斯平滑, 磨掉轮廓栅格锯齿')
    ap.add_argument('--poly-smooth', type=int, default=3,
                    help='外轮廓多边形 Taubin 平滑次数(0=关闭, 默认3), '
                         '从源头去掉 cell 级楼梯')
    ap.add_argument('--band-smooth', type=float, default=0.035,
                    help='贴合带平滑阈值(×spacing, 0=关闭; 磨掉陡坡棱面)')
    ap.add_argument('--mtl', default=None,
                    help='覆盖 mtllib 行的材料文件名(默认自动找本目录 *.mtl)')
    args = ap.parse_args()

    if args.obj is None:
        args.obj = default_obj_path()
        if args.obj is None:
            raise RuntimeError('未指定 --obj, 且本目录没有可自动探测的模型文件')
    obj_path = resolve(args.obj)
    obj_dir = os.path.dirname(obj_path)
    if args.index is None:
        args.index = os.path.join(obj_dir, 'index.txt')
    if args.out is None:
        stem = os.path.splitext(os.path.basename(obj_path))[0]
        args.out = os.path.join(obj_dir, stem + '_wrap_result.obj')
    idx_path, out_path = resolve(args.index), resolve(args.out)

    print('[1/6] 读取 %s ...' % os.path.basename(obj_path))
    verts, faces = load_obj(obj_path)
    print('    %d 顶点, %d 三角面' % (len(verts), len(faces)))
    gidx = load_index(idx_path, len(verts))
    print('    index.txt: %d 个眼镜顶点' % len(gidx))

    # ---- 眼镜子网格 & 正面轴向
    in_g = np.zeros(len(verts), bool)
    in_g[gidx] = True
    sel = in_g[faces].all(axis=1)
    sub = faces[sel]
    if len(sub) < 10:
        raise RuntimeError('眼镜子网格面太少(%d), 请检查 index.txt' % len(sub))
    print('    眼镜子网格: %d 面' % len(sub))
    # 整 mesh 体积符号 → 面法线是否朝外
    p0, p1, p2 = verts[faces[:, 0]], verts[faces[:, 1]], verts[faces[:, 2]]
    vol = np.einsum('ij,ij->i', p0, np.cross(p1, p2)).sum() / 6.0
    wind = 1.0 if vol >= 0 else -1.0
    coarse = max(0.15, args.cell * 1.5)
    print('    判定包装视线轴(封闭镂空最大者):')
    k, sgn, ax_u, ax_v = pick_front_axis(verts, sub, wind, coarse)
    print('    正面视线轴 = %s, u=轴%d(%s), v=轴%d(%s)'
          % ('xyz'[k] + ('+' if sgn > 0 else '-'),
             ax_u, 'xyz'[ax_u], ax_v, 'xyz'[ax_v]))

    tri_uv = verts[sub][:, :, [ax_u, ax_v]]
    tri_w = sgn * verts[sub][:, :, k]

    # ---- 朝前面掩码(法线与视线轴余弦) → 只包装镜框正面, 排除向后弯折的镜腿
    sp0, sp1, sp2 = verts[sub[:, 0]], verts[sub[:, 1]], verts[sub[:, 2]]
    nrm = np.cross(sp1 - sp0, sp2 - sp0)
    ln = np.linalg.norm(nrm, axis=1)
    front = (sgn * nrm[:, k] / np.where(ln > 1e-12, ln, 1.0)) > args.front
    print('    朝前面: %d/%d (front>%.2f)' % (front.sum(), len(sub), args.front))

    # ---- 栅格化 → 外轮廓(闭运算补缝 → 只取含镜片镂空的镜框核心组件,
    #      内部封闭孔洞填满; 半径自适应有限加大, 不会封死鼻梁凹口)
    print('[2/6] 栅格化朝前表面 (cell=%.2f) ...' % args.cell)
    cover, topz, (u0, v0) = rasterize(tri_uv, tri_w, args.cell, front)
    cover_all, topz_all, (ua, va), cell_tris = rasterize(
        tri_uv, tri_w, args.cell, collect=True)
    assert (ua, va) == (u0, v0)               # 两次栅格化共用同一网格
    print('    朝前表面覆盖 %d/%d 格' % (cover.sum(), cover.size))
    r = args.closing
    cap = max(args.closing + 4, int(round(1.0 / args.cell)))
    cover_area = float(cover.sum()) * args.cell ** 2
    while True:
        # 镜片镂空改以"全表面轮缘环"的封闭孔判定:
        # 轮缘正面覆盖稀疏/有缺口时, 正面环闭不上会漏填整块镜片
        # (632587 右镜片曾缺失 1/3); 全表面环更完整, 而 开向边界的
        # 鼻梁凹口不属于封闭孔, 不会被误填
        ring = binary_closing(cover_all, r)
        holes = fill_enclosed(ring) & ~ring
        foot = largest_filled_component(binary_closing(cover, r) | holes)
        # 上/下缘的掠射面不在正面覆盖内但属于镜架:
        # 正面 footprint 附近 0.2 内的轮缘格也纳入, 包装才铺满镜框外缘
        near = foot.copy()
        for _ in range(4):
            near = dilate(near)
        foot = largest_filled_component(foot | (ring & near))
        cycle, n_sig = trace_outer_cycle(foot)
        cyc_area = abs(signed_area(cycle)) * args.cell ** 2
        foot_area = float(foot.sum()) * args.cell ** 2
        # 判据一: 外轮廓单环完整(追踪不破);
        # 判据二: 镜片镂空被填上(footprint 达到朝前覆盖的 2 倍以上,
        #         覆盖稀疏/有缺口的数据需要更大的闭运算半径)
        if (cyc_area >= 0.9 * foot_area and foot_area >= 2.0 * cover_area) \
                or r >= cap:
            break
        print('    闭运算半径 %d → 外轮廓覆盖 footprint 的 %.0f%%, '
              '镂空填充率 %.2f, 继续加大...'
              % (r, 100 * cyc_area / foot_area, foot_area / cover_area))
        r += 2
    if foot_area < 1.5 * cover_area:
        print('    注意: 镜片镂空可能未被填满(填充率 %.2f), '
              '可尝试加大 --closing' % (foot_area / cover_area))
    print('    镜框正面 footprint: %d 格(面积 %.1f), 闭运算半径 %d'
          % (foot.sum(), foot_area, r))
    poly = np.array([[u0 + i * args.cell, v0 + j * args.cell] for (i, j) in cycle],
                    dtype=np.float64)
    print('    外轮廓多边形: %d 个角点' % len(poly))
    if args.poly_smooth > 0:
        poly = smooth_closed_curve(poly, passes=args.poly_smooth)
        print('    外轮廓 Taubin 平滑 %d 次(消除 cell 级栅格楼梯, 不收缩轮廓)'
              % args.poly_smooth)

    # ---- 薄膜网格点
    print('[3/6] 铺薄膜网格 (spacing=%.2f) ...' % args.spacing)
    s = args.spacing
    gu0 = poly[:, 0].min() - s; gu1 = poly[:, 0].max() + s
    gv0 = poly[:, 1].min() - s; gv1 = poly[:, 1].max() + s
    ni = int(np.ceil((gu1 - gu0) / s)) + 1
    nj = int(np.ceil((gv1 - gv0) / s)) + 1
    GU, GV = np.meshgrid(gu0 + np.arange(ni) * s,
                         gv0 + np.arange(nj) * s, indexing='ij')
    pts = np.stack([GU.ravel(), GV.ravel()], axis=1)
    ins = points_in_poly(pts, poly).reshape(ni, nj)
    print('    轮廓内网格点: %d' % ins.sum())
    if ins.sum() < 4:
        raise RuntimeError('轮廓内网格点太少, 请增大 --spacing 或减小 --cell')

    # ---- 按镜架整体投影剪裁: 闭运算/轮缘带会把 footprint 向外胀出
    #      0.1~0.2 的薄裙边(镜架上缘/下缘可见"超出"), 以镜架外轮廓
    #      (全表面投影 + 补采样缝 + 填充内部镂空)为准裁掉
    frame_sil = fill_enclosed(binary_closing(cover_all, 2))
    ic = np.clip(np.round((GU - u0) / args.cell).astype(int),
                 0, frame_sil.shape[0] - 1)
    jc = np.clip(np.round((GV - v0) / args.cell).astype(int),
                 0, frame_sil.shape[1] - 1)
    n_clip = int((ins & ~frame_sil[ic, jc]).sum())
    if n_clip:
        ins &= frame_sil[ic, jc]
        print('    按镜架投影裁掉外轮廓之外的裙边: %d 个网格点' % n_clip)
    # 趋势拟合的锚点只取正面表面(cover 膨胀 3 格)上的贴合高度:
    # 外缘陡边/掠射面处薄膜本该顺侧影下潜, 高度不能拿来外推镜片平面
    cov_d = cover.copy()
    for _ in range(3):
        cov_d = dilate(cov_d)
    front_sel = cov_d[ic, jc]

    # ---- 顶面高度: 贴表面 + 镂空处与镜框平齐
    print('[4/6] 计算薄膜高度(贴眼镜表面, 镜片处与镜框平齐) ...')
    Z = exact_film_heights(pts, tri_uv, tri_w, cell_tris,
                           u0, v0, args.cell).reshape(ni, nj)
    known = ~np.isnan(Z) & ins
    n_hole = int((ins & ~known).sum())
    print('    %d 点直接贴到表面, %d 点位于镂空处(趋势平齐填充)'
          % (known.sum(), n_hole))
    Z = trend_fill(Z, ins, select=front_sel, blend_passes=args.hole_blend)

    # ---- 贴合带平滑: 镜框侧弯坡度大处, 高度场曲率被坡度放大成棱面
    #      (相邻面法线跳变 30-60 度), 受限松弛磨掉棱面并软化接缝;
    #      浅坡处与均匀斜面本来就光滑, 不会被改动
    if args.band_smooth > 0:
        thr = args.band_smooth * args.spacing
        Zs, w = smooth_adherence(Z, known, ins, thr, args.spacing)
        n_up = int((np.abs(Zs - Z) > 1e-6).sum())
        if n_up:
            # 松弛后的点(含软化接缝带)作为新边界, 镂空处重新平齐外推
            Z = trend_fill(np.where(w > 0, Zs, np.nan), ins,
                           select=front_sel, blend_passes=args.hole_blend)
            print('    贴合带平滑: %d 个点被松弛(阈值 %.4f), '
                  '镂空处重新平齐' % (n_up, thr))

    # ---- 组装闭合 mesh
    print('[5/6] 组装水密包装 mesh (thickness=%.2f) ...' % args.thickness)
    # 底面膜深度钳制: 不超过镜架最后缘(留 0.1 薄膜厚度),
    # 防止薄膜沿侧影线下潜时背面膜探出镜架范围
    w_floor = float((sgn * verts[sub][:, :, k]).min()) - 0.1
    wverts, wfaces, ftags = build_wrap_mesh(GU, GV, Z, ins, args.thickness, w_floor)
    if args.rim_smooth > 0:
        # 周边收边: 顶膜边缘沿边界曲线 1D 拉普拉斯, 磨掉栅格楼梯,
        # 底膜边缘同步位移(侧壁随缘平移), 薄膜内部不动。
        # 注意: 不做整 mesh 的拉普拉斯松弛——侧壁细长三角形会把
        # 伞形算子变成"跨 3~4 个单位距离的强混合", 边缘反而出锯齿。
        print('    周边收边: 边缘曲线拉普拉斯平滑 %d 次 ...' % args.rim_smooth)
        wverts = smooth_rim_curve(wverts, wfaces, ftags, passes=args.rim_smooth)
    closed, chi = check_closed(wfaces, len(wverts))
    print('    %d 顶点, %d 面; 水密=%s, 欧拉示性数=%d'
          % (len(wverts), len(wfaces), closed, chi))
    if not closed:
        raise RuntimeError('生成的网格不是水密闭合网格, 请检查参数')

    # 变换回世界坐标: (u, v, w) -> 原坐标
    W = sgn * wverts[:, 2]
    verts_world = np.zeros_like(wverts)
    verts_world[:, ax_u] = wverts[:, 0]
    verts_world[:, ax_v] = wverts[:, 1]
    verts_world[:, k] = W

    # 顶点法线(面积加权)
    p0, p1, p2 = verts_world[wfaces[:, 0]], verts_world[wfaces[:, 1]], verts_world[wfaces[:, 2]]
    fn = np.cross(p1 - p0, p2 - p0)
    vn = np.zeros_like(verts_world)
    for c in range(3):
        np.add.at(vn, wfaces[:, c], fn)
    ln = np.linalg.norm(vn, axis=1, keepdims=True)
    vn = vn / np.where(ln > 1e-12, ln, 1.0)

    # ---- 写出 obj
    print('[6/6] 写出 %s ...' % os.path.basename(out_path))
    mtl = args.mtl
    if mtl is None:                      # 先找模型同目录的 mtl, 再退回脚本目录
        found = sorted(glob.glob(os.path.join(obj_dir, '*.mtl'))) or \
            sorted(glob.glob(os.path.join(SCRIPT_DIR, '*.mtl')))
        mtl = os.path.basename(found[0]) if found else None
    n_v = n_vn = 0
    with open(out_path, 'w') as out, open(obj_path, 'r', errors='ignore') as src:
        for line in src:                       # 原模型逐行原样拷贝
            if line.startswith('v '):
                n_v += 1
            elif line.startswith('vn '):
                n_vn += 1
            elif line.startswith('mtllib ') and mtl:
                line = 'mtllib %s\n' % mtl      # 指向本目录实际存在的 mtl
            elif line.startswith('usemtl ') and line.strip().endswith('.001'):
                line = 'usemtl %s\n' % line.split()[1][:-4]   # Blender 重命名残留
            out.write(line)
        # 追加包装 mesh
        out.write('o Transparent\n')
        out.write('s 0\n')
        out.write('usemtl Transparent\n')
        for p in verts_world:
            out.write('v %.6f %.6f %.6f\n' % (p[0], p[1], p[2]))
        base_v, base_n = n_v, n_vn
        for nn in vn:
            out.write('vn %.6f %.6f %.6f\n' % (nn[0], nn[1], nn[2]))
        for f in wfaces:
            out.write('f %d//%d %d//%d %d//%d\n' % (
                f[0] + 1 + base_v, f[0] + 1 + base_n,
                f[1] + 1 + base_v, f[1] + 1 + base_n,
                f[2] + 1 + base_v, f[2] + 1 + base_n))

    lo, hi = W.min(), W.max()
    if sgn < 0:
        lo, hi = -hi, -lo
    print('完成! 输出: %s' % out_path)
    print('  新增 Transparent mesh: %d 顶点 / %d 面 (参考结果为 9999/19994)'
          % (len(wverts), len(wfaces)))
    print('  沿视线轴 %s 的包装范围: [%.2f, %.2f]' % ('xyz'[k], lo, hi))
    print('  包装 xy 外接盒: x[%.2f, %.2f], y[%.2f, %.2f]'
          % (verts_world[:, ax_u].min(), verts_world[:, ax_u].max(),
             verts_world[:, ax_v].min(), verts_world[:, ax_v].max()))


if __name__ == '__main__':
    main()