OpenGL 计算着色器

1.简介

计算着色器(Compute Shader)是一个特殊类型的着色器，其独立于OpenGL的图像渲染管道之外。他们被设计出用于充分利用图像处理器GPU的大规模并发计算能力，因此它们的作用并不仅仅限于进行图像渲染。需要注意⚠️的是计算着色器是OpenGL 4.3及其之上的特性。

接下来将会从如下三个点开始计算着色器来详细介绍计算着色器。

• 如何创建编译和分发计算着色器
• 如何在计算着色器之间传递数据
• 如何同步运行计算着色器，确保它们之间的有序性

2.使用计算着色器

和普通着色器一样，我们使用OpenGL的API来创建，填充源码，编译和附着计算作色器，并在完成这些操作后链接OpenGL程序。

// 1. 创建计算着色器
GLuint shader = glCreateShader(GL_FRAGMENT_SHADER);
// 2. 为计算着色器填充源码
glShaderSource(shader, 1, &shaderSource, nil);
// 3. 编译计算着色器
glCompileShader(shader);
// 4. 附着计算着色器
glAttachShader(self.glProgram, shader);
// 5. 链接OpenGL程序
glLinkProgram(self.glProgram);

计算着色器不能够和其他类型的着色器混用，也就是说在一个OpenGL的程序中不能同时附着计算着色器和其他类型的着色器。包含计算着色器的OpenGL程序也称为计算程序(compute program)，相对应的是包含其他类型着色器的程序称为图像程序(graphics program)。

除了一个专属的变量和标示符，计算着色器的写法和其他类型着色器类似，都是使用GLSL编写其源码，下面是一个简单的计算着色器源码。

#version 430 core
layout (local_size_x = 32, local_size_y = 32, local_size_z = 32) in;
void main() {
    // Do nothing
}

这里先要提到一个工作区(work group)的概念，这是一个抽象的概率，和GPU的处理单元无必然联系。工作组中的最小单位称为工作项(working item)，它们和一次计算着色器的调用相关。计算着色器都运行在GPU的特殊的处理单元上，一个处理单元同一时间只能够执行一次计算着色器的调用。并且这些处理单元有时间切片的概念，也就是说一个处理单元可能执行某个着色器内部分指令后再去执行其他着色器的指令。工作组是一个三维的抽象模型，计算程序运行时需要指定该工作组的大小，我们使用x、y、z三个分量来描述。

上面的代码中，第二行指定了工作组的大小为32✖️32✖️32，在着色器中也可以声明1维和2维的工作组，当x、y、z某个分量未指定时，系统OpenGL默认其值为1。下面声明了一个512大小的一维工作组。

layout (local_size_x = 512) in;

和其他类型的OpenGL着色器一样，计算着色器内部也能使用统一变量、统一闭包、着色器闭包等数据类型。为计算着色器内部这些变量赋值的方式也和普通的着色器一样，先调用函数glUseProgram()设置当前使用的程序，再调用类似函数glUniform4fv()就可以为对应的变量赋值。

2.1 运行计算着色器

指定OpenGL计算程序使用的工作组大小也称为工作组的调度，通过如下两个函数完成。

void glDispatchCompute(GLuint num_groups_x, 
                       GLuint num_groups_y,
                       GLuint num_groups_z);

void glDispatchComputeIndirect(GLintptr indirect);

第一个函数直接指定了全局工作组中包含的本地工作组数量，这在后面会详细介绍。而第二个函数使得我们可以通过一个缓存对象来指定这些参数，而参数indirect表示的是在缓存对象中内存偏移量。通常我们使用如下结构体的方式定义一套工作区参数。

typedef struct {
    GLuint num_groups_x; 
    GLuint num_groups_y; 
    GLuint num_groups_z;
    } DispatchIndirectCommand;

2.1.1 全局和局部工作组

计算着色器被运行在工作组(work group)中，每次调用函数glDispatchCompute和glDispatchComputeIndirectOpenGL内部都会先获得一个全局工作组 (global work group)，然后再将其细分为多个局部工作组(local work group)，局部工作组包含的处理单元数量由这两个函数只能够的参数指定，也就是x、y、z三个分量的乘积。局部工作组可以理解为具有三维立体性质。

单个局部工作组的三维尺寸可以通过函数glGetIntegeri_v和参数GL_MAX_COMPUTE_
WORK_GROUP_SIZE获得，其中参数target传入1、2、3分别对应获取x、y、z轴的最大处理单元数量，它们的值一定大于1024。通过指定参数GL_MAX_COMPUTE_WORK_
GROUP_INVOCATIONS查到的是整个工作组的最大大小，也可以理解为工作组的体积，也就是x、y、z三个分量乘积的最大数字，获取到的值也一定大于1024。

当计算着色器被链接到OpenGL程序后可以调用方法glGetProgramiv()并传入参数GL_COMPUTE_WORK_GROUP_SIZE获取其包含的计算着色器工作组的大小。示例代码如下。

int size[3];
glGetProgramiv(program, GL_COMPUTE_WORKGROUP_SIZE, size);
printf("Work group size is %d x %d % xd items.\n", size[0], size[1], size[2]);

一旦我们确定了计算着色器所需要的工作区大小，我们就能够通过函数glDispatchCompute和glDispatchComputeIndirect为计算程序调度工作组。在使用后面的函数时，相关的数据写入和读取使用到的缓存对象需要和GL_DISPATCH_INDIRECT_BUFFER这种类型的缓存绑定点关联。前面提到这两个函数调用后我们得到的是全局工作组。

需要注意的是，全局工作组的大小可以大于局部工作组，也就是说我们甚至可以使用一个3维的全局工作组来执行2维，甚至是1维的局部工作组，OpenGL内部会对全局工作组进行细分。

2.1.2 计算着色器的输入和输出

首先需要注意，计算着色器内部不包含任何内建输出变量，也不能像其他OpenGL着色器一样随意定义输出变量，它的输出变量定义必须遵循特有的规则。但是它仍然包含一些内建输入变量给出了其在局部工作组和全局工作组的位置以及一些其他必要信息。

