高斯模糊算法的全面优化过程分享。

Posted 2020-09-01 只(挚)爱图像处理

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了高斯模糊算法的全面优化过程分享。相关的知识，希望对你有一定的参考价值。

在高斯模糊算法的全面优化过程分享（一）一文中我们已经给出了一种相当高性能的高斯模糊过程，但是优化没有终点，经过上一个星期的发愤图强和测试，对该算法的效率提升又有了一个新的高度，这里把优化过程中的一些心得和收获用文字的形式记录下来。

第一个尝试直接使用内联汇编替代intrinsics代码（无效）

我在某篇博客里看到说intrinsics语法虽然简化了SSE编程的难度，但是他无法直接控制XMM0-XMM7寄存器，很多指令中间都会用内存做中转，所以我就想我如果直接用汇编写效率肯定还能有进一步的提高，于是我首先尝试把GaussBlurFromLeftToRight_SSE优化，仔细观察这个函数，如果弄得好，确实能有效的利用这些寄存器，有关代码如下：

void GaussBlurFromLeftToRight_SSE(float *Data, int Width, int Height, float B0, float B1, float B2, float B3)
{
    float *MemB3 = (float *)_mm_malloc(4 * sizeof(float), 16);
    MemB3[0] = MemB3[1] = MemB3[2] = MemB3[3] = B3;
    int Stride = Width * 4 * sizeof(float);
    _asm
    {
        mov     ecx, Height
        movss   xmm0, B0
        shufps  xmm0, xmm0, 0            //    xmm0 = B0
        movss   xmm1, B1
        shufps  xmm1, xmm1, 0            //    xmm1 = B1
        movss   xmm2, B2
        shufps  xmm2, xmm2, 0            //    xmm2 = B2
        mov     edi, MemB3
    LoopH24 :
        mov     esi, ecx
        imul    esi, Stride
        add     esi, Data                //    LinePD = Data + Y * Width * 4
        mov     eax, Width
        movaps  xmm3, [esi]              //    xmm3 = V1
        movaps  xmm4, xmm3                //  xmm4 = V2 = V1
        movaps  xmm5, xmm3                //     xmm5 = V3 = V1
    LoopW24 :
    movaps  xmm6, [esi]                //    xmm6 = V0
        movaps  xmm7, xmm3                //    xmm7 = V1
        mulps   xmm5, [edi]                //    xmm5 = V3 * B3
        mulps   xmm7, xmm1                //    xmm7 = V1 * B1
        mulps   xmm6, xmm0                //    xmm6 = V0 * B0
        addps   xmm6, xmm7                //    xmm6 = V0 * B0 + V1 * B1
        movaps  xmm7, xmm4                //    xmm7 = V2
        mulps   xmm7, xmm2                //    xmm7 = V2 * B2
        addps   xmm5, xmm7                //    xmm5 = V3 * B3 + V2 * B2
        addps   xmm6, xmm5                //    xmm6 = V0 * B0 + V1 * B1 + V3 * B3 + V2 * B2
        movaps  xmm5, xmm4                //    V3 = V2            
        movaps  xmm4, xmm3                //    V2 = V1
        movaps [esi], xmm6
        movaps  xmm3, xmm6                //    V1 = V0
        add     esi, 16
        dec     eax
        jnz     LoopW24
        dec     ecx
        jnz     LoopH24
    }
    _mm_free(MemB3);
}

　　看上面的代码，基本上把XMM0-XMM7这8个寄存器都充分利用了，在我的预想中应该能有速度的提升的，可是一执行，真的好悲剧，和原先相比速度毫无变化，这是怎么回事呢。

后来我反编译intrinsics的相关代码，发现编译器真的很厉害，他的汇编代码和我上面的基本一致，只是寄存器的利用顺序有所不同而已，后面又看了其他的几个函数，发现编译器的汇编码都写的非常高效，基本上我们是超不过他了，而且编译器还能充分调整指令执行的顺序，使得有关指令还能实现指令层次的并行，而如果我们自己写ASM，这个对功力的要求就更高了，所以说网络上的说法也不可以完全相信，而如果不是有特别强的汇编能力，也不要去挑战编译器。

