新算法

.gitignore

@@ -78,3 +78,36 @@ color_change/__pycache__/calib_io.cpython-310.pyc
 .gitignore
 .gitignore
 color_change/__pycache__/calibration.cpython-310.pyc
+fill2sc/632723.jpg
+fill2sc/632723.mtl
+fill2sc/fill_preview.png
+fill2sc/fill_report.txt
+.gitignore
+fill2sc/index.txt
+fill2sc/data/20261102/20261102_glasses_vertices.txt
+fill2sc/data/20261102/20261102_glasses.blend
+fill2sc/data/20261102/20261102.mtl
+fill2sc/data/20261102/20261102.zip
+fill2sc/data/20261102/20261102Tex1.jpg
+fill2sc/data/20261102/fill_20261102_preview.png
+fill2sc/data/20261102/index.txt
+fill2sc/data/631230/631230_glasses_vertices.txt
+fill2sc/data/631230/631230_glasses.blend
+fill2sc/data/631230/631230.jpg
+fill2sc/data/631230/631230.mtl
+fill2sc/data/631230/631230.zip
+fill2sc/data/631230/fill_631230_preview.png
+fill2sc/data/631230/index.txt
+fill2sc/data/632723/632723.jpg
+fill2sc/data/632723/632723.mtl
+fill2sc/data/632723/fill_632723_preview.png
+fill2sc/data/632723/fill_preview.png
+fill2sc/data/632723/fill_report.txt
+fill2sc/data/632723/index.txt
+fill2sc/data/636394/636394_glasses_vertices.txt
+fill2sc/data/636394/636394_glasses.blend
+fill2sc/data/636394/636394.jpg
+fill2sc/data/636394/636394.mtl
+fill2sc/data/636394/636394.zip
+fill2sc/data/636394/fill_636394_preview.png
+fill2sc/data/636394/index.txt

fill2sc/make_shell.py

