【译】A GPU Approach to Particle Physics

想到原文中有提到参考的教程，就去看了下，发现对一些逻辑的理解很有帮助，顺便翻译记录一下。

我的 GPGPU 系列的下一个项目是一个粒子物理引擎，它在 GPU 上计算整个物理模拟。粒子受重力影响，会与场景几何体产生反弹。这个 WebGL 演示使用了着色器功能，并不需要严格按照 OpenGL ES 2.0 规范要求，因此它可能在某些平台上无法工作，尤其是在移动设备上。这将在本文后面讨论。

它是可交互的。鼠标光标是一个让粒子反弹的圆形障碍物，单击将在模拟中放置永久性障碍物。你可以绘制粒子可以流过的结构。

这是示例的 HTML5 视频展示，出于必要，它以每秒 60 帧的高比特率录制，所以它相当大。视频编解码器不能很好地处理全屏所有粒子，较低的帧率也不能很好地捕捉效果。我还添加了一些在实际演示中听不到的声音。

• 视频播放地址：https://nullprogram.s3.amazon...

在现代 GPU 上，它可以以每秒 60 帧的速度模拟并且绘制超过 4 百万个粒子。请记住，这是一个 JavaScript 应用程序，我没有真正花时间优化着色器，它受 WebGL 的约束，而不是像 OpenCL 或至少桌面 OpenGL 这样更适合一般计算的东西。

粒子状态编码为颜色

就像 Game of Life 和 path finding 项目一样，模拟状态存储在成对的纹理中，大部分工作是在片元着色器中通过它们之间逐像素映射完成。我不会重复这个设置细节，所以如果你需要了解它是如何工作的，请参考 Game of Life 一文。

对于这个模拟，这些纹理中有四个而不是两个：一对位置纹理和一对速度纹理。为什么是成对的纹理？有 4 个通道，因此其中的每一个部分（x、y、dx、dy）都可以打包到自己的颜色通道中。这似乎是最简单的解决方案。

这个方案的问题是缺乏精确性。对于 R8G8B8A8 内部纹理格式，每个通道为一个字节。总共有 256 个可能的值。显示区域为 800×600 像素，因此显示区域上的每个位置不能都显示。幸运的是，两个字节（总计 65536 个值）对于我们来说已经足够了。

下一个问题是如何跨这两个通道编码值。它需要覆盖负值（负速度），并应尽量充分利用动态范围，比如尝试使用所有 65536 个范围内的值。

要对一个值编码，将该值乘以一个标量，将其扩展到编码的动态范围。选择标量时，所需的最高值（显示的尺寸）是编码的最高值。

接下来，将动态范围的一半添加到缩放值。这会将所有负值转换为正值，0 表示最小值。这种表示法称为 Excess-K 。其缺点是用透明黑色清除纹理（glClearColor）不能将解码值设置为 0 。

最后，将每个通道视为基数为 256 的数字。OpenGL ES 2.0 着色器语言没有按位运算符，因此这是使用普通的除法和模来完成的。我用 JavaScript 和 GLSL 制作了一个编码器和解码器。JavaScript 需要它来写入初始值，并且出于调试目的，它可以读回粒子位置。

vec2 encode(float value) {
    value = value * scale + OFFSET;
    float x = mod(value, BASE);
    float y = floor(value / BASE);
    return vec2(x, y) / BASE;
}

float decode(vec2 channels) {
    return (dot(channels, vec2(BASE, BASE * BASE)) - OFFSET) / scale;
}

JavaScript 与上面的标准化 GLSL 值（0.0-1.0）不同，这会生成一个字节的整数（0-255），用于打包到类型化数组中。

function encode(value, scale) {
    var b = Particles.BASE;
    value = value * scale + b * b / 2;
    var pair = [
        Math.floor((value % b) / b * 255),
        Math.floor(Math.floor(value / b) / b * 255)
    ];
    return pair;
}

