Marching Cubes：3D 形状是怎么从数学公式里"长出来"的？（C++ 从零实现）

前言：一个奇怪的问题

如果有人给你一个数学公式，比如：

$$f(x, y, z) = \sqrt{x^2 + y^2 + z^2} - 1$$

你能在 3D 空间里把它”画”出来吗？

答案是可以的。$f(x,y,z) = 0$ 的所有点，恰好就是一个半径为 1 的球面。这种”由方程定义的曲面”叫做等值面（Isosurface）。

问题是：计算机显卡只会画三角形，不会直接画”方程”。我们需要一个算法，把等值面转换成三角形。

这就是 Marching Cubes（行进立方体）的工作。

今天我们从零实现这个算法，得到下面这张图：

这是 5 个球用”平滑融合”连接在一起的有机形状，由 15,196 个三角形组成，生成耗时 0.078 秒。

---

一、等值面是什么？用 SDF 定义三维形状

隐式表示 vs 显式表示

3D 形状有两种定义方式：

显式：直接列出表面上的点（三角网格、点云）。
隐式：用一个函数 $f(x,y,z)$ 描述空间中每个点和形状的关系。

隐式表示里最常用的是 SDF（Signed Distance Field，有符号距离场）：

• $f(x,y,z) < 0$：点在形状内部
• $f(x,y,z) = 0$：点在形状表面（这就是等值面）
• $f(x,y,z) > 0$：点在形状外部

球的 SDF 最简单——就是点到球心的距离减去半径：

float sdSphere(Vec3 p, Vec3 center, float radius) {
    // 点 p 到球心 center 的距离，减去 radius
    // < 0：在球内部
    // = 0：在球面上
    // > 0：在球外部
    float dx = p.x - center.x;
    float dy = p.y - center.y;
    float dz = p.z - center.z;
    return sqrtf(dx*dx + dy*dy + dz*dz) - radius;
}

直觉：SDF 告诉你”你离最近的表面有多远，以及你在内部还是外部”。

平滑并集（Smooth Union）

多个球可以融合在一起，形成有机形状。普通的并集（取最小值）会产生硬边缘：

// 普通并集（硬边缘）
float d = min(sdSphere(p, center1, r1), sdSphere(p, center2, r2));

Smooth Union（平滑并集） 在交界处增加平滑过渡：

float smoothMin(float a, float b, float k) {
    // h 计算两个距离有多接近（相差在 k 范围内）
    float h = std::max(k - std::abs(a - b), 0.0f) / k;
    // 返回插值结果，在交界处产生平滑过渡
    return std::min(a, b) - h * h * k * 0.25f;
    //              ↑普通最小值  ↑平滑修正项（越接近交界越大）
}

参数 k 控制过渡区域的大小：

• k = 0：等同于硬并集，没有过渡
• k = 0.5：0.5 单位范围内平滑融合（本项目使用的值）
• k = 2.0：很大的融合区域，形状变得很”黏”

把 5 个球融合在一起：

float scalarField(Vec3 p) {
    float d = sdSphere(p, {0, 0, 0},      1.2f);  // 中心大球（半径1.2）
    d = smoothMin(d, sdSphere(p, { 1.4f,  0.5f, 0},    0.7f), 0.5f);  // 右卫星
    d = smoothMin(d, sdSphere(p, {-1.4f,  0.5f, 0},    0.7f), 0.5f);  // 左卫星
    d = smoothMin(d, sdSphere(p, {0,     -1.4f, 0.5f}, 0.6f), 0.5f);  // 下卫星
    d = smoothMin(d, sdSphere(p, {0,      0.5f,-1.4f}, 0.65f),0.5f);  // 后卫星
    return d;
}

这个函数就是我们的”标量场”——Marching Cubes 接下来要把它变成三角形。

---

二、Marching Cubes 核心算法

大思路：把空间切成小方块

Marching Cubes 的核心想法：把 3D 空间划分成很多小立方体（体素），逐个处理每个体素。

对于每个体素，看它的 8 个角点，判断每个角点是在等值面内部还是外部：

体素的 8 个角：

      7──────6
     /|      /|
    / |     / |    z
   4──────5  |    ↑
   |  3── | ─2    |
   | /    | /     +──→ x
   |/     |/     /
   0──────1    y↗

对每个角点：

• 如果 scalarField(角点) < 0（在等值面内部）→ 标记为 1
• 如果 scalarField(角点) >= 0（在等值面外部）→ 标记为 0

8 个角点，每个 0 或 1，总共有 $2^8 = 256$ 种可能的组合。

关键洞察：256 种情况

每种 256 种组合，等值面穿过这个体素的方式是不同的。但有规律可循：

比如：