内建变量gl_LocalInvocationID是赋值着色器运行的硬件单元在工作组中的索引，它是一个uvec3类型的变量，每个元素的取值范围分布从0到之前在着色器中声明的本地工作组大小local_size_x-1、local_size_y-1和local_size_z-1。变量gl_WorkGroupSize为工作组的大小，它也是一个uvec3类型的变量，其值为(local_size_x, local_size_y, local_size_z )。需要注意的是当你使用了1维或者2维的本地工作组，gl_LocalInvocationID在未使用的维度上为0，gl_WorkGroupSize则为1。

变量gl_NumWorkGroups表示本地工作组的数量，gl_WorkGroupID表示当前的本地工作组在全局工作组中的索引值。这几个变量的的含义入下图所示，图中描述了一个包含3✖️4✖️8个局部工作组的全局工作组。其中每个局部工作组包含了6✖️4✖️1个处理单元。

内建变量示意图

OpenGL也提供变量gl_GlobalInvocationID是我们能够快速定位当前着色器运行的硬件在全局工作组中的索引。它内部的计算过程如下。

gl_GlobalInvocationID = gl_WorkGroupID * gl_WorkGroupSize + gl_LocalInvocationID;

变量gl_LocalInvocationIndex可以理解为是gl_LocalInvocationID的扁平模式，它内部的计算公式如下。也就是将3维的索引映射到1维中索引中。

gl_LocalInvocationIndex =
    gl_LocalInvocationID.z * gl_WorkGroupSize.x * gl_WorkGroupSize.y + 
    gl_LocalInvocationID.y * gl_WorkGroupSize.x + 
    gl_LocalInvocationID.x;

当前着色器运行在3维工作组中的索引值可以被用于数组、纹理坐标、随机种子等场景中。后面会详细介绍。

计算着色器产生的数据必须直接写入到显存中。例如在计算着色器中，你可以将数据写入到着色器闭包中，使用如imageStore相关的图像操作函数，或者递增和递减原子计数等。这些操作有副作用，这意味着这些操作直接更新了显存中的值将会被检测到，或者有明显的外部影响。

在下面的计算着色器中，从一张图片中读取了像素值，并取每个像素的补色，再将得到的结果写入到另外一张图片中。

#version 430 core
layout (local_size_x = 32, local_size_y = 32) in;
layout (binding = 0, rgba32f) uniform image2D img_input;
layout (binding = 1) uniform image2D img_output;

void main() {
    vec4 texel;
    ivec2 p = ivec2(gl_GlobalInvocationID.xy);
    texel = imageLoad(img_input, p); 
    texel = vec4(1.0) - texel; 
    imageStore(img_output, p, texel);
}

在执行计算上面的着色器之前，我们需要为着色器中声明的两个统一图片变量赋值。这里需要注意我们的工作组大小为32✖️32，因此我们的图片大小最好是它的整数倍。然后我们再调用函数glDispatchCompute分配硬件资源。

// Bind input image
glBindImageTexture(0, tex_input, 0, GL_FALSE,
                   0, GL_READ_ONLY, GL_RGBA32F);
// Bind output image
glBindImageTexture(1, tex_output, 0, GL_FALSE,
                   0, GL_WRITE_ONLY, GL_RGBA32F);
// Dispatch the compute shader
glDispatchCompute(IMAGE_WIDTH / 32, IMAGE_HEIGHT / 32, 1);

2.2 计算着色器通信

计算着色器在多个工作组中每个工作单元上运行方式，和曲面细分控制着色器在多个批次中的每个控制点执行逻辑类似，工作组和控制点的批次都是根据调用组来创建的。对于曲面细分控制着色器，在单个批次控制点的执行逻辑内，通过使用关键字patch声明变量，并且在它们被正确同步执行的情况下，这些着色器能够正确的写入和读取这些patch关键字声明的变量。一个批次内的曲目细分着色器调用之间有有限种通信模式，这里不详细讨论。但是这种通信方式有很大的限制，例如能够使用patch修饰的变量数量受到严格限制，单个批次内的控制点数量也很小。

计算着色器有着类似的机制，但是其更加灵活和强大。你可以使用shared关键字修饰一个变量，使其在同一个局部工作组调用内的所有计算着色器中成为共享变量(shared variables)。访问共享变量的速度要明显高于通过图片或者存储闭包访问主显存中的数据。因此当多个着色器调用需要访问到同一份数据时，你可以先将这份数据从显存中拷贝到共享变量中，然后在需要的地方访问这些变量，甚至你可以在合适的地方更新它们，最后在将处理好的结果写回到显存中。

2.2.1 同步计算着色器

尽管我们说一个工作组内的计算着色器调用是并发执行的，但其内部仍然有着时间切片结构设计，也就是出它们不能真正的彻底并发执行。处理器通常将每个局部工作组分为更小的块，常见的块大小为16、32或者64，每个块内的调用具有一致性(lockstep)，也就是一个块内的计算着色器都会同时调用某一个指令。OpenGL会为这些块分配计算资源，块与块之间的调用是完全无序的，但这些块之间可以进行通信。

OpenGL提供了一种称为屏障(barrier)的同步机制，也就是函数barrier()，该函数建立了一个流动的控制屏障(flow control barrier)。当在着色器内部调用函数barrier()后，它会阻塞当前的逻辑直到同一个局部工作组内的所有计算着色器都执行到该函数。在前面讲到曲面细分控制着色器时已经简单聊过了这个点。在一个时间切片架构中，执行函数barrier()意味着着色器以及其存在的工作块将放弃它们的时间片段，使得其他着色器在它运行到该函数之前仍能执行。一旦一个局部工作组内所有的调用执行到了这个函数，它们会将不再等待，会继续执行剩下的逻辑。

当使用共享显存特性时，流控制屏障十分重要，他们使你知道其他着色器调用内执行的逻辑到达了和当前着色器的同一个逻辑点。如果当前着色器需要向共享显存中写入数据，那么你必须确保其他着色器也完成数据的写入然后，这样在执行后续的逻辑从这个共享显存中读取数据才是安全的。下面是一个使用屏障的计算着色器源码。

#version 430 core
layout (local_size_x = 1024) in;
layout (binding = 0, r32ui) uniform uimageBuffer image_in; 
layout (binding = 1) uniform uimageBuffer image_out;
shared uint temp_storage[1024];

