第1章：3D Mesh基础

      v2
      /|\
     / | \
    /  |  \
   v3--v--v1
    \  |  /
     \ | /
      \|/
      v4

struct HalfEdge {
    Vertex* origin;      // 起始顶点
    HalfEdge* twin;      // 孪生半边（反向）
    HalfEdge* next;      // 同一面内的下一条半边
    HalfEdge* prev;      // 同一面内的上一条半边
    Face* face;          // 所属的面
};

struct Vertex {
    Vector3 position;    // 空间坐标
    HalfEdge* halfedge;  // 任意一条以此为起点的半边
};

struct Face {
    HalfEdge* halfedge;  // 面的任意一条半边
};

        v2
        /|\
       / | \
     h1  |  h0
     /   |   \
    v0---h2---v1

def vertex_one_ring(v):
    neighbors = []
    h = v.halfedge
    h_start = h

    do:
        neighbors.append(h.twin.origin)
        h = h.twin.next
    while h != h_start

    return neighbors

def vertex_valence(v):
    count = 0
    h = v.halfedge
    h_start = h

    do:
        count += 1
        h = h.twin.next
    while h != h_start

    return count

# 注释
v  x y z [w]          # 顶点坐标，w默认为1.0
vt u v [w]            # 纹理坐标，w默认为0
vn x y z              # 顶点法线（已归一化）
f  v/vt/vn ...        # 面定义（顶点/纹理/法线索引）

g group_name          # 组定义
o object_name         # 对象定义
mtllib material.mtl   # 材质库引用
usemtl material_name  # 使用材质
s smoothing_group     # 平滑组（0或off表示不平滑）

newmtl material_name
Ka r g b              # 环境光反射
Kd r g b              # 漫反射
Ks r g b              # 镜面反射
Ns exponent           # 镜面高光指数
d alpha               # 透明度（或Tr）
illum model           # 光照模型
map_Kd texture.jpg    # 漫反射贴图

ply
format ascii 1.0      # 或 binary_little_endian / binary_big_endian
comment 扫描数据来源
element vertex 8
property float x
property float y
property float z
property uchar red    # 顶点颜色
property uchar green
property uchar blue
element face 6
property list uchar int vertex_indices
end_header
[数据部分]

property float nx     # 法线
property float ny
property float nz
property float confidence  # 扫描置信度
property float intensity   # 激光雷达强度
property float curvature   # 局部曲率

solid name
  facet normal nx ny nz
    outer loop
      vertex x1 y1 z1
      vertex x2 y2 z2
      vertex x3 y3 z3
    endloop
  endfacet
  ...
endsolid name

[80字节头部]
[4字节无符号整数 - 三角形数量]
[每个三角形50字节：
  12字节 - 法线向量 (3个float)
  36字节 - 三个顶点 (9个float)
  2字节  - 属性字节计数（通常为0）]

# 转换管线示例
Source → Import → Normalize → Process → Export → Target
         ↓         ↓           ↓         ↓
      验证拓扑  统一坐标系  优化网格  验证输出

    非流形边（T型）        非流形顶点（沙漏）
        /|\                    /|\
       / | \                  / | \
      /  |  \                /  |  \
     /___|___\              /   *   \
         |                 /___|_____\

         |                 /___|_____\
         |                     |

def is_manifold(mesh):
    # 1. 检查边的面数
    edge_face_count = {}
    for face in mesh.faces:
        for edge in face.edges():
            key = tuple(sorted(edge))
            edge_face_count[key] = edge_face_count.get(key, 0) + 1

    # 非流形边检测
    for edge, count in edge_face_count.items():
        if count > 2:  # T型边
            return False

    # 2. 检查顶点的扇形结构
    for vertex in mesh.vertices:
        if not is_fan_structure(vertex):
            return False

    return True

def verify_topology(mesh, expected_genus):
    chi = mesh.num_vertices - mesh.num_edges + mesh.num_faces
    if mesh.is_closed:
        computed_genus = (2 - chi) // 2
        return computed_genus == expected_genus
    return False