第二个尝试水平方向的模糊一次执行二行（15%提速）

这个尝试纯粹是随意而为，谁知道居然非常有效果，具体而言就是在GaussBlurFromLeftToRight_SSE和GaussBlurFromRightToLeft_SSE函数的Y循环内部，一次性处理二行代码，我们以LeftToRight为例，示意代码如下：

    __m128 CofB0 = _mm_set_ps(0, B0, B0, B0);
    __m128 CofB1 = _mm_set_ps(0, B1, B1, B1);
    __m128 CofB2 = _mm_set_ps(0, B2, B2, B2);
    __m128 CofB3 = _mm_set_ps(0, B3, B3, B3);
    __m128 V1 = _mm_load_ps(LineP1);                //    起点重复数据
    __m128 W1 = _mm_load_ps(LineP2);
    __m128 V2 = V1, V3 = V1;
    __m128 W2 = W1, W3 = W1;
    for (int X = 0; X < Length; X++, LineP1 += 4, LineP2 += 4)            
    {
        __m128 V0 = _mm_load_ps(LineP1);
        __m128 W0 = _mm_load_ps(LineP2);
        __m128 V01 = _mm_add_ps(_mm_mul_ps(CofB0, V0), _mm_mul_ps(CofB1, V1));
        __m128 W01 = _mm_add_ps(_mm_mul_ps(CofB0, W0), _mm_mul_ps(CofB1, W1));
        __m128 V23 = _mm_add_ps(_mm_mul_ps(CofB2, V2), _mm_mul_ps(CofB3, V3));
        __m128 W23 = _mm_add_ps(_mm_mul_ps(CofB2, W2), _mm_mul_ps(CofB3, W3));
        __m128 V = _mm_add_ps(V01, V23);
        __m128 W = _mm_add_ps(W01, W23);
        V3 = V2;    V2 = V1;    V1 = V;
        W3 = W2;    W2 = W1;    W1 = W;
        _mm_store_ps(LineP1, V);
        _mm_store_ps(LineP2, W);
    }

　　就是把原来的代码复制一份，在稍微调整一下，当然注意这个时候Y分量一次要递增2行了，还有如果Height是奇数，还要对最后一行做处理。这些活都是细活，稍微注意就不会出错了。

就这么样的简单的一个调整，经过测试性能居然能有15%的提升，真是意想不到，分析具体的原因，我觉得Y循环变量的计数耗时的减少在这里是无关紧要的，核心可能还是这些intrinsics内部寄存器的一些调整，是的更多的指令能并行执行。

但是，在垂直方向的SSE代码用类似的方式调整似乎没有性能的提升，还会到底代码的可读性较差。

　　第三种尝试：不使用中间内存实现的近似效果（80%提速）

以前我在写高斯模糊时考虑到内存占用问题，采用了一种近似的方式，在水平方向计算时，只需要分配一行大小的浮点数据，然后每一行都利用这一行数据做缓存，当一行数据的水平模糊计算完后，就把这些数据转换为字节数据保存到结果图像中，当水平方向都计算完成后，在进行列方向的处理。列方向也是只分配高度大小的一列中间浮点缓存数据，然后进行高度方向处理，每列处理完后，把浮点的结果转换成字节数据。

可见，上述过程存在的一定的精度损失，因为在行方向的处理完成后的浮点到字节数据的转换丢失了部分数据。但是考虑到是模糊，这种丢失对于结果在视觉上是基本察觉不到的。因此，是可以接受的，测试表明，纯C版本的这种做法和纯C版本的标准做法在速度上基本相当。

我们考虑这种做法的SSE优化，第一，是水平方向的处理，想想很简单，核心的代码和之前的是没有区别的，当然我们也应该带上我们的两行一次性处理这种诀窍的。

但是垂直方向呢，如果按照上述方式处理，就无法利用到SSE的优势了，因为上述方式要求每次都是隔行取样，Cache miss的可能性太高，那么还能不能利用我们在高斯模糊算法的全面优化过程分享（一）提高的那种方式呢。