void main(void) {
    // Load from the input image
    uint n = imageLoad(image_in, gl_LocalInvocationID.x).x; 

    // Store into shared storage
    temp_storage[gl_LocalInvocationID.x] = n;

    // Uncomment this to avoid the race condition
    // barrier();
    // memoryBarrierShared();

    // Read the data written by the invocation ‘‘to the left’’
    n = temp_storage[(gl_LocalInvocationID.x - 1) & 1023];
    // Write new data into the buffer
    imageStore(image_out, gl_LocalInvocationID.x, n); 
}

该着色器从一个图片缓存对象中读取数据到一个数组类型的共享变量中，写入使用的数组索引和其调用索引一致。然后再从该数组类型的共享变量中读取其存入位置左侧是数据，并将得到的结果写入到一个图像缓存对象中。总代来说这个着色器的逻辑可以理解为将输入的图片缓存中的数据向左移动了一位。下图颜色了这个着色器的执行逻辑。

计算着色器中的竞争

图中颜色了工作组中连续的ABCD4个计算着色器的调用竞争一个GPU处理单元的情形，t0到t12是该处理单元的12个时间切片。我们可以明显看到在t5时间片段上，计算着色器的调用C获取到了硬件资源，它执行了从共享变量中读取B调用计算出数据的指令，然而此时B调用并没有执行到写入数据的指令，这样一定会导致数据错误。这就是计算着色器的资源竞争(race condition).

将上面程序清单中的函数barrier()解注释，此时这4个计算着色器调用如下图所示。

使用流控制屏障处理计算着色器竞争

图中t0时间切片着色器调用A执行到了barrier()，此后该调用将会放弃获得到的时间片段。直至在t4时间段D调用也执行到了屏障函数，此时其余的调用就会收到通知使得其内部的逻辑能够继续执行。这样我们最后就能够得到正确的数据。

需要注意的是尽管上面的两张图片中仅仅演示了4个着色器调用去竞争一个计算单元。但是实际上GPU在运行时通常有着几百个线程去竞争几十个计算资源。

3 示例

接下来我们将用几个实际的例子来颜色计算着色器如何使用。在第一个例子并行前缀求和(parallel prefix sum)中，我们将颜色如果高效并行的执行算数计算，尽管这个计算第一眼看上去像是串行的过程。第二个例子演示了群集算法(flocking algorithm)。

3.1 计算着色器并行前缀求和

前缀求和操作(prefix sum operation)，已知一个数组A作为输入变量，计算出一个新的数组B，数组B中的每个元素的值是数组A中从索引0的元素累加到当前索引元素的和。包含当前索引元素的前缀求和操作被称为是包含前缀求和(inclusive prefix sum)，不包含当前元素的称为排除(exclusive prefix sum)。下面是常见的c++代码执行的逻辑。

void prefix_sum(const float_* in_array, 
                float * out array, 
                int elements, 
                bool inclusive)  {

    float f = 0.0f;
    int I;
    if (inclusive) {
        for (i = 0; i < elements; i++) {
            f += in_array[I];
            out_array[i] = f;
        }
    } else {
        for (i = 0; i < elements; i++) {
            out_array[i] = f;
            f += in_array[I];
        }
    }
}

一个包含前缀求和的实际计算结果如下图。

包含前缀求和的实例

表面上看输出数组B中的每个元素都依赖于它前一个元素的值和数组A中当前索引值的和，这看上去似乎只能是一个串行的计算逻辑，但是实际上它可以拆解成大量的并行计算。假定我们输入数组A包含元素I0到I3，输出数组B包含元素O0到O3，则输出数组的计算公式可以表示如下。

并发的关键是将这些计算任务拆解为更小相互独立的子任务，从而使得它们可以并发执行。我们可以选择如下的方式来执行前缀求和算法，第一步如下。

然后我们再执行如下的操作。

可以看到尽管两个大阶段之间的计算是有依赖关系的，但是每个阶段内的多个计算之间都是相互独立的，因此能够并行计算。实际上我们能够将任意尺寸的输入数组拆分为多个块，直到所有的块包含的元素都不超过两个。拆分过程如下图所示。

前缀求和块拆分示意图

通常使用块拆分方法实际执行的加法运算次数比直接计算的方式更多。上面的例子中如果使用常规的顺序执行算法需要15次加法运算，而使用块拆分方式需要32次加法运算。但是这并不意味着块拆分方案的总耗时更长，加入我们能够有8个资源并行计算，那么世界上我们只需要相当于四个单次计算的时间就能完成整个前序求和，这要比常规的串行计算几乎快出4倍。可以想象随着输入变量数组的容量增加，这种块拆分的计算方式优势将更明显。

下面是该方法的着色器实现源码。这里并不会详细讲道拆分算法的原理，可以自己看下面例子中相关部分代码，我们更关注的是和OpenGL相关的知识。其中函数barrier和memoryBarrier的区别在于前者是指令控制屏障，这在前面已经讲过，后者是显存写入屏障，这意味着它会阻塞向显存中的数据写入，确保该指令后一定能够正确读取前面写入的数据。而memoryBarrierShared是更精细的控制，它只对声明为shared的变量写入才会有内存屏障作用。

#version 430 core
layout (local_size_x = 1024) in;
layout (binding = 0) coherent buffer block1 {
    float input_data[gl_WorkGroupSize.x];
};
layout (binding = 1) coherent buffer block2 {
    float output_data[gl_WorkGroupSize.x];
};
shared float shared_data[gl_WorkGroupSize.x * 2];

void main(void) {
    uint id = gl_LocalInvocationID.x; 
    uint rd_id;
    uint wr_id;
    uint mask;
    // The number of steps is the log base 2 of the
    // work group size, which should be a power of 2
    const uint steps = uint(log2(gl_WorkGroupSize.x)) + 1; uint step = 0;
    // Each invocation is responsible for the content of
    // two elements of the output array
    shared_data[id * 2] = input_data[id * 2];
    shared_data[id * 2 + 1] = input_data[id * 2 + 1];
    // Synchronize to make sure that everyone has initialized
    // their elements of shared_data[] with data loaded from
    // the input arrays
    barrier();
    memoryBarrierShared();
    // For each step...
    for (step = 0; step < steps; step++) {
        // Calculate the read and write index in the
        // shared array
        mask = (1 << step) - 1;
        rd_id = ((id >> step) << (step + 1)) + mask;
        wr_id = rd_id + 1 + (id & mask);
        // Accumulate the read data into our element
        shared_data[wr_id] += shared_data[rd_id];
        // Synchronize again to make sure that everyone
        // has caught up with us
        barrier();
        memoryBarrierShared();
    }
    // Finally write our data back to the output image
    output_data[id * 2] = shared_data[id * 2];
    output_data[id * 2 + 1] = shared_data[id * 2 + 1];
}