• 全是 0（都在外部）→ 等值面不经过这个体素
• 全是 1（都在内部）→ 等值面不经过这个体素
• 角 0 是 1，其余是 0 → 等值面在角 0 附近切一个小三角形

Lorensen & Cline（1987）的论文预先计算出了所有 256 种情况的三角形方案，存成一个查找表：

// triTable[256][16]
// 每行最多描述 5 个三角形（15 条边索引 + -1 结束符）
// 边索引 0-11 对应体素的 12 条边
static const int triTable[256][16] = {
    {-1,...},         // case 0: 全外部，无三角形
    {0, 8, 3, -1,...},// case 1: 角0在内部，1个三角形
    {0, 1, 9, -1,...},// case 2: 角1在内部
    // ... 共256行
};

每行的数字是”边索引”——三角形的顶点在哪条边上（等值面和体素的边的交点）。

代码实现：逐体素处理

for (int iz = 0; iz < gridRes-1; iz++) {
for (int iy = 0; iy < gridRes-1; iy++) {
for (int ix = 0; ix < gridRes-1; ix++) {

    // 步骤1：计算体素8个角的标量值
    float values[8];
    Vec3  corners[8];
    for (int k = 0; k < 8; k++) {
        // cornerOffset[k] 给出第 k 个角的偏移量 (0或1)
        int cx = ix + cornerOffset[k][0];
        int cy = iy + cornerOffset[k][1];
        int cz = iz + cornerOffset[k][2];
        corners[k] = gridOrigin + Vec3(cx, cy, cz) * voxelSize;
        values[k] = grid[cx][cy][cz];  // 预先采样好的标量值
    }

    // 步骤2：计算这个体素的"案例索引"
    int cubeIndex = 0;
    for (int k = 0; k < 8; k++) {
        if (values[k] < 0.0f) {
            cubeIndex |= (1 << k);  // 第k位设为1
        }
    }
    // cubeIndex 现在是 0-255 之间的整数

    // 步骤3：查表得到三角形
    const int* tris = triTable[cubeIndex];
    for (int t = 0; tris[t] != -1; t += 3) {
        // tris[t], tris[t+1], tris[t+2] 是三角形的三条边索引
        Triangle tri;
        for (int v = 0; v < 3; v++) {
            int edgeIdx = tris[t + v];
            // 这条边连接哪两个角？
            int v0 = edgeVertices[edgeIdx][0];
            int v1 = edgeVertices[edgeIdx][1];
            // 线性插值找到等值面和这条边的交点
            tri.v[v] = interpolate(corners[v0], corners[v1],
                                   values[v0], values[v1]);
        }
        triangles.push_back(tri);
    }
}}}

精确找交点：线性插值

等值面（value = 0）和体素某条边的交点，通过线性插值求得：

角点A（value=-0.8，在内部）     角点B（value=0.4，在外部）
         A●─────────────────────●B
              ↑等值面在这里

交点在 A 向 B 走 t 的位置，其中：

$$t = \frac{0 - v_A}{v_B - v_A} = \frac{-(-0.8)}{0.4 - (-0.8)} = \frac{0.8}{1.2} \approx 0.667$$

Vec3 interpolate(Vec3 p0, Vec3 p1, float v0, float v1) {
    // 找到 value=0 的交点
    if (std::abs(v1 - v0) < 1e-6f) return (p0 + p1) * 0.5f;

    // t = (目标值 - v0) / (v1 - v0) = (0 - v0) / (v1 - v0)
    float t = -v0 / (v1 - v0);

    // 线性插值位置
    return {
        p0.x + t * (p1.x - p0.x),
        p0.y + t * (p1.y - p0.y),
        p0.z + t * (p1.z - p0.z)
    };
}

这确保了三角形的顶点精确地落在等值面上（标量值 = 0 的地方），而不是随意放在体素边的中点。

---

三、计算法线：数值梯度

有了三角形顶点，还需要法线才能做光照计算。

法线垂直于等值面，恰好就是标量场的梯度方向：

$$\vec{n} = \nabla f = \left(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}, \frac{\partial f}{\partial z}\right)$$

SDF 的梯度等于”指向外部的方向”——这正是法线！

在代码里用中心差分数值近似梯度：

Vec3 computeNormal(Vec3 p) {
    const float eps = 0.001f;  // 微小扰动（差分步长）

    // 在 x 方向各偏移一点点，用差值近似偏导数
    float dx = scalarField({p.x + eps, p.y,       p.z})
             - scalarField({p.x - eps, p.y,       p.z});
    float dy = scalarField({p.x,       p.y + eps, p.z})
             - scalarField({p.x,       p.y - eps, p.z});
    float dz = scalarField({p.x,       p.y,       p.z + eps})
             - scalarField({p.x,       p.y,       p.z - eps});

    return Vec3(dx, dy, dz).normalized();
    //                        ↑ 归一化为单位向量
}