@@ -0,0 +1,482 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+"""
+make_shell.py — 为每个眼镜实例的顶点连通区域生成"最小壳" (minimal tight shell)。
+
+输入支持两种方式:
+  1. data.json (推荐, 支持多眼镜实例, 每个实例单独处理并各导出一个壳):
+       {"input": "xxx.obj",
+        "instances": [
+          {"id": 0, "type": "ring-pair", "score": 0.91, "count": 1234,
+           "indices": [12, 13, 14, ...]},
+          {"id": 1, "type": "ring-pair", "score": 0.87, "count": 1180,
+           "indices": [5012, 5013, ...]}
+        ]}
+     "input" 相对 data.json 所在目录解析 (解析不到时按同名 obj 兜底); 也可用 --obj 显式指定 (优先级更高)。
+     各实例独立计算, 壳最后合并导出为单个 obj (默认 shell.obj)。
+  2. index.txt (兼容旧的单实例流程): 顶点索引, 每行一个, 输出单个 shell.obj。
+
+不是 AABB, 也不是带 padding 的膨胀包裹, 而是 alpha shape (3D 凹包):
+  - 壳面由区域最外侧的原始采样点构成, 紧贴外侧轮廓 (含凹陷, 如眼镜框);
+  - alpha 自动取"使壳保持单一连通的最小值" (二分搜索), 即最小壳;
+  - 桥接只发生在顶点漏采的空洞处, 桥的粗细由最小 alpha 自然决定;
+  - 最后做可选的 Taubin 平滑 (不收缩), 兼顾"最小"与"光滑"。
+
+流程 (每个实例独立执行一遍):
+  1. 读 obj (全局一次) 与实例顶点索引 (默认 0-based), 取区域顶点 +
+     "三顶点均在区域内"的三角形重心/边中点加密 -> 采样点;
+  2. 3D Delaunay 四面体化, 计算每个四面体的外接球半径;
+  3. 二分搜索最小 alpha: alpha 复合体 (外接球半径 <= alpha 的四面体之并)
+     满足 (a) 连通 (b) 覆盖全部采样点;
+  4. 若最小 alpha 下边界非流形/不水密, alpha 向上微调至水密;
+  5. 抽复合体的边界三角形 (外向绕序), 得到紧贴的壳;
+  6. 可选 Taubin 平滑 + 数值校验 (点偏移, 体积, 水密性);
+  7. 全部实例的壳合并 (同一坐标系直接拼接) 导出单个 obj。
+
+用法:
+  python make_shell.py                          # 读同目录 data.json, 各实例合并导出 shell.obj
+  python make_shell.py --data other.json        # 指定 data.json
+  python make_shell.py --index index.txt        # 兼容旧单索引, 导出 shell.obj
+  python make_shell.py --alpha 3.0              # 手动指定 alpha (更大=更圆润)
+  python make_shell.py --smooth 0               # 关闭平滑, 严格贴点
+"""
+import argparse
+import glob
+import json
+import os
+
+import numpy as np
+from scipy.spatial import Delaunay, cKDTree
+from scipy.sparse import coo_matrix
+from scipy.sparse.csgraph import connected_components
+from scipy import ndimage
+from skimage import measure
+import trimesh
+from trimesh.smoothing import filter_taubin
+
+
+# ---------------------------------------------------------------- IO
+def load_obj(path):
+    verts, faces = [], []
+    with open(path) as f:
+        for line in f:
+            if line.startswith("v "):
+                verts.append([float(x) for x in line.split()[1:4]])
+            elif line.startswith("f "):
+                faces.append([int(t.split("/")[0]) - 1 for t in line.split()[1:4]])
+    return np.asarray(verts, np.float64), np.asarray(faces, np.int64)
+
+
+def load_data(args):
+    """解析输入, 返回 (instances, input_obj)。
+
+    instances: [(实例名, 顶点索引 ndarray), ...]
+    input_obj: data.json 声明的 obj 路径 (未解析相对路径, 可能为 None)
+
+    优先级: --index 单实例 txt > --data 多实例 json。
+    """
+    if args.index:
+        idx = np.loadtxt(args.index, dtype=np.int64).ravel()
+        if args.one_based:
+            idx = idx - 1
+        return [("index", idx)], None
+
+    if not os.path.exists(args.data):
+        raise SystemExit(f"找不到 {args.data}, 请用 --data 指定 data.json, "
+                         f"或用 --index 指定单实例索引文件")
+    with open(args.data) as f:
+        data = json.load(f)
+
+    insts = []
+    for inst in data.get("instances", []):
+        name = str(inst.get("id", len(insts)))
+        idx = np.asarray(inst.get("indices", []), dtype=np.int64).ravel()
+        if args.one_based:
+            idx = idx - 1
+        cnt = inst.get("count")
+        if cnt is not None and cnt != len(idx):
+            print(f"[warn] 实例 {name}: count={cnt} 与 indices 数量 {len(idx)} 不一致")
+        if len(idx) == 0:
+            print(f"[warn] 实例 {name}: 索引为空, 跳过")
+            continue
+        insts.append((name, idx))
+    if not insts:
+        raise SystemExit(f"{args.data} 中没有可用实例 (instances 为空)")
+    return insts, data.get("input")
+
+
+def resolve_obj(args, input_obj):
+    """确定 obj 路径: --obj 显式指定 > data.json 的 input 字段 > data 目录内唯一 obj。
+
+    input 字段可能是相对其他根目录的路径 (如 "datasets/xxx/xxx.obj"), 本地解析不到时
+    按 basename 在 data 目录内找同名 obj 兜底。
+    """
+    if args.obj:
+        return args.obj
+    base = os.path.dirname(os.path.abspath(args.data))
+    if input_obj:
+        p = input_obj if os.path.isabs(input_obj) else os.path.join(base, input_obj)
+        if os.path.exists(p):
+            return p
+        name = os.path.basename(input_obj)     # 按文件名在 data 目录内找同名 obj
+        cands = sorted(glob.glob(os.path.join(base, name)))
+        if len(cands) == 1:
+            print(f"[warn] data.json 的 input '{input_obj}' 本地不存在, "
+                  f"使用同名 obj: {cands[0]}")
+            return cands[0]
+        print(f"[warn] data.json 的 input '{input_obj}' 本地不存在, 尝试目录内查找 obj")
+    cands = sorted(glob.glob(os.path.join(base, "*.obj")))
+    if len(cands) == 1:
+        return cands[0]
+    if not cands:
+        raise SystemExit(f"未找到 obj: 请用 --obj 指定, 或在 {base} 放置 .obj 文件")
+    raise SystemExit(f"{base} 下有多个 obj, 请用 --obj 指定其一: "
+                     f"{', '.join(os.path.basename(c) for c in cands)}")
+
+
+def densify_points(verts, faces, region_mask):
+    """区域顶点 + 区域内部三角形的重心/边中点采样点。"""
+    pts = [verts[region_mask]]
+    inner = faces[np.all(region_mask[faces], axis=1)]
+    if len(inner):
+        tris = verts[inner]                                   # (F,3,3)
+        for a, b, c in ((1/3, 1/3, 1/3),
+                        (0.5, 0.5, 0.0), (0.5, 0.0, 0.5), (0.0, 0.5, 0.5)):
+            pts.append(a * tris[:, 0] + b * tris[:, 1] + c * tris[:, 2])
+    pts = np.concatenate(pts, axis=0)
+    return np.unique(pts, axis=0), inner                      # 去重, 避免退化四面体
+
+
+def largest_cluster(pts, r):
+    """半径 r 欧氏聚类, 返回最大簇的布尔掩码 (剔除分割噪声)。"""
+    tree = cKDTree(pts)
+    pairs = np.array(list(tree.query_pairs(r=r)), dtype=np.int64).reshape(-1, 2)
+    n = len(pts)
+    if len(pairs):
+        m = coo_matrix((np.ones(len(pairs)), (pairs[:, 0], pairs[:, 1])), shape=(n, n))
+        _, lab = connected_components(m, directed=False)
+        sizes = np.bincount(lab)
+        return lab == int(np.argmax(sizes)), int(sizes.max())
+    return np.ones(n, bool), n
+
+
+# ------------------------------------------------------ Delaunay 预处理
+class AlphaStructure:
+    """Delaunay 四面体化 + 面邻接/外接球半径, 供 alpha 二分搜索复用。"""
+
+    def __init__(self, pts):
+        self.pts = pts
+        tri = Delaunay(pts)
+        self.tri = tri                                      # find_simplex 用于凸包内部栅格化
+        tets = tri.simplices.copy()
+        a = pts[tets[:, 0]]
+        m = np.stack([pts[tets[:, 1]] - a, pts[tets[:, 2]] - a, pts[tets[:, 3]] - a], axis=1)
+        det = np.linalg.det(m)
+        scale = float(np.abs(det).max()) ** (1 / 3) + 1e-30
+        good = np.abs(det) > 1e-10 * scale                    # 剔除退化四面体
+        tets, det, a = tets[good], det[good], a[good]
+        tets[det < 0] = tets[det < 0][:, [0, 1, 3, 2]]        # 统一正体积绕序
+        self.tets = tets
+        self.T = len(tets)
+
+        # 外接球半径: 解 2(b-a).x = |b|^2-|a|^2 线性方程组
+        e = [np.einsum("ij,ij->i", pts[tets[:, k]], pts[tets[:, k]]) for k in range(4)]
+        A = 2.0 * np.stack([pts[tets[:, 1]] - a, pts[tets[:, 2]] - a, pts[tets[:, 3]] - a], axis=1)
+        rhs = np.stack([e[1] - e[0], e[2] - e[0], e[3] - e[0]], axis=1)
+        center = np.linalg.solve(A, rhs)
+        self.circum = np.linalg.norm(center - a, axis=1)
+
+        # 每个四面体的 4 个面 (tet 优先排列: 线性位置 t*4+f), 外向绕序 + 面共享关系
+        faces = np.concatenate([tets[:, [1, 2, 3]], tets[:, [0, 3, 2]],
+                                tets[:, [0, 1, 3]], tets[:, [0, 2, 1]]],
+                               axis=1).reshape(self.T * 4, 3)
+        fcent = pts[faces].mean(axis=1)
+        tcent = np.repeat(pts[tets].mean(axis=1), 4, axis=0)
+        n = np.cross(pts[faces[:, 1]] - pts[faces[:, 0]], pts[faces[:, 2]] - pts[faces[:, 0]])
+        flip = np.einsum("ij,ij->i", n, fcent - tcent) < 0
+        faces[flip] = faces[flip][:, [0, 2, 1]]                # 法向朝外
+        self.faces = faces
+
+        keys = np.sort(self.faces, axis=1)
+        _, inv = np.unique(keys, axis=0, return_inverse=True)
+        self.face_id = np.asarray(inv).reshape(-1)
+        cnt = np.bincount(self.face_id)
+        shared = np.flatnonzero(cnt == 2)
+        # 每个共享面 -> 一对四面体 (面在线性数组中的两个出现位置)
+        order = np.argsort(self.face_id, kind="stable")
+        inv_sorted = self.face_id[order]
+        starts = np.searchsorted(inv_sorted, shared, side="left")
+        self.pairs = np.stack([order[starts], order[starts + 1]], axis=1) // 4
+        self.n_pts = len(pts)
+
+    # ---- 给定 alpha 的复合体状态 --------------------------------------
+    def status(self, alpha):
+        keep = self.circum <= alpha
+        used = np.zeros(self.n_pts, bool)
+        used[self.tets[keep].ravel()] = True
+        ok_used = bool(used.all())
+        pe = keep[self.pairs[:, 0]] & keep[self.pairs[:, 1]]
+        g = coo_matrix((np.ones(int(pe.sum()), bool),
+                        (self.pairs[pe, 0], self.pairs[pe, 1])), shape=(self.T, self.T))
+        nc, lab = connected_components(g, directed=False)
+        sizes = np.bincount(lab, weights=keep.astype(np.int64))
+        ok_conn = nc - int((~keep).sum()) == 1 and int(sizes.max()) == int(keep.sum())
+        return keep, ok_used, ok_conn
+
+    def boundary_faces(self, keep):
+        """复合体的边界三角形 (恰属于一个保留四面体的面, 外向绕序)。"""
+        mask = np.repeat(keep, 4)                             # faces 按 [tet 的 4 面] 排列
+        fid = self.face_id[mask]
+        cnt = np.bincount(fid)                                # 该面在保留四面体中的出现次数
+        return self.faces[mask][cnt[fid] == 1]
+
+
+# ------------------------------------------------------ 眼镜盒式外壳
+def build_case(pts, tri, res, sigma, faces_target):
+    """凸包(最外侧点连成的最小凸外壳) -> 内部栅格化 -> 高斯平滑 -> MC。
+
+    壳贴着最外侧点(误差 ~ 1 voxel), 边角以 sigma*voxel 尺度圆润, 均匀光滑, 水密。
+    """
+    hull = trimesh.Trimesh(pts).convex_hull
+    lo, hi = pts.min(axis=0), pts.max(axis=0)
+    voxel = (hi - lo).max() / res
+    pad = 6 * sigma * voxel + 0.5
+    lo2 = lo - pad
+    dims = np.maximum(np.ceil((hi + pad - lo2) / voxel).astype(int) + 1, 2)
+
+    # 体素中心是否在凸包内 (Delaunay find_simplex, C 实现)
+    grids = np.meshgrid(*[lo2[i] + (np.arange(dims[i]) + 0.5) * voxel for i in range(3)], indexing="ij")
+    centers = np.stack([g.ravel() for g in grids], axis=1)
+    occ = (tri.find_simplex(centers) >= 0).reshape(dims)   # 体素中心在凸包内
+
+    field = ndimage.gaussian_filter(occ.astype(np.float32), sigma=sigma)
+    v, f, _, _ = measure.marching_cubes(field, level=0.5)
+    v = (v + 0.5) * voxel + lo2                            # 网格索引 -> 世界坐标
+    m = trimesh.Trimesh(vertices=v, faces=f[:, ::-1], process=True)
+    m = max(m.split(only_watertight=False), key=lambda x: x.volume)
+    filter_taubin(m, lamb=0.5, nu=0.5, iterations=5)
+
+    if faces_target and len(m.faces) > faces_target:       # pymeshlab 简化 + 轻平滑
+        import pymeshlab as pml
+        ms = pml.MeshSet()
+        ms.add_mesh(pml.Mesh(vertex_matrix=m.vertices, face_matrix=m.faces.astype(np.int64)))
+        ms.meshing_decimation_quadric_edge_collapse(targetfacenum=int(faces_target),
+                                                    preservenormal=True)
+        out = ms.current_mesh()
+        m = trimesh.Trimesh(vertices=np.array(out.vertex_matrix()),
+                            faces=np.array(out.face_matrix()), process=False)
+        filter_taubin(m, lamb=0.5, nu=0.5, iterations=3)
+    return m, hull
+
+
+# ------------------------------------------------------ 单实例流水线
+def run_instance(verts, faces, idx, args):
+    """处理一个眼镜实例 (一组顶点索引) 的完整流程, 返回该实例的壳 mesh (由调用方合并导出)。"""
+    assert idx.min() >= 0 and idx.max() < len(verts), \
+        "index 超出顶点范围, 请确认索引含义 (0/1-based, --one-based 可切换)"
+    region_mask = np.zeros(len(verts), bool)
+    region_mask[idx] = True                                     # Blender 0-based 顶点索引
+
+    # ---- 噪声过滤: 只保留最大空间簇 (index 中混有 SAM3 多视角假阳性) ----
+    if args.filter == "main":
+        ridx = np.flatnonzero(region_mask)
+        mkeep, nkeep = largest_cluster(verts[ridx], args.cluster_r)
+        region_mask[ridx[~mkeep]] = False
+        print(f"[filter] 聚类(r={args.cluster_r}) 保留最大簇 {nkeep}/{len(ridx)} 顶点, "
+              f"剔除 {len(ridx) - nkeep} 个噪声点")
+
+    pts, inner = densify_points(verts, faces, region_mask)
+    print(f"[region] 顶点 {int(region_mask.sum())}, 内部面 {len(inner)}, "
+          f"采样点 {len(pts)}, bbox {np.round(pts.max(0) - pts.min(0), 3)}")
+
+    print("[delaunay] 四面体化中 ...")
+    st = AlphaStructure(pts)
+    print(f"[delaunay] 四面体 {st.T}, 采样密度决定的最小可行 alpha ~ "
+          f"{np.percentile(st.circum, 50):.3f}")
+
+    if args.mode == "case":
+        # ---- 眼镜盒模式: 最外侧点连成凸外壳 + 栅格化高斯隐式面 + MC -----
+        mesh, hull = build_case(pts, st.tri, args.res, args.sigma, args.case_faces)
+        kd = cKDTree(mesh.vertices)
+        dd, _ = kd.query(pts, k=1)
+        print(f"[case] 凸包体积 {hull.volume:.1f} -> 壳体积 {mesh.volume:.1f} "
+              f"({mesh.volume / hull.volume * 100:.1f}%)")
+        print(f"[shell] 顶点 {len(mesh.vertices)}, 面 {len(mesh.faces)}, "
+              f"水密 {mesh.is_watertight}, 欧拉数 {mesh.euler_number}, "
+              f"bbox {np.round(mesh.extents, 3)} (区域 bbox {np.round(pts.max(0) - pts.min(0), 3)})")
+        print(f"[fit] 采样点到壳最近距离: mean={dd.mean():.4f}, "
+              f"p99={np.percentile(dd, 99):.4f}, max={dd.max():.4f}")
+        return mesh
+
+    # ---- 最小 alpha 二分搜索 ------------------------------------------
+    if args.alpha is not None:
+        alpha = args.alpha
+        keep, ok_used, ok_conn = st.status(alpha)
+        print(f"[alpha] 手动 alpha={alpha:.4f}, 连通={ok_conn}, 覆盖={ok_used}")
+    else:
+        cand = np.unique(st.circum)
+        lo, hi = 0, len(cand) - 1
+        keep, ok_used, ok_conn = st.status(cand[hi])
+        assert ok_conn and ok_used, "最大 alpha 仍不连通, 数据异常"
+        while lo < hi:
+            mid = (lo + hi) // 2
+            k, u, c = st.status(cand[mid])
+            if u and c:
+                hi, keep, ok_used, ok_conn = mid, k, u, c
+            else:
+                lo = mid + 1
+        alpha = float(cand[lo])
+        print(f"[alpha] 最小连通 alpha = {alpha:.4f}")
+
+    # ---- 边界抽取 + 掐点拓扑修复 (仅局部补四面体, 不牺牲最小性) ---------
+    def try_extract(keep):
+        """抽边界; 若有非 2 重边(掐点), 把该边周围的四面体扇补进复合体再试。"""
+        keep = keep.copy()
+        for _ in range(12):
+            bf = st.boundary_faces(keep)
+            m = trimesh.Trimesh(vertices=pts, faces=bf, process=False)
+            uniq, cnt = np.unique(m.edges_sorted, axis=0, return_counts=True)
+            bad = uniq[cnt != 2]
+            if len(bad) == 0:
+                return keep, m, m.is_watertight
+            added = 0
+            for a, b in bad:
+                fan = (st.tets == a).any(axis=1) & (st.tets == b).any(axis=1) & ~keep
+                keep |= fan
+                added += int(fan.sum())
+            if added == 0:                    # 无可补的四面体, 修复失败
+                return keep, m, False
+        return keep, m, m.is_watertight
+
+    keep, mesh, ok = try_extract(keep)
+    if not ok:
+        for _ in range(6):                    # 兜底: 逐步放大 alpha 重试
+            alpha = alpha * 1.1
+            k1, _, _ = st.status(alpha)
+            keep, mesh, ok = try_extract(k1)
+            print(f"[alpha] 局部修复不足, alpha 上调 10% -> {alpha:.4f}")
+            if ok:
+                break
+    if not mesh.is_watertight:
+        print("[warn] 仍未水密, 直接导出当前结果")
+
+    # ---- 外壳模式: 搜最小 alpha 使边界为单一 genus-0 闭合外壳 ----------
+    if args.mode == "outer" and mesh.euler_number != 2:
+        def pred(a):
+            """alpha=a 时边界(修复掐点后)是否已是单一 genus-0 外壳。"""
+            k, _, _ = st.status(a)
+            _, m, _ = try_extract(k)
+            parts = m.split(only_watertight=False)
+            outer = [p for p in parts if p.volume > 0]
+            if len(outer) != 1:
+                return False, None
+            m = outer[0]
+            if m.is_watertight and m.euler_number == 2:
+                return True, m
+            return False, None
+
+        a_lo, a_hi, found = alpha, None, None          # a_lo: 当前非外壳
+        a = alpha
+        for _ in range(20):                            # 指数上扫找可行 alpha
+            a = a * 1.3
+            okf, mf = pred(a)
+            if okf:
+                a_hi, found = a, mf
+                break
+        if a_hi is None:
+            print("[warn] 未找到 genus-0 外壳, 保持紧凑壳")
+        else:
+            for _ in range(8):                         # 二分细化到最小外壳 alpha
+                mid = (a_lo + a_hi) / 2
+                okf, mf = pred(mid)
+                if okf:
+                    a_hi, found = mid, mf
+                else:
+                    a_lo = mid
+            alpha, mesh = a_hi, found
+            print(f"[alpha] 外壳模式: 最小 genus-0 外壳 alpha = {alpha:.4f}")
+    vol0 = mesh.volume
+
+    # 去掉内腔壳 (负体积分量), 只留外壳
+    parts = mesh.split(only_watertight=False)
+    outer = [p for p in parts if p.volume > 0]
+    if len(outer) != len(parts):
+        print(f"[shell] 去除 {len(parts) - len(outer)} 个内腔面片")
+    mesh = max(outer, key=lambda m: m.volume) if len(outer) > 1 else outer[0]
+
+    # ---- Taubin 平滑 (不收缩) ------------------------------------------
+    if args.smooth > 0:
+        before = mesh.vertices.copy()
+        filter_taubin(mesh, lamb=0.5, nu=0.5, iterations=args.smooth)
+        drift = np.linalg.norm(mesh.vertices - before, axis=1)
+        print(f"[smooth] Taubin x{args.smooth}, 顶点漂移 mean={drift.mean():.4f} max={drift.max():.4f}")
+
+    # ---- 数值校验 ------------------------------------------------------
+    prox = trimesh.proximity.ProximityQuery(mesh)
+    d = np.abs(prox.signed_distance(pts))
+    print(f"[shell] 顶点 {len(mesh.vertices)}, 面 {len(mesh.faces)}, "
+          f"水密 {mesh.is_watertight}, 欧拉数 {mesh.euler_number}")
+    print(f"[shell] 体积 {mesh.volume:.1f} (平滑前 {vol0:.1f}), "
+          f"bbox {np.round(mesh.extents, 3)}")
+    print(f"[fit] 采样点到壳距离 |d|: mean={d.mean():.4f}, "
+          f"p99={np.percentile(d, 99):.4f}, max={d.max():.4f}")
+    return mesh
+
+
+# ------------------------------------------------------ 主流程
+def main():
+    ap = argparse.ArgumentParser(description="生成每个眼镜实例的最小壳 (alpha shape)")
+    here = os.path.dirname(os.path.abspath(__file__))
+    ap.add_argument("--obj", default=None,
+                    help="obj 路径; 缺省依次取 data.json 的 input 字段 / data 目录内唯一 obj")
+    ap.add_argument("--data", default=os.path.join(here, "data.json"),
+                    help="多实例 data.json (instances[].indices), 每实例单独处理")
+    ap.add_argument("--index", default=None,
+                    help="兼容旧的单实例 index.txt 顶点索引文件; 指定时优先于 --data")
+    ap.add_argument("--out", default=os.path.join(here, "shell.obj"),
+                    help="输出路径, 所有实例的壳合并导出为这一个文件 (默认 shell.obj)")
+    ap.add_argument("--mode", choices=["case", "outer", "tight"], default="case",
+                    help="case=眼镜盒式光滑凸外壳(最外侧连成一体), outer=跟随外形的 genus-0 外壳, "
+                         "tight=贴表面的紧身壳 (默认 case)")
+    ap.add_argument("--res", type=int, default=128, help="case 模式体素分辨率 (默认 128)")
+    ap.add_argument("--sigma", type=float, default=2.5, help="case 模式高斯平滑 sigma, 体素单位 (默认 2.5)")
+    ap.add_argument("--case-faces", type=int, default=25000, help="case 模式简化目标面数, 0=不简化 (默认 25000)")
+    ap.add_argument("--alpha", type=float, default=None,
+                    help="outer/tight 模式手动指定 alpha; 缺省自动搜索")
+    ap.add_argument("--smooth", type=int, default=10,
+                    help="Taubin 平滑迭代数, 0=关闭取严格最小壳 (默认 10)")
+    ap.add_argument("--filter", choices=["main", "all"], default="main",
+                    help="main=只保留最大空间簇(去 SAM3 噪声), all=用全部 index 点 (默认 main)")
+    ap.add_argument("--cluster-r", type=float, default=2.5,
+                    help="噪声过滤的聚类半径, 世界单位 (默认 2.5)")
+    ap.add_argument("--one-based", action="store_true",
+                    help="索引为 1-based (默认按 Blender 导出的 0-based 读取)")
+    args = ap.parse_args()
+
+    instances, input_obj = load_data(args)
+    obj_path = resolve_obj(args, input_obj)
+    verts, faces = load_obj(obj_path)
+    print(f"[load] {obj_path}: 顶点 {len(verts)}, 面 {len(faces)}, "
+          f"待处理实例 {len(instances)} 个")
+
+    meshes, n_ok = [], 0
+    for name, idx in instances:
+        print(f"\n===== 实例 {name}: {len(idx)} 顶点 =====")
+        try:
+            meshes.append(run_instance(verts, faces, idx, args))
+            n_ok += 1
+        except Exception as e:
+            print(f"[error] 实例 {name} 处理失败, 跳过: {e}")
+    if not meshes:
+        raise SystemExit("[error] 所有实例处理失败, 未导出任何壳")
+
+    # ---- 合并所有实例的壳 (同一坐标系, 直接拼接), 导出单个 obj ----
+    combined = meshes[0] if len(meshes) == 1 else trimesh.util.concatenate(meshes)
+    combined.export(args.out)
+    print(f"\n[summary] 完成 {n_ok}/{len(instances)} 个实例, 合并导出: {args.out} "
+          f"(顶点 {len(combined.vertices)}, 面 {len(combined.faces)}, "
+          f"连通分量 {len(combined.split())})")
+
+
+if __name__ == "__main__":
+    main()