def classify_mesh(mesh):
    chi = compute_euler_characteristic(mesh)
    if chi == 2:
        return "球面拓扑"
    elif chi == 0:
        return "环面拓扑"
    elif chi < 0:
        g = (2 - chi) // 2
        return f"{g}-亏格曲面"

def reduce_genus(mesh, target_genus):
    current_genus = compute_genus(mesh)

    while current_genus > target_genus:
        # 找到一个非分离循环
        loop = find_handle_loop(mesh)

        # 切开并填充
        mesh = cut_and_fill(mesh, loop)

        current_genus -= 1

    return mesh

def find_boundaries(mesh):
    edge_count = {}

    # 统计每条边的面数
    for face in mesh.faces:
        for edge in face.edges():
            key = tuple(sorted(edge))
            edge_count[key] = edge_count.get(key, 0) + 1

    # 找出边界边
    boundary_edges = [e for e, c in edge_count.items() if c == 1]

    # 组织成边界圈
    boundary_loops = organize_into_loops(boundary_edges)

    return boundary_loops

def smooth_boundary(mesh, iterations=10):
    boundary_vertices = find_boundary_vertices(mesh)

    for _ in range(iterations):
        new_positions = {}

        for v in boundary_vertices:
            # 只在边界上移动
            neighbors = get_boundary_neighbors(v)
            if len(neighbors) == 2:
                # 线性平滑
                new_pos = (neighbors[0].pos + neighbors[1].pos) / 2
                # 投影到切平面
                new_positions[v] = project_to_tangent(v, new_pos)

        # 更新位置
        for v, pos in new_positions.items():
            v.position = pos

def triangle_quality(v0, v1, v2):
    # 计算边长
    e0 = np.linalg.norm(v1 - v0)
    e1 = np.linalg.norm(v2 - v1)
    e2 = np.linalg.norm(v0 - v2)

    # 计算面积
    area = 0.5 * np.linalg.norm(np.cross(v1 - v0, v2 - v0))

    # 计算圈度
    perimeter = e0 + e1 + e2
    roundness = 4 * np.pi * area / (perimeter ** 2)

    # 计算最小角
    angles = compute_angles(e0, e1, e2)
    min_angle = np.min(angles)

    # 综合评分 (0-1)
    quality = roundness * (min_angle / 60.0)

    return quality

def valence_distribution(mesh):
    valences = {}
    for v in mesh.vertices:
        val = v.valence()
        valences[val] = valences.get(val, 0) + 1

    # 计算度数6的顶点比例
    regular_ratio = valences.get(6, 0) / len(mesh.vertices)

    return valences, regular_ratio

def compute_curvature_tensor(vertex):
    # 构造协方差矩阵
    C = np.zeros((3, 3))

    for edge in vertex.edges():
        e = edge.vector()
        kappa = edge_curvature(edge)
        C += kappa * np.outer(e, e)

    # 特征值分解
    eigenvalues, eigenvectors = np.linalg.eigh(C)

    # 主曲率和主方向
    k1, k2 = eigenvalues[2], eigenvalues[1]
    d1, d2 = eigenvectors[:, 2], eigenvectors[:, 1]

    return k1, k2, d1, d2