为什么用中心差分而不是前向差分？

中心差分（两边各取一点）精度更高（二阶精度），前向差分（只往一边走）精度低（一阶精度）。每次多调用一次 scalarField，但法线更准确。

在 Marching Cubes 里，法线计算这样集成：

// 先计算每个角点的法线（在网格采样阶段）
Vec3 cornerNormals[8];
for (int k = 0; k < 8; k++) {
    cornerNormals[k] = computeNormal(corners[k]);
}

// 三角形顶点的法线：对两个角点的法线做同样的插值
tri.n[v] = interpolate_normal(cornerNormals[v0], cornerNormals[v1], t);
tri.n[v] = tri.n[v].normalized();

这样每个三角形顶点有自己的法线，Phong 着色时可以在三角形面内插值法线，得到平滑的光照效果（而不是每个三角形面都是同一个颜色）。

---

四、正交投影和光栅化

有了三角形网格，接下来把它渲染成 2D 图像。

正交投影

我们用正交投影（也叫平行投影），而不是透视投影。区别：

• 透视投影：近大远小，模拟真实摄像机，适合游戏画面
• 正交投影：物体大小不随距离变化，适合技术展示、工程图纸

透视投影：                正交投影：
    ↓ 近                     ↓ 所有光线平行
   大物体                    相同大小
    ↓ 远
   小物体

代码：

bool project(Vec3 p, float& sx, float& sy, float& depth) {
    // 把世界坐标转换到相机坐标系
    Vec3 rel = p - camPos;
    float x = rel.dot(camRight);   // 相机右方向分量
    float y = rel.dot(camUp);      // 相机上方向分量
    float z = rel.dot(camFwd);     // 相机前方向分量（深度）

    // 正交投影：直接用 x/y 映射到屏幕（不除以 z！）
    sx = (x / orthoWidth + 0.5f) * W;   // 映射到 [0, W]
    sy = (0.5f - y / orthoHeight) * H;  // 注意 y 轴翻转（屏幕 y 向下）
    depth = z;  // 深度用于 Z-Buffer
    return true;
}

Z-Buffer：解决遮挡

多个三角形重叠时，近的应该遮住远的。Z-Buffer 为每个像素记录”当前最近的三角形的深度”：

// 每个像素一个深度值，初始化为无穷远
std::vector<float> zbuffer(W * H, 1e30f);

// 光栅化时，对每个像素检查深度
float z = w0*depth[0] + w1*depth[1] + w2*depth[2];  // 插值深度
int pidx = y * W + x;

if (z >= zbuffer[pidx]) continue;  // 已经有更近的三角形，跳过
zbuffer[pidx] = z;                  // 更新深度记录
// 计算颜色并写入图像

重心坐标：判断像素是否在三角形内

对屏幕上的每个像素，需要判断它是否在当前三角形内。用重心坐标：

对于点 $P$，如果 $P = w_0 \cdot A + w_1 \cdot B + w_2 \cdot C$，其中 $w_0 + w_1 + w_2 = 1$，那么：

• $w_0, w_1, w_2$ 全都 $\geq 0$：点在三角形内
• 有任何一个 $< 0$：点在三角形外

重心坐标还有另一个用途：插值三角形内部的属性（法线、颜色等）。

// 2D 重心坐标计算
float area = (sx[1]-sx[0])*(sy[2]-sy[0]) - (sx[2]-sx[0])*(sy[1]-sy[0]);

float w0 = ((sx[1]-sx[2])*(py-sy[2]) + (sy[2]-sy[1])*(px-sx[2])) / area;
float w1 = ((sx[2]-sx[0])*(py-sy[0]) + (sy[0]-sy[2])*(px-sx[0])) / area;
float w2 = 1.0f - w0 - w1;

if (w0 < 0 || w1 < 0 || w2 < 0) continue;  // 在三角形外，跳过

// 利用重心坐标插值法线
Vec3 N = tri.n[0]*w0 + tri.n[1]*w1 + tri.n[2]*w2;
N = N.normalized();

---

五、Phong 光照

有了每个像素的法线，就可以计算光照颜色。用经典的 Phong 光照模型：

$$\text{颜色} = \text{环境光} + \text{漫反射} + \text{高光}$$

环境光（Ambient）：均匀的基础亮度，避免完全黑暗的暗面
漫反射（Diffuse）：取决于法线和光线的角度（正对光源最亮）
高光（Specular）：镜面反射，产生亮点

