再谈快速的高斯模糊算法（使用多次均值滤波逼近和扩展的二项式滤波滤波器）及其优化。

关于高斯模糊，我在我早期的博客里也有两篇文章予以描述：

　　一个是递归的IIR滤波器，一个Deriche滤波器，他们的速度都已经是顶级的了，而且都能够使用SIMD指令优化，其中有讲到《Recursive implementation of the Gaussian filter》这个方法在半径较大的时候会出现一定的瑕疵，核心原因是大半径会导致其中的某些系数特别小，因此造成浮点精度的丢失，因此，要保证效果就必须在计算过程中使用double数据类型，而使用了double，普通的sse指令集的增速效果就不是很明显了，因此，为了速度可能需要使用AVX或者更高的AVX512。

　　那么这两个都是从离散角度来说比较精确的算法，因为有了SIMD指令，使得他们在PC上即使有了大量的浮点计算，计算的速度也比较不错。在一些特殊的情况或者特殊的硬件中，还是存在浮点计算非常慢的情况（比如FPGA），因此，还是有不上资料和文章提出了一些近似的算法来减少计算量，在celerychen分享的文章快速高斯滤镜算法中，除了上述两个算法外，还提到了均值滤波逼近高斯滤波以及扩展二项式滤波逼近高斯滤波两个方法。

一、Binomial Filter 二项式滤波滤波器

多年前我也看过这个文章，那个时候也没有怎么在意，最近在研究halcon的一些滤波器时，偶尔翻到其binomial_filter函数的说明时，突然看到其有如下的一些描述：

The binomial filter is a very good approximation of a Gaussian filter that can be implemented extremely efficiently using only integer operations. Hence, binomial_filter is very fast.

　他说这个binomial filter只有整形的计算，因此速度非常快。

是不是这样的呢，于是我在网络上搜索这方面的资料，大部分都指向 extended binomial filter for Fast gaussian Blur这篇文章，可以参考道客巴巴里的这篇文章：https://www.doc88.com/p-1896883955048.html

　　这个文章写得很详细，有着非常丰富的数学推导和公式，还有很多漂亮的曲线，比如下面这个描述了不同的半径r和不同的n阶毕竟的精度。当n=3或者n=4时，也已经相当的准确了。特别是半径比价大时（模糊程度）。

　　这个文章的末尾提供了相关的参考代码，但是那个参考代码有点凌乱，而且不知道是啥语言，好像是可以为PS写插件的脚本一样，不过大概还是能猜到是什么意思。

celerychen在他的文章说里面几个重要的初始参数不知道如何确定，主要可能是指的scaleFactor，这个其实因该他做界面预览时的一些缩放值，我们做算法时，他就是1。

另外，那个文章是比较早的了，作者提供的那个网站里有提供最新的一批代码，详细见：Efficient Gaussian Blur Algorithm

这个文章中提供的一个重要的计算公式即为：

　　其中为用户输入的标准差，n为离散的阶数，r为计算得到半径，由上式我们可以推到出：

　　　　　　　　　　　　　　　　r = sqrtf(12.0 * Sigma * Sigma / n + 1)

　　注意，这里得到的r应该是浮点数，但是浮点数，我们无法进行均值模糊，所以一般需要四舍五入取整。当然，如果要求精度的，那就要去上下两个半径值分别做处理后，在对结果进行插值。

　　这个公式在均值滤波逼近高斯滤波的文章里也有提到。

我们从Efficient Gaussian Blur Algorithm 这个网站的相关链接里可以有如下描述：

gauss2.ffp: The program for the first degree approximation is a great improvement to the box blur. It is as easy to understand and has only slight approximation error.
gauss3.ffp: The program for the second degree approximation is accurate enough for most applications.
gauss4.ffp: The third degree approximation has nearly no difference to the Gaussian blur for 8-bit images apart of rounding errors

我想说的是gauss2.ffp其实是Third degree approximation对应的代码（n = 3）,gauss3.ffp和gauss4.ffp都是Fourth degree approximation对应的代码(n=4)，但是他们里面的叠加顺序写法又有所不同（里面的变量sum, deri1, deri2, diff叠加的顺序），那么真正的Second degree approximation文章里没有给链接，但是可以用在链接 http://members.chello.at/easyfilter/gauss1.ffp 里打开。感觉像是作者故意隐藏了一样。

