移动立方体算法 —— 三维等值面提取的经典方法以及其扩展和改进

移动立方体算法（Marching Cubes Algorithm）是由 Lorensen 和 Cline 在1987年提出的经典算法，用于从三维标量场数据中提取等值面，构建三角网格模型。该算法在医学影像（CT、MRI）、科学可视化、地质建模等领域有着广泛应用。其核心思想是将三维空间划分为立方体网格，通过检查每个立方体的顶点标量值，确定等值面与立方体的交线，进而生成三角面片。

一、算法基本概念

1.1 问题描述

给定一个三维标量场函数 $f(x, y, z)$ 和阈值 $\tau$，我们需要提取满足 $f(x, y, z) = \tau$ 的所有点的集合，即等值面（isosurface），并用三角网格表示。

在实际应用中，$f(x, y, z)$ 通常以离散采样点的形式给出，即在一个规则的三维网格上，每个网格点 $(i, j, k)$ 对应一个标量值 $V_{i,j,k}$。

1.2 基本思想

移动立方体算法的基本思想是：

空间划分：将三维空间划分为规则的立方体网格
顶点分类：对每个立方体的8个顶点，根据其标量值与阈值的关系，标记为”内部”（小于阈值）或”外部”（大于等于阈值）
查找表：根据8个顶点的标记状态，查找对应的三角化模式
插值计算：在等值面与立方体边的交点处进行线性插值，确定顶点位置
三角化：根据查找表的模式，生成三角面片

二、算法流程详解

2.1 立方体的结构

对于每个立方体(体素)，都共有八个顶点(体元)。立方体的体素和体元

2.2 立方体顶点分类

对于每个立方体，有8个顶点，每个顶点有两种状态：

内部点：标量值 $< \tau$（$记为0）
外部点：标量值 $\geq \tau$（记为1）

因此，8个顶点共有 $2^8 = 256$ 种可能的组合状态。由于对称性，实际只需要考虑15种不同的情况。

移动立方体算法15种基本模式

2.3 等值面与立方体的交点

等值面与立方体的交点必然位于立方体的边上。对于每条边，如果两个端点的标量值一个小于阈值、一个大于等于阈值，则等值面与该边有且仅有一个交点。

1.线性插值：交点位置可以通过线性插值计算：

\[P = P_1 + \frac{\tau - V_1}{V_2 - V_1}(P_2 - P_1)\]

其中：

$P_1, P_2$ 是边的两个端点
$V_1, V_2$ 是对应的标量值
$\tau$ 是阈值
$P$ 是交点位置

2.光线穿透：通过光线穿透穿透原始物体表面找到真是体素对应体元间的交点。

移动立方体算法流程

graph TD
    A["开始：输入标量场数据<br/>网格分辨率 (nx, ny, nz)<br/>阈值 τ"] --> B["初始化顶点列表<br/>和三角面片列表"]
    B --> C["遍历所有立方体<br/>for i=0 to nx-2<br/>for j=0 to ny-2<br/>for k=0 to nz-2"]
    C --> D["获取当前立方体<br/>8个顶点的标量值"]
    D --> E["对8个顶点分类<br/>V < τ: 标记为0(内部)<br/>V ≥ τ: 标记为1(外部)"]
    E --> F["计算立方体索引<br/>index = Σ(vertex[i] << i)"]
    F --> G{"索引是否为0或255?<br/>(全内部或全外部)"}
    G -->|是| H["跳过该立方体<br/>无等值面穿越"]
    G -->|否| I["查找边表<br/>edges = edgeTable[index]"]
    I --> J["计算每条边与<br/>等值面的交点"]
    J --> K["根据查找表获取<br/>三角化模式"]
    K --> L["生成三角面片<br/>添加到面片列表"]
    L --> M{"是否遍历完<br/>所有立方体?"}
    M -->|否| C
    M -->|是| N["输出三角网格<br/>顶点列表 + 面片列表"]
    N --> O["结束"]
    
    style A fill:#e1f5fe
    style O fill:#c8e6c9
    style G fill:#fff3e0
    style H fill:#ffcdd2
    style I fill:#dcedc8

算法特点：

🔵 蓝色区域：算法初始化和输入数据
🟡 黄色区域：决策节点（判断是否需要处理）
🔴 红色区域：跳过无交点的立方体
🟢 绿色区域：核心处理流程（交点计算和三角化）

立方体索引计算

