c# - 正确实现 2 pass 高斯模糊

Question

我正在尝试利用高斯核是可分离的这一事实来实现高性能的高斯模糊，即。 e.您可以将二维卷积表示为两个一维卷积的组合。

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

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

给定一个 3 的 sigma，这会在每个方向产生以下结果。

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

完美地加起来为 1，所以我走在正确的道路上。

将内核应用于图像被证明是困难的，因为出了点问题，虽然我不确定是什么。

我正在使用以下代码来运行 2-pass 算法，该算法在对甚至在两个方向上的卷积核(如用于边缘检测的 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);
                }
            }
        });
}

我离目标越来越近了，因为现在出现了模糊效果。不幸的是，它不正确。

如果仔细观察，您会发现垂直方向有双 strip ，水平方向没有足够的模糊。当我将 sigma 提高到 10 时，下图清楚地表明了这一点。

任何助手都会很棒。

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 西格玛的代码的输出。