Cloth Simulation — 布料模拟

一件衬衫被风吹起的样子，一块丝绸落在桌面上的纹理，游戏里的旗帜随风飘动——这些看似简单的视觉效果，背后隐藏着一个有趣的物理问题：如何用计算机模拟柔软的织物？

今天实现了经典的质点弹簧布料模拟，这是游戏引擎和影视特效中处理布料最常见的基础方案。

最终效果

主渲染图（含风力的最终状态）

布料模拟最终状态

红色布料从水平展开的初始状态，在重力拉扯下垂落，包裹住下方的球体，随后被侧风吹起一角。Phong 着色区分正面（深红）和背面（暗红），形成明显的布料折叠感。

物理演化序列（0 → 150 → 300 → 600 帧）

布料垂落序列

从左上到右下：布料从展开状态，逐渐在重力作用下垂落，绕过球体弯曲，最终稳定在平衡态。

第一章：为什么布料模拟难？

布料的本质是大量纤维交织的网。真正从纤维层面模拟是不现实的——一件 T 恤大概有几十亿根纤维。

工程上的解法是离散化：把布料简化成有限个质点，再用弹簧连接它们，让弹簧产生的弹力近似真实布料的弹性。这就是 Mass-Spring System（质点弹簧系统）。

直觉模型：想象一个渔网。每个网结是一个质点（有质量，能运动），绳子是弹簧（被拉伸时产生回弹力）。把渔网竖起来，重力拉着每个结往下走，弹簧弹力阻止它们走得太远，最终网会挂成一个自然的弧形。

这套方案的优点是简单、可控；缺点是参数多（弹簧刚度、阻尼、步长），调不好会”数值爆炸”（弹簧力越来越大，最终粒子飞出去）。

第二章：质点网格与三类弹簧

布料被离散为 25×25 的质点网格（625个质点），质点之间通过三类弹簧连接。

为什么需要三类弹簧？

只用一类弹簧（比如只连接横纵邻居）的话，布料会出现一个严重问题：过度剪切——布料沿对角线方向几乎没有阻力，会像橡皮泥一样被拉成菱形，完全失去织物的感觉。

弹簧类型	连接方式	刚度	物理含义
结构弹簧	横向/纵向邻居（间距1）	800 N/m	抵抗拉伸，保持整体形状
剪切弹簧	对角线邻居（间距√2）	400 N/m	抵抗剪切变形，防止”菱形塌陷”
弯曲弹簧	跨一个质点（间距2）	200 N/m	抵抗弯曲，模拟布料的刚度

// 结构弹簧：横向和纵向直接相邻的质点
for (int r = 0; r < rows; r++) {
    for (int c = 0; c < cols; c++) {
        int idx = r * cols + c;
        if (c + 1 < cols) addSpring(idx, idx + 1,     struct_k);  // 横向
        if (r + 1 < rows) addSpring(idx, idx + cols,  struct_k);  // 纵向
    }
}

// 剪切弹簧：对角线方向
for (int r = 0; r < rows - 1; r++) {
    for (int c = 0; c < cols - 1; c++) {
        int idx = r * cols + c;
        addSpring(idx,     idx + cols + 1, shear_k); // 右下
        addSpring(idx + 1, idx + cols,     shear_k); // 左下
    }
}

// 弯曲弹簧：跨一个质点，提供抗弯刚度
for (int r = 0; r < rows; r++) {
    for (int c = 0; c < cols - 2; c++) {
        int idx = r * cols + c;
        addSpring(idx, idx + 2,        bend_k);  // 横向跨一格
    }
}
for (int r = 0; r < rows - 2; r++) {
    for (int c = 0; c < cols; c++) {
        int idx = r * cols + c;
        addSpring(idx, idx + cols * 2, bend_k);  // 纵向跨一格
    }
}

刚度比例 4:2:1 的直觉：结构弹簧最硬（布料不容易被拉长），剪切次之（可以有一定形变），弯曲最软（布料可以弯折）。这个比例是经典配置，来自实践经验。

第三章：弹簧力计算——胡克定律

每根弹簧按照胡克定律产生力：拉伸越多，力越大；压缩越多，反方向力越大。

$$F = k \cdot (|r| - L_0) \cdot \hat{r}$$

其中：

$k$ 是弹簧刚度（越大越硬）
$|r| - L_0$ 是形变量（当前长度 - 原始静止长度）
$\hat{r}$ 是从 A 指向 B 的单位方向向量

