使用计算着色器（Compute Shader）模拟粒子效果【OpenGL】【GLSL】

个人感觉计算着色器很像 CUDA，都是利用显卡的强大计算能力来加速，只不过 CUDA 仅适用于 N 卡，而计算着色器具有跨平台的能力(Shader Model 5.0以上才支持)

效果如图：







关键代码及注释如下：

C++ 代码
void initialize()
{
// 计算着色器
GLuint compute_shader = buildShader(
&compute_shader_source,
GL_COMPUTE_SHADER,
"Error in compiling the compute shader\n");

compute_prog = buildComputeProg(compute_shader);

setupData();
// ---------------------------------------------
// 渲染着色器
GLuint vs = buildShader(
&render_vs,
GL_VERTEX_SHADER,
"Error in compiling the vertex shader\n");

GLuint fs = buildShader(
&render_fs,
GL_FRAGMENT_SHADER,
"Error in compiling the fragment shader\n");

render_prog = buildRenderProg(vs, fs);

// ---------------------------------------------
}

C++代码：
// 激活顶点/纹理属性
void setupData()
{
dt_location = glGetUniformLocation(compute_prog, "dt");

// 创建 VAO
glGenVertexArrays(1, &render_vao);
glBindVertexArray(render_vao);

// 创建 VBO
glGenBuffers(2, buffers);

//glEnableVertexAttribArray(1); // 启用索引为 1 的顶点属性——粒子位置
//glVertexAttribPointer(0, 4, GL_FLOAT, GL_FALSE, 0, NULL); // 表明 buffer 内数据的格式

glBindBuffer(GL_ARRAY_BUFFER, position_buffer);
glBufferData(GL_ARRAY_BUFFER, PARTICLE_COUNT * sizeof(glm::vec4), NULL, GL_DYNAMIC_COPY); // 由于数据是不断变化的，所以需要动态拷贝

glm::vec4 * positions = (glm::vec4 *)glMapBufferRange(GL_ARRAY_BUFFER,
0,
PARTICLE_COUNT * sizeof(glm::vec4),
GL_MAP_WRITE_BIT | GL_MAP_INVALIDATE_BUFFER_BIT);

for (int i = 0; i < PARTICLE_COUNT; i++)
{
positions[i] = glm::vec4(randomVector(-10.0f, 10.0f), randomFloat());
}

glUnmapBuffer(GL_ARRAY_BUFFER);

glEnableVertexAttribArray(0); // 启用索引为 0 的顶点属性——粒子速度
glVertexAttribPointer(0, 4, GL_FLOAT, GL_FALSE, 0, NULL); // 表明 buffer 内数据的格式

glBindBuffer(GL_ARRAY_BUFFER, velocity_buffer); // 表明使用的 buffer 数据来源
glBufferData(GL_ARRAY_BUFFER, PARTICLE_COUNT * sizeof(glm::vec4), NULL, GL_DYNAMIC_COPY); // 由于数据是不断变化的，所以需要动态拷贝

glm::vec4 * velocities = (glm::vec4 *)glMapBufferRange(GL_ARRAY_BUFFER,
0,
PARTICLE_COUNT * sizeof(glm::vec4),
GL_MAP_WRITE_BIT | GL_MAP_INVALIDATE_BUFFER_BIT);

for (int i = 0; i < PARTICLE_COUNT; i++)
{
velocities[i] = glm::vec4(randomVector(-0.1f, 0.1f), 0.0f);
}

glUnmapBuffer(GL_ARRAY_BUFFER);

// ---------------------
// 创建 TBO
glGenTextures(2, tbos);

for (int i = 0; i < 2; i++)
{
glBindTexture(GL_TEXTURE_BUFFER, tbos[i]);
glTexBuffer(GL_TEXTURE_BUFFER, GL_RGBA32F, buffers[i]); // ☆ VBO 将从 TBO 中获取数据，注意：buffers 是 position_buffer/velocity_buffer 的别名
}

// ---------------------
// 创建 UBO
glGenBuffers(1, &attractor_buffer);
glBindBuffer(GL_UNIFORM_BUFFER, attractor_buffer);
glBufferData(GL_UNIFORM_BUFFER, 32 * sizeof(glm::vec4), NULL, GL_STATIC_DRAW); // 仅初始化，静态拷贝

// attractor 的最后一位，表示粒子质量
for (int i = 0; i < MAX_ATTRACTORS; i++)
{
attractor_masses[i] = 0.5f + randomFloat() * 0.5f;
}

glBindBufferBase(GL_UNIFORM_BUFFER, 0, attractor_buffer);

// --------------------------------------
}