仔细看看（一）中的过程，很明显他一次性只会利用到4行的数据，同时，相邻两行的处理数据有3行是重叠的，那么这就为我们的低内存占用同时又能高效的利用SSE提供了可能性，我们只需要分配4行的浮点缓存区，然后每次交换行行之间的指针，对垂直方向的处理就能利用相同的SIMD优化算法了。

但是这样做又会带来另外一个小问题，就是在Top到Bottom的处理过程中，每一行处理完后又会有一个浮点到字节数据的精度丢失，这种丢失经过测试也是可以接受的。

还有一个问题就是，这样做会增加很多次自己数据到浮点数据的转换，这种转换的耗时是否对最后的结果又重要的影响呢，只有实测才知道。我们待会再分析，这里贴出这种近似的优化的有关代码：

 void GaussBlur_FastAndLowMemory_SSE(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, float Radius)
 {
     float B0, B1, B2, B3;                            
     float *Line0, *Line1, *Line2, *Line3, *Temp;
     int Y = 0;
     CalcGaussCof(Radius, B0, B1, B2, B3);

     float *Buffer = (float *)_mm_malloc(Width * 4 * 4 * sizeof(float), 16);                //    最多需要4行缓冲区

     //    行方向的优化，这个是没有啥精度损失的
     for (; Y < Height - 1; Y += 2)                                                            //    两行执行的代码比单行快
     {
         ConvertBGR8U2BGRAF_Line_SSE(Src + (Y + 0) * Stride, Buffer, Width);
         ConvertBGR8U2BGRAF_Line_SSE(Src + (Y + 1) * Stride, Buffer + Width * 4, Width);            //    读取两行数据
         GaussBlurLeftToRight_TwoLine_SSE(Buffer, Width, B0, B1, B2, B3);                    //    分开来执行速度比写在一起有要快些
         GaussBlurRightToLeft_TwoLine_SSE(Buffer, Width, B0, B1, B2, B3);
         ConvertBGRAF2BGR8U_Line_SSE(Buffer, Dest + (Y + 0) * Stride, Width);                    //    浮点转换为字节数据
         ConvertBGRAF2BGR8U_Line_SSE(Buffer, Dest + (Y + 1) * Stride, Width);
     }
     for (; Y < Height; Y++)                                                                //    执行剩下的单行
     {
         ConvertBGR8U2BGRAF_Line_SSE(Src + Y * Stride, Buffer, Width);
         GaussBlurLeftToRight_OneLine_SSE(Buffer, Width, B0, B1, B2, B3);
         GaussBlurRightToLeft_OneLine_SSE(Buffer, Width, B0, B1, B2, B3);
         ConvertBGRAF2BGR8U_Line_SSE(Buffer, Dest + Y * Stride, Width);
     }

     //    列方向考虑优化，多了一次浮点到字节类型的转换，有精度损失
     ConvertBGR8U2BGRAF_Line_SSE(Dest, Buffer + 3 * Width * 4, Width);
     memcpy(Buffer + 0 * Width * 4, Buffer + 3 * Width * 4, Width * 4 * sizeof(float));            //    起始值取边界的值
     memcpy(Buffer + 1 * Width * 4, Buffer + 3 * Width * 4, Width * 4 * sizeof(float));
     memcpy(Buffer + 2 * Width * 4, Buffer + 3 * Width * 4, Width * 4 * sizeof(float));

     Line0 = Buffer + 0 * Width * 4;    Line1 = Buffer + 1 * Width * 4;
     Line2 = Buffer + 2 * Width * 4;    Line3 = Buffer + 3 * Width * 4;
     for (Y = 0; Y < Height; Y++)
     {
         ConvertBGR8U2BGRAF_Line_SSE(Dest + Y * Stride, Line3, Width);                                //    转换当前行到浮点缓存
         GaussBlurTopToBottom_LowMemory_SSE(Line0, Line1, Line2, Line3, Width, B0, B1, B2, B3);    //    垂直方向处理
         ConvertBGRAF2BGR8U_Line_SSE(Line3, Dest + Y * Stride, Width);                                //    又再次转换为字节数据
         Temp = Line0;    Line0 = Line1;    Line1 = Line2;    Line2 = Line3;    Line3 = Temp;            //    交换行缓存
     }    

     ConvertBGR8U2BGRAF_Line_SSE(Dest + (Height - 1) * Stride, Buffer + 3 * Width * 4, Width);        //    重复边缘像素
     memcpy(Buffer + 0 * Width * 4, Buffer + 3 * Width * 4, Width * 4 * sizeof(float));
     memcpy(Buffer + 1 * Width * 4, Buffer + 3 * Width * 4, Width * 4 * sizeof(float));
     memcpy(Buffer + 2 * Width * 4, Buffer + 3 * Width * 4, Width * 4 * sizeof(float));

     Line0 = Buffer + 0 * Width * 4;    Line1 = Buffer + 1 * Width * 4;
     Line2 = Buffer + 2 * Width * 4;    Line3 = Buffer + 3 * Width * 4;
     for (Y = Height - 1; Y > 0; Y--)                                                            //    垂直向上处理
     {
         ConvertBGR8U2BGRAF_Line_SSE(Dest + Y * Stride, Line0, Width);
         GaussBlurBottomToTop_LowMemory_SSE(Line0, Line1, Line2, Line3, Width, B0, B1, B2, B3);
         ConvertBGRAF2BGR8U_Line_SSE(Line0, Dest + Y * Stride, Width);
         Temp = Line3;    Line3 = Line2;    Line2 = Line1;    Line1 = Line0;    Line0 = Temp;
     }
     _mm_free(Buffer);
 }