function decode(pair, scale) {
    var b = Particles.BASE;
    return (((pair[0] / 255) * b +
             (pair[1] / 255) * b * b) - b * b / 2) / scale;
}

更新每个粒子的片元着色器在该粒子的“索引”处对位置和速度纹理进行采样，解码它们的值，对它们进行操作，然后将它们编码回一种颜色，以便写入输出纹理。因为我使用的是 WebGL ，它缺少多个渲染目标（尽管支持 gl_FragData ），所以片元着色器只能输出一种颜色。位置在一个过程中更新，速度在另一个过程中更新为两个单独的绘图。缓冲区在两个过程完成后才会交换，因此速度着色器（有意）不会使用更新的位置值。

最大纹理大小有一个限制，通常为 8192 或 4096 ，因此纹理不是以一维纹理排列，而是保持方形。粒子由二维坐标索引。

看到直接绘制到屏幕上而不是正常显示的位置或速度纹理非常有趣。这是观看模拟的另一个领域，它甚至帮助我发现了一些其它方面很难看到的问题。输出是一组闪烁的颜色，但有明确的模式，展示了系统的许多状态（或不在其中的状态）。我想分享一段视频，但编码比普通显示更不切实际。以下是截图：位置，然后是速度。这里没有捕捉到阿尔法分量。

在 GPU 上运行这样的模拟最大的挑战之一是缺少随机值。着色器语言中没有 rand() 函数，因此默认情况下整个过程都是确定性的。所有的状态都来自 CPU 填充的初始纹理状态。当粒子聚集并匹配状态时，可能是一起流过一个障碍物，很难将它们重新分离，因为模拟会以相同的方式处理它们。

为了缓解这个问题，第一条规则是尽可能地保持状态。当一个粒子离开显示区域底部时，会将其移回顶部进行“重置”。如果通过将粒子的 Y 值设置为 0 来完成此操作，那么信息将被销毁。这必须避免！尽管在相同的迭代过程中退出，但显示底部边缘以下粒子的 Y 值往往略有不同。代替重置为 0 的做法是添加一个常量：显示区域的高度。Y 值仍然不同，因此这些粒子在碰撞障碍物时更有可能遵循不同的路线。

我使用的下一种技术是通过 uniform 为每次迭代提供一个新随机值，该值被添加到重置粒子的位置和速度。相同的值用于该特定迭代的所有粒子，因此这无助于重叠粒子，但有助于分离“流”。这些都是清晰可见的粒子线，都沿着相同的路径。每一个都会在不同的迭代中退出显示的底部，因此随机值会将它们稍微分开。最终，这会在每次迭代的模拟中加入一些新的状态。

或者，可以向着色器提供包含随机值的纹理。CPU 必须经常填充和上传纹理，另外还有选择纹理采样位置的问题，纹理本身需要一个随机值。

最后，为了处理完全重叠的粒子，在重置时缩放粒子的唯一二维索引，并将其添加到位置和速度中，将它们分开。随机值的符号乘以索引，以避免在任何特定方向上出现偏差。

为了在演示中看到所有这些，制作一个大圆形来捕获所有粒子，让它们流入一个点。这将去除系统中的所有状态。现在清除障碍。它们都会变成一个紧密的团。在顶部重置时，它仍然会有一些聚集，但你会看到它们稍微分开（添加了粒子索引）。它们将在稍有不同的时间离开底部，因此随机值发挥了作用，使它们更加分离。几轮之后，粒子应该会再次均匀分布。

状态的最后一个来源是你的鼠标。在场景中移动它时，会干扰粒子，并向模拟引入一些噪声。

纹理作为顶点属性缓冲区

在阅读 OpenGL ES 着色器语言规范（PDF）时，我产生了这个项目的想法。我一直想做一个粒子系统，但我卡在如何绘制粒子的问题上。表示位置的纹理数据需要以某种方式作为顶点反馈到管道中。通常，缓冲区纹理——由数组缓冲区支持的纹理——或像素缓冲区对象——异步纹理数据复制——可用于此操作，但 WebGL 没有这些功能。从 GPU 中提取纹理数据，并将其作为每帧上的数组缓冲区重新加载是不可能的。

