Correctly implement a 2 pass Gaussian blur

Question

I'm attempting to implement a performant Gaussian blur taking advantage of the fact that a Gaussian kernel is separable, i. e. you can express a 2D convolution as a combination of two 1D convolutions.

I'm able to generate two kernels which I believe are correct using the following code.

/// <summary>
/// Create a 1 dimensional Gaussian kernel using the Gaussian G(x) function
/// </summary>
/// <param name="horizontal">Whether to calculate a horizontal kernel.</param>
/// <returns>The <see cref="T:float[,]"/></returns>
private float[,] CreateGaussianKernel(bool horizontal)
{
    int size = this.kernelSize;
    float[,] kernel = horizontal ? new float[1, size] : new float[size, 1];
    float sum = 0.0f;

    float midpoint = (size - 1) / 2f;
    for (int i = 0; i < size; i++)
    {
        float x = i - midpoint;
        float gx = this.Gaussian(x);
        sum += gx;
        if (horizontal)
        {
            kernel[0, i] = gx;
        }
        else
        {
            kernel[i, 0] = gx;
        }
    }

    // Normalise kernel so that the sum of all weights equals 1
    if (horizontal)
    {
        for (int i = 0; i < size; i++)
        {
            kernel[0, i] = kernel[0, i] / sum;
        }
    }
    else
    {
        for (int i = 0; i < size; i++)
        {
            kernel[i, 0] = kernel[i, 0] / sum;
        }
    }

    return kernel;
}

/// <summary>
/// Implementation of 1D Gaussian G(x) function
/// </summary>
/// <param name="x">The x provided to G(x)</param>
/// <returns>The Gaussian G(x)</returns>
private float Gaussian(float x)
{
    const float Numerator = 1.0f;
    float deviation = this.sigma;
    float denominator = (float)(Math.Sqrt(2 * Math.PI) * deviation);

    float exponentNumerator = -x * x;
    float exponentDenominator = (float)(2 * Math.Pow(deviation, 2));

    float left = Numerator / denominator;
    float right = (float)Math.Exp(exponentNumerator / exponentDenominator);

    return left * right;
}

The kernel size is calculated from the sigma as follows.

this.kernelSize = ((int)Math.Ceiling(sigma) * 2) + 1;
this.sigma = sigma;

Given a sigma of 3 This yields the following in each direction.

0.106288522,
0.140321344,
0.165770069,
0.175240144,
0.165770069,
0.140321344,
0.106288522

Which sums up to 1 perfectly so I'm in the right path.

Applying the kernels against an image is proving to be difficult as something is going wrong though I'm not sure what.

I'm using the following code which to run a 2-pass algorithm which works perfectly when convolving kernels that are even in both directions like Sobel or Prewitt for edge detection.

protected override void Apply(ImageBase target, 
ImageBase source, 
Rectangle targetRectangle, 
Rectangle sourceRectangle, 
int startY, 
int endY)
{
    float[,] kernelX = this.KernelX;
    float[,] kernelY = this.KernelY;
    int kernelYHeight = kernelY.GetLength(0);
    int kernelYWidth = kernelY.GetLength(1);
    int kernelXHeight = kernelX.GetLength(0);
    int kernelXWidth = kernelX.GetLength(1);
    int radiusY = kernelYHeight >> 1;
    int radiusX = kernelXWidth >> 1;

    int sourceY = sourceRectangle.Y;
    int sourceBottom = sourceRectangle.Bottom;
    int startX = sourceRectangle.X;
    int endX = sourceRectangle.Right;
    int maxY = sourceBottom - 1;
    int maxX = endX - 1;

    Parallel.For(
        startY,
        endY,
        y =>
        {
            if (y >= sourceY && y < sourceBottom)
            {
                for (int x = startX; x < endX; x++)
                {
                    float rX = 0;
                    float gX = 0;
                    float bX = 0;
                    float rY = 0;
                    float gY = 0;
                    float bY = 0;

                    // Apply each matrix multiplier to the
                    // color components for each pixel.
                    for (int fy = 0; fy < kernelYHeight; fy++)
                    {
                        int fyr = fy - radiusY;
                        int offsetY = y + fyr;

                        offsetY = offsetY.Clamp(0, maxY);

                        for (int fx = 0; fx < kernelXWidth; fx++)
                        {
                            int fxr = fx - radiusX;
                            int offsetX = x + fxr;

                            offsetX = offsetX.Clamp(0, maxX);

                            Color currentColor = source[offsetX, offsetY];
                            float r = currentColor.R;
                            float g = currentColor.G;
                            float b = currentColor.B;

                            if (fy < kernelXHeight)
                            {
                                rX += kernelX[fy, fx] * r;
                                gX += kernelX[fy, fx] * g;
                                bX += kernelX[fy, fx] * b;
                            }

                            if (fx < kernelYWidth)
                            {
                                rY += kernelY[fy, fx] * r;
                                gY += kernelY[fy, fx] * g;
                                bY += kernelY[fy, fx] * b;
                            }
                        }
                    }

                    float red = (float)Math.Sqrt((rX * rX) + (rY * rY));
                    float green = (float)Math.Sqrt((gX * gX) + (gY * gY));
                    float blue = (float)Math.Sqrt((bX * bX) + (bY * bY));

                    Color targetColor = target[x, y];
                    target[x, y] = new Color(red, 
                                             green, blue, targetColor.A);
                }
            }
        });
}