C++ 代码：
void display()
{
static const GLuint start_ticks = ::GetTickCount() - 100000;
GLuint current_ticks = ::GetTickCount();
static GLuint last_ticks = current_ticks;
float time = ((start_ticks - current_ticks) & 0xFFFFF) / float(0xFFFFF);
float delta_time = (float)(current_ticks - last_ticks) * 0.075f;

// ----------------------------------------------------
// 映射 UBO
glm::vec4 * attractors = (glm::vec4 *)glMapBufferRange(GL_UNIFORM_BUFFER,
0,
32 * sizeof(glm::vec4),
GL_MAP_WRITE_BIT | GL_MAP_INVALIDATE_BUFFER_BIT);

int i;

for (i = 0; i < 32; i++)
{
attractors[i] = glm::vec4(
sinf(time * (float)(i + 4) * 7.5f * 20.0f) * 50.0f,
cosf(time * (float)(i + 7) * 3.9f * 20.0f) * 50.0f,
sinf(time * (float)(i + 3) * 5.3f * 20.0f) * cosf(time * (float)(i + 5) * 9.1f) * 100.0f,
attractor_masses[i]);
//std::cout << "attractors = " << attractors[i][0] << std::endl;
}

glUnmapBuffer(GL_UNIFORM_BUFFER);

// ----------------------------------------------------
// 激活 计算着色器，并绑定到 TBO（存储粒子的位置和速度）
glUseProgram(compute_prog);
glBindImageTexture(0, velocity_tbo, 0, GL_FALSE, 0, GL_READ_WRITE, GL_RGBA32F);// 既可读又可写，注意可以省去 glActiveTexture 的调用
glBindImageTexture(1, position_tbo, 0, GL_FALSE, 0, GL_READ_WRITE, GL_RGBA32F);
// 为 计算着色器传参——delta time
// If dt is too large, the system could explode, so cap it to
// some maximum allowed value
if (delta_time >= 2.0f)
{
delta_time = 2.0f;
}
glUniform1f(dt_location, delta_time);
// 分发计算任务 num_groups_x * num_groups_y * num_groups_z
glDispatchCompute(PARTICLE_GROUP_COUNT, 1, 1);

// 同步 计算着色器的 Invocations
glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);

// 实时改变 MVP 矩阵
glm::mat4 projection = glm::perspective(45.0f, aspect_ratio, 0.1f, 1000.0f);
glm::mat4 view = glm::mat4(1.f);
view = glm::translate(view, glm::vec3(0.0f, 0.0f, -60.0f));
view = glm::rotate(view, time * 100.0f, glm::vec3(0.0f, 1.0f, 0.0f));

glm::mat4 mvp = projection * view;

// ----------------------------------------------------
// 清屏，并切换到 渲染着色器
glClearColor(0., 0., 0., 0.);
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glDisable(GL_DEPTH_TEST);

glUseProgram(render_prog);
glUniformMatrix4fv(0, 1, GL_FALSE, (const GLfloat *)&mvp[0][0]); // 传入 MVP 矩阵
glBindVertexArray(render_vao); // 绑定 VAO

glEnable(GL_BLEND);
glBlendFunc(GL_ONE, GL_ONE);
glPointSize(2.0f);
glDrawArrays(GL_POINTS, 0, PARTICLE_COUNT); // 图元为 Points

last_ticks = current_ticks;

glutSwapBuffers();
}

计算着色器 GLSL代码：
#define STRINGIZE(a) #a