算法内部输入变量input_data仅包含1024个元素，而计算过程中取值input_data[id * 2]没有触发数组越界需要验证。

该计算着色器声明了一个大小为1024的局部工作区，这意味着共享变量数组shared_data会有2048个元素，多一倍的设计是为了算法的实现，这里不具体分析。每次调用的开始都从输入的数据中读取了两个原始元素然后立即调用了指令屏障和共享内存屏障函数，是的每个计算着色器的调用在其内部的循环开始之前共享变量shared_data中已经填充完初始数据。

在每次迭代内部都先计算出待做加法的两个数组下标，然后在求和，紧接着调用了指令屏障和共享内存屏障函数，保证所有的计算着色器调用都执行完同一步后再同时进入到下一次的循环中。最后将计算得到的结果写入到输出缓存对象中。

前缀计算算法能够以一种拆分的方式(separable manner)应用到多维到数据集中，如二维的图像。在前面讲到辉光滤镜时，我们使用到的高斯滤镜就是一个拆分算法应用的实例。对二维图像执行拆分前缀算法时我们首先对每行像素执行前缀求和计算得到一张新的图像，然后我们再对每列像素执行前缀求和计算得到最终的图像。最终得到的图像中每个像素的值就是从原点到该像素点构成的矩阵中所有像素点值的和。该过程的示意图如下。

2维前缀求和示意图

上图中a为原始图像，对每行像素前缀求和后得到了图b，再对每列像素前缀求和后得到了图c，这被称为求和空间表(summed area table)，在很多计算机图像程序中都是一个十分重要的概念。

我们对前面的计算着色器源码进行简单的修改，使用2维图像纹理替代着色器存储闭包。修改后的代码如下。需要注意的是，着色器读取的是行像素却写入到了列像素中，着意味着输出图像将是输入图像的转置。我们巧妙的执行这个着色器两次，两次转置后我们最终得到的图像和原始的图像就是同一个方向的。当然你也可以选择不使用转置的复杂方法。

#version 430 core
layout (local_size_x = 1024) in;
shared float shared_data[gl_WorkGroupSize.x * 2];
layout (binding = 0, r32f) readonly uniform image2D input_image; 
layout (binding = 1, r32f) writeonly uniform image2D output_image;

void main(void) {
    uint id = gl_LocalInvocationID.x; 
    uint rd_id;
    uint wr_id;
    uint mask;
    ivec2 P = ivec2(id * 2, gl_WorkGroupID.x);

    const uint steps = uint(log2(gl_WorkGroupSize.x)) + 1;
    uint step = 0;

    shared_data[id * 2] = imageLoad(input_image, P).r; 
    shared_data[id * 2 + 1] = imageLoad(input_image, P + ivec2(1, 0)).r;
    barrier();
    memoryBarrierShared();

    for (step = 0; step < steps; step++) {
        mask = (1 << step) - 1;
        rd_id = ((id >> step) << (step + 1)) + mask;
        wr_id = rd_id + 1 + (id & mask);
        shared_data[wr_id] += shared_data[rd_id];
        barrier();
        memoryBarrierShared();
    }

    imageStore(output_image, P.yx, vec4(shared_data[id * 2]));     
    imageStore(output_image, P.yx + ivec2(0, 1), vec4(shared_data[id * 2 + 1]));
}

在上面的代码中局部工作组仍是一维的，在第一次处理图像时，我们调度一个局部工作组数量和图像行数相等的全局工作组，同样的在第二次处理图像时，我们调度的全局工作组包含的局部工作数量和原始图像的列数相同。这样就可以使用本地工作组调用索引和本地工作组在全局工作组中的索引来计算出对应的图像坐标。

给定一个求和区域表(summed area table)，只需要4个特定的点就可以计算出原始图像中任意一个矩形范围内的所有值的和。其计算过程如下。

在求和区域矩形中计算原始图像中任意一个矩形内的所有像素和

当我们能够计算原始图像任意矩形内所有像素值和时，我们再除以像素的数量，就能很快的计算改矩形内像素值平均值。尽管这种求平均值的方式很粗暴，但是在特定的程序中还是很有用的，这种计算方式的一种应用就是盒式滤镜(box filter)。更有意义的是，当我们在计算以任意一个像素点为中心，特定大小的矩形内像素平均值创建可变尺寸滤镜(variable sized filter)时，我们还能对每个像素选择不同的矩形大小。

下面是这种可变尺寸滤镜的一个例子。同种从左向右滤镜效果逐渐加强，我们可以明显看到右边的图像比左边度图像更加模糊。

可变尺寸滤镜处理后的图像

基于同样的原理我们可以做更多有趣的效果，如景深效果(depth of field effect)。对于相机而言和这个特效相关的两个概念是焦距(focal distance)和焦深(focal depth)。这里仅对上述两个概念做简单介绍。

摄像机的镜头是一系列凸透镜和凹透镜的组合。整体上镜头符合凸透镜成像规律，即1 /f = 1/a + 1/b，其中f指的是焦距，简单的说就是镜头透镜组的光学中心到焦点的距离；a指的是物距，也就是被拍摄物体到镜头光学中心的距离；b指的是相距，即被拍摄物体向各个方向发散的光经过透镜汇聚后的实像平面和焦点的距离。

在确定到镜头某个距离平面图像符合凸透镜成像公式，能够完全清晰显示时。从该平面靠近或者远离镜头仍有一定的距离空间内的景色能够清晰成像，这个距离长度称为景深。在这个区间内最近和最远平面的实像面之间的距离称为焦深。

焦距越长，透镜边缘的光线被弯折的角度越小，其两个弥散圆到实像之间的距离越大，也就是焦深越大，即景深越大。反之则景深越小。