然而，我想出了一个很酷的技巧，比这两个都好。着色器函数 texture2D 用于对纹理中的像素进行采样。通常情况下，片元着色器将其用作计算一个像素颜色过程的一部分。但是着色器语言规范提到，texture2D 也可以在顶点着色器中使用。就在那时，一个点子击中了我。顶点着色器本身可以执行从纹理到顶点的转换。

它的工作原理是将前面提到的二维粒子索引作为顶点属性传递，使用它们从顶点着色器中查找粒子位置。着色器将以 GL_POINTS 模式运行，发射点粒子。这是简略的版本：

attribute vec2 index;

uniform sampler2D positions;
uniform vec2 statesize;
uniform vec2 worldsize;
uniform float size;

// float decode(vec2) { ...

void main() {
    vec4 psample = texture2D(positions, index / statesize);
    vec2 p = vec2(decode(psample.rg), decode(psample.ba));
    gl_Position = vec4(p / worldsize * 2.0 - 1.0, 0, 1);
    gl_PointSize = size;
}

真实版本也会对速度进行采样，因为它会调节颜色（缓慢移动的粒子比快速移动的粒子更亮）。

然而，有一个潜在的问题：允许实现将顶点着色器纹理绑定的数量限制为 0（GL_MAX_vertex_texture_IMAGE_UNITS）。所以从技术上讲，顶点着色器必须始终支持 texture2D ，但它们不需要支持实际的纹理。这有点像飞机上不载客的餐饮服务。有些平台不支持这种技术。到目前为止，我只在一些移动设备上遇到过这个问题。

除了缺乏一些平台的支持之外，这允许模拟的每个部分都留在 GPU 上，并为纯 GPU 粒子系统铺平道路。

一个重要的观察结果是粒子之间不相互作用。这不是一个 n 体模拟。然而，它们确实与世界其它地方互动：它们直观地从这些静止的圆圈上反弹。该环境由另一个纹理表示，该纹理在正常迭代期间不会更新。我称之为障碍物纹理。

障碍物纹理上的颜色是曲面法线。也就是说，每个像素都有一个指向它的方向，一个将粒子导向某个方向的流。空隙有一个特殊的常值（0，0）。这不是单位向量（长度不为1），因此它是对粒子没有影响的带外值。

粒子只需对障碍物纹理进行采样即可检查碰撞。如果在其位置找到法线，则使用着色器函数 reflect 更改其速度。该函数通常用于在 3D 场景中反射光线，但对于缓慢移动的粒子同样适用。其效果是粒子以自然的方式从圆上反弹。

有时粒子以低速或零速度落在障碍物上或落在障碍物中。为了把它们从障碍物上移开，会将它们朝着正常的方向轻轻推一推。你会在斜坡上看到这一点，在那里，缓慢的粒子像跳跃的豆子一样摇动着向下自由移动。

为了使障碍物纹理对用户友好，实际的几何图形在 JavaScript 的 CPU 端进行维护。这些圆会保留在一个列表中，并在更新时，这个列表中重新绘制障碍物纹理。例如，每次在屏幕上移动鼠标时，都会出现这种情况，从而产生移动障碍物。纹理提供了对几何体的着色器友好访问。两种表现有两个目的。

当我开始编写程序的这一部分时，我设想除了圆以外其它可以放置的形状。例如，实心矩形：法线看起来像这样。

到目前为止，这些尚未得到实施。

我还没试过，但我想知道粒子是否也可以通过将自己绘制到障碍物纹理上来相互作用。附近的两个粒子会相互反弹。也许整个 liquid demo 可以像这样在 GPU 上运行。如果我猜想正确，粒子会增大体积，形成碗状的障碍物会填满，而不是将粒子集中到一个点上。

我认为这个项目还有一些需要探索的地方。

• A GPU Approach to Particle Physics