Here's the input image:

And here's the image with an attempted blur using 3 sigma.

As you can see something is amiss. It's like I'm sampling from the wrong point or something.

Any ideas? I appreciate that this is a long winded question.

Update

So on the suggestion of Nico Schertler I have split the algorithm into two passes as follows:

protected override void Apply(
    ImageBase target,
    ImageBase source,
    Rectangle targetRectangle,
    Rectangle sourceRectangle,
    int startY,
    int endY)
{
    float[,] kernelX = this.KernelX;
    float[,] kernelY = this.KernelY;
    int kernelXHeight = kernelX.GetLength(0);
    int kernelXWidth = kernelX.GetLength(1);
    int kernelYHeight = kernelY.GetLength(0);
    int kernelYWidth = kernelY.GetLength(1);
    int radiusXy = kernelXHeight >> 1;
    int radiusXx = kernelXWidth >> 1;
    int radiusYy = kernelYHeight >> 1;
    int radiusYx = kernelYWidth >> 1;

    int sourceY = sourceRectangle.Y;
    int sourceBottom = sourceRectangle.Bottom;
    int startX = sourceRectangle.X;
    int endX = sourceRectangle.Right;
    int maxY = sourceBottom - 1;
    int maxX = endX - 1;

    // Horizontal blur
    Parallel.For(
        startY,
        endY,
        y =>
        {
            if (y >= sourceY && y < sourceBottom)
            {
                for (int x = startX; x < endX; x++)
                {
                    float rX = 0;
                    float gX = 0;
                    float bX = 0;

                    // Apply each matrix multiplier to the color 
                    // components for each pixel.
                    for (int fy = 0; fy < kernelXHeight; fy++)
                    {
                        int fyr = fy - radiusXy;
                        int offsetY = y + fyr;

                        offsetY = offsetY.Clamp(0, maxY);

                        for (int fx = 0; fx < kernelXWidth; fx++)
                        {
                            int fxr = fx - radiusXx;
                            int offsetX = x + fxr;

                            offsetX = offsetX.Clamp(0, maxX);

                            Color currentColor = source[offsetX, offsetY];
                            float r = currentColor.R;
                            float g = currentColor.G;
                            float b = currentColor.B;

                            rX += kernelX[fy, fx] * r;
                            gX += kernelX[fy, fx] * g;
                            bX += kernelX[fy, fx] * b;
                        }
                    }

                    float red = rX;
                    float green = gX;
                    float blue = bX;

                    Color targetColor = target[x, y];
                    target[x, y] = new Color(red, green, blue, targetColor.A);
                }
            }
        });

    // Vertical blur
    Parallel.For(
        startY,
        endY,
        y =>
        {
            if (y >= sourceY && y < sourceBottom)
            {
                for (int x = startX; x < endX; x++)
                {
                    float rY = 0;
                    float gY = 0;
                    float bY = 0;

                    // Apply each matrix multiplier to the 
                    // color components for each pixel.
                    for (int fy = 0; fy < kernelYHeight; fy++)
                    {
                        int fyr = fy - radiusYy;
                        int offsetY = y + fyr;

                        offsetY = offsetY.Clamp(0, maxY);

                        for (int fx = 0; fx < kernelYWidth; fx++)
                        {
                            int fxr = fx - radiusYx;
                            int offsetX = x + fxr;

                            offsetX = offsetX.Clamp(0, maxX);

                            Color currentColor = source[offsetX, offsetY];
                            float r = currentColor.R;
                            float g = currentColor.G;
                            float b = currentColor.B;

                            rY += kernelY[fy, fx] * r;
                            gY += kernelY[fy, fx] * g;
                            bY += kernelY[fy, fx] * b;
                        }
                    }

                    float red = rY;
                    float green = gY;
                    float blue = bY;

                    Color targetColor = target[x, y];
                    target[x, y] = new Color(red, green, blue, targetColor.A);
                }
            }
        });
}