class TopologyChecker:
    def __init__(self, mesh):
        self.mesh = mesh
        self.issues = []

    def check_manifold(self):
        """\u68c0\u67e5\u975e\u6d41\u5f62\u7ed3\u6784"""
        edge_faces = {}

        for face in self.mesh.faces:
            for edge in face.edges():
                key = tuple(sorted(edge))
                edge_faces.setdefault(key, []).append(face)

        for edge, faces in edge_faces.items():
            if len(faces) > 2:
                self.issues.append({
                    'type': 'non-manifold-edge',
                    'edge': edge,
                    'faces': faces
                })

    def check_intersections(self):
        """\u68c0\u67e5\u81ea\u76f8\u4ea4"""
        bvh = BVH(self.mesh)  # \u4f7f\u7528BVH\u52a0\u901f

        for i, face1 in enumerate(self.mesh.faces):
            candidates = bvh.query(face1.bbox())

            for j in candidates:
                if j <= i:
                    continue

                face2 = self.mesh.faces[j]
                if triangles_intersect(face1, face2):
                    self.issues.append({
                        'type': 'face-intersection',
                        'faces': (i, j)
                    })

    def check_degeneracies(self):
        """\u68c0\u67e5\u9000\u5316\u5143\u7d20"""
        # \u96f6\u9762\u79ef\u4e09\u89d2\u5f62
        for i, face in enumerate(self.mesh.faces):
            area = face.area()
            if area < 1e-10:
                self.issues.append({
                    'type': 'zero-area-face',
                    'face': i,
                    'area': area
                })

        # \u91cd\u590d\u9876\u70b9
        vertex_map = {}
        for i, v in enumerate(self.mesh.vertices):
            key = tuple(np.round(v.position, 6))
            if key in vertex_map:
                self.issues.append({
                    'type': 'duplicate-vertex',
                    'vertices': (vertex_map[key], i)
                })
            else:
                vertex_map[key] = i

    def check_orientation(self):
        """\u68c0\u67e5\u65b9\u5411\u4e00\u81f4\u6027"""
        visited = set()

        def propagate_orientation(face, expected_orientation):
            if face in visited:
                return

            visited.add(face)

            for neighbor in face.neighbors():
                shared_edge = face.shared_edge(neighbor)

                # \u68c0\u67e5\u5171\u4eab\u8fb9\u7684\u65b9\u5411
                face_dir = face.edge_direction(shared_edge)
                neighbor_dir = neighbor.edge_direction(shared_edge)

                if face_dir == neighbor_dir:
                    self.issues.append({
                        'type': 'inconsistent-orientation',
                        'faces': (face, neighbor)
                    })

                propagate_orientation(neighbor, -expected_orientation)

        # \u4ece\u7b2c\u4e00\u4e2a\u9762\u5f00\u59cb\u4f20\u64ad
        if self.mesh.faces:
            propagate_orientation(self.mesh.faces[0], 1)

def quad_advancing_front(boundary):
    front = boundary.copy()
    quads = []

    while len(front) > 4:
        # 找最适合形成四边形的位置
        best_score = float('inf')
        best_quad = None

        for i in range(len(front)):
            # 尝试用连续4个顶点形成四边形
            v0, v1, v2, v3 = front[i:i+4]
            score = evaluate_quad_quality(v0, v1, v2, v3)

            if score < best_score:
                best_score = score
                best_quad = (v0, v1, v2, v3)
                best_idx = i

        # 添加四边形并更新前沿
        quads.append(best_quad)
        front = front[:best_idx] + front[best_idx+4:]

    return quads

struct PolygonFace {
    std::vector<Vertex*> vertices;
    Vector3 normal;

    bool isConvex() const {
        // 检查所有内角是否小于180度
        for (int i = 0; i < vertices.size(); i++) {
            Vector3 v0 = vertices[i]->position;
            Vector3 v1 = vertices[(i+1)%n]->position;
            Vector3 v2 = vertices[(i+2)%n]->position;

            if (cross(v1-v0, v2-v1).dot(normal) < 0)
                return false;
        }
        return true;
    }
};

class HybridMesh:
    def __init__(self):
        self.triangles = []
        self.quads = []
        self.polygons = []

    def add_face(self, vertices):
        n = len(vertices)
        if n == 3:
            self.triangles.append(vertices)
        elif n == 4:
            self.quads.append(vertices)
        else:
            self.polygons.append(vertices)

    def to_triangle_mesh(self):
        """转换为纯三角形网格"""
        tri_mesh = []

        # 三角形直接添加
        tri_mesh.extend(self.triangles)

        # 四边形分割为两个三角形
        for quad in self.quads:
            # 选择更好的对角线
            diag1_len = distance(quad[0], quad[2])
            diag2_len = distance(quad[1], quad[3])

            if diag1_len < diag2_len:
                tri_mesh.append([quad[0], quad[1], quad[2]])
                tri_mesh.append([quad[0], quad[2], quad[3]])
            else:
                tri_mesh.append([quad[0], quad[1], quad[3]])
                tri_mesh.append([quad[1], quad[2], quad[3]])

        # 多边形三角化
        for poly in self.polygons:
            tri_mesh.extend(triangulate_polygon(poly))

        return tri_mesh