　　上述代码中的ConvertBGR8U2BGRAF_Line_SSE和ConvertBGRAF2BGR8U_Line_SSE是之前的相关函数的单行版。

经过测试，上述改进后的算法在同样配置的电脑上，针对3000*2000的彩色图像耗时约为86ms，和之前的145ms相比，提速了近一倍，而基本不占用额外的内存，可是为什么呢，似乎代码中还增加了很多浮点到字节和字节到浮点数据的转换代码，总的计算量应该是增加的啊。按照我的分析，我认为这是这里分配的辅助内存很小，被分配到一级缓存或者二级缓存或其他更靠近CPU的位置的内尺寸区域的可能性更大，而第一版本的内存由于过大，只可能分配堆栈中，同时我们算法里有着大量访问内存的地方，这样虽然总的转换量增加了，但是内存访问节省的时间已经超越了转换增加的时间了。

　　第四种尝试：列方向直接使用BGR而不是BGRA的SSE优化（100%提速）

在高斯模糊算法的全面优化过程分享（一）中，为了解决水平方向上的SSE优化问题，我们将BGR数据转换为了BGRA格式的浮点数后再进行处理，这样在列方向处理时同样需要处理A的数据，但是在经过第三种尝试后，在垂直方向的处理我们还有必要处理这个多余的A吗，当然没有必要，这样垂直方向整体上又可以减少约25%的时间，耗时只有75ms左右了，实现了约100%的提速。

第五种尝试：算法稳定性的考虑和最终的妥协

　在上一篇文章中，我们提到了由于float类型的精度问题，当模糊的半径较大时，算法的结果会出现很大的瑕疵，一种方式就是用double类型来解决，还有一种方式就是可以用Deriche滤波器来解决，为了完美解决这个问题，我还是恨着头皮用SSE实现了Deriche滤波器，这里简要说明如下：

　　Deriche滤波器和高斯滤波器有很多类似的地方：The Deriche filter is a smoothing filter (low-pass) which was designed to optimally detect, along with a derivation operator, the contours in an image (Canny criteria optimization). Besides, as this filter is very similar to a gaussian filter, but much simpler to implement (based on simple first order IIR filters), it is also much used for general image filtering.

按照维基的解释，Deriche滤波器可按照如下的步骤执行：详见https://en.wikipedia.org/wiki/Deriche_edge_detector。