graph LR
    subgraph "立方体顶点编号"
        A["v0: (0,0,0)<br/>位0"] 
        B["v1: (1,0,0)<br/>位1"]
        C["v2: (1,1,0)<br/>位2"]
        D["v3: (0,1,0)<br/>位3"]
        E["v4: (0,0,1)<br/>位4"]
        F["v5: (1,0,1)<br/>位5"]
        G["v6: (1,1,1)<br/>位6"]
        H["v7: (0,1,1)<br/>位7"]
    end
    
    I["顶点分类<br/>V < τ: 0<br/>V ≥ τ: 1"] --> J["计算索引<br/>index = v0×1 + v1×2 +<br/>v2×4 + v3×8 +<br/>v4×16 + v5×32 +<br/>v6×64 + v7×128"]
    J --> K["查找边表<br/>edgeTable[index]"]
    K --> L["查找三角表<br/>triTable[index]"]
    L --> M["生成三角面片"]
    
    style A fill:#e1f5fe
    style B fill:#e1f5fe
    style C fill:#e1f5fe
    style D fill:#e1f5fe
    style E fill:#e1f5fe
    style F fill:#e1f5fe
    style G fill:#e1f5fe
    style H fill:#e1f5fe
    style M fill:#c8e6c9

关键理解：

顶点编号：每个立方体的8个顶点按特定顺序编号（0-7）
位编码：每个顶点用一个位表示（0或1），8个顶点组成一个8位整数
索引范围：0-255，对应256种可能的顶点状态组合
对称性：利用互补对称性（如果所有顶点状态取反，等值面不变），可将256种情况减少到15种基本模式

以二维举例描述编码规则： MC二维Square编码规则

三、查找表构造

3.1 边表（Edge Table）

边表用于快速确定哪些边与等值面相交。对于立方体的12条边，使用12位整数表示：

\[edgeMask = \sum_{i=0}^{11} edge_i \times 2^i\]

其中 $edge_i = 1$ 表示第 $i$ 条边与等值面相交，否则为0。

3.2 三角表（Triangle Table）

三角表存储每个立方体索引对应的三角化模式。每个模式由一系列顶点索引组成，每3个连续索引构成一个三角形。

3.3 15种基本模式

由于立方体的对称性和互补性，256种情况可以归结为15种基本模式：

情况0：全内部或全外部，无等值面
情况1-5：单顶点在外或在内
情况6-10：两个相对顶点在外或在内
情况11-14：更复杂的相交模式

四、算法实现细节

4.1 顶点插值计算

对于每条与等值面相交的边，需要计算交点的精确位置：

// 伪代码：线性插值计算交点
Point3D interpolate(Point3D p1, Point3D p2, float v1, float v2, float threshold) {
    if (abs(v1 - v2) < epsilon) return p1;  // 避免除零
    
    float t = (threshold - v1) / (v2 - v1);
    return p1 + t * (p2 - p1);
}

数学表达式：

\[t = \frac{\tau - V_1}{V_2 - V_1}\] \[P = (1-t) \cdot P_1 + t \cdot P_2\]

4.2 法向量计算

为了渲染高质量的三维模型，需要计算每个顶点的法向量。法向量可以通过梯度估计得到：

\[\nabla f = \left(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}, \frac{\partial f}{\partial z}\right)\]

在离散网格中，可以使用中心差分法近似：

\[g_x = \frac{V_{i+1,j,k} - V_{i-1,j,k}}{2h}\] \[g_y = \frac{V_{i,j+1,k} - V_{i,j-1,k}}{2h}\] \[g_z = \frac{V_{i,j,k+1} - V_{i,j,k-1}}{2h}\]

其中 $h$ 是网格间距。归一化后的法向量为：

\[\mathbf{n} = \frac{\nabla f}{|\nabla f|}\]

五、应用示例

示例1：简单球体提取

问题描述： 从一个距离场中提取半径为 $R$ 的球面。

解决方案：

定义标量场： $f(x, y, z) = \sqrt{x^2 + y^2 + z^2}$
设置阈值： $\tau = R$
应用算法：
- 对于每个立方体，计算8个顶点到原点的距离
- 根据距离与 $R$ 的关系标记顶点
- 生成等值面三角网格

结果： 得到一个由三角面片组成的球面网格。

六、算法优化与改进

6.1 计算复杂度

时间复杂度： $O(n^3)$，其中 $n$ 是每个维度的网格分辨率
空间复杂度： $O(n^3)$ 用于存储标量场数据

6.2 优化策略

空间分割： 使用八叉树等数据结构，只处理与等值面相交的区域
并行计算： 每个立方体可以独立处理，适合并行化
自适应细化： 在等值面附近使用更细的网格，远离等值面的区域使用粗网格
缓存优化： 重用相邻立方体共享的顶点和边

6.3 歧义性问题

原始移动立方体算法在某些配置下存在歧义性，可能导致网格出现空洞。解决方案包括：