下图是一个用专业相机拍出的一个景深较低的照片，图中最前面的杯子是完全对焦的，向远处图像开始变得越来越模糊。我们能够使用基本的盒式滤镜模拟出这种实际的光学效果。

使用数码相机拍出的具有景深效果图片

首先我们渲染出一张无景深效果的图片，但是需要存储每个片段的深度(也就是它们和相机的距离)。当深度值和模拟出的相机焦距相同时，片段被认为和正常的光学相机下的完全对焦类似。当片段的深度值偏离焦距时，我们增加该片段的模糊度。

下面演示使用该原理的一个例子。在示例程序中，对前文的计算着色器稍作修改，将计算的每个值由float类型改变为vec3类型，从而计算出图片的求和区域表。另外我们将片段的深度值存储在alpha通道中，然后通过如下的片段着色器获取每个片段的深度值。通过深度值计算出应用模糊效果的像素块尺寸(即对宽高为多少个像素的区域进行模糊)，以当前片段为中心，从求和区域表中计算出模糊效果区域的所有像素在RGB各个通道的和。最后再计算出平均值，即当前像素的颜色值，从而模拟景深效果。

#version 430 core
layout (binding = 0) uniform sampler2D input_image;
layout (location = 0) out vec4 color;

uniform float focal_distance = 50.0; 
uniform float focal_depth = 30.0;

void main(void) {
    // s will be used to scale our texture coordinates before looking up data in 
    // our SAT image.
    vec2 s = 1.0 / textureSize(input_image, 0);
    // C is the center of the filter
    vec2 C = gl_FragCoord.xy;
    // First, retrieve the value of the SAT at the center of the filter. The last 
    // channel of this value stores the view-space depth of the pixel.
    vec4 v = texelFetch(input_image, ivec2(gl_FragCoord.xy), 0).rgba;
    // M will be the radius of our filter kernel
    float m;
    // For this application, we clear our depth image to zero before rendering to 
    // it, so if it’s still zero, we haven’t rendered to the image here. Thus, we 
    // set our radius to 0.5 (i.e., a diameter of 1.0) and move on.
    if (v.w == 0.0) {
        m = 0.5; 
    } else {
        // Calculate a circle of confusion
        m = abs(v.w - focal_distance);
        // Simple smoothstep scale and bias. Minimum radius is 0.5 (diameter 1.0), 
        // maximum is 8.0. Box filter kernels greater than about 16 pixels don’t 
        // look good at all.
        m = 0.5 + smoothstep(0.0, focal_depth, m) * 7.5;
    }

    // Calculate the positions of the four corners of our area to sample from.
    vec2 P0 = vec2(C * 1.0) + vec2(-m, -m);
    vec2 P1 = vec2(C * 1.0) + vec2(-m, m);
    vec2 P2 = vec2(C * 1.0) + vec2(m, -m); 
    vec2 P3 = vec2(C * 1.0) + vec2(m, m);
    // Scale our coordinates.
    P0 *= s;
    P1 *= s;
    P2 *= s;
    P3 *= s;
    // Fetch the values of the SAT at the four corners
    vec3 a = textureLod(input_image, P0, 0).rgb;
    vec3 b = textureLod(input_image, P1, 0).rgb; 
    vec3 c = textureLod(input_image, P2, 0).rgb; 
    vec3 d = textureLod(input_image, P3, 0).rgb;
    // Calculate the sum of all pixels inside the kernel.
    vec3 f = a - b - c + d;
    // Scale radius -> diameter.
    m *= 2;
    // Divide through by area
    f /= float(m * m);
    // Output final color
    color = vec4(f, 1.0); }

上面的片段着色器读取了一个图像纹理，该纹理的RGB通道为像素颜色数据，alpha通道为预先计算好的片段景深数据，此外还接收外部输入的模拟相机参数。通过片段的深度和相机的景深差值的绝对值计算出滤镜矩形尺寸（该片段需要应用模糊效果的像素块尺寸）。然后从求和区域表中读取矩形的四个顶点对应的值，仔计算出该区域内的像素的平均颜色，并将其写入到该片段的颜色缓存中。下图是程序运行后得到的图片。

上图中，我们设定向屏幕深处的第2个龙雕刻屏幕为模拟的合焦面(当前的传感器能够接收到完全聚焦清晰的景物平面)，该平面上的片段最清晰，远离这个平面逐渐模糊，超过景深值后模糊程度不再增加。上图中的第二个龙雕塑最清晰，其余的雕塑都是模糊的，这样就在这一排龙雕塑场景中模拟出了景深效果。

通过改变设定的焦距合景深，可以模拟出更多有趣的场景，下图的左侧图片中，向屏幕深处第1个龙雕塑位于模拟合焦面上。在中间的图片中，最远处的龙雕塑位于模拟合焦面上，并且景深设置很浅，使得其他的雕塑都很模糊。在右侧的图中，设定的景深值较大，大部分雕塑看上去都比较清晰。

上面的例子中为了简化程序使其更容易理解，输入纹理的数据格式被指定为32位浮点型，这样造成两个后果。第一个很好理解，那就是纹理会占据更多的显存空间，第二个后果是会有精度丢失。

在了解精度丢失之前先回顾一下浮点型数据的特点，当其尾数越大，也就是其占据更多的有效位时，其精度越低。而上文提到的求和区域表最后一个元素的值是整个图片前面所有像素颜色值的和，可以想象这个值一定会很大。而我们在计算单个片段的具体颜色时，使用了求和区域表中的4个坐标计算模糊滤镜生效像素块的平均值，这样损失的精度就会被转移到最终计算出的结果中，导致图片出现噪声。

为了解决这个问题，我们可以通过以下三部来实现。

• 使用16位浮点型数据来保存最初渲染得到的图像。
• 将片段的深度数据存储在一个单一的纹理中(或者将它们存入深度缓存中)。
• 在渲染最初的图像时，将每个像素的值都减去0.5，这样求和区域表中的值就会在0附近，而不会出现很大的值，这样能提高精度。

3.2 计算着色器群

下面是一个利用计算着色器实现集群算法的例子。集群算法(flocking algorithms)通过单独更新一个组内某个元素的属性实现一些突发行为。在自然界中这种突发行为很常见，如蜂群、鸟群和鱼群的移动，尽管其中的每个成员互不交流，但是它们的移动轨迹表面上看仍然很协调。也就是说个体的决定仅仅取决于它对周围其他同伴的感知，目前的观察表明鱼群中并不存在指挥者。因为每个成员都是独立的，因此它们自己的属性可以并行计算，这就能够充分的利用GPU的性能。