Vec3 shade(Vec3 pos, Vec3 normal, Vec3 baseColor) {
    Vec3 ambient = {0.08f, 0.08f, 0.12f};  // 深蓝色环境光（让暗部不死黑）

    auto computeLight = [&](Vec3 lightPos, Vec3 lightColor, float intensity) {
        Vec3 L = (lightPos - pos).normalized();  // 指向光源的方向
        Vec3 V = (camPos - pos).normalized();    // 指向摄像机的方向
        Vec3 H = (L + V).normalized();           // 半角向量（L和V的中间方向）

        // 漫反射：法线和光线方向夹角的余弦值
        // 正对光源（夹角=0°）cos=1 最亮，侧面（90°）cos=0 无光，背面截断为0
        float diff = max(0.0f, normal.dot(L));

        // 高光：法线和半角向量夹角的高次幂
        // 指数越高，高光点越小越亮（越"金属感"）
        float spec = pow(max(0.0f, normal.dot(H)), 64.0f);

        return (baseColor * diff + lightColor * spec * 0.6f) * intensity;
    };

    // 主光源：暖白色，右上前方
    Vec3 light1 = computeLight({3.0f, 5.0f, 4.0f}, {1.0f, 0.95f, 0.9f}, 1.0f);
    // 补光：冷蓝色，左后方（增加立体感）
    Vec3 light2 = computeLight({-4.0f, 2.0f, -2.0f}, {0.5f, 0.6f, 0.8f}, 0.4f);

    return ambient + light1 + light2;  // 三者叠加
}

为什么用半角向量 H 而不是反射向量 R？
传统 Phong 用反射向量，Blinn-Phong 改用半角向量。两者效果相似，但半角向量计算更简单（不需要 reflect 函数），而且在大角度时有更好的表现。

---

六、完整渲染流程

把上面所有模块串在一起：

1. 构建标量场
   → scalarField(p) = smoothUnion(5个球的SDF)

2. 网格采样（64×64×64 = 262,144 个体素）
   → 每个角点采样 scalarField 和法线

3. Marching Cubes
   → 逐体素 → 查 triTable → 线性插值交点 → 生成三角形

4. 软件渲染
   → 正交投影 → 重心坐标光栅化 → Z-Buffer → Phong 着色

5. 输出 PNG
   → 自制 zlib + PNG 编码（无第三方依赖）

运行结果：

• 三角形数：15,196
• 网格分辨率：64×64×64
• 总耗时：0.078 秒

---

七、自制 PNG 输出

这个项目没有用任何第三方库，连 PNG 输出都是自己写的。简单讲讲原理：

PNG 文件格式 = 文件头 + IHDR块（图像信息）+ IDAT块（压缩像素数据）+ IEND块

像素数据用 zlib/deflate 压缩。我们用最简单的”存储块”（不压缩，只打包）：

// deflate 存储块格式（BTYPE=00，不压缩）
out.push_back(is_last ? 0x01 : 0x00);   // BFINAL | BTYPE
out.push_back(block_len & 0xFF);          // LEN 低字节
out.push_back(block_len >> 8);            // LEN 高字节
out.push_back(~block_len & 0xFF);         // NLEN（LEN取反，校验用）
out.push_back(~block_len >> 8);
// 紧接着是原始数据

还需要 CRC32（块校验）和 Adler-32（整体校验），代码里都有完整实现。

---

八、参数和调优

参数	当前值	调大效果	调小效果
网格分辨率	64	更多三角形，更平滑	更快但有锯齿
smoothMin k	0.5	融合过渡更大	趋近硬并集
等值面阈值	0.0	形状缩小	形状变大
中心差分 eps	0.001	法线误差大	太小会数值不稳定

---

九、完整代码

编译运行：

g++ -O2 -std=c++17 main.cpp -o marching_cubes
./marching_cubes
# 输出: marching_cubes_output.png, marching_cubes_side.png

• GitHub：03-07-marching-cubes

---

总结

Marching Cubes 把一个”数学函数定义的形状”变成三角网格，核心步骤：

1. 划分体素：把空间切成小方块
2. 判断角点内外：8 个角各是 0 或 1
3. 查表：8 位索引 → 256 种情况 → 预定义的三角形方案
4. 插值交点：等值面在边上的精确位置
5. 计算法线：数值梯度，中心差分

配合 Smooth Union SDF，可以生成有机融合的形状——游戏里的地形、角色体积建模、医学影像的 CT 重建，都是 Marching Cubes 的实际应用场景。

---

参考资料

• Lorensen & Cline (1987) - Marching Cubes: A High Resolution 3D Surface Construction Algorithm
• Paul Bourke - Polygonising a Scalar Field
• Inigo Quilez - Smooth Minimum