const char* compute_shader_source =
STRINGIZE(
#version 430 core\n
// Uniform Block
layout(std140, binding = 0) uniform attractor_block
{
vec4 attractor[64]; // xyz = position, w = mass
};

layout(local_size_x = 128) in;

layout(rgba32f, binding = 0) uniform imageBuffer velocity_buffer;
layout(rgba32f, binding = 1) uniform imageBuffer position_buffer;

uniform float dt = 1.0;

void main(void)
{
// 从 TBO 中取出数据，imageLoad 和 texelFetch 相似，因为省去了 filtering 的过程所以更高效
vec4 vel = imageLoad(velocity_buffer, int(gl_GlobalInvocationID.x));
vec4 pos = imageLoad(position_buffer, int(gl_GlobalInvocationID.x));

int i;

pos.xyz += vel.xyz * dt;
pos.w -= 0.0001 * dt;

for (i = 0; i < 4; i++)
{
vec3 dist = (attractor[i].xyz - pos.xyz);
vel.xyz += dt * dt * attractor[i].w * normalize(dist) / (dot(dist, dist) + 10.0);
}

if (pos.w <= 0.0)
{
pos.xyz = -pos.xyz * 0.01;
vel.xyz *= 0.01;
pos.w += 1.0f;
}
// 经过计算以后再写回 TBO
imageStore(position_buffer, int(gl_GlobalInvocationID.x), pos); // 类似的还有 imageSize
imageStore(velocity_buffer, int(gl_GlobalInvocationID.x), vel);
}
);
注：



1）
gl_WorkGroupSize：存储 local workgroup 的大小（三维）

gl_NumWorkGroups：存储了组在三个维度上的个数

gl_LocalInvocationID：当前 Invocation 在 local workgroup 中的位置（三维）
范围在 [uvec3(0), gl_WorkGroupSize - uvec3(1)] 之间

gl_LocalInvocationIndex：意义同上，区别在于它是一维的
它相当于 gl_LocalInvocationID.z * gl_WorkGroupSize.x * gl_WorkGroupSize.y + gl_LocalInvocationID.y * gl_WorkGroupSize.x + gl_LocalInvocationID.x

gl_GlobalInvocationID：当前 Invocation 在 global workgroup 中的位置（三维）
它相当于 gl_WorkGroupID * gl_WorkGroupSize + gl_LocalInvocationID

gl_WorkGroupID：当前 local workgroup 在 global workgroup 中的位置

范围在 [uvec3(0), gl_NumWorkGroups - uvec3(1)] 之间

2）local_size_x, local_size_y, local_size_z 声明了 local workgroup 的大小；

3）可以通过 glGetProgramiv() 搭配 GL_MAX_COMPUTE_WORK_GROUP_SIZE 查询 Local workgroup 的大小；

4）shared 类型的变量，意味着位于同一 local workgroup 中的 Invocation 共享该变量，通常访问共享变量的性能要优于访问 image 和 shader storage buffer；

5）Invocation 的同步
barrier：同步同一个 local workgroup 的 Invocation，确保所有 Invocation 都到达 barrier 之后才能往下执行；
memoryBarrier：确保所有内存的写入操作在此之前都已完成（没有数据驻留在缓存或者计划放入缓存）；

memoryBarrierAtomicCounter：等待所有更新原子计数器的操作完毕，才会继续执行；

memoryBarrierBuffer/memoryBarrierImage：等待所有写入 buffer 或 image 变量的操作都完成

memoryBarrierShared：等待所有更新 shared 变量的操作完毕，才会继续执行

❤ 但是以上这些函数并不能确保其他 invocation 都能到达这个点，所以仍然需要 barrier

groupMemoryBarrier 是高效版的 memoryBarrier，但是 groupMemoryBarrier 仅对 local workgroup 有效，而其他的 memoryBarrer 是全局的，即 global workgroup 的写入请求都已完成。

6）image类型（比如 imagebuffer ）可用于一般化的数据存储。image 类型与 sampler 类型相似，有两点区别：其一，image 类型仅表示单层的纹理，并没有完整的 mipmap 链；其二，image 类型并不支持一些 sampler 操作，比如滤波（filtering）和深度比较。（注：image 类型在声明时还需要带上格式布局修饰符——format layout qualifier，[例如 rgba32f，具体可以参考 OpenGL 红宝书第八版第11章 P566 ] 需要和 glBindImageTexture 中的类型一致）

相关资料：
【OpenGL】向Shader中传递数据

GLAPI/glTexBuffer