libs/MVS/SceneTexture.cpp

@@ -795,6 +795,16 @@ public:
 		const Mesh::TexCoordArr& existingTexcoords,  // 添加已有UV参数
 		const Mesh::TexIndexArr& existingTexindices   // 添加已有纹理索引参数
 	);
+
+	// 在 MeshTexture 类定义中添加
+	struct VertexColorInfo {
+		bool visibleInAnyView = false;  // 至少在1个视图中可见
+		cv::Vec3f color;                // 最终颜色（扩散后）
+		cv::Vec3f mainColor;            // 主视图颜色
+		std::vector<int> neighborVertices;  // 邻接顶点
+	};
+	std::vector<VertexColorInfo> m_vertexColors;
+
 	std::vector<uint32_t> faceToPatchID;  // faceID → 新 patchID 映射
 	void MergeSameViewPatches();
 	void BuildSeamEdgesFromFaceToPatchID();
@@ -16410,12 +16420,13 @@ bool MeshTexture::RasterizeVirtualFaces(
     int textureSize = ComputeOptimalTextureSizeAdaptive(virtualFaceMap, virtualFaceViews, nTextureSizeMultiple);
     if (textureSize < 1024)  textureSize = 1024;
     if (textureSize > 16384) textureSize = 16384;
-    textureSize = std::max(textureSize, 4096);
+    textureSize = std::max(textureSize, 8192);
 
     // 2. atlas + score
     outTextures.emplace_back(textureSize, textureSize);
     Image8U3& atlas = outTextures.back();
     atlas.setTo(cv::Scalar(colEmpty.b, colEmpty.g, colEmpty.r));
+
     m_texelScores.assign(textureSize * textureSize, TexelScore{});
     m_texelPatchID.assign(textureSize * textureSize, NO_ID);
 
@@ -16427,14 +16438,40 @@ bool MeshTexture::RasterizeVirtualFaces(
     }
     currentTextureSize = textureSize;
 
-	VERBOSE("[Raster] Propagating main views for consistency...");
-	PropagateMainViews(virtualFaceMap, virtualFaceViews, virtualFaceViewWeights, faceViewsHighFreq);
+    VERBOSE("[Raster] Propagating main views for consistency...");
+    PropagateMainViews(virtualFaceMap, virtualFaceViews, virtualFaceViewWeights, faceViewsHighFreq);
 
     if (!m_gainsEstimated) EstimateGlobalPhotometricCorrection();
 
+    // ★★★ 提前初始化顶点颜色结构 ★★★
+    m_vertexColors.resize(scene.mesh.vertices.size());
+    for (auto& vc : m_vertexColors) {
+        vc.visibleInAnyView = false;
+        vc.color = cv::Vec3f(0, 0, 0);
+        vc.mainColor = cv::Vec3f(0, 0, 0);
+        vc.neighborVertices.clear();
+    }
+    
+    // 建立顶点邻接关系
+    for (FIndex fi = 0; fi < faces.size(); ++fi) {
+        const Face& f = faces[fi];
+        for (int j = 0; j < 3; ++j) {
+            int v0 = f[j];
+            int v1 = f[(j + 1) % 3];
+            
+            auto& nb0 = m_vertexColors[v0].neighborVertices;
+            if (std::find(nb0.begin(), nb0.end(), v1) == nb0.end())
+                nb0.push_back(v1);
+            
+            auto& nb1 = m_vertexColors[v1].neighborVertices;
+            if (std::find(nb1.begin(), nb1.end(), v0) == nb1.end())
+                nb1.push_back(v0);
+        }
+    }
+
     VERBOSE("[Raster] Step 2: starting patch loop, total faces=%zu", virtualFaceMap.size());
 
-    // ========== 新增：预加载所有图像到浮点缓存（加速多次采样）==========
+    // ========== 预加载所有图像到浮点缓存 ==========
     std::vector<cv::Mat> imagesFloat(images.size());
     for (size_t vi = 0; vi < images.size(); vi++) {
         if (!images[vi].image.empty()) {
@@ -16442,26 +16479,10 @@ bool MeshTexture::RasterizeVirtualFaces(
         }
     }
 
-	// ★★★ 在循环外面预先计算全局 homography ★★★
-    cv::Mat globalH;
-    IIndex globalView = NO_ID;
-    for (int ii = 0; ii < (int)virtualFaceMap.size(); ++ii) {
-        const VirtualFaceGeometry& g = m_virtualFaceGeometries[ii];
-        if (g.isValid && !virtualFaceViews[ii].empty()) {
-            globalView = virtualFaceViews[ii][0];
-            auto it = g.viewHomographies.find(globalView);
-            if (it != g.viewHomographies.end()) {
-                globalH = it->second.clone();
-                break;
-            }
-        }
-    }
-    DEBUG_EXTRA("[TEST] Using global view %d for all faces", globalView);
-
     // 3. 逐三角形光栅化（分频带合成版）
-#ifdef _USE_OPENMP
-#pragma omp parallel for schedule(dynamic)
-#endif
+    #ifdef _USE_OPENMP
+    #pragma omp parallel for schedule(dynamic)
+    #endif
     for (int i = 0; i < (int)virtualFaceMap.size(); ++i) {
 
         if (i % 10000 == 0)
@@ -16472,13 +16493,12 @@ bool MeshTexture::RasterizeVirtualFaces(
         if (!geom.isValid || vf.faces.empty() || virtualFaceViews[i].empty()) continue;
 
         // ========== 分频带视图获取 ==========
-        IIndex mainView = faceViewsHighFreq[i];  // 高频视图（主视图）
-        const auto& lowFreqViews = faceViewsLowFreq[i];  // 低频视图列表
+        IIndex mainView = faceViewsHighFreq[i];
+        const auto& lowFreqViews = faceViewsLowFreq[i];
         
         if (mainView == NO_ID || mainView >= (IIndex)images.size()) continue;
         if (images[mainView].image.empty()) continue;
         
-        // 确保低频视图也有效
         bool hasLowFreqViews = false;
         for (const auto& lfv : lowFreqViews) {
             if (lfv.viewId < (IIndex)images.size() && !images[lfv.viewId].image.empty()) {
@@ -16489,7 +16509,7 @@ bool MeshTexture::RasterizeVirtualFaces(
 
         float currentScore = virtualFaceViewWeights[i].empty() ? -1.0f : virtualFaceViewWeights[i][0];
 
-		int pad = 8;
+        int pad = 8;
         int minX = std::max(0, (int)floor(geom.uvBounds.ptMin.x() * textureSize - pad));
         int maxX = std::min(textureSize - 1, (int)ceil(geom.uvBounds.ptMax.x() * textureSize + pad));
         int minY = std::max(0, (int)floor(geom.uvBounds.ptMin.y() * textureSize - pad));
@@ -16503,62 +16523,52 @@ bool MeshTexture::RasterizeVirtualFaces(
             continue;
         }
 
-// ========== 用重心坐标 + 3D 投影替代 homography ==========
+        // ========== 用重心坐标 + 3D 投影 ==========
         cv::Mat mapX(patchH, patchW, CV_32FC1), mapY(patchH, patchW, CV_32FC1);
         mapX.setTo(-1);
         mapY.setTo(-1);
 
-        // 获取该虚拟面对应的原始三角面
         if (vf.faces.empty()) continue;
         FIndex faceIdx = vf.faces[0];
         const Face& face = scene.mesh.faces[faceIdx];
 
-        // 获取相机参数
+        Point3f v[3];
+        TexCoord uv[3];
+        for (int j = 0; j < 3; ++j) {
+            v[j] = scene.mesh.vertices[face[j]];
+            uv[j] = scene.mesh.faceTexcoords[faceIdx * 3 + j];
+        }
+
         const Camera& cam = images[mainView].camera;
         bool camValid = (cam.K(0,0) > 0 && cam.K(1,1) > 0);
         if (!camValid) {
-            // fallback 到 homography
             const float* H = geom.homography.ptr<float>();
             for (int y = minY; y <= maxY; ++y) {
                 for (int x = minX; x <= maxX; ++x) {
-                    float u = (float)x / (float)textureSize, v = (float)y / (float)textureSize;
-                    float w = H[6]*u + H[7]*v + H[8];
+                    float u = (float)x / (float)textureSize, v_coord = (float)y / (float)textureSize;
+                    float w = H[6]*u + H[7]*v_coord + H[8];
                     if (std::abs(w) < 1e-12f) { mapX.at<float>(y-minY,x-minX)=-1; mapY.at<float>(y-minY,x-minX)=-1; continue; }
-                    mapX.at<float>(y-minY,x-minX) = (H[0]*u+H[1]*v+H[2])/w;
-                    mapY.at<float>(y-minY,x-minX) = (H[3]*u+H[4]*v+H[5])/w;
+                    mapX.at<float>(y-minY,x-minX) = (H[0]*u+H[1]*v_coord+H[2])/w;
+                    mapY.at<float>(y-minY,x-minX) = (H[3]*u+H[4]*v_coord+H[5])/w;
                 }
             }
         } else {
-			// 提取 3D 顶点和 UV
-			Point3f v[3];
-			TexCoord uv[3];
-			for (int j = 0; j < 3; ++j) {
-				v[j] = scene.mesh.vertices[face[j]];
-				uv[j] = scene.mesh.faceTexcoords[faceIdx * 3 + j];  // 注意这里的索引方式
-			}
-
-            // 对 patch 内每个像素做重心坐标投影
             for (int py = minY; py <= maxY; ++py) {
                 for (int px = minX; px <= maxX; ++px) {
                     float u = (float)px / (float)textureSize;
                     float v_coord = (float)py / (float)textureSize;
                     
-                    // 计算重心坐标
                     float b0, b1, b2;
-                    if (!ComputeBarycentric(u, v_coord, uv[0], uv[1], uv[2], b0, b1, b2)) {
-                        continue;
-                    }
+                    bool inside = ComputeBarycentric(u, v_coord, uv[0], uv[1], uv[2], b0, b1, b2);
+                    if (!inside) continue;
                     
-                    // 3D 插值点
                     Point3f P = v[0] * b0 + v[1] * b1 + v[2] * b2;
                     
-                    // 投影到图像
-					Point2f imgPt = cam.ProjectPointP(P);
-					// 检查是否在图像有效范围内
-					if (imgPt.x < 0 || imgPt.x >= images[mainView].image.cols ||
-						imgPt.y < 0 || imgPt.y >= images[mainView].image.rows) {
-						continue;
-					}
+                    Point2f imgPt = cam.ProjectPointP(P);
+                    if (imgPt.x < 0 || imgPt.x >= images[mainView].image.cols ||
+                        imgPt.y < 0 || imgPt.y >= images[mainView].image.rows) {
+                        continue;
+                    }
                     
                     mapX.at<float>(py - minY, px - minX) = imgPt.x;
                     mapY.at<float>(py - minY, px - minX) = imgPt.y;
@@ -16566,94 +16576,217 @@ bool MeshTexture::RasterizeVirtualFaces(
             }
         }
 
-        // ========== 从主视图采样（全分辨率）==========
+        // ========== 从主视图采样 ==========
         cv::Mat mainPatch;
         cv::remap(images[mainView].image, mainPatch, mapX, mapY, cv::INTER_LINEAR, cv::BORDER_CONSTANT, cv::Scalar(0,0,0));
 
-		// ========== 从辅助视图采样（全部用主视图的映射）==========
-		cv::Mat lowFreqPatch(patchH, patchW, CV_32FC3, cv::Scalar(0,0,0));
-		if (hasLowFreqViews) {
-			float totalWeight = 0;
-			for (const auto& lfv : lowFreqViews) {
-				if (lfv.viewId >= (IIndex)images.size() || images[lfv.viewId].image.empty()) continue;
-				
-				// ★★★ 所有视图都用主视图的 mapX/mapY，不再各自构建 homography ★★★
-				cv::Mat auxPatch;
-				cv::remap(images[lfv.viewId].image, auxPatch, mapX, mapY, 
-						cv::INTER_LINEAR, cv::BORDER_CONSTANT, cv::Scalar(0,0,0));
-				
-				cv::Mat auxFloat;
-				auxPatch.convertTo(auxFloat, CV_32FC3, 1.0f / 255.0f);
-				
-				lowFreqPatch += auxFloat * lfv.weight;
-				totalWeight += lfv.weight;
-			}
-			if (totalWeight > 0) {
-				lowFreqPatch /= totalWeight;
-			} else {
-				cv::Mat mainFloat;
-				mainPatch.convertTo(mainFloat, CV_32FC3, 1.0f / 255.0f);
-				lowFreqPatch = mainFloat.clone();
-			}
-		} else {
-			cv::Mat mainFloat;
-			mainPatch.convertTo(mainFloat, CV_32FC3, 1.0f / 255.0f);
-			lowFreqPatch = mainFloat.clone();
-		}
-		
-		// ========== 多视图高频混合 ==========
-		cv::Mat mainFloat;
-		mainPatch.convertTo(mainFloat, CV_32FC3, 1.0f / 255.0f);
-
-		// 主视图高频
-		cv::Mat mainLowFreq;
-		cv::GaussianBlur(mainFloat, mainLowFreq, cv::Size(31, 31), 10.0);
-		cv::Mat highFreqAccum = (mainFloat - mainLowFreq) * 0.6f;
-		float highFreqTotalWeight = 0.6f;
-
-		// 辅助视图也贡献高频（也用主视图的 mapX/mapY）
-		if (hasLowFreqViews) {
-			for (const auto& lfv : lowFreqViews) {
-				if (lfv.viewId == mainView) continue;
-				if (lfv.viewId >= (IIndex)images.size() || images[lfv.viewId].image.empty()) continue;
-				
-				// ★★★ 用主视图的 mapX/mapY（3D投影的那个），不是各视图自己的homography ★★★
-				cv::Mat auxPatch;
-				cv::remap(images[lfv.viewId].image, auxPatch, mapX, mapY, 
-						cv::INTER_LINEAR, cv::BORDER_CONSTANT, cv::Scalar(0,0,0));
-				
-				cv::Mat auxFloat;
-				auxPatch.convertTo(auxFloat, CV_32FC3, 1.0f / 255.0f);
-				
-				cv::Mat auxLowFreq;
-				cv::GaussianBlur(auxFloat, auxLowFreq, cv::Size(31, 31), 10.0);
-				
-				float auxWeight = lfv.weight * 0.4f;
-				highFreqAccum += (auxFloat - auxLowFreq) * auxWeight;
-				highFreqTotalWeight += auxWeight;
-			}
-		}
-
-		if (highFreqTotalWeight > 0) {
-			highFreqAccum /= highFreqTotalWeight;
-		}
-
-		cv::Mat highFreq = highFreqAccum;
-
-		// 在合成之前，对 lowFreqPatch 再做一次轻微平滑
-		// cv::GaussianBlur(lowFreqPatch, lowFreqPatch, cv::Size(5, 5), 1.0);
-
-        // ========== 合成：低频用多视图融合，高频用主视图细节 ==========
-        // cv::Mat finalPatch = lowFreqPatch + highFreq;
-		// cv::Mat finalPatch = lowFreqPatch;
-		// cv::Mat finalPatch = highFreq.clone();
-		// finalPatch = cv::Mat::zeros(finalPatch.size(), finalPatch.type());
-
-		cv::Mat finalPatch = mainFloat;  // 只用主视图
-		// cv::Mat finalPatch = lowFreqPatch * 0.3f + highFreq * 0.7f;
-		// 低频只做主视图，不做多视图平均（避免发灰）
-		// cv::Mat finalPatch = mainFloat * 0.7f + highFreq * 0.3f;
+        // ========== 从辅助视图采样 ==========
+        cv::Mat lowFreqPatch(patchH, patchW, CV_32FC3, cv::Scalar(0,0,0));
+        if (hasLowFreqViews) {
+            float totalWeight = 0;
+            for (const auto& lfv : lowFreqViews) {
+                if (lfv.viewId >= (IIndex)images.size() || images[lfv.viewId].image.empty()) continue;
+                
+                const Camera& auxCam = images[lfv.viewId].camera;
+                bool auxCamValid = (auxCam.K(0,0) > 0 && auxCam.K(1,1) > 0);
+                
+                cv::Mat auxMapX(patchH, patchW, CV_32FC1), auxMapY(patchH, patchW, CV_32FC1);
+                auxMapX.setTo(-1);
+                auxMapY.setTo(-1);
+                
+                if (!auxCamValid) {
+                    const auto& geomRef = m_virtualFaceGeometries[i];
+                    auto it = geomRef.viewHomographies.find(lfv.viewId);
+                    if (it == geomRef.viewHomographies.end()) continue;
+                    const float* H = it->second.ptr<float>();
+                    for (int y = minY; y <= maxY; ++y) {
+                        for (int x = minX; x <= maxX; ++x) {
+                            float u = (float)x / (float)textureSize, v_coord = (float)y / (float)textureSize;
+                            float w = H[6]*u + H[7]*v_coord + H[8];
+                            if (std::abs(w) < 1e-12f) { 
+                                auxMapX.at<float>(y-minY,x-minX) = -1; 
+                                auxMapY.at<float>(y-minY,x-minX) = -1; 
+                                continue; 
+                            }
+                            auxMapX.at<float>(y-minY,x-minX) = (H[0]*u+H[1]*v_coord+H[2])/w;
+                            auxMapY.at<float>(y-minY,x-minX) = (H[3]*u+H[4]*v_coord+H[5])/w;
+                        }
+                    }
+                } else {
+                    for (int py = minY; py <= maxY; ++py) {
+                        for (int px = minX; px <= maxX; ++px) {
+                            float u = (float)px / (float)textureSize;
+                            float v_coord = (float)py / (float)textureSize;
+                            
+                            float b0, b1, b2;
+                            bool inside = ComputeBarycentric(u, v_coord, uv[0], uv[1], uv[2], b0, b1, b2);
+                            if (!inside) continue;
+                            
+                            Point3f P = v[0] * b0 + v[1] * b1 + v[2] * b2;
+                            
+                            Point2f imgPt = auxCam.ProjectPointP(P);
+                            if (imgPt.x < 0 || imgPt.x >= images[lfv.viewId].image.cols ||
+                                imgPt.y < 0 || imgPt.y >= images[lfv.viewId].image.rows) {
+                                continue;
+                            }
+                            
+                            auxMapX.at<float>(py - minY, px - minX) = imgPt.x;
+                            auxMapY.at<float>(py - minY, px - minX) = imgPt.y;
+                        }
+                    }
+                }
+                
+                cv::Mat auxPatch;
+                cv::remap(images[lfv.viewId].image, auxPatch, auxMapX, auxMapY, 
+                        cv::INTER_LINEAR, cv::BORDER_CONSTANT, cv::Scalar(0,0,0));
+                
+                cv::Mat auxFloat;
+                auxPatch.convertTo(auxFloat, CV_32FC3, 1.0f / 255.0f);
+                
+                lowFreqPatch += auxFloat * lfv.weight;
+                totalWeight += lfv.weight;
+            }
+            if (totalWeight > 0) {
+                lowFreqPatch /= totalWeight;
+            } else {
+                cv::Mat mainFloat;
+                mainPatch.convertTo(mainFloat, CV_32FC3, 1.0f / 255.0f);
+                lowFreqPatch = mainFloat.clone();
+            }
+        } else {
+            cv::Mat mainFloat;
+            mainPatch.convertTo(mainFloat, CV_32FC3, 1.0f / 255.0f);
+            lowFreqPatch = mainFloat.clone();
+        }
         
+        // ========== 多视图高频混合 ==========
+        cv::Mat mainFloat;
+        mainPatch.convertTo(mainFloat, CV_32FC3, 1.0f / 255.0f);
+
+        cv::Mat mainLowFreq;
+       	cv::GaussianBlur(mainFloat, mainLowFreq, cv::Size(5, 5), 1.5);
+        cv::Mat highFreqAccum = (mainFloat - mainLowFreq) * 0.9f;
+        float highFreqTotalWeight = 0.6f;
+
+        if (hasLowFreqViews) {
+            for (const auto& lfv : lowFreqViews) {
+                if (lfv.viewId == mainView) continue;
+                if (lfv.viewId >= (IIndex)images.size() || images[lfv.viewId].image.empty()) continue;
+                
+                cv::Mat auxPatch;
+                cv::remap(images[lfv.viewId].image, auxPatch, mapX, mapY, 
+                        cv::INTER_LINEAR, cv::BORDER_CONSTANT, cv::Scalar(0,0,0));
+                
+                cv::Mat auxFloat;
+                auxPatch.convertTo(auxFloat, CV_32FC3, 1.0f / 255.0f);
+                
+                cv::Mat auxLowFreq;
+                cv::GaussianBlur(auxFloat, auxLowFreq, cv::Size(31, 31), 10.0);
+                
+                float auxWeight = lfv.weight * 0.4f;
+                highFreqAccum += (auxFloat - auxLowFreq) * auxWeight;
+                highFreqTotalWeight += auxWeight;
+            }
+        }
+
+        if (highFreqTotalWeight > 0) {
+            highFreqAccum /= highFreqTotalWeight;
+        }
+
+        cv::Mat highFreq = highFreqAccum;
+
+		// ========== 每个顶点独立判断可见性（增强版）==========
+		bool vertexVisible[3] = {false, false, false};
+		cv::Vec3f vertexBestColor[3] = {cv::Vec3f(0,0,0), cv::Vec3f(0,0,0), cv::Vec3f(0,0,0)};
+
+		// 检查主视图
+		for (int j = 0; j < 3; ++j) {
+			int vx = (int)(uv[j].x * textureSize) - minX;
+			int vy = (int)(uv[j].y * textureSize) - minY;
+			if (vx >= 0 && vx < patchW && vy >= 0 && vy < patchH) {
+				cv::Vec3f color = mainFloat.at<cv::Vec3f>(vy, vx);
+				if (color[0] > 0.0f || color[1] > 0.0f || color[2] > 0.0f) {
+					vertexVisible[j] = true;
+					vertexBestColor[j] = color;
+				}
+			}
+		}
+
+		// 如果主视图不可见，检查辅助视图
+		if (hasLowFreqViews && (!vertexVisible[0] || !vertexVisible[1] || !vertexVisible[2])) {
+			for (const auto& lfv : lowFreqViews) {
+				if (lfv.viewId >= (IIndex)images.size() || images[lfv.viewId].image.empty()) continue;
+				
+				// 对该辅助视图重新采样该顶点位置
+				const Camera& auxCam = images[lfv.viewId].camera;
+				bool auxCamValid = (auxCam.K(0,0) > 0 && auxCam.K(1,1) > 0);
+				
+				for (int j = 0; j < 3; ++j) {
+					if (vertexVisible[j]) continue; // 已经可见就跳过
+					
+					// 用3D投影得到该顶点在辅助视图中的像素位置
+					Point2f imgPt;
+					if (auxCamValid) {
+						imgPt = auxCam.ProjectPointP(v[j]);
+					} else {
+						// 用 homography fallback
+						const auto& geomRef = m_virtualFaceGeometries[i];
+						auto it = geomRef.viewHomographies.find(lfv.viewId);
+						if (it == geomRef.viewHomographies.end()) continue;
+						
+						const float* H = it->second.ptr<float>();
+						float u = uv[j].x;
+						float v_coord = uv[j].y;
+						float w = H[6]*u + H[7]*v_coord + H[8];
+						if (std::abs(w) < 1e-12f) continue;
+						
+						float projX = (H[0]*u + H[1]*v_coord + H[2]) / w;
+						float projY = (H[3]*u + H[4]*v_coord + H[5]) / w;
+						imgPt = Point2f(projX, projY);
+					}
+					
+					if (imgPt.x >= 0 && imgPt.x < images[lfv.viewId].image.cols &&
+						imgPt.y >= 0 && imgPt.y < images[lfv.viewId].image.rows) {
+						// 采样该像素颜色
+						cv::Vec3f auxColor = imagesFloat[lfv.viewId].at<cv::Vec3f>((int)imgPt.y, (int)imgPt.x);
+						if (auxColor[0] > 0.0f || auxColor[1] > 0.0f || auxColor[2] > 0.0f) {
+							vertexVisible[j] = true;
+							vertexBestColor[j] = auxColor;
+						}
+					}
+				}
+				
+				// 如果三个顶点都已可见，提前退出
+				if (vertexVisible[0] && vertexVisible[1] && vertexVisible[2])
+					break;
+			}
+		}
+
+		#pragma omp critical (vertex_visibility)
+		{
+			for (int j = 0; j < 3; ++j) {
+				VIndex vi = face[j];
+				if (vi < m_vertexColors.size() && vertexVisible[j]) {
+					m_vertexColors[vi].visibleInAnyView = true;
+					m_vertexColors[vi].mainColor = vertexBestColor[j];
+				}
+			}
+		}
+
+        // ========== 空洞填充与合成 ==========
+        cv::Mat filledLowFreq = lowFreqPatch.clone();
+        for (int y = 0; y < patchH; ++y) {
+            for (int x = 0; x < patchW; ++x) {
+                cv::Vec3f color = filledLowFreq.at<cv::Vec3f>(y, x);
+                if (color[0] == 0.0f && color[1] == 0.0f && color[2] == 0.0f) {
+                    filledLowFreq.at<cv::Vec3f>(y, x) = mainFloat.at<cv::Vec3f>(y, x);
+                }
+            }
+        }
+        
+        cv::Mat finalPatch = filledLowFreq * 0.8f + mainFloat * 0.2f;
+        // cv::Mat finalPatch = mainFloat;
+
         // 钳制到 [0, 1]
         cv::Mat clampedPatch;
         cv::max(finalPatch, 0.0f, clampedPatch);
@@ -16663,50 +16796,255 @@ bool MeshTexture::RasterizeVirtualFaces(
         cv::Mat finalPatch8U;
         clampedPatch.convertTo(finalPatch8U, CV_8UC3, 255.0f);
 
-		// 在写入 atlas 时，对边界像素做衰减混合
-		#pragma omp critical (atlas_write)
-		{
-			for (int y = 0; y < patchH; ++y) {
-				for (int x = 0; x < patchW; ++x) {
-					cv::Vec3b color = finalPatch8U.at<cv::Vec3b>(y, x);
-					
-					if (color[0] == 0 && color[1] == 0 && color[2] == 0) continue;
+        // ========== 写入 atlas ==========
+        #pragma omp critical (atlas_write)
+        {
+            for (int y = 0; y < patchH; ++y) {
+                for (int x = 0; x < patchW; ++x) {
+                    cv::Vec3b color = finalPatch8U.at<cv::Vec3b>(y, x);
 
-					int atlasX = x + minX, atlasY = y + minY;
-					if (atlasX<0||atlasX>=textureSize||atlasY<0||atlasY>=textureSize) continue;
-					size_t idx = atlasY * textureSize + atlasX;
+                    if (color[0] == 0 && color[1] == 0 && color[2] == 0) continue;
 
-					// ===== 新增：边界衰减权重 =====
-					float weight = 1.0f;
-					float distToEdgeX = std::min((float)x, (float)(patchW - 1 - x));
-					float distToEdgeY = std::min((float)y, (float)(patchH - 1 - y));
-					float minDist = std::min(distToEdgeX, distToEdgeY);
-					
-					if (minDist < pad) {
-						// 边界 8 像素内做线性衰减
-						weight = minDist / (float)pad;
-						weight = weight * weight;  // 平方衰减更平滑
-					}
-					
-					// 如果已有更高分的 texel，跳过（保持高分优先）
-					if (currentScore > m_texelScores[idx].score) {
-						// 如果有旧值，做加权混合
-						if (m_texelScores[idx].score > 0 && weight < 1.0f) {
-							cv::Vec3b oldColor = atlas.at<cv::Vec3b>(atlasY, atlasX);
-							float oldWeight = 1.0f - weight;
-							color[0] = (uchar)(color[0] * weight + oldColor[0] * oldWeight);
-							color[1] = (uchar)(color[1] * weight + oldColor[1] * oldWeight);
-							color[2] = (uchar)(color[2] * weight + oldColor[2] * oldWeight);
-						}
-						
-						atlas.at<cv::Vec3b>(atlasY, atlasX) = color;
-						m_texelScores[idx].score = currentScore;
-						m_texelScores[idx].viewID = mainView;
-						m_texelPatchID[idx] = i;
-					}
-				}
-			}
-		}
+                    int atlasX = x + minX, atlasY = y + minY;
+                    if (atlasX<0||atlasX>=textureSize||atlasY<0||atlasY>=textureSize) continue;
+                    size_t idx = atlasY * textureSize + atlasX;
+
+                    // 边界衰减权重
+                    float weight = 1.0f;
+                    float distToEdgeX = std::min((float)x, (float)(patchW - 1 - x));
+                    float distToEdgeY = std::min((float)y, (float)(patchH - 1 - y));
+                    float minDist = std::min(distToEdgeX, distToEdgeY);
+                    
+                    if (minDist < pad) {
+                        weight = minDist / (float)pad;
+                        weight = weight * weight;
+                    }
+                    
+                    if (currentScore > m_texelScores[idx].score) {
+                        atlas.at<cv::Vec3b>(atlasY, atlasX) = color;
+                        m_texelScores[idx].score = currentScore;
+                        m_texelScores[idx].viewID = mainView;
+                        m_texelPatchID[idx] = i;
+                    }
+                }
+            }
+        }
+    }
+
+    // ========== OpenMP 循环结束，执行顶点级别颜色扩散 ==========
+
+    // BFS扩散
+    std::queue<VIndex> bfsQueue;
+    std::vector<bool> visited(m_vertexColors.size(), false);
+
+    for (VIndex vi = 0; vi < m_vertexColors.size(); ++vi) {
+        if (m_vertexColors[vi].visibleInAnyView) {
+            m_vertexColors[vi].color = m_vertexColors[vi].mainColor;
+            visited[vi] = true;
+            bfsQueue.push(vi);
+        }
+    }
+
+    while (!bfsQueue.empty()) {
+        VIndex cur = bfsQueue.front();
+        bfsQueue.pop();
+        
+        for (int nb : m_vertexColors[cur].neighborVertices) {
+            if (!visited[nb]) {
+                m_vertexColors[nb].color = m_vertexColors[cur].color;
+                visited[nb] = true;
+                bfsQueue.push(nb);
+            }
+        }
+    }
+
+    // 处理孤立顶点
+    for (VIndex vi = 0; vi < m_vertexColors.size(); ++vi) {
+        if (!visited[vi]) {
+            cv::Vec3f avgColor(0, 0, 0);
+            int count = 0;
+            for (int nb : m_vertexColors[vi].neighborVertices) {
+                if (visited[nb]) {
+                    avgColor += m_vertexColors[nb].color;
+                    count++;
+                }
+            }
+            if (count > 0) {
+                m_vertexColors[vi].color = avgColor / (float)count;
+            } else {
+                m_vertexColors[vi].color = cv::Vec3f(0.35f, 0.65f, 0.95f);
+            }
+        }
+    }
+    // ★★★ 在这里插入：强制给仍然为黑色的顶点赋色 ★★★
+    for (VIndex vi = 0; vi < m_vertexColors.size(); ++vi) {
+        if (m_vertexColors[vi].color[0] == 0.0f && 
+            m_vertexColors[vi].color[1] == 0.0f && 
+            m_vertexColors[vi].color[2] == 0.0f) {
+            
+            cv::Vec3f avgColor(0, 0, 0);
+            int count = 0;
+            for (int nb : m_vertexColors[vi].neighborVertices) {
+                if (m_vertexColors[nb].color[0] > 0.0f || 
+                    m_vertexColors[nb].color[1] > 0.0f || 
+                    m_vertexColors[nb].color[2] > 0.0f) {
+                    avgColor += m_vertexColors[nb].color;
+                    count++;
+                }
+            }
+            if (count > 0) {
+                m_vertexColors[vi].color = avgColor / (float)count;
+            } else {
+                // 实在找不到，用主视图平均色
+                m_vertexColors[vi].color = cv::Vec3f(0.35f, 0.55f, 0.75f); // 灰色保底
+            }
+        }
+    }
+
+    DEBUG_EXTRA("Vertex color diffusion: %zu vertices, %zu visible, %zu invisible",
+                m_vertexColors.size(), 
+                std::count(visited.begin(), visited.end(), true),
+                std::count(visited.begin(), visited.end(), false));
+
+    // ========== 用顶点颜色修复空洞 ==========
+    #ifdef _USE_OPENMP
+    #pragma omp parallel for schedule(dynamic)
+    #endif
+    for (int atlasY = 0; atlasY < textureSize; ++atlasY) {
+        for (int atlasX = 0; atlasX < textureSize; ++atlasX) {
+            size_t idx = atlasY * textureSize + atlasX;
+            
+            // 只处理空洞像素
+            Pixel8U px = atlas(atlasY, atlasX);
+            if (px[0] != colEmpty.b || px[1] != colEmpty.g || px[2] != colEmpty.r)
+                continue;
+            
+            // 找到这个像素属于哪个面
+            int patchID = m_texelPatchID[idx];
+if (patchID >= (int)virtualFaceMap.size())
+    continue;
+
+if (patchID == NO_ID) {
+    // 尝试找最近的已填充像素
+    int bestDist = INT_MAX;
+    int bestPatchID = NO_ID;
+    int searchRad = 16;
+    for (int dy = -searchRad; dy <= searchRad; ++dy) {
+        for (int dx = -searchRad; dx <= searchRad; ++dx) {
+            int sx = atlasX + dx, sy = atlasY + dy;
+            if (sx < 0 || sx >= textureSize || sy < 0 || sy >= textureSize) continue;
+            size_t sidx = sy * textureSize + sx;
+            if (m_texelPatchID[sidx] != NO_ID) {
+                int dist = dx*dx + dy*dy;
+                if (dist < bestDist) {
+                    bestDist = dist;
+                    bestPatchID = m_texelPatchID[sidx];
+                }
+            }
+        }
+    }
+    if (bestPatchID != NO_ID) {
+        patchID = bestPatchID;
+    } else {
+        continue; // 实在找不到归属面，跳过
+    }
+}
+            
+            const VirtualFace& vf = virtualFaceMap[patchID];
+            if (vf.faces.empty()) continue;
+            
+            FIndex faceIdx = vf.faces[0];
+            const Face& face = scene.mesh.faces[faceIdx];
+            
+            // 获取顶点颜色
+            cv::Vec3f vertexColors[3];
+            for (int j = 0; j < 3; ++j) {
+                VIndex vi = face[j];
+                vertexColors[j] = (vi < m_vertexColors.size()) ? 
+                    m_vertexColors[vi].color : cv::Vec3f(0.5f, 0.5f, 0.5f);
+            }
+            
+            // 计算 UV 坐标
+            float u = (float)atlasX / (float)textureSize;
+            float v_coord = (float)atlasY / (float)textureSize;
+            
+            // 获取该面的 UV
+            TexCoord uv[3];
+            for (int j = 0; j < 3; ++j) {
+                uv[j] = scene.mesh.faceTexcoords[faceIdx * 3 + j];
+            }
+            
+            // 重心坐标插值
+            float b0, b1, b2;
+            bool inside = ComputeBarycentric(u, v_coord, uv[0], uv[1], uv[2], b0, b1, b2);
+            
+            cv::Vec3f fillColor;
+            if (inside) {
+                fillColor = vertexColors[0] * b0 + vertexColors[1] * b1 + vertexColors[2] * b2;
+            } else {
+                // 最近顶点
+                float minDist = FLT_MAX;
+                int bestVertex = 0;
+                for (int j = 0; j < 3; ++j) {
+                    float du = u - uv[j].x;
+                    float dv = v_coord - uv[j].y;
+                    float dist = du * du + dv * dv;
+                    if (dist < minDist) { minDist = dist; bestVertex = j; }
+                }
+                fillColor = vertexColors[bestVertex];
+            }
+            
+// ★★★ 新增：检查 fillColor 是否有效（接近黑色）★★★
+            bool fillIsBlack = (fillColor[0] < 0.05f && fillColor[1] < 0.05f && fillColor[2] < 0.05f);
+            if (fillIsBlack) {
+                // Fallback：从邻近非空洞像素采样
+                cv::Vec3f sampledColor(0, 0, 0);
+                int sampleCount = 0;
+                int sampleRad = 5;
+                
+                for (int dy = -sampleRad; dy <= sampleRad && sampleCount < 16; ++dy) {
+                    for (int dx = -sampleRad; dx <= sampleRad && sampleCount < 16; ++dx) {
+                        int sx = atlasX + dx, sy = atlasY + dy;
+                        if (sx < 0 || sx >= textureSize || sy < 0 || sy >= textureSize) continue;
+                        
+                        Pixel8U spx = atlas(sy, sx);
+                        // 检查是否是非空洞像素（与 colEmpty 比较）
+                        if (spx[0] != colEmpty.b || spx[1] != colEmpty.g || spx[2] != colEmpty.r) {
+                            sampledColor[0] += spx[2] / 255.0f; // R
+                            sampledColor[1] += spx[1] / 255.0f; // G
+                            sampledColor[2] += spx[0] / 255.0f; // B
+                            sampleCount++;
+                        }
+                    }
+                }
+                
+                if (sampleCount > 0) {
+                    fillColor = sampledColor / (float)sampleCount;
+                } else {
+                    // 实在不行，用该面的平均色
+                    if (patchID < (int)patchAvgColor.size()) {
+                        fillColor = cv::Vec3f(
+                            patchAvgColor[patchID][0] / 255.0f,
+                            patchAvgColor[patchID][1] / 255.0f,
+                            patchAvgColor[patchID][2] / 255.0f
+                        );
+                    } else {
+                        fillColor = cv::Vec3f(0.50f, 0.55f, 0.65f); // 最终保底
+                    }
+                }
+            }
+            
+            // 钳制并写入
+            fillColor[0] = std::max(0.0f, std::min(1.0f, fillColor[0]));
+            fillColor[1] = std::max(0.0f, std::min(1.0f, fillColor[1]));
+            fillColor[2] = std::max(0.0f, std::min(1.0f, fillColor[2]));
+            
+            atlas(atlasY, atlasX) = Pixel8U(
+                (uchar)(fillColor[2] * 255.0f),
+                (uchar)(fillColor[1] * 255.0f),
+                (uchar)(fillColor[0] * 255.0f)
+            );
+        }
     }
 
     VERBOSE("[Raster] Step 3: rasterization done, building RC patches...");
@@ -16732,7 +17070,7 @@ bool MeshTexture::RasterizeVirtualFaces(
     }
     DEBUG_EXTRA("Created %zu RC patches", rcPatches.size());
 
-    // ★ 统计每个 patch 的平均色
+    // 统计每个 patch 的平均色
     const int NP = (int)rcPatches.size();
     patchAvgColor.assign(NP, Color(0,0,0));
     std::vector<int> pixelCount(NP, 0);
@@ -16772,13 +17110,130 @@ bool MeshTexture::RasterizeVirtualFaces(
     int validPixels = 0;
     for (size_t i = 0; i < m_texelPatchID.size(); ++i)
         if (m_texelPatchID[i] != NO_ID) validPixels++;
-    DEBUG_EXTRA("m_texelPatchID: %d valid pixels out of %zu", validPixels, m_texelPatchID.size());
+		DEBUG_EXTRA("m_texelPatchID: %d valid pixels out of %zu", validPixels, m_texelPatchID.size());
+		
+	// ★★★ 全局空洞终极填充（优化版）★★★
+	cv::Mat holeMask(textureSize, textureSize, CV_8UC1);
+	for (int y = 0; y < textureSize; ++y) {
+		for (int x = 0; x < textureSize; ++x) {
+			Pixel8U px = atlas(y, x);
+			holeMask.at<uchar>(y, x) = (px[0] == colEmpty.b && px[1] == colEmpty.g && px[2] == colEmpty.r) ? 255 : 0;
+		}
+	}
 
-	// cv::imwrite("/home/algo/Documents/openMVS/data/546893/out.546893/debug_atlas_before_return.png", atlas);
+	// 多次膨胀，逐步向内填充
+	cv::Mat kernel = cv::getStructuringElement(cv::MORPH_ELLIPSE, cv::Size(5, 5));
+	cv::Mat dilatedMask;
+	cv::Mat prevMask = holeMask.clone();
 
-    return true;
+	for (int iter = 0; iter < 30; ++iter) {  // 最多膨胀30次（约150像素半径）
+		cv::dilate(prevMask, dilatedMask, kernel);
+		
+		// 找出新覆盖的区域
+		cv::Mat newPixels = dilatedMask & ~prevMask;
+		if (cv::countNonZero(newPixels) == 0) break;  // 没有新像素了
+		
+		// 对新像素赋值为最近的非空洞像素颜色
+		for (int y = 0; y < textureSize; ++y) {
+			for (int x = 0; x < textureSize; ++x) {
+				if (newPixels.at<uchar>(y, x)) {
+					// 在原始非空洞区域找最近的颜色
+					int bestDist = INT_MAX;
+					Pixel8U bestColor(128, 128, 128);
+					int rad = 3;
+					for (int dy = -rad; dy <= rad && bestDist > 1; ++dy) {
+						for (int dx = -rad; dx <= rad && bestDist > 1; ++dx) {
+							int sx = x + dx, sy = y + dy;
+							if (sx < 0 || sx >= textureSize || sy < 0 || sy >= textureSize) continue;
+							if (!holeMask.at<uchar>(sy, sx)) {  // 原始非空洞像素
+								Pixel8U spx = atlas(sy, sx);
+								int dist = dx*dx + dy*dy;
+								if (dist < bestDist) {
+									bestDist = dist;
+									bestColor = spx;
+								}
+							}
+						}
+					}
+					atlas(y, x) = bestColor;
+				}
+			}
+		}
+		
+		prevMask = dilatedMask;
+	}
+
+	// 如果还有剩余空洞，用整体平均色填充
+	cv::Scalar avgColor = cv::mean(atlas, holeMask == 0);
+	Pixel8U avgPx((uchar)avgColor[0], (uchar)avgColor[1], (uchar)avgColor[2]);
+	for (int y = 0; y < textureSize; ++y) {
+		for (int x = 0; x < textureSize; ++x) {
+			if (holeMask.at<uchar>(y, x)) {
+				atlas(y, x) = avgPx;
+			}
+		}
+	}
+
+	// ========== 内置 AO 计算与暗部提亮 ==========
+	// 先释放不再需要的 RC patches 大内存
+	{
+		std::vector<RCPatch>().swap(rcPatches);
+		std::vector<unsigned int>().swap(m_texelPatchID);
+		DEBUG_EXTRA("[AO] Freed RC patches and texelPatchID memory");
+	}
+
+	if (false)
+	{
+		DEBUG_EXTRA("[AO] Loading pre-baked AO from Blender...");
+		
+		std::string aoPath = "/home/algo/Documents/openMVS/data/603454/out.603454/aon.jpg";
+		cv::Mat aoMapSmall = cv::imread(aoPath, cv::IMREAD_GRAYSCALE);
+		if (aoMapSmall.empty()) {
+			DEBUG_EXTRA("[AO] Failed to load AO map from %s, skipping", aoPath.c_str());
+		} else {
+			cv::resize(aoMapSmall, aoMapSmall, cv::Size(1024, 1024), 0, 0, cv::INTER_LINEAR);
+			aoMapSmall.convertTo(aoMapSmall, CV_32FC1, 1.0f / 255.0f);
+			
+			cv::Mat aoMapFull;
+			cv::resize(aoMapSmall, aoMapFull, cv::Size(textureSize, textureSize), 0, 0, cv::INTER_LINEAR);
+			
+			cv::Mat atlasFloat;
+			atlas.convertTo(atlasFloat, CV_32FC3, 1.0f / 255.0f);
+					
+			float strength = 0.38f;
+			float aoDarkThresh = 0.62f;
+			float exponent = 0.92f;
+
+			for (int y = 0; y < textureSize; ++y) {
+				for (int x = 0; x < textureSize; ++x) {
+					float aoVal = aoMapFull.at<float>(y, x);
+					if (aoVal >= aoDarkThresh) continue;
+					
+					cv::Vec3f color = atlasFloat.at<cv::Vec3f>(y, x);
+					
+					float lum = 0.299f * color[2] + 0.587f * color[1] + 0.114f * color[0];
+					if (lum < 0.04f) continue;  // 只保护最暗的材质
+					
+					// ★★★ 线性映射，不加 cap ★★★
+					float factor = (aoDarkThresh - aoVal) / aoDarkThresh;
+					factor = std::pow(factor, exponent);
+					float add = factor * strength;
+					
+					color[0] = std::min(1.0f, color[0] + add);
+					color[1] = std::min(1.0f, color[1] + add);
+					color[2] = std::min(1.0f, color[2] + add);
+					atlasFloat.at<cv::Vec3f>(y, x) = color;
+				}
+			}
+			
+			atlasFloat.convertTo(atlas, CV_8UC3, 255.0f);
+			DEBUG_EXTRA("[AO] AO-based lightening applied using Blender AO map");
+		}
+	}
+
+	return true;
 }
-
+ 
 bool MeshTexture::BuildSeamVertsFromEdges()
 {
     const size_t NP = rcPatches.size();