void step(double dt, bool apply_wind) {
    // 清空所有力
    for (auto& p : particles) p.force = {0, 0, 0};

    // 1. 重力：每个自由质点都受 mg 向下
    for (auto& p : particles) {
        if (!p.pinned)
            p.force += {0, gravity * p.mass, 0};  // gravity = -9.8
    }

    // 2. 风力：沿 x/z 方向施加小扰动
    if (apply_wind) {
        for (auto& p : particles)
            if (!p.pinned) p.force += wind * 0.1;
    }

    // 3. 弹簧力：遍历所有弹簧，对两端质点施加等大反向的力
    for (auto& s : springs) {
        Vec3 diff = particles[s.b].pos - particles[s.a].pos;
        double len = diff.length();
        if (len < 1e-12) continue;

        double extension = len - s.rest_length;  // 正数=拉伸，负数=压缩
        Vec3 force_dir   = diff.normalized();    // 从 a 指向 b
        Vec3 f = force_dir * (s.stiffness * extension);

        // 牛顿第三定律：a 受 +f，b 受 -f
        if (!particles[s.a].pinned) particles[s.a].force += f;
        if (!particles[s.b].pinned) particles[s.b].force -= f;
    }
    // ...（后续积分）
}

一个关键细节：顶边一行的质点标记为 pinned=true，模拟被夹住悬挂的效果。固定点不参与力的计算和位置更新，只是充当”锚点”。

第四章：Verlet 积分——稳定的时间推进

有了合力，下一步是更新每个质点的速度和位置。最直觉的方法是欧拉积分：

$$\mathbf{v}_{t+dt} = \mathbf{v}t + \mathbf{a} \cdot dt$$
$$\mathbf{x}{t+dt} = \mathbf{x}t + \mathbf{v}{t+dt} \cdot dt$$

但欧拉积分在布料模拟中有个大问题：能量不守恒。它会不断地向系统注入微小的能量，导致弹簧越震越猛，最终爆炸。

Verlet 积分的解法很巧妙——不显式存储速度，而是从当前位置和上一帧位置的差推算速度：

$$\mathbf{v}_t \approx \frac{\mathbf{x}t - \mathbf{x}{t-dt}}{dt}$$

$$\mathbf{x}_{t+dt} = \mathbf{x}_t + (\mathbf{x}t - \mathbf{x}{t-dt}) \cdot \text{damping} + \mathbf{a} \cdot dt^2$$

// Verlet 积分
for (auto& p : particles) {
    if (p.pinned) continue;
    Vec3 acceleration = p.force / p.mass;

    // vel ≈ (pos - prev_pos)，再乘 damping 耗散能量
    Vec3 vel = (p.pos - p.prev_pos) * damping;  // damping = 0.99

    Vec3 new_pos = p.pos + vel + acceleration * (dt * dt);
    p.prev_pos = p.pos;   // 保存当前帧位置作为"上一帧"
    p.pos = new_pos;
}

阻尼系数 0.99 的作用：每帧让速度乘以 0.99，相当于有微小的空气阻力。没有阻尼，布料会永远震荡不停；阻尼太大（比如 0.9）则布料运动过于迟缓。0.99 是让布料”自然减速”的经验值。

Verlet 的优点：

数值更稳定（时间复杂度 O(N)）
不需要单独存储速度向量（省内存）
阻尼直接乘在速度项上，实现简洁

第五章：球体碰撞检测

布料落下时会遇到球体，需要防止质点穿入球体内部。

碰撞检测：法线推出法（Normal Pushout）

原理非常简单：检查质点是否在球体内部，如果是，就沿着球心到质点的方向，把质点推到球面外：

// 球体碰撞处理
for (auto& p : particles) {
    if (p.pinned) continue;

    Vec3 diff = p.pos - sphere_center;      // 球心到质点的向量
    double dist = diff.length();             // 当前距离

    // epsilon = 0.02，额外留一点间隙，避免抖动
    if (dist < sphere_radius + 0.02) {
        // 将质点推到球面上
        p.pos = sphere_center + diff.normalized() * (sphere_radius + 0.02);
        // 注意：不修改 prev_pos，这样隐式速度会指向外侧（粒子"贴着"球面滑动）
    }
}

为什么不修改 prev_pos？ 因为 Verlet 积分中，速度由 pos - prev_pos 决定。如果同时修改 prev_pos，质点的隐式速度会变成零（完全粘住）。只修改 pos 的话，质点在球面上有切向速度，可以自然地沿球面滑动——这正是我们想要的效果。

第六章：渲染——从物理到像素

自制软光栅化器

本项目没有使用 OpenGL，完全用 CPU 自制光栅化：

投影：将 3D 坐标通过透视投影变换到屏幕坐标：

std::array<double, 3> project(const Vec3& p) const {
    Vec3 forward = (target - eye).normalized();
    Vec3 right   = forward.cross(up).normalized();
    Vec3 up_vec  = right.cross(forward);

    Vec3 dir      = p - eye;
    double d_fwd  = dir.dot(forward);   // 深度（距相机距离）
    double d_right= dir.dot(right);
    double d_up   = dir.dot(up_vec);

    // 透视除法：近处大，远处小
    double tan_fov = std::tan(fov * 0.5 * M_PI / 180.0);
    double sx = d_right / (d_fwd * tan_fov * aspect);
    double sy = d_up    / (d_fwd * tan_fov);

    double px = (sx + 1.0) * 0.5 * w;     // [-1,1] → [0, width]
    double py = (1.0 - (sy + 1.0) * 0.5) * h;  // Y 轴翻转
    return {px, py, d_fwd};
}

