Separable Gaussian Blur in GLSL Compute Shaders

A separable Gaussian blur optimizes image processing by breaking a two-dimensional convolution kernel into two sequential one-dimensional passes: a horizontal pass and a vertical pass. This mathematical separation reduces computational complexity from \(\mathcal{O}(K^2)\) to \(\mathcal{O}(2K)\) operations per pixel, where \(K\) is the kernel diameter. When implemented in OpenGL Shading Language (GLSL) compute shaders, these passes are executed sequentially using a ping-pong texture approach, leveraging compute workgroups and thread-group shared memory to cache texels and drastically reduce texture bandwidth consumption.

The Principle of Separability

The 2D Gaussian function is defined as:

\[G(x, y) = \frac{1}{2\pi\sigma^2} e^{-\frac{x^2 + y^2}{2\sigma^2}}\]

Because the exponential term decomposes into independent factors (\(e^{-\frac{x^2 + y^2}{2\sigma^2}} = e^{-\frac{x^2}{2\sigma^2}} \cdot e^{-\frac{y^2}{2\sigma^2}}\)), the operation can be factored into two 1D operations:

\[G(x, y) = G(x) \times G(y)\]

Instead of gathering an \(11 \times 11\) grid (121 texture samples per pixel), a separable filter performs an 11-tap horizontal blur followed by an 11-tap vertical blur, totaling only 22 texture samples per pixel.

Architecture of Two-Pass Execution

Executing the filter in GLSL compute pipelines requires two distinct dispatches managed by a ping-pong buffer pattern:

  1. Horizontal Pass: Reads the original source texture via an image sampler or image2D binding, applies the 1D Gaussian kernel along the X-axis, and writes the output to an intermediate texture.
  2. Pipeline Barrier: The host application invokes an execution barrier (e.g., glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT)) to ensure all writes to the intermediate texture are finalized before the next pass reads them.
  3. Vertical Pass: Reads from the intermediate texture, applies the identical 1D Gaussian weights along the Y-axis, and writes the final blurred result to the destination texture.

Workgroup Shared Memory Optimization

A naive compute shader implementation fetches neighboring texels individually per thread. Because neighboring pixels share most of the same samples within the filter kernel, GLSL compute shaders can preload image data into fast shared memory (L1/tile memory on the GPU).

Horizontal Compute Pass

In the horizontal pass, a 1D or 2D compute workgroup caches the target pixels plus the apron (padding equal to the filter radius) along the X-axis.

#version 450 core

layout (local_size_x = 128, local_size_y = 1) in;

layout (binding = 0, rgba8) uniform readonly image2D inputImage;
layout (binding = 1, rgba8) uniform writeonly image2D outputImage;

const int RADIUS = 5;
const int TILE_SIZE = 128 + 2 * RADIUS;

shared vec4 sharedTile[TILE_SIZE];

// Precomputed 1D Gaussian weights (sum = 1.0)
const float weights[6] = float[](0.227027, 0.1945946, 0.1216216, 0.054054, 0.016216, 0.003126);

void main() {
    ivec2 coord = ivec2(gl_GlobalInvocationID.xy);
    int localX = int(gl_LocalInvocationID.x);

    // Populate interior pixels
    sharedTile[localX + RADIUS] = imageLoad(inputImage, coord);

    // Populate the left and right aprons
    if (localX < RADIUS) {
        sharedTile[localX] = imageLoad(inputImage, coord - ivec2(RADIUS, 0));
        sharedTile[localX + 128 + RADIUS] = imageLoad(inputImage, coord + ivec2(128, 0));
    }

    barrier(); // Wait until shared memory is fully loaded

    // Apply the 1D horizontal filter
    vec4 result = sharedTile[localX + RADIUS] * weights[0];
    for (int i = 1; i <= RADIUS; ++i) {
        result += sharedTile[localX + RADIUS + i] * weights[i];
        result += sharedTile[localX + RADIUS - i] * weights[i];
    }

    imageStore(outputImage, coord, result);
}

Vertical Compute Pass

The vertical pass uses identical weights, but rotates the tile caching and sampling logic across the Y-axis.

#version 450 core

layout (local_size_x = 1, local_size_y = 128) in;

layout (binding = 0, rgba8) uniform readonly image2D intermediateImage;
layout (binding = 1, rgba8) uniform writeonly image2D finalImage;

const int RADIUS = 5;
const int TILE_SIZE = 128 + 2 * RADIUS;

shared vec4 sharedTile[TILE_SIZE];

const float weights[6] = float[](0.227027, 0.1945946, 0.1216216, 0.054054, 0.016216, 0.003126);

void main() {
    ivec2 coord = ivec2(gl_GlobalInvocationID.xy);
    int localY = int(gl_LocalInvocationID.y);

    sharedTile[localY + RADIUS] = imageLoad(intermediateImage, coord);

    if (localY < RADIUS) {
        sharedTile[localY] = imageLoad(intermediateImage, coord - ivec2(0, RADIUS));
        sharedTile[localY + 128 + RADIUS] = imageLoad(intermediateImage, coord + ivec2(0, 128));
    }

    barrier();

    vec4 result = sharedTile[localY + RADIUS] * weights[0];
    for (int i = 1; i <= RADIUS; ++i) {
        result += sharedTile[localY + RADIUS + i] * weights[i];
        result += sharedTile[localY + RADIUS - i] * weights[i];
    }

    imageStore(finalImage, coord, result);
}

Host Execution Flow

On the CPU side, executing the separable blur involves:

  1. Binding the First Pass: Bind the source texture to unit 0 and the intermediate texture to unit 1.
  2. Dispatching Horizontal Compute: Call glDispatchCompute(imageWidth / 128, imageHeight, 1).
  3. Inserting a Barrier: Invoke glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT) so subsequent texture fetch operations wait until all image write operations complete.
  4. Binding the Second Pass: Bind the intermediate texture to unit 0 and the target display/framebuffer texture to unit 1.
  5. Dispatching Vertical Compute: Call glDispatchCompute(imageWidth, imageHeight / 128, 1).
  6. Final Barrier: Call glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT | GL_TEXTURE_FETCH_BARRIER_BIT) if the target texture will be consumed immediately by downstream rendering pipelines.