I'm getting closer to my goal as there is now a blurring effect. Unfortunately it is not correct.

If you look closely you will see that there is double banding in the vertical direction and not enough blur in the horizontal. The following images demonstrate that clearly when I bump up the sigma to 10.

Any assistant would be great.

You can't separate the convolutions like this. You have to run two separate passes. The second convolution should use the results of the first one. — Nico Schertler, Nov 06 '15 at 15:17
Basically, call your method twice. In the first call, omit the vertical convolution. In the second call, omit the horizontal convolution. — Nico Schertler, Nov 06 '15 at 15:56
You will want to keep this artisctic cross-hatch filter as well ;-) — TaW, Nov 06 '15 at 16:29
@TaW I'll definitely be playing about with that kind of stuff once I've got the guts of my framework running. — James South, Nov 09 '15 at 04:23

James South · Accepted Answer · 2015-11-07T09:49:54.260

Ok so with that last update I was being a little silly and not creating an interim image to convolve against for the second pass.

Here's a complete working sample of the convolution algorithm. The original gaussian kernel creation code was correct:

/// <inheritdoc/>
protected override void Apply(
    ImageBase target,
    ImageBase source,
    Rectangle targetRectangle,
    Rectangle sourceRectangle,
    int startY,
    int endY)
{
    float[,] kernelX = this.KernelX;
    float[,] kernelY = this.KernelY;

    ImageBase firstPass = new Image(source.Width, source.Height);
    this.ApplyConvolution(firstPass, source, sourceRectangle, startY, endY, kernelX);
    this.ApplyConvolution(target, firstPass, sourceRectangle, startY, endY, kernelY);
}

/// <summary>
/// Applies the process to the specified portion of the specified <see cref="ImageBase"/> at the specified location
/// and with the specified size.
/// </summary>
/// <param name="target">Target image to apply the process to.</param>
/// <param name="source">The source image. Cannot be null.</param>
/// <param name="sourceRectangle">
/// The <see cref="Rectangle"/> structure that specifies the portion of the image object to draw.
/// </param>
/// <param name="startY">The index of the row within the source image to start processing.</param>
/// <param name="endY">The index of the row within the source image to end processing.</param>
/// <param name="kernel">The kernel operator.</param>
private void ApplyConvolution(
    ImageBase target,
    ImageBase source,
    Rectangle sourceRectangle,
    int startY,
    int endY,
    float[,] kernel)
{
    int kernelHeight = kernel.GetLength(0);
    int kernelWidth = kernel.GetLength(1);
    int radiusY = kernelHeight >> 1;
    int radiusX = kernelWidth >> 1;

    int sourceY = sourceRectangle.Y;
    int sourceBottom = sourceRectangle.Bottom;
    int startX = sourceRectangle.X;
    int endX = sourceRectangle.Right;
    int maxY = sourceBottom - 1;
    int maxX = endX - 1;

    Parallel.For(
        startY,
        endY,
        y =>
        {
            if (y >= sourceY && y < sourceBottom)
            {
                for (int x = startX; x < endX; x++)
                {
                    float red = 0;
                    float green = 0;
                    float blue = 0;
                    float alpha = 0;

                    // Apply each matrix multiplier to the color components for each pixel.
                    for (int fy = 0; fy < kernelHeight; fy++)
                    {
                        int fyr = fy - radiusY;
                        int offsetY = y + fyr;

                        offsetY = offsetY.Clamp(0, maxY);

                        for (int fx = 0; fx < kernelWidth; fx++)
                        {
                            int fxr = fx - radiusX;
                            int offsetX = x + fxr;

                            offsetX = offsetX.Clamp(0, maxX);

                            Color currentColor = source[offsetX, offsetY];

                            red += kernel[fy, fx] * currentColor.R;
                            green += kernel[fy, fx] * currentColor.G;
                            blue += kernel[fy, fx] * currentColor.B;
                            alpha += kernel[fy, fx] * currentColor.A;
                        }
                    }

                    target[x, y] = new Color(red, green, blue, alpha);
                }
            }
        });
}

Here's the output of the code with 10 sigma.

Correctly implement a 2 pass Gaussian blur

Update

1 Answers1