正确实现双通道高斯模糊处理

Question

我正在尝试实现一个高性能的高斯模糊，利用了高斯核是可分离的事实，即可以将2D卷积表示为两个1D卷积的组合。

我使用以下代码能够生成我认为是正确的两个核。

/// <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;
}

核大小是根据 sigma 计算的，具体如下。

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

给定sigma为3，这会导致每个方向如下。

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

这些加起来完全等于1，所以我走在正确的道路上。

将卷积核应用于图像正在变得困难，因为出现了一些问题，尽管我不确定是什么问题。

我正在使用以下代码运行一个双通道算法，当卷积内核在两个方向上都是偶数时，如Sobel或Prewitt用于边缘检测时，它可以完美地工作。

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);
                }
            }
        });
}

这是输入图像：

这是一张尝试使用3个标准差模糊的图片。

正如您所看到的，有些地方不对劲。就像我从错误的点进行采样一样。

有什么想法吗？我知道这是一个冗长的问题，感谢您的帮助。

---

更新

根据Nico Schertler的建议，我将算法拆分为两个步骤，如下：

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);
                }
            }
        });
}

由于出现了模糊效果，我离我的目标越来越近。不幸的是，这并不正确。

如果你仔细观察，就会发现垂直方向有双重条纹，而水平方向的模糊不足。当我将sigma增加到10时，以下图像清楚地展示了这一点。

任何助手都很好。

- James South

---

你不能这样分离卷积。你必须运行两个单独的传递。第二个卷积应该使用第一个卷积的结果。 - Nico Schertler

你能否提供一个代码示例来进一步说明吗？ - James South

基本上，调用你的方法两次。在第一次调用中，省略垂直卷积。在第二次调用中，省略水平卷积。 - Nico Schertler

您也会想要保留这个艺术性的交叉阴影滤镜;-) - TaW

@TaW 一旦我把框架的核心部分运行起来，我肯定会尝试着玩一下那种东西。 - James South

Answer 1

好的，所以在上一次更新中，我有点儿傻，没有创建一个过渡图像来进行第二次卷积。

这是卷积算法的完整工作示例。原始的高斯核创建代码是正确的：

/// <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);
                }
            }
        });
}

这是10 sigma的代码输出。

- James South