另外还要注意在论文里的sum, deri1, deri2, diff和链接里的顺序也不太一样。

　　为了方便大家测试，我从链接里整理了一个可以运行的代码，很丑陋，但是确实效果和速度都还不错。

int IM_BinomialFilter4_PureC_Gray_Raw(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Sigma)
{
    int Channel = Stride / Width;
    if ((Src == NULL) || (Dest == NULL) || (Dest == NULL))        return IM_STATUS_NULLREFRENCE;
    if ((Width <= 0) || (Height <= 0))                            return IM_STATUS_INVALIDPARAMETER;
    if (Channel != 1)                                            return IM_STATUS_NOTSUPPORTED;
    int Status = IM_STATUS_OK;

    Sigma = IM_Min(Sigma, 256);            //    如果全部用整形的数据，半径不能大于30， Sum和Deri用浮点半径可以取得较大

    int Radius = (int)(sqrtf(3.0 * Sigma  * Sigma + 1) + 0.4999999f);
    int Weight = Radius * Radius * Radius * Radius;
    float InvWeight = 1.0f / Weight;
    unsigned char *Buffer = (unsigned char *)malloc((IM_Max(Width, Height) + 4 * Radius) * sizeof(unsigned char));
    unsigned char *Temp = (unsigned char *)malloc(Height * Stride * sizeof(unsigned char));
    int *RowOffset = (int *)malloc((Width + 4 * Radius) * sizeof(int));
    int *ColOffset = (int *)malloc((Height + 4 * Radius) * sizeof(int));

    if ((Buffer == NULL) || (Temp == NULL) || (RowOffset == NULL) || (ColOffset == NULL))
    {
        Status = IM_STATUS_OUTOFMEMORY;
        goto FreeMemory;
    }
    //    注意起点位置是偏移了一个1个像素的
    for (int X = 0; X < Width + 4 * Radius; X++)
        RowOffset[X] = IM_GetMirrorPos(Width, X - 2 * Radius + 1);

    for (int Y = 0; Y < Height + 4 * Radius; Y++)
        ColOffset[Y] = IM_GetMirrorPos(Height, Y - 2 * Radius + 1);

    for (int Y = 0; Y < Height; Y++)     // horizontal blur...
    {
        unsigned char *LinePS = Src + Y * Stride;
        unsigned char *LinePD = Temp + Y * Stride;
        for (int X = 0; X < 2 * Radius - 1; X++)
            Buffer[X] = LinePS[RowOffset[X]];
        memcpy(Buffer + 2 * Radius - 1, LinePS, Width);
        for (int X = 2 * Radius - 1 + Width; X < Width + 4 * Radius; X++)
            Buffer[X] = LinePS[RowOffset[X]];

        float Sum = 0, Deri1 = 0;
        int Deri2 = 0, Diff = 0;
        for (int X = -4 * Radius; X < 0; X++)
        {
            Sum += Deri1;
            Deri1 += Deri2;
            Deri2 += Diff;  // accumulate pixel blur
            Diff += Buffer[X + 4 * Radius];
            if (X + 3 * Radius >= 0) Diff -= 4 * Buffer[X + 3 * Radius];    //  -4,
            if (X + 2 * Radius >= 0) Diff += 6 * Buffer[X + 2 * Radius];    //      +6,
            if (X + Radius >= 0) Diff -= 4 * Buffer[X + 1 * Radius];        //          -4, (+1)}
        }
        for (int X = 0; X < Width; X++)
        {
            Sum += Deri1;
            Deri1 += Deri2;
            Deri2 += Diff;  
            Diff += Buffer[X] - 4 * Buffer[X + Radius] + 6 * Buffer[X + 2 * Radius] - 4 * Buffer[X + 3 * Radius] + Buffer[X + 4 * Radius];        //{+1,    -4,    +6,    -4,    +1}
            LinePD[X] = (int)(Sum * InvWeight + 0.499999f);            
            //LinePD[X] = (Sum + Weight / 2) / Weight;
        }
    }
    for (int X = 0; X < Width; X++)　 //　　Vertical Blur
    {
        unsigned char *LinePS = Temp + X;
        unsigned char *LinePD = Dest + X;
        for (int Y = 0; Y < 2 * Radius - 1; Y++)
            Buffer[Y] = LinePS[ColOffset[Y] * Stride];
        for (int Y = 0; Y < Height; Y++)
        {
            Buffer[Y + 2 * Radius - 1] = LinePS[Y * Stride];
        }
        for (int Y = 2 * Radius - 1 + Height; Y < Height + 4 * Radius; Y++)
            Buffer[Y] = LinePS[ColOffset[Y] * Stride];

        float Sum = 0, Deri1 = 0;
        int Deri2 = 0, Diff = 0;
        for (int Y = -4 * Radius; Y < 0; Y++)
        {
            Sum += Deri1;
            Deri1 += Deri2;
            Deri2 += Diff;  // accumulate pixel blur
            Diff += Buffer[Y + 4 * Radius];
            if (Y + 3 * Radius >= 0) Diff -= 4 * Buffer[Y + 3 * Radius];    //  -4,
            if (Y + 2 * Radius >= 0) Diff += 6 * Buffer[Y + 2 * Radius];    //      +6,
            if (Y + Radius >= 0) Diff -= 4 * Buffer[Y + 1 * Radius];        //          -4, (+1)}
        }
        for (int Y = 0; Y < Height; Y++)
        {
            Sum += Deri1;
            Deri1 += Deri2;
            Deri2 += Diff;
            Diff += Buffer[Y] - 4 * Buffer[Y + Radius] + 6 * Buffer[Y + 2 * Radius] - 4 * Buffer[Y + 3 * Radius] + Buffer[Y + 4 * Radius];        //{+1,    -4,    +6,    -4,    +1}
            LinePD[0] = (int)(Sum * InvWeight + 0.499999f);
            //LinePD[0] = (Sum + Weight / 2) / Weight;
            LinePD += Stride;
        }
    }
FreeMemory:

    if (Buffer != NULL)        free(Buffer);
    if (Temp != NULL)        free(Temp);
    if (RowOffset != NULL)    free(RowOffset);
    if (ColOffset != NULL)    free(ColOffset);
    return Status;
}