我们使用了一个计算着色器实现集群算法。集群中的每个成员都独立的存储于一个着色器闭包缓存对象中。计算着色器从一个缓存对象中读取数据，计算出某个成员的位置和速度属性后写入到另一个缓存对象中。该缓存对象接下来和一个顶点缓存绑定，再使用绑定后的缓存作为顶点着色器的输入。集群中的每个元素都通过一次调用绘制命令调用来渲染。顶点着色器负责将模型(在这个例子中就是纸飞机)变换到正确的位置和方向。接着算法开始循环，进入下一次的迭代逻辑，再次执行计算着色器，复用上一次迭代计算得到的位置和速度数据。需要注意的是在这整个过程中，数据都只存在GPU的显存中，CPU并未参与任何计算。

我们使用一对缓存对象来存储集群中每个元素当前的位置，统一的我们使用两个顶点数组对象(VAO)来管理循环中交替的两个阶段的顶点数组状态。这些顶点数组对象都还持有模型的顶点数据。这种双缓存的设计模式是为了避免我们正在使用某个顶点数组对象渲染场景的同时去更新这些缓存的部分内容。下图简要说明了集群算法实现的两个交替阶段。

图片的左上角我们更新了偶数帧。索引为0的缓存对象作为计算着色器的输入，其中包含位置和速度数据，计算着色器计算出新的位置和速度值后将它们写入到索引为1的缓存对象中。接下来执行渲染逻辑，即图中右上角部分，和计算着色器一样使用0号索引的缓存作为顶点着色器的输入。需要注意的是这里我们在更新阶段和渲染阶段使用了同一个作为输入，这样使得它们之间并无依赖。这意味着OpenGL能够在更新阶段完成前开始渲染阶段的工作。另外渲染阶段的顶点着色器还从额外的几何缓存中读取模型的顶点数据。

在图片的左下角中，我们接着处理下一帧。首先我们交换的缓存，现在1号索引的缓存被作为输入，而0号索引的缓存被作为输出。接着来到上图右下角部分，我们渲染下一帧的画面。同样的1号索引缓存被错位顶点着色器的输入。需要注意的是，基偶阶段的顶点数组对象是用来同一个几何缓存作为模型的顶点数据，因为模型的数据不会变换，没有必要再拷贝一份。

前文描述的逻辑代码实现如下。

glGenBuffers(2, flock_buffer); 
glBindBuffer(GL_SHADER_STORAGE_BUFFER, flock_buffer[0]); 
glBufferData(GL_SHADER_STORAGE_BUFFER, 
             FLOCK_SIZE * sizeof(flock_member), 
             NULL, 
             GL_DYNAMIC_COPY);
glBindBuffer(GL_SHADER_STORAGE_BUFFER, flock_buffer[1]); 
glBufferData(GL_SHADER_STORAGE_BUFFER,
             FLOCK_SIZE * sizeof(flock_member), 
             NULL,
             GL_DYNAMIC_COPY);

glGenBuffers(1, &geometry_buffer);
glBindBuffer(GL_ARRAY_BUFFER, geometry_buffer); 
glBufferData(GL_ARRAY_BUFFER, sizeof(geometry), geometry, GL_STATIC_DRAW);
    
glGenVertexArrays(2, flock_render_vao);

for (i = 0; i < 2; i++) {
    glBindVertexArray(flock_render_vao[i]); 
    glBindBuffer(GL_ARRAY_BUFFER, geometry_buffer); 
    glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 0, NULL); 
    glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, 0, 
                          (void *)(8 * sizeof(vmath::vec3)));
    glBindBuffer(GL_ARRAY_BUFFER, flock_buffer[i]); 
    glVertexAttribPointer(2, 3, GL_FLOAT, GL_FALSE, sizeof(flock_member), NULL); 
     // 书中第一个组件在顶点属性数组的第一个元素中的偏移量指定为vec4，个人认为应该为vec3
    glVertexAttribPointer(3, 3, GL_FLOAT, GL_FALSE, sizeof(flock_member),
                          (void *)sizeof(vmath::vec4)); 
    // 开启2、3号顶点属性多实例特性，每1次实例调用改变一次
    glVertexAttribDivisor(2, 1);
    glVertexAttribDivisor(3, 1);

    glEnableVertexAttribArray(0);
    glEnableVertexAttribArray(1);
    glEnableVertexAttribArray(2);
    glEnableVertexAttribArray(3);
}

在运行上面的代码之前，我们需要为集群中的每个成员指定初始的位置和速度数据，我们可以生存一个随机向量表示其方向，并将所有成员的速度设置为0。

现在我们需要一个渲染循环来完成集群成员的属性更新，以及这些成员的渲染，渲染循环的代码实现如下。我们能够清晰的看到其中存在两个阶段不断循环。首先将OpenGL程序update_program设置为当前活跃程序，并使用改程序更新集群成员的位置和速度数据。将绑定至0号和1号GL_SHADER_STORAGE_BUFFER索引的着色器存储缓存分别用于读取和写入数据，接着开始调用计算着色器开始计算。

接下来清空屏幕，激活渲染程序，更新仿射矩阵，绑定顶点数组对象，渲染场景。渲染的实例数量和模拟的集群成员数量相同，顶点的数量是描述单个纸飞机模型的顶点数量。

glUseProgram(flock_update_program);
vmath::vec3 goal = vmath::vec3(sinf(t * 0.34f),
                               cosf(t * 0.29f),
                               sinf(t * 0.12f) * cosf(t * 0.5f));
goal = goal * vmath::vec3(15.0f, 15.0f, 180.0f);

glUniform3fv(uniforms.update.goal, 1, goal);

glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, flock_buffer[frame_index]); 
glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, flock_buffer[frame_index ^ 1]);
 
glDispatchCompute(NUM_WORKGROUPS, 1, 1);

glViewport(0, 0, info.windowWidth, info.windowHeight); 
glClearBufferfv(GL_COLOR, 0, black); 
glClearBufferfv(GL_DEPTH, 0, &one);
    