标准T-junction：        虚拟T-junction：        复杂T-junction：
    A                       A                       A───B
    |                       |                         \ |
    |                       |                          \|
B───*───C               B───+···C                   C───*───D
    |                   (虚拟顶点)                      |

    |                   (虚拟顶点)                      |
    |                                                   |

    D                                                   E

def detect_t_junctions(mesh, tolerance=1e-6):
    t_junctions = []

    # 构建边的空间索引
    edge_tree = KDTree()
    for edge in mesh.edges:
        edge_tree.insert(edge)

    # 检查每个顶点
    for vertex in mesh.vertices:
        # 查找附近的边
        nearby_edges = edge_tree.query_radius(vertex.position, tolerance)

        for edge in nearby_edges:
            # 跳过包含该顶点的边
            if vertex in edge.vertices:
                continue

            # 计算点到线段的距离
            t = point_to_segment_parameter(vertex.position, edge)

            if 0 < t < 1:  # 在边内部
                dist = point_to_segment_distance(vertex.position, edge)

                if dist < tolerance:
                    t_junctions.append({
                        'vertex': vertex,
                        'edge': edge,
                        'parameter': t,
                        'distance': dist
                    })

    return t_junctions

def fix_by_vertex_insertion(mesh, t_junction):
    vertex = t_junction['vertex']
    edge = t_junction['edge']
    t = t_junction['parameter']

    # 在边上插入新顶点
    new_vertex = mesh.split_edge(edge, t)

    # 合并重合顶点
    mesh.merge_vertices(vertex, new_vertex, tolerance)

    # 更新相关面的拓扑
    affected_faces = get_affected_faces(edge)
    for face in affected_faces:
        face.update_topology()

def fix_by_edge_split(mesh, t_junctions):
    # 按边分组T-junctions
    edge_splits = {}
    for tj in t_junctions:
        edge = tj['edge']
        if edge not in edge_splits:
            edge_splits[edge] = []
        edge_splits[edge].append(tj['parameter'])

    # 处理每条边
    for edge, parameters in edge_splits.items():
        # 排序分割点
        parameters.sort()

        # 从后向前分割（保持索引有效）
        for t in reversed(parameters):
            mesh.split_edge_at_parameter(edge, t)

def fix_by_retriangulation(mesh, t_junction):
    # 收集影响区域的面
    affected_faces = []
    for face in mesh.faces:
        if t_junction['vertex'] in face.vertices:
            affected_faces.append(face)
        elif any(v in t_junction['edge'].vertices for v in face.vertices):
            affected_faces.append(face)

    # 提取边界
    boundary = extract_boundary(affected_faces)

    # 删除原始面
    for face in affected_faces:
        mesh.remove_face(face)

    # 重新三角化
    new_faces = triangulate_polygon(boundary)
    for face in new_faces:
        mesh.add_face(face)

class QuadtreeMesh:
    def __init__(self, domain, max_level):
        self.root = QuadNode(domain, level=0)
        self.max_level = max_level

    def adapt(self, criterion_func):
        """根据准则函数自适应细化"""
        queue = [self.root]

        while queue:
            node = queue.pop(0)

            if node.level >= self.max_level:
                continue

            # 评估是否需要细化
            if criterion_func(node):
                # 细分为四个子节点
                node.subdivide()
                queue.extend(node.children)

        # 确保2:1平衡
        self.balance()

    def balance(self):
        """确保相邻单元级别差不超过1"""
        changed = True
        while changed:
            changed = False

            for node in self.get_all_nodes():
                if node.is_leaf():
                    # 检查邻居
                    for neighbor in node.get_neighbors():
                        if neighbor.level > node.level + 1:
                            node.subdivide()
                            changed = True
                            break