It\'s possible to separate the process of obtaining the value of a two-dimensional Deriche filter into two parts. In first part, image array is passed in the horizontal direction from left to right according to the following formula:

y_{{ij}}^{1}=a_{1}x_{{ij}}+a_{2}x_{{ij-1}}+b_{1}y_{{ij-1}}^{1}+b_{2}y_{{ij-2}}^{1}

and from right to left according to the formula:

y_{{ij}}^{2}=a_{3}x_{{ij+1}}+a_{4}x_{{ij+2}}+b_{1}y_{{ij+1}}^{2}+b_{2}y_{{ij+2}}^{2}

The result of the computation is then stored into temporary two-dimensional array:

\\theta _{{ij}}=c_{1}(y_{{ij}}^{1}+y_{{ij}}^{2})

The second step of the algorithm is very similar to the first one. The two-dimensional array from the previous step is used as the input. It is then passed in the vertical direction from top to bottom and bottom-up according to the following formulas:

y_{{ij}}^{1}=a_{5}\\theta _{{ij}}+a_{6}\\theta _{{i-1j}}+b_{1}y_{{i-1j}}^{1}+b_{2}y_{{i-2j}}^{1}

y_{{ij}}^{2}=a_{7}\\theta _{{i+1j}}+a_{8}\\theta _{{i+2j}}+b_{1}y_{{i+1j}}^{2}+b_{2}y_{{i+2j}}^{2}

\\Theta _{{ij}}=c_{2}(y_{{ij}}^{1}+y_{{ij}}^{2})

　　可见他们也是行列可分离的算法。

同样为了节省内存，我们也使用了类似上述第三种和第四重尝试的方式，但是考虑到Deriche的特殊性（主要是 $\\Theta _{{ij}}=c_{2}(y_{{ij}}^{1}+y_{{ij}}^{2})$ 这里），他还是需要一份中间内存的，为了效率和内存，我们再次以牺牲精度为准备，中间使用了一份和图像一样的字节数据内存。

由于计算量较原先的高斯有所增加，这里最后的优化完成的耗时约为100ms。

第六：多线程

在水平方向计算时，各行之间的计算时独立的，因此是可以并行处理的，但是垂直方向由于是前后依赖的，是无法并行的。比如用openmp做2个线程的并行，大概速度能提高到（高斯）到60ms，但是这个东西在不是本文这里的重点。

　第七：比较

同标准的高斯滤波相比，Deriche滤波器由于其特性，无法支持In-Place操作，也就是说Src和Dest不能相同，如果一定要相同，就只有通过一个中间内存来过渡了，而标准高斯是可以的。第二就是高斯是可以不占用太多额外的内存就可以实现的，而Deriche需要一份同样大小的内存。第三就是标准高斯速度还是要快一点。第四Deriche滤波器的精度在float类型时精度要比标准高斯高。综合选择，我觉得还是以后选择Deriche代替标准的高斯模糊。

总结：有心就有成绩

同opencv的cvsmooth相比，同样的机器上同样的3000*2000大小的彩图，Ksize我取100时，需要1200多ms，比我这里慢了10倍，但是cvsmooth似乎对ksize参数敏感，他并不是与核大小无关的，ksize较小时还会很快的，不过除了一些特效外，在很多场合我们其实需要大半径的高斯的（比如图像增强、锐化上）。

做完了在回头看看优化的过程，觉得和看书是一个道理，先是越看越厚，通了就像一张薄纸一样。

最后总结下，就是一件事情，只要你有时间和信心，就能有进步，坚持是胜利的必要条件。

提供一个测试的Demo: http://files.cnblogs.com/files/Imageshop/FastGaussBlur.rar

由测试Demo可以测试出，当选择低内存近似版本或者准确版本时，当半径较大时，如果连续的拖动滚动条，图像会出现闪烁，而如果选择Deriche时，则图像变换很平缓，而当半径特别大时，如果选择低内存近似版本或者准确版本，则图像有可能会出现线条或者色块，只有Deriche滤波的效果是完美的。