glUseProgram(flock_render_program);
vmath::mat4 mv_matrix = vmath::lookat(vmath::vec3(0.0f, 0.0f, -400.0f),
                        vmath::vec3(0.0f, 0.0f, 0.0f),
                        vmath::vec3(0.0f, 1.0f, 0.0f));
vmath::mat4 proj_matrix = 
    vmath::perspective(60.0f, 
                       (float)info.windowWidth / (float)info.windowHeight,
                       0.1f, 3000.0f);
vmath::mat4 mvp = proj_matrix * mv_matrix;

glUniformMatrix4fv(uniforms.render.mvp, 1, GL_FALSE, mvp); 
glBindVertexArray(flock_render_vao[frame_index]); 
glDrawArraysInstanced(GL_TRIANGLE_STRIP, 0, 8, FLOCK_SIZE); 
frame_index ^= 1;

我们继续分析着色器内部的代码，这是这个程序最有趣的部分。集群算法依赖于一套规则，这个规则决定其内部每个成员的运动方向。每一条规则都认为单个成员的属性更新都是独立的。大多数规则都要求当前成员变量能够读取其他成员变量的信息，因此在程序update_program中使用了共享存储缓存来解决这个问题。下面是计算着色器的前部分代码，其中声明了将要使用到的统一变量，成员属性结构，两个缓存对象分别负责读取和写入数据，以及一个各个成员在更新属性时所需要使用到的共享变量。

#version 430 core
layout (local_size_x = 256) in;

uniform float closest_allowed_dist2 = 50.0; 
uniform float rule1_weight = 0.18;
uniform float rule2_weight = 0.05;
uniform float rule3_weight = 0.17;
uniform float rule4_weight = 0.02; 
// 飞行目的地
uniform vec3 goal = vec3(0.0); 
uniform float timestep = 0.5;

struct flock_member {
    vec3 position;
    vec3 velocity; 
};

layout (std430, binding = 0) buffer members_in {
    flock_member member[];
} input_data;

layout (std430, binding = 1) buffer members_out {
    flock_member member[];
} output_data;

shared flock_member shared_member[gl_WorkGroupSize.x];

在我们的实例程序中，我们用于确定每个成员变量属性的一套规则如下。

• 成员之间不能相互碰撞，任何时候它们之间的距离都应该大于一个固定值。
• 每个成员都会尝试调整自己的飞行方向，使其与周围的成员一致。
• 所有成员都尝试飞向同一个目标。
• 每个成员都尝试和整个集群的飞行速度保持一致，它们都会尝试将自己保持在集群的中心。

前两个规则是成员内部规则，也就是说每个成员之间的相互影响被认为是独立的，也就是说它们不会依赖于整个集群。下面是第一条规则的代码实现。如果某个成员变量和另外一个成员变量过近，将它移开即可。

vec3 rule1(vec3 my_position, 
           vec3 my_velocity, 
           vec3 their_position, 
           vec3 their_velocity) {

    vec3 d = my_position - their_position; 
    if (dot(d, d) < closest_allowed_dist2) {
        return d; 
    } else {
        return vec3(0.0);
    }
}

第二条规则等代码如下，它返回了当前成员的速度向量。首先计算出两个成员之间的速度向量差，然后使用距离的平方加调整常量的倒数作为权重，得到最终的速度向量。这里使用10作为调整常量是为了增大分母，从而得到的速度更小，使得集群效果模拟更稳定。

vec3 rule2(vec3 my_position, 
           vec3 my_velocity,
           vec3 their_position, 
           vec3 their_velocity) {

    vec3 d = their_position - my_position; 
    vec3 dv = their_velocity - my_velocity; 
    return dv / (dot(d, d) + 10.0);
}

第三条和第四条规则对于每个成员只会应用一次，另外第四条规则在计算时还需要使用到整个集群成员的平均位置。这部分的代码将在计算着色器的主体代码片段中涉及。

在计算着色器的主体代码中，集群被分为了多个组，每个组可以认为其对应一个本地工作组组，其成员量为256。由于集群的每个成员都要和其他所有成员进行交互，该算法的时间复杂度是O(N^2)。 另外使用了一个共享存储缓存来保存成员的属性，使得一个本地工作组中的每次计算着色器调用可以直接从中获取其他成员的属性，而不需要从着色器存储闭包中读取，这样可以提高效率。

对于每一个集群成员都对应了计算着色器的调用，我们循环遍历所有的局部工作组，将其在每个局部工作组中对应位置的成员属性存入共享变量shared_member中。然后调用内存屏障和指令屏障函数，确保同一个局部工作组内256次调用都完成共享数据拷贝，然后再执行后面的指令。接下来我们遍历整个共享数组，对其中每个非自身成员(即取到的成员不是自己)应用一次规则1和规则2，将得到的结果进行累加。然后再次调用内存屏障和指令屏障函数，确保单个局部工作组内所有着色器调用都执行到这条指令后再开启下一次的循环，好再次向共享变量中写入新的数据。计算着色器的主体代码如下。

void main(void) {
    uint i, j;
    int global_id = int(gl_GlobalInvocationID.x); 
    int local_id = int(gl_LocalInvocationID.x);

    flock_member me = input_data.member[global_id]; 
    flock_member new_me;
    vec3 acceleration = vec3(0.0);
    vec3 flock_center = vec3(0.0);

    for (i = 0; i < gl_NumWorkGroups.x; i++) {
        flock_member them = input_data.member[i * gl_WorkGroupSize.x + local_id];
        shared_member[local_id] = them; memoryBarrierShared();
        barrier();

        for (j = 0; j < gl_WorkGroupSize.x; j++) {
            them = shared_member[j];
            flock_center += them.position;
            if (i * gl_WorkGroupSize.x + j != global_id) {
                acceleration += rule1(me.position,
                                      me.velocity,
                                      them.position,
                                      them.velocity) * rule1_weight;
                acceleration += rule2(me.position,
                                      me.velocity,
                                      them.position,
                                      them.velocity) * rule2_weight;
            }
        }
        barrier(); 
    }

    flock_center /= float(gl_NumWorkGroups.x * gl_WorkGroupSize.x); 
    new_me.position = me.position + me.velocity * timestep;
    acceleration += normalize(goal - me.position) * rule3_weight; 
    acceleration += normalize(flock_center - me.position) * rule4_weight; 
    new_me.velocity = me.velocity + acceleration * timestep;
    if (length(new_me.velocity) > 10.0) {
        new_me.velocity = normalize(new_me.velocity) * 10.0;
    }
    new_me.velocity = mix(me.velocity, new_me.velocity, 0.4);
    output_data.member[global_id] = new_me;
}