双线性插值方法： 使用更复杂的插值函数
歧义性检测： 检测并特殊处理歧义情况
改进算法： 如Marching Tetrahedra等变体

6.4 多属性mc算法

不同的体元节点有不同的属性，在体元节点判断内外的时候的时候只跟自己同属性的判断。因此带标签的体元输入时，结果可能为多个体。简单描述就是不同体带不同属性进入MC构造一次。但是由于多层进入MC初始化加重新遍历，n个标签就要操作n次，时间复杂度高。因此，根据论文方案采取了利用细分最小基本体素单元，更换体素构型的方案(体素重新划分8个块根据8个角点的不同属性构造)。

多次构造：
MC(sdf_tooth) → Mesh_tooth
MC(sdf_bone) → Mesh_bone
MC(sdf_gum) → Mesh_gum

单次构造：
MCTT(sdf_bone_tooth_gum) -> Mesh_tooth,Mesh_bone,Mesh_gum

对具有两个以上面中心的立方体中的线段进行三角剖分。

6.5 跨体素算法

（算法是一种将不同分辨率体素数据生成的相邻三角网格无缝拼接在一起的方法，从而使细节层次（LOD）能够应用于大型体素数据集）

Transvoxel算法的工作原理是在以一种分辨率采样的体素数据和以该分辨率一半采样的体素数据之间的边界上，在常规单元之间插入特殊的过渡单元。我们并非考虑过渡边界处全分辨率和半分辨率数据的所有可能的体素状态组合（这将需要处理大约120万个案例），而是仅考虑高分辨率数据的9个样本。这使得我们需要处理的案例数量大大减少，仅为512个，而这些案例恰好属于图2所示的73个等价类。每个过渡单元都填充有这些三角形图案之一，以便完美填充不同分辨率网格之间出现的接缝、裂缝和孔洞，正如图1中洞穴的处理方式一样。

📄 跨体素算法的网站

Transvoxel算法效果

相关参考资料：
📄 Octree-Based Decimation of Marching Cubes Surface - 基于多材质MC行进立方体算法

6.5 OCtree加速构建算法

传统 MC 的问题

分辨率固定 → 细节要靠整体加密
内存爆炸（$N^3$）
量空区域也在跑 Marching Cube
牙齿 / 骨骼这种 局部高曲率、局部细节 非常浪费

Octree + MC 的优势”哪里复杂就细，哪里简单就粗” Octree 的思想其实非常直观： 1.把这件事和 MC 结合起来，等价于： 2.不在整个体素网格上跑 MC，只在“真正需要的地方”跑

1. Octree 的细分判定（核心逻辑）

对于一个 Octree cell，我只做了一件事：

在 8 个角点采样 SDF

如果全部 > 0 或全部 < 0
→ 这个 cell 完全在内 / 外，不可能有表面
如果正负混合
→ 等值面穿过该 cell

是否继续细分，取决于：

cell 的尺寸是否大于最小阈值
SDF / 法向变化是否足够大
否则，该 cell 就作为 叶节点，交给 MC 处理。

2. 在 Octree 叶节点上跑 MC

这一阶段非常“朴素”：

每个叶节点就是一个小立方体
使用标准 MC 的 edge table / tri table
在这个 cell 内生成三角形

也就是说：

Octree 负责“剪枝”，MC 负责“成型”。

相关参考资料：
📄 Octree-Based Decimation of Marching Cubes Surface - 基于八叉树的移动立方体表面简化方法

具体实现和代码： 📄 参考meshlab中对possion重建中的MC算法优化和实现

七、算法步骤总结以及代码

移动立方体算法的完整步骤：

初始化： 输入标量场数据 $V_{i,j,k}$，设置阈值 $\tau$，初始化顶点列表和面片列表
遍历立方体： 对每个立方体 $(i, j, k)$，其中 $i \in [0, nx-2]$，$j \in [0, ny-2]$，$k \in [0, nz-2]$
顶点分类： 获取8个顶点的标量值，根据与阈值的关系标记为0（内部）或1（外部）
计算索引： $index = \sum_{i=0}^{7} vertex_i \times 2^i$
判断： 如果 $index = 0$ 或 $index = 255$，跳过该立方体（无等值面穿越）
查找边表： 根据 $index$ 查找 $edgeTable[index]$，确定哪些边与等值面相交
计算交点： 对每条相交的边，使用线性插值计算交点位置： $P = P_1 + \frac{\tau - V_1}{V_2 - V_1}(P_2 - P_1)$
三角化： 根据 $triTable[index]$ 获取三角化模式，生成三角面片
计算法向量： （可选）使用梯度估计计算顶点法向量
输出结果： 返回顶点列表和三角面片索引列表