def error_based_refinement(mesh, target_error):
    while True:
        # 计算每个元素的误差
        errors = compute_element_errors(mesh)

        # 全局误差
        global_error = np.sqrt(np.sum(errors**2))

        if global_error < target_error:
            break

        # 标记需要细化的元素
        threshold = np.percentile(errors, 80)
        marked = errors > threshold

        # 细化标记的元素
        for i, elem in enumerate(mesh.elements):
            if marked[i]:
                mesh.refine_element(elem)

class RedGreenRefinement:
    def __init__(self, mesh):
        self.mesh = mesh

    def red_refinement(self, triangle):
        """红色细化：将三角形分为4个相似子三角形"""
        v0, v1, v2 = triangle.vertices

        # 创建边中点
        m01 = self.mesh.add_vertex((v0.pos + v1.pos) / 2)
        m12 = self.mesh.add_vertex((v1.pos + v2.pos) / 2)
        m20 = self.mesh.add_vertex((v2.pos + v0.pos) / 2)

        # 创建4个新三角形
        self.mesh.add_triangle(v0, m01, m20)
        self.mesh.add_triangle(m01, v1, m12)
        self.mesh.add_triangle(m20, m12, v2)
        self.mesh.add_triangle(m01, m12, m20)

        # 删除原三角形
        self.mesh.remove_triangle(triangle)

    def green_refinement(self, triangle, edge):
        """绿色细化：二分三角形以消除悬挂节点"""
        v0, v1, v2 = triangle.vertices

        # 找到需要细分的边
        if edge == (v0, v1):
            m = self.mesh.add_vertex((v0.pos + v1.pos) / 2)
            self.mesh.add_triangle(v0, m, v2)
            self.mesh.add_triangle(m, v1, v2)
        # ... 其他边类似

        self.mesh.remove_triangle(triangle)

def parallel_adaptive_refinement(mesh, num_threads):
    """并行自适应细化"""
    import threading
    from queue import Queue

    # 1. 标记阶段（并行）
    marked_elements = parallel_mark_elements(mesh, num_threads)

    # 2. 着色分组（避免冲突）
    color_groups = mesh_coloring(marked_elements)

    # 3. 并行细化每个颜色组
    for color in color_groups:
        threads = []
        for i in range(num_threads):
            # 分配元素子集
            subset = color[i::num_threads]
            t = threading.Thread(
                target=refine_elements,
                args=(mesh, subset)
            )
            threads.append(t)
            t.start()

        # 等待当前颜色组完成
        for t in threads:
            t.join()

    # 4. 修复和平衡
    mesh.fix_hanging_nodes()
    mesh.balance_2to1()

# OBJ文件读取时的常见错误
vertex_index = int(token) - 1  # 必须-1！
# 防御性编程
assert 0 <= vertex_index < len(vertices)

def check_normal_consistency(mesh):
    for edge in mesh.interior_edges():
        face1, face2 = edge.adjacent_faces()
        # 共享边在两个面中方向应该相反
        if edge.direction_in(face1) == edge.direction_in(face2):
            return False
    return True

def vertices_equal(v1, v2, epsilon=1e-6):
    return np.linalg.norm(v1 - v2) < epsilon

# 使用空间哈希加速
def spatial_hash(vertex, cell_size):
    return tuple(np.floor(vertex / cell_size).astype(int))

def is_degenerate(v0, v1, v2, area_threshold=1e-10):
    area = 0.5 * np.linalg.norm(np.cross(v1-v0, v2-v0))
    if area < area_threshold:
        return True

    # 检查极细长三角形
    edges = [np.linalg.norm(v1-v0), 
             np.linalg.norm(v2-v1), 
             np.linalg.norm(v0-v2)]
    if max(edges) / min(edges) > 1000:
        return True

    return False

def parse_obj_face(line):
    vertices = []
    for token in line.split()[1:]:
        parts = token.split('/')
        v_idx = int(parts[0]) - 1
        vt_idx = int(parts[1]) - 1 if len(parts) > 1 and parts[1] else None
        vn_idx = int(parts[2]) - 1 if len(parts) > 2 and parts[2] else None
        vertices.append((v_idx, vt_idx, vn_idx))
    return vertices