每个成员处理和其他所有成员比较并应用规则1和规则2外，计算出的位移量还需要进行调整，是的每个成员还需要向着共同的目标，以及集群中心移动，我们使用规则3和规则4来完善这个模拟算法。此外，如果计算出的最新速度太高，我们需要将其截断到允许的最大值。另外，我们并不会直接输出计算出的新速度，相反我们会取新旧速度值中间的一个值。这样可以防止成员变量过快或者的加速或者减速，更重要的是它们的运动方向改变会更加自然。

执行完整个计算着色器的代码后更新阶段就已经完成了，剩下的就是渲染工作了。渲染程序使用了更新阶段计算着色器得到的每个成员的位置和速度数据作为一个顶点数组的属性，从而更新单个模型的最新位置，顶点着色器部分代码如下。

#version 430 core
layout (location = 0) in vec3 position;
layout (location = 1) in vec3 normal;
layout (location = 2) in vec3 bird_position;
layout (location = 3) in vec3 bird_velocity;

out VS_OUT {
    flat vec3 color; 
} vs_out;
uniform mat4 mvp;

在上面的着色器中，位置和法向量输入变量是一般的顶点属性，它们从几何缓存中读取，该缓存包含了单个纸飞机的模型数据。输入变量bird_position和bird_velocity是实例顶点属性，也就是对应一个模型内所有的顶点，这两个变量的值相同，其值由计算着色器运行得到。实例顶点属性通过函数glVertexAttribDivisor()向顶点着色器中传入数据。顶点着色器的主体代码如下，其中使用了集群成员的速度属性计算处理一个仿射矩阵，使纸飞机的机头总是朝向屏幕深处。

mat4 make_lookat(vec3 forward, vec3 up) {
    vec3 side = cross(forward, up);
    vec3 u_frame = cross(side, forward);
    // 原书中的光照矩阵计算算法有待商榷，x轴单位向量的计算方式应该为cross(up, forward);
    return mat4(vec4(side, 0.0), 
                vec4(u_frame, 0.0), 
                vec4(forward, 0.0),
                vec4(0.0, 0.0, 0.0, 1.0));
}

vec3 choose_color(float f) {
    float R = sin(f * 6.2831853);
    float G = sin((f + 0.3333) * 6.2831853); 
    float B = sin((f + 0.6666) * 6.2831853);
    return vec3(R, G, B) * 0.25 + vec3(0.75); 
}
 
void main(void) {
    // 原书中的光照矩阵计算算法有待商榷，x轴单位向量的计算方式应该为
    mat4 lookat = make_lookat(normalize(bird_velocity), vec3(0.0, 1.0, 0.0)); 

    vec4 obj_coord = lookat * vec4(position.xyz, 1.0);
    gl_Position = mvp * (obj_coord + vec4(bird_position, 0.0)); 

    vec3 N = mat3(lookat) * normal;
    vec3 C = choose_color(fract(float(gl_InstanceID / float(1237.0))));     

    vs_out.color = mix(C * 0.2, C, smoothstep(0.0, 0.8, abs(N).z));
}

这里构建lookat矩阵的原理和前面章节3D图像学数学基础中描述的lookat矩阵类似，也就是说可以理解为我们构建了一个模型-视口矩阵(model-view matrix)。如下图左侧，默认的飞机模型加载出来后机头正对着屏幕外部，机翼平面和Z轴平行，如果我们不指定一个lookat矩阵对模型进行仿射变换，那么我们看到的总是昨天那样单元的线条。下图右侧演示了计算出的“视图坐标系”和世界坐标系的关系(其中蓝色的为世界坐标系，红色的为“视图”坐标系)。需要注意的是这里的视图坐标系和前面章节3D图像学数学基础中不太一样，这里对X轴做了水平镜像处理，这样便能只改变纸飞机机头的方向而不改变其机翼的左右关系。

接下来根据当前顶点着色器的实例调用索引计算出每个模型的颜色值，再使用构建的lookat矩阵对模型顶点的法向量进行了变换，这样我们就能够正确的模拟光照效果。将最后得到的颜色值写入到顶点着色器的颜色输出变量中。片段着色器十分简单，直接将得到的颜色输出即可，这里就不再展示其源码。最后运行程序得到的渲染效果如下。

该示例程序的一个优化点事通过计算着色器来构建lookat矩阵，在上面的例子中我们在顶点着色器中计算该矩阵，这样会有重复的计算工作。这里因为每个模型的顶点数并不多，因此不会有太多的性能问题，但是对于单个模型网格数大的情况最好在计算着色器中构建这个矩阵，并将其作为额外的顶点实例属性，这样能极大的提升程序效率。我们也能在计算着色器中使用更多的规则来模拟真实的物理效果，而不仅仅是例子中那几条简单的规则。例如，我们可以模拟重力，使得向下飞比向上飞更容易，或者我们可以允许飞机碰撞或者弹开其他的模型。在本文中我们关注的是计算着色器的高并发效率的特性，而不是一些复杂的物理学原理，因此不再继续展开，但是需要记住，计算着色器是实现你想法的一个不错的工具。

4 总结

在本文中，我们深入了解了计算着色器，它可以构成独立于OpenGL图像渲染管道之外的计算管道，它可以是的我们在图像渲染之外的领域内充分利用现代图像处理单元(GPU)的强大并发计算能力。

首先介绍了计算着色器的执行模式，对工作组，同步以及本地工作组内的通信都做了详细说明。接下来使用了两个例子详细的讲解了计算着色器程序如何使用。第一个例子侧重图像处理，利用计算着色器实现光线相机的景深效果。第二个例子侧重科学计算，以及如何利用计算着色器实现物理仿真，实现集群算法(flocking algorithm)。

其实巧妙的运用几何着色器我们能完成很多有趣的工作，如人工智能，预处理和后处理，甚至是音频程序。