Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Correctly implement a 2 pass Gaussian blur

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:

input image

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

Fail blurred image

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.

Incorrect blur

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.

original image

Double banded blur

Any assistant would be great.

like image 439
James South Avatar asked Nov 06 '15 14:11

James South


People also ask

How do you implement a Gaussian blur?

TLDR: A Gaussian blur is applied by convolving the image with a Gaussian function. In English, this means that we'll take the Gaussian function and we'll generate an n x m matrix. Using this matrix and the height of the Gaussian distribution at that pixel location, we'll compute new RGB values for the blurred image.

What is sigmaX and sigmaY in Gaussian blur?

sigmax - standard deviation in X direction; if 0, calculated from kernel size. sigmay - standard deviation in Y direction; if sigmaY is None, sigmaY is taken to equal sigmaX.

Is Gaussian blur a low-pass filter?

Gaussian blur is a low-pass filter, attenuating high frequency signals.

Which tool is used for Gaussian blur window?

Now, before you apply the Gaussian Blur, we suggest you use the Refine Edges tool. This way, you will ensure that any sharp edges left due to the blur effect do not spoil your image. Once you are satisfied, click OK, and start applying the blur.

How to blur images using Gaussian blurring?

Gaussian Blurring makes use of a function called Gaussian Blur () function to reduce the clarity of images or to make the images distinct or to remove the noise from the images or to reduce the details from the images. The image that is to be blurred is read using imread () function.

Is Gaussian blur a low pass filter?

Gaussian blur is a low-pass filter, attenuating high frequency signals. [3] Its amplitude Bode plot (the log scale in the frequency domain) is a parabola . How much does a Gaussian filter with standard deviation smooth the picture? In other words, how much does it reduce the standard deviation of pixel values in the picture?

Why do the values in a Gaussian kernel add up to one?

Usually, the values in the kernel add up to one. This is to make sure no energy is added or removed from the image after the operation. Specifically, a Gaussian kernel (used for Gaussian blur) is a square array of pixels where the pixel values correspond to the values of a Gaussian curve (in 2D).

What is Gaussian smoothing in image processing?

Gaussian smoothing is also used as a pre-processing stage in computer vision algorithms in order to enhance image structures at different scales—see scale space representation and scale space implementation. Mathematically, applying a Gaussian blur to an image is the same as convolving the image with a Gaussian function.


1 Answers

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.

blurred china image

blurred balls

like image 175
James South Avatar answered Sep 19 '22 13:09

James South