　　注意这个是翻译自gauss3.ffp的代码并做了适当的调整，我们看到程序还是使用了很多的浮点数的，否则当Sigma稍微大一点，结果就溢出了，当然，代码里也给出了如果SIgma不大于30，则可以全部用整形实现。

　　我们处理了一副3000*2000的灰度图，出乎意料的是，就这个很随意的代码能获得82毫秒的成绩。

　　优化：

乍一看起来，Vertical Blur那一部分的代码访问内存的方式很不友好，有着大量的跨行访问内存的地方，确实原始的代码是有这个问题，但是这种模糊都有这个特性，在高斯模糊算法的全面优化过程分享（一）中我们探讨过垂直方向处理算法一般不宜直接写，而应该用一个临时的行缓存进行处理，这里同样可以用这个方法来处理，而且处理的代码还能通用各种不同的通道，这里不多赘述。

Horizontal blur部分，看代码知道是存在着严重的像素前后依赖的，这种特性非常不利于指令集优化，处理的方式呢，我在博客里也多次提到，把多行（比如四行）合并到一行后，进行处理，这样合并的一行里就有连续4个像素的处理是毫无关系的，就可以借助于SSE之类的指令予以处理了，这部分方式可以参考我的博文： SSE图像算法优化系列五：超高速指数模糊算法的实现和优化（10000*10000在100ms左右实现）。

经过一番折腾呢，这个算法可以优化到13ms左右（4阶），当然，这个代码明显有个特性，随着Sigma的增大，程序里的边缘部分的计算量也会增加，因此，算法的耗时和参数有一定的关系，但是不是特别明显。

如果是使用的2阶算法，则内部计算要更加的少，同样的图片大约在11ms左右可以处理完成，而原始的Deriche算法嘛，耗时大概在18ms上下。

2阶的有个比较好的特性，就是他的全部的计算过程里的累加量都在整形的有效范围呢，因此确实可以用整形计算实现效果。

