Correctly implement 2-way Gaussian blur

Question

I am trying to implement strong Gaussian blur using the fact that the Gaussian core is separable, i.e. E. You can express a two-dimensional convolution as a combination of two one-dimensional convolutions.

I can generate two kernels that 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 sigma as follows.

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

Given sigma 3, we get the following in each direction.

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

Which sums up to 1, so I'm on the right track.

Applying kernels to an image is difficult because something goes wrong, although I'm not sure what.

I use the following code to run a two-pass algorithm that works great when collapsing kernels that even in both directions, such as Sobel or Prewitt to detect the edge.

 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 is an image with an attempt to blur using 3 sigma.

As you can see, something is wrong. It looks like I'm taking a sample from the wrong point or something else.

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

---

Update

So, at the suggestion of Nico Shertler, I divided 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 blur effect. Unfortunately, this is not true.

If you look closely, you will see that in the vertical direction there is a double band and not enough blur in the horizontal plane. The following images demonstrate this clearly when I increase sigma to 10.

Any assistant would be great.

Answer 1

Ok, so with this latest update, I was a little dumb and didn't create an intermediate image to curl up against the second pass.

Here is a complete working example of a convolution algorithm. The source code for creating a Gaussian kernel 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 code output with 10 sigma.