Z-Buffer：记录每个像素已渲染的最近深度，只更新更近的片元：

bool test(int x, int y, double z) {
    if (z < buf[y*w + x]) {
        buf[y*w + x] = z;  // 更新深度
        return true;        // 允许渲染
    }
    return false;           // 被遮挡，跳过
}

双面渲染 + Phong 着色

布料是双面的——正面和背面颜色不同。通过法线与视线方向的点积来判断：

Vec3 n = edge1.cross(edge2).normalized();  // 面法线
Vec3 view = eye - face_center;              // 视线方向

// 点积 > 0 表示法线朝向相机（正面），否则是背面
Vec3 color = (n.dot(view) > 0) ? cloth_color_front : cloth_color_back;

Phong 光照模型：环境光 + 漫反射 + 高光：

Vec3 phong(const Vec3& pos, const Vec3& normal, ...) {
    Vec3 N = normal.normalized();
    Vec3 L = (light_pos - pos).normalized();  // 光线方向
    Vec3 V = (eye - pos).normalized();         // 视线方向
    Vec3 R = 2.0 * N.dot(L) * N - L;         // 反射方向

    double ambient  = 0.15;
    double diffuse  = std::max(0.0, N.dot(L));              // 朗伯余弦律
    double specular = std::pow(std::max(0.0, R.dot(V)), 32); // 高光指数 32

    return cloth_color * (ambient + diffuse * 0.75)
         + light_color * specular * 0.3;
}

ACES 色调映射

最后通过 ACES 曲线把 HDR 线性颜色压缩到 [0,1] 范围，避免过曝：

Vec3 aces(Vec3 c) {
    // ACES Film Tone Mapping Curve (近似)
    // 阴影细节保留 + 高光柔和过渡 + 中性色偏
    double a=2.51, b=0.03, d=2.43, e=0.59, f=0.14;
    return { (c.x*(a*c.x+b))/(c.x*(d*c.x+e)+f), ... };
}

第七章：物理参数与稳定性

关键参数

const int    ROWS = 25, COLS = 25;   // 质点网格大小
const double CLOTH_SIZE = 2.4;        // 布料物理尺寸（单位：米）
const double DT   = 0.008;            // 时间步长（每步8毫秒）
const double struct_k = 800.0;        // 结构弹簧刚度
const double shear_k  = 400.0;        // 剪切弹簧刚度
const double bend_k   = 200.0;        // 弯曲弹簧刚度
const double damping  = 0.99;         // Verlet 阻尼系数

数值稳定性：弹簧刚度与步长的约束

一个经验法则：要保持 Verlet 积分稳定，需要满足：

$$k \cdot dt^2 < m$$

代入我们的参数：$800 \times 0.008^2 = 0.0512 < 1.0$（每个质点质量为 1.0 kg），满足条件。

如果把 struct_k 提高到 8000，或者把 dt 增大到 0.08，系统就会不稳定——弹簧力在每步内推动质点超过弹簧原长，然后下一步又更大的力推回来，振幅越来越大，最终”爆炸”。

模拟阶段

阶段1（600步，无风）：布料在重力下自然垂落，包裹球体
阶段2（200步，有风）：侧风让布料一角飘起
阶段3（100步，无风）：稳定至平衡态
总物理时间：900 × 0.008s = 7.2秒

量化验证结果

粒子最低点  y = -2.500（到达地面限制 ✅）
粒子最高点  y = -0.203（顶边固定点附近 ✅）
粒子平均 y  = -1.744（布料充分垂落 ✅，初始 y=0）

布料像素覆盖率：7.2%（34,330 / 480,000像素）
布料平均颜色：RGB(151, 51, 25) — 红色调正确 ✅
运行时间：0.077秒 ✅（全部CPU，无GPU加速）

第八章：进阶方向

本实现的局限

无摩擦约束：质点在球面滑动时没有摩擦力，布料会滑落
无撕裂：弹簧永远不会断裂（真实布料会撕破）
自碰撞缺失：布料折叠时不会检测布料与自身的碰撞，会穿插
Verlet 精度：对于大变形或快速运动，精度不如 Runge-Kutta 4（RK4）

工业级方案

PBD（Position-Based Dynamics）：Nvidia FleX 和 Unreal Engine 的布料引擎都用 PBD，通过位置约束而非力约束来保证稳定性，允许更大的时间步长
FEM（有限元法）：精度最高，用于影视特效（如 Pixar 的布料模拟），但计算量大
自碰撞：通常用 BVH 加速碰撞对检测，复杂度从 O(N²) 降至 O(N log N)