完成了这个算法后，回头在去看Halcon的binomial_filter ，感觉味道又变了，本来以为他的那个MaskWidth和MaskHeight就是对应的水平和垂直方向的Sigma的，结果测试不是，他的结果在MaskWidth和MaskHeight很大时，模糊的程度也比较小。目前还搞不清楚是怎么回事。

不过有一点，Halcon的这个算子的速度那真叫一个牛逼，3000*2000的图都只要几豪秒（随MaskWidth和MaskHegiht的增加耗时也在增加）。

　　二、使用多次均值模糊模拟高斯模糊

　　这个算法的参考文献的正式名字应该是Fast Almost-Gaussian Filtering，而不是celerychen文章里提的 Arbitrary Gaussian Filtering with 25 Addtions and 5 Multiplications per Pixel，那个文章里根本没看到说25次加法和5次乘法。

这个文章的理论基础其实也是这个公式：

英文的原文为：

　　他这里给出了用均值模糊模拟高斯模糊的半径的计算方式，即上图中公式5，当给定均方差，给定均值模拟的次数，就可以计算出m。

　　在https://blog.ivank.net/fastest-gaussian-blur.html这里作者给出了相关的简易的示例代码，稍作整理如下所示：

//    https://github.com/amazingyyc/fasted_gauss_blur
//    https://www.cnblogs.com/amazingyyc/p/4530634.html
//    https://www.doc88.com/p-1896883955048.html
//    使用 N 个 BoxBlur 拟合 Sigma的高斯函数
//    参考：Arbitrary Gaussian Filtering with 25 Addtions and 5 Multiplications per Pixel.pdf
//    Fast Almost - Gaussian Filtering.pdf
void IM_GetBoxRadiusForGauss(float Sigma, int *Size, int N)
{
    float IdealW = sqrt(12.0 * Sigma * Sigma / N + 1.0);            //    论文公式1
    int LowerW = floor(IdealW);                                        //    向下取整
    if (LowerW % 2 == 0)    LowerW--;                                //    the first odd integer that is less than IdealW
    int UpperW = LowerW + 2;                                        //    the other being being the first odd integer greater than IdealW
    int M = round((12.0 * Sigma * Sigma - N * LowerW * LowerW - 4 * N * LowerW - 3 * N) / (-4 * LowerW - 4));    //    论文公式2
    for (int Y = 0; Y < N; Y++)                                        //    Apply an averaging filter of width LowerW M times.
        Size[Y] = (Y < M ? LowerW : UpperW);                        //    Apply an averaging filter of width UpperW N - M times.
}

int IM_FastGaussBlur(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Radius)
{
    int Channel = Stride / Width;
    if ((Src == NULL) || (Dest == NULL) || (Dest == NULL))        return IM_STATUS_NULLREFRENCE;
    if ((Width <= 0) || (Height <= 0))                            return IM_STATUS_INVALIDPARAMETER;
    if ((Channel != 1) && (Channel != 3) && (Channel != 4))        return IM_STATUS_NOTSUPPORTED;
    int Status = IM_STATUS_OK;

    unsigned char *Temp = (unsigned char *)malloc(Height * Stride * sizeof(unsigned char));
    if (Temp == NULL)
    {
        Status = IM_STATUS_OUTOFMEMORY;
        goto FreeMemory;
    }
    int BoxRadius[3];                                //    使用3次均值模糊代替高斯
    IM_GetBoxRadiusForGauss(Radius, BoxRadius, 3);        //    计算拟合的半径
    Status = IM_BoxBlur(Src, Dest, Width, Height, Stride, (BoxRadius[0] - 1) / 2);
    if (Status != IM_STATUS_OK)    goto FreeMemory;
    Status = IM_BoxBlur(Dest, Temp, Width, Height, Stride, (BoxRadius[1] - 1) / 2);
    if (Status != IM_STATUS_OK)    goto FreeMemory;
    Status = IM_BoxBlur(Temp, Dest, Width, Height, Stride, (BoxRadius[2] - 1) / 2);
    if (Status != IM_STATUS_OK)    goto FreeMemory;
FreeMemory:
    if (Temp != NULL)    free(Temp);
    return Status;
}