一文中大家已经交给了一种异常高性能的高斯模糊过程,这里把优化过程中的一些体会和取得用文字的样式记录下来

     
相关链接: 高斯模糊算法的系数优化过程分享(一)

     
相关链接: 高斯模糊算法的两全优化过程分享(一)

   
 这里的高斯模糊使用的是舆论《Recursive implementation of the Gaussian
filter》里描述的递归算法。

   
 在高斯模糊算法的一揽子优化过程分享(一)一文中我们已经付出了一种异常高性能的高斯模糊过程,不过优化没有终点,经过上一个星期的奋斗和测试,对该算法的功效提高又有了一个新的莫大,这里把优化过程中的一些经验和得到用文字的款式记录下来。

   
 在高斯模糊算法的两全优化过程分享(一)一文中大家已经付出了一种分外高性能的高斯模糊过程,然则优化没有终点,经过上一个星期的辛劳奋斗和测试,对该算法的效能提高又有了一个新的惊人,这里把优化过程中的一些经验和取得用文字的样式记录下来。

           
  图片 1

     第一个尝试   直接使用内联汇编替代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,那一个对功夫的渴求就更高了,所以说网络上的传道也不可以完全信任,而如若不是有特意强的汇编能力,也毫无去挑衅编译器。

     第一个尝试   直接选拔内联汇编替代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,这一个对功夫的渴求就更高了,所以说网络上的传教也不可以完全倚重,而只要不是有专门强的汇编能力,也毫不去挑衅编译器。

     
仔细察看和清楚上述公式,在forward过程中,n是递增的,由此,即便在展开forward此前,把in数据先整体的赋值给w,然后式子(9a)就可以改为:

    第二个尝试   水平方向的歪曲两遍实施二行(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%的涨价。

    第二个尝试   水平方向的模糊两遍施行二行(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%的涨价。

      w[n] = B w[n] +
(b1 w[n-1] + b2 w[n-2] +
b3 w[n-3]) / b0;          ———>    
(1a)

      第五种尝试:算法稳定性的考虑和末段的折衷

  在上一篇著作中,我们提到了是因为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:

图片 2

and from right to left according to the formula:

图片 3

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

图片 4 
                                                

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:

图片 5

图片 6

图片 7

  可见他们也是行列可分其它算法。

   
 同样为了节省内存,大家也运用了类似上述第两种和第四重尝试的主意,然而考虑到Deriche的特殊性(重如若图片 8此地),他要么需要一份中间内存的,为了功能和内存,大家重新以献身精度为准备,中间使用了一份和图像一样的字节数据内存。

   
由于总结量较原先的高斯有所增多,这里最终的优化完成的耗时约为100ms。

      第五种尝试:算法稳定性的设想和最终的投降

  在上一篇散文中,我们提到了是因为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:

图片 9

and from right to left according to the formula:

图片 10

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

图片 11 
                                                

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:

图片 12

图片 13

图片 14

  可见他们也是行列可分另外算法。

   
 同样为了省去内存,我们也采用了近似上述第二种和第四重尝试的主意,不过考虑到Deriche的特殊性(首借使图片 15这里),他要么需要一份中间内存的,为了效能和内存,大家再度以牺牲精度为准备,中间使用了一份和图像一样的字节数据内存。

   
由于总计量较原先的高斯有所增多,这里最终的优化完成的耗时约为100ms。

 
  在backward过程中,n是递减的,因此在backward前,把w的数目完整的拷贝到out中,则式子(9b)变为:

     第六:多线程

   
 在档次方向总结时,各行之间的统计时独立的,因而是可以并行处理的,但是垂直方向由于是前后看重的,是不能并行的。比如用openmp做2个线程的相互,大概速度能加强到(高斯)到60ms,不过这一个事物在不是本文这里的关键。

     第六:多线程

   
 在档次方向总结时,各行之间的乘除时独自的,由此是足以并行处理的,可是垂直方向由于是左右依赖的,是力不从心并行的。比如用openmp做2个线程的相互,大概速度能增强到(高斯)到60ms,不过这一个东西在不是本文这里的第一。

     out[n] = B out[n] +
(b1 
out[n+1] + b2 out[n+2] +
b3 
out[n+3]) / b0 ;     <———     (1b)

  第七:比较

   
同专业的高斯滤波比较,Deriche滤波器由于其特点,不能够支撑In-Place操作,也就是说Src和Dest不可以平等,如果一定要一如既往,就只有因此一个中等内存来过渡了,而正规高斯是可以的。第二就是高斯是足以不占用太多额外的内存就可以实现的,而Deriche需要一份同样大小的内存。第三就是规范高斯速度仍旧要快一些。第四Deriche滤波器的精度在float类型时精度要比标准高斯高。综合取舍,我觉得仍然之后选取Deriche代替标准的高斯模糊。

  第七:比较

   
同专业的高斯滤波相比较,Deriche滤波器由于其性状,不可能支撑In-Place操作,也就是说Src和Dest不可以平等,假如一定要一致,就只有由此一个中等内存来过渡了,而正规高斯是可以的。第二就是高斯是足以不占用太多额外的内存就可以实现的,而Deriche需要一份同样大小的内存。第三就是规范高斯速度仍旧要快一些。第四Deriche滤波器的精度在float类型时精度要比正规高斯高。综合取舍,我认为如故之后采纳Deriche代替标准的高斯模糊。

   
 从编程角度看来,backward中的拷贝是一心没有必要的,因而 (1b)可以直接写为:

     总计:有心就有成绩

   
同opencv的cvsmooth相比较,同样的机器上同样的3000*2000大小的彩图,Ksize我取100时,需要1200多ms,比自己这边慢了10倍,可是cvsmooth似乎对ksize参数敏感,他并不是与核大小无关的,ksize较刻钟还会很快的,但是除了有些特效外,在广大场子我们实际上需要大半径的高斯的(比如图像增强、锐化上)。

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

   
最后总括下,就是一件工作,只要您有时光和信念,就能有发展,坚定不移是战胜的必要条件。

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

   
由测试Demo可以测试出,当采取低内存近似版本或者纯粹版本时,当半径较大时,如若老是的拖动滚动条,图像会并发闪烁,而只要采取Deriche时,则图像变换很温和,而当半径特别大时,尽管采用低内存近似版本或者纯粹版本,则图像有可能会出现线条依然色块,唯有Deriche滤波的效益是完善的。

图片 16

  高斯模糊的优化到此截止,假如有何人有用GPU实现的,还请报告我下大约的耗时。

     拒绝无脑索取代码。

图片 17

 

     统计:有心就有成就

   
同opencv的cvsmooth相相比,同样的机械上平等的3000*2000轻重的彩图,Ksize我取100时,需要1200多ms,比我那边慢了10倍,但是cvsmooth似乎对ksize参数敏感,他并不是与核大小无关的,ksize较时辰还会快速的,不过除了有的神效外,在众多场所大家实在需要大半径的高斯的(比如图像增强、锐化上)。

   
做完了在自查自纠看看优化的经过,觉得和看书是一个道理,先是越看越厚,通了就像一张薄纸一样。

   
最终统计下,就是一件事情,只要你有时间和信心,就能有开拓进取,坚定不移是胜利的必要条件。

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

   
由测试Demo可以测试出,当选拔低内存近似版本或者纯粹版本时,当半径较大时,即便老是的拖动滚动条,图像会现出闪烁,而只要采纳Deriche时,则图像变换很温和,而当半径特别大时,假设采纳低内存近似版本或者纯粹版本,则图像有可能会合世线条仍然色块,唯有Deriche滤波的效益是一揽子的。

图片 18

  高斯模糊的优化到此截止,如果有何人有用GPU实现的,还请告知自己下大约的耗时。

     拒绝无脑索取代码。

图片 19

 

           w[n] = B w[n] +
(b1 w[n+1] + b2 w[n+2] +
b3 w[n+3]) / b0 ;               <———    
(1c)

   
 
从进度上考虑,浮点除法很慢,因而,大家再一次定义b1 = b1 / b0, b2 = b2 /
b0, b3 = b3 / b0, 最终得到大家应用的递归公式:

           w[n] = B w[n] +
b1 w[n-1] + b2 w[n-2] +
b3 w[n-3];          ———>     (a)

        w[n] = B w[n] +
b1 w[n+1] + b2 w[n+2] +
b3 w[n+3] ;             <———      (b)

   
上述公式是一维的高斯模糊总计方法,针对二维的图像,正确的做法就是先对每个图像行举行模糊处理得到中间结果,再对中等结果的每列举行模糊操作,最后得到二维的模糊结果,当然也足以采纳先列后行这样的做法。

     
另外注意点就是,边缘像素的拍卖,大家看来在公式(a)或者(b)中,w[n]的结果个别凭借于前五个或者后六个元素的值,而对于边缘地方的值,这个都是不在合理范围内的,平常的做法是镜像数据或者另行边缘数据,实践注明,由于是递归的算法,起始值的不等会将结果直接延续下去,因而,不同的办法对边缘部分的结果要么有必然的震慑的,这里大家使用的简要的再一次边缘像素值的法门。

     
由于地方公式中的周详均为浮点类型,因而,总计一般也是对浮点举办的,也就是说一般需要先把图像数据转换为浮点,然后开展高斯模糊,在将结果转换为字节类型的图像,因而,上述过程可以分级用一下多少个函数完成:

               
CalcGaussCof          
//  总结高斯模糊中行使到的周全
      ConvertBGR8U2BGRAF      //  将字节数据转换为浮点数据 
      GaussBlurFromLeftToRight    //  水平方向的前向传播
      GaussBlurFromRightToLeft    //  水平方向的反向传播
      GaussBlurFromTopToBottom  
 //   垂直方向的前向传播
      GaussBlurFromBottomToTop  
 //   垂直方向的反向传播
      ConvertBGRAF2BGR8U      
 //   将结果转换为字节数据

   大家逐一分析之。

     
 CalcGaussCof,这么些很简短,就遵照杂谈中提议的总括公式依据用户输入的参数来统计,当然结合下方面我们提到的除法变乘法的优化,注意,为了持续的片段标号的问题,我讲上述公式(a)和(b)中的全面B写为b0了。

void CalcGaussCof(float Radius, float &B0, float &B1, float &B2, float &B3)
{
    float Q, B;
    if (Radius >= 2.5)
        Q = (double)(0.98711 * Radius - 0.96330);                            //    对应论文公式11b
    else if ((Radius >= 0.5) && (Radius < 2.5))
        Q = (double)(3.97156 - 4.14554 * sqrt(1 - 0.26891 * Radius));
    else
        Q = (double)0.1147705018520355224609375;

    B = 1.57825 + 2.44413 * Q + 1.4281 * Q * Q + 0.422205 * Q * Q * Q;        //    对应论文公式8c
    B1 = 2.44413 * Q + 2.85619 * Q * Q + 1.26661 * Q * Q * Q;
    B2 = -1.4281 * Q * Q - 1.26661 * Q * Q * Q;
    B3 = 0.422205 * Q * Q * Q;

    B0 = 1.0 - (B1 + B2 + B3) / B;
    B1 = B1 / B;
    B2 = B2 / B;
    B3 = B3 / B;
}

  由上述代码可见,B0+B1+B2+B3=1,即是归一化的周到,这也是算法可以递归举行的重要之一。

   
 接着为了便于中间经过,我们先将字节数据转换为浮点数据,这部分代码也很粗略:

void ConvertBGR8U2BGRAF(unsigned char *Src, float *Dest, int Width, int Height, int Stride)
{
    //#pragma omp parallel for
    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePS = Src + Y * Stride;
        float *LinePD = Dest + Y * Width * 3;
        for (int X = 0; X < Width; X++, LinePS += 3, LinePD += 3)
        {
            LinePD[0] = LinePS[0];    LinePD[1] = LinePS[1];    LinePD[2] = LinePS[2];
        }
    }
}

  大家可以尝尝自己把内部的X循环再拓展探访效果。

   
 水平方向的前向传播依据公式去写也是很简短的,但是一向动用公式里面有为数不少拜访数组的进程(不必然就慢),我稍稍改造下写成如下格局:

void GaussBlurFromLeftToRight(float *Data, int Width, int Height, float B0, float B1, float B2, float B3)
{
    //#pragma omp parallel for
    for (int Y = 0; Y < Height; Y++)
    {
        float *LinePD = Data + Y * Width * 3;
        float BS1 = LinePD[0], BS2 = LinePD[0], BS3 = LinePD[0];          //  边缘处使用重复像素的方案
        float GS1 = LinePD[1], GS2 = LinePD[1], GS3 = LinePD[1];
        float RS1 = LinePD[2], RS2 = LinePD[2], RS3 = LinePD[2];
        for (int X = 0; X < Width; X++, LinePD += 3)
        {
            LinePD[0] = LinePD[0] * B0 + BS1 * B1 + BS2 * B2 + BS3 * B3;
            LinePD[1] = LinePD[1] * B0 + GS1 * B1 + GS2 * B2 + GS3 * B3;         // 进行顺向迭代
            LinePD[2] = LinePD[2] * B0 + RS1 * B1 + RS2 * B2 + RS3 * B3;
            BS3 = BS2, BS2 = BS1, BS1 = LinePD[0];
            GS3 = GS2, GS2 = GS1, GS1 = LinePD[1];
            RS3 = RS2, RS2 = RS1, RS1 = LinePD[2];
        }
    }
}

  不多描述,请大家把上述代码和公式(a)对应一下就精通了。

     
有了GaussBlurFromLeftToRight的参阅代码,GaussBlurFromRightToLeft的代码就不会有怎样大的孤苦了,为了制止有些懒人不看随笔不想想,这里自己不贴这段的代码。

     
那么垂直方向上粗略的做只需要变更下循环的样子,以及历次指针扩充量更改为Width
* 3即可。

     
那么我们来设想下垂直方向能不可能不这么处理呢,指针每一次扩充Width *
3会促成惨重的Cache
miss,特别是对此宽度稍微大点的图像,我们精心考察垂直方向,会发觉仍旧可以按照Y
 — X那样的大循环形式也是可以的,代码如下:

void GaussBlurFromTopToBottom(float *Data, int Width, int Height, float B0, float B1, float B2, float B3)
{
    for (int Y = 0; Y < Height; Y++)
    {
        float *LinePD3 = Data + (Y + 0) * Width * 3;
        float *LinePD2 = Data + (Y + 1) * Width * 3;
        float *LinePD1 = Data + (Y + 2) * Width * 3;
        float *LinePD0 = Data + (Y + 3) * Width * 3;
        for (int X = 0; X < Width; X++, LinePD0 += 3, LinePD1 += 3, LinePD2 += 3, LinePD3 += 3)
        {
            LinePD0[0] = LinePD0[0] * B0 + LinePD1[0] * B1 + LinePD2[0] * B2 + LinePD3[0] * B3;
            LinePD0[1] = LinePD0[1] * B0 + LinePD1[1] * B1 + LinePD2[1] * B2 + LinePD3[1] * B3;
            LinePD0[2] = LinePD0[2] * B0 + LinePD1[2] * B1 + LinePD2[2] * B2 + LinePD3[2] * B3;
        }
    }
}

  就是说我们不是一念之差就把列方向的前向传播举行完,而是每一次只举行一个像素的扩散,当一行所有像素都进行完了列方向的前向传播后,在切换来下一行,这样处理就制止了Cache
miss出现的状态。

   
 注意到列方向的边缘地方,为了便于代码的拍卖,大家在中度方向上上下分别增添了3个行的像素,当举行完中间的实用行的行方向前向和反向传来后,遵照前述的重复边缘像素的办法填充上下那空出的三行数据。

   
 同样的道理,GaussBlurFromBottomToTop的代码可由读者自己补充进去。

   
 最后的ConvertBGRAF2BGR8U也很简单,就是一个for循环:

void ConvertBGRAF2BGR8U(float *Src, unsigned char *Dest, int Width, int Height, int Stride)
{
    //#pragma omp parallel for
    for (int Y = 0; Y < Height; Y++)
    {
        float *LinePS = Src + Y * Width * 3;
        unsigned char *LinePD = Dest + Y * Stride;
        for (int X = 0; X < Width; X++, LinePS += 3, LinePD += 3)
        {
            LinePD[0] = LinePS[0];    LinePD[1] = LinePS[1];    LinePD[2] = LinePS[2];
        }
    }
}

   在有效的界定内,上述浮点总括的结果不会超出byte所能表达的限制,由此也不需要特地的进展Clamp操作。

     
 最终就是局部内存分配和刑释解教的代码了,如下所示:

void GaussBlur(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, float Radius)
{
    float B0, B1, B2, B3;
    float *Buffer = (float *)malloc(Width * (Height + 6) * sizeof(float) * 3);

    CalcGaussCof(Radius, B0, B1, B2, B3);
    ConvertBGR8U2BGRAF(Src, Buffer + 3 * Width * 3, Width, Height, Stride);

    GaussBlurFromLeftToRight(Buffer + 3 * Width * 3, Width, Height, B0, B1, B2, B3);
    GaussBlurFromRightToLeft(Buffer + 3 * Width * 3, Width, Height, B0, B1, B2, B3);        //    如果启用多线程,建议把这个函数写到GaussBlurFromLeftToRight的for X循环里,因为这样就可以减少线程并发时的阻力

    memcpy(Buffer + 0 * Width * 3, Buffer + 3 * Width * 3, Width * 3 * sizeof(float));
    memcpy(Buffer + 1 * Width * 3, Buffer + 3 * Width * 3, Width * 3 * sizeof(float));
    memcpy(Buffer + 2 * Width * 3, Buffer + 3 * Width * 3, Width * 3 * sizeof(float));

    GaussBlurFromTopToBottom(Buffer, Width, Height, B0, B1, B2, B3);

    memcpy(Buffer + (Height + 3) * Width * 3, Buffer + (Height + 2) * Width * 3, Width * 3 * sizeof(float));
    memcpy(Buffer + (Height + 4) * Width * 3, Buffer + (Height + 2) * Width * 3, Width * 3 * sizeof(float));
    memcpy(Buffer + (Height + 5) * Width * 3, Buffer + (Height + 2) * Width * 3, Width * 3 * sizeof(float));

    GaussBlurFromBottomToTop(Buffer, Width, Height, B0, B1, B2, B3);

    ConvertBGRAF2BGR8U(Buffer + 3 * Width * 3, Dest, Width, Height, Stride);

    free(Buffer);
}

  正如上所述,分配了Height +
6行的内存区域,首假如为着便于垂直方向的拍卖,在执行GaussBlurFromTopToBottom从前遵照重复边缘的规范复制3行,然后在GaussBlurFromBottomToTop此前在复制底部边缘的3行像素。

     
至此,一个简单而又急迅的高斯模糊就主旨到位了,代码数量也不多,也从没啥难度,不了解为何许六个人搞不定。

     
依照自己的测试,上述措施代码在一台I5-6300HQ
2.30GHZ的笔记本上针对一副3000*2000的24位图像的拍卖时间大体需要370ms,假设在C++的编译选项的代码生成选项里的启用增强指令集拔取–>流式处理SIMD增添2(/arch:sse2),则编译后的主次大概需要220ms的时间。

     
我们尝试使用系统的部分资源进一步提升速度,首先大家想到了SSE优化,关于这上边英特尔有一篇官方的著作和代码:

        IIR Gaussian Blur Filter
Implementation using Intel® Advanced Vector
Extensions
 [PDF
513KB]

     source
code: gaussian_blur.cpp [36KB]

     
我只是简短的看了下这篇随笔,感觉他其中用到的总括公式和Deriche滤波器的很像,和本文参考的Recursive
implementation
不太雷同,并且其SSE代码对能处理的图还有众多范围,对本身这边的参阅意义不大。

     
大家先看下核心的臆想的SSE优化,注意到  GaussBlurFromLeftToRight
的代码中,
其主旨的计量部分是多少个乘法,可是他只有3个乘法总计,虽然可以凑成四行,那么就足以充裕利用SSE的批量划算效能了,也就是一旦能扩张一个坦途,使得GaussBlurFromLeftToRight变为如下形式:

void GaussBlurFromLeftToRight(float *Data, int Width, int Height, float B0, float B1, float B2, float B3)
{
    //#pragma omp parallel for
    for (int Y = 0; Y < Height; Y++)
    {
        float *LinePD = Data + Y * Width * 4;
        float BS1 = LinePD[0], BS2 = LinePD[0], BS3 = LinePD[0];                //  边缘处使用重复像素的方案
        float GS1 = LinePD[1], GS2 = LinePD[1], GS3 = LinePD[1];
        float RS1 = LinePD[2], RS2 = LinePD[2], RS3 = LinePD[2];
        float AS1 = LinePD[3], AS2 = LinePD[3], AS3 = LinePD[3];
        for (int X = 0; X < Width; X++, LinePD += 4)
        {
            LinePD[0] = LinePD[0] * B0 + BS1 * B1 + BS2 * B2 + BS3 * B3;
            LinePD[1] = LinePD[1] * B0 + GS1 * B1 + GS2 * B2 + GS3 * B3;         // 进行顺向迭代
            LinePD[2] = LinePD[2] * B0 + RS1 * B1 + RS2 * B2 + RS3 * B3;
            LinePD[3] = LinePD[3] * B0 + AS1 * B1 + AS2 * B2 + AS3 * B3;
            BS3 = BS2, BS2 = BS1, BS1 = LinePD[0];
            GS3 = GS2, GS2 = GS1, GS1 = LinePD[1];
            RS3 = RS2, RS2 = RS1, RS1 = LinePD[2];
            AS3 = AS2, AS2 = AS1, AS1 = LinePD[3];
        }
    }
}

  则很容易就把上述代码转换成SSE可以规范处理的代码了。

  而对于Y方向的代码,你精心观看会发现,无论是及通道的图,天然的就可以使用SSE进行拍卖,详见下边的代码。

  好,我们依然一个一个的来分析:

   第一个函数
CalcGaussCof 无需举行其他的优化。

     
第二个函数 ConvertBGR8U2BGRAF依据上述分析需要再行写,因为急需追加一个大路,新的康庄大道的值填0或者其他值都得以,但提出填0,这对有些SSE函数很有用,我把这一个函数的SSE实现共享一下:

void ConvertBGR8U2BGRAF_SSE(unsigned char *Src, float *Dest, int Width, int Height, int Stride)
{
    const int BlockSize = 4;
    int Block = (Width - 2) / BlockSize;
    __m128i Mask = _mm_setr_epi8(0, 1, 2, -1, 3, 4, 5, -1, 6, 7, 8, -1, 9, 10, 11, -1);            //    Mask为-1的地方会自动设置数据为0
    __m128i Zero = _mm_setzero_si128();
    //#pragma omp parallel for
    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePS = Src + Y * Stride;
        float *LinePD = Dest + Y * Width * 4;
        int X = 0;
        for (; X < Block * BlockSize; X += BlockSize, LinePS += BlockSize * 3, LinePD += BlockSize * 4)
        {
            __m128i SrcV = _mm_shuffle_epi8(_mm_loadu_si128((const __m128i*)LinePS), Mask);        //    提取了16个字节,但是因为24位的,我们要将其变为32位的,所以只能提取出其中的前12个字节
            __m128i Src16L = _mm_unpacklo_epi8(SrcV, Zero);
            __m128i Src16H = _mm_unpackhi_epi8(SrcV, Zero);
            _mm_store_ps(LinePD + 0, _mm_cvtepi32_ps(_mm_unpacklo_epi16(Src16L, Zero)));        //    分配内存时已经是16字节对齐了,然后每行又是4的倍数的浮点数,因此,每行都是16字节对齐的
            _mm_store_ps(LinePD + 4, _mm_cvtepi32_ps(_mm_unpackhi_epi16(Src16L, Zero)));        //    _mm_stream_ps是否快点?
            _mm_store_ps(LinePD + 8, _mm_cvtepi32_ps(_mm_unpacklo_epi16(Src16H, Zero)));
            _mm_store_ps(LinePD + 12, _mm_cvtepi32_ps(_mm_unpackhi_epi16(Src16H, Zero)));
        }
        for (; X < Width; X++, LinePS += 3, LinePD += 4)
        {
            LinePD[0] = LinePS[0];    LinePD[1] = LinePS[1];    LinePD[2] = LinePS[2];    LinePD[3] = 0;        //    Alpha最好设置为0,虽然在下面的CofB0等SSE常量中通过设置ALPHA对应的系数为0,可以使得计算结果也为0,但是不是最合理的
        }
    }
}

  稍作解释,_mm_loadu_si128几遍性加载16个字节的数目到SSE寄存器中,对于24位图像,16个字节里富含了5

  • 1 /
    3个像素的消息,而我们的目标是把这个多少转换为4通道的音信,由此,大家不得不一回性的领到处16/4=4个像素的值,这借助于_mm_shuffle_epi8函数和适度的Mask来贯彻,_mm_unpacklo_epi8/_mm_unpackhi_epi8分别领取了SrcV的高8位和低8位的8个字节数据并将它们转换为8个16进制数保存到Src16L和Src16H中,
    而_mm_unpacklo_epi16/_mm_unpackhi_epi16则更为把16位数据增添到32位整形数据,最后通过_mm_cvtepi32_ps函数把整形数据转换为浮点型。

  可能有人注意到了 int Block = (Width – 2) /
BlockSize;
这一行代码,为何要-2操作呢,这也是本身在频繁测试发现先后连接出现内存错误时才发现到的,因为_mm_loadu_si128五回性加载了5

  • 1 / 3个像素的消息,当在拍卖最终一行像素时(其他行不会),假诺Block
    取Width/BlockSize,
    则很有可能访问了超出像素范围内的内存,而-2不是-1就是因为非常额外的1/3像素起的坚守。

  接着看档次方向的前向传来的SSE方案:

void GaussBlurFromLeftToRight_SSE(float *Data, int Width, int Height, float B0, float B1, float B2, float B3)
{
    const  __m128 CofB0 = _mm_set_ps(0, B0, B0, B0);
    const  __m128 CofB1 = _mm_set_ps(0, B1, B1, B1);
    const  __m128 CofB2 = _mm_set_ps(0, B2, B2, B2);
    const  __m128 CofB3 = _mm_set_ps(0, B3, B3, B3);
    //#pragma omp parallel for
    for (int Y = 0; Y < Height; Y++)
    {
        float *LinePD = Data + Y * Width * 4;
        __m128 V1 = _mm_set_ps(LinePD[3], LinePD[2], LinePD[1], LinePD[0]);
        __m128 V2 = V1, V3 = V1;
        for (int X = 0; X < Width; X++, LinePD += 4)            //    还有一种写法不需要这种V3/V2/V1递归赋值,一次性计算3个值,详见 D:\程序设计\正在研究\Core\ShockFilter\Convert里的高斯模糊,但速度上没啥区别
        {
            __m128 V0 = _mm_load_ps(LinePD);
            __m128 V01 = _mm_add_ps(_mm_mul_ps(CofB0, V0), _mm_mul_ps(CofB1, V1));
            __m128 V23 = _mm_add_ps(_mm_mul_ps(CofB2, V2), _mm_mul_ps(CofB3, V3));
            __m128 V = _mm_add_ps(V01, V23);
            V3 = V2;    V2 = V1;    V1 = V;
            _mm_store_ps(LinePD, V);
        }
    }
}

  和上面的4通道的GaussBlurFromLeftToRight_SSE相比较,你会发现差不多语法上未曾此外的区别,实在是太简单了,注意自己从未用_mm_storeu_ps,而是直接运用_mm_store_ps,这是因为,第一,分配Data内存时,我使用了_mm_malloc分配了16字节对齐的内存,而Data每行元素的个数又都是4的倍数,由此,每行先河地方处的内存也是16字节对齐的,所以从来_mm_store_ps完全是足以的。

   
 同理,在笔直方向的前向传来的SSE优化代码就更直接了:

void GaussBlurFromTopToBottom_SSE(float *Data, int Width, int Height, float B0, float B1, float B2, float B3)
{
    const  __m128 CofB0 = _mm_set_ps(0, B0, B0, B0);
    const  __m128 CofB1 = _mm_set_ps(0, B1, B1, B1);
    const  __m128 CofB2 = _mm_set_ps(0, B2, B2, B2);
    const  __m128 CofB3 = _mm_set_ps(0, B3, B3, B3);
    for (int Y = 0; Y < Height; Y++)
    {
        float *LinePS3 = Data + (Y + 0) * Width * 4;
        float *LinePS2 = Data + (Y + 1) * Width * 4;
        float *LinePS1 = Data + (Y + 2) * Width * 4;
        float *LinePS0 = Data + (Y + 3) * Width * 4;
        for (int X = 0; X < Width * 4; X += 4)
        {
            __m128 V3 = _mm_load_ps(LinePS3 + X);
            __m128 V2 = _mm_load_ps(LinePS2 + X);
            __m128 V1 = _mm_load_ps(LinePS1 + X);
            __m128 V0 = _mm_load_ps(LinePS0 + X);
            __m128 V01 = _mm_add_ps(_mm_mul_ps(CofB0, V0), _mm_mul_ps(CofB1, V1));
            __m128 V23 = _mm_add_ps(_mm_mul_ps(CofB2, V2), _mm_mul_ps(CofB3, V3));
            _mm_store_ps(LinePS0 + X, _mm_add_ps(V01, V23));
        }
    }
}

  对下边的代码不想做其他解释了。

  最有难度的应有算是ConvertBGRAF2BGR8U的SSE版本了,由于一些原因,我不太情愿分享那一个函数的代码,有趣味的敌人可以参见opencv的关于实现。

   
 最终的SSE版本高斯模糊的严重性代码如下:

void GaussBlur_SSE(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, float Radius)
{
    float B0, B1, B2, B3;
    float *Buffer = (float *)_mm_malloc(Width * (Height + 6) * sizeof(float) * 4, 16);

    CalcGaussCof(Radius, B0, B1, B2, B3);
    ConvertBGR8U2BGRAF_SSE(Src, Buffer + 3 * Width * 4, Width, Height, Stride);

    GaussBlurFromLeftToRight_SSE(Buffer + 3 * Width * 4, Width, Height, B0, B1, B2, B3);        //    在SSE版本中,这两个函数占用的时间比下面两个要多,不过C语言版本也是一样的
    GaussBlurFromRightToLeft_SSE(Buffer + 3 * Width * 4, Width, Height, B0, B1, B2, B3);        //    如果启用多线程,建议把这个函数写到GaussBlurFromLeftToRight的for X循环里,因为这样就可以减少线程并发时的阻力

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

    GaussBlurFromTopToBottom_SSE(Buffer, Width, Height, B0, B1, B2, B3);

    memcpy(Buffer + (Height + 3) * Width * 4, Buffer + (Height + 2) * Width * 4, Width * 4 * sizeof(float));
    memcpy(Buffer + (Height + 4) * Width * 4, Buffer + (Height + 2) * Width * 4, Width * 4 * sizeof(float));
    memcpy(Buffer + (Height + 5) * Width * 4, Buffer + (Height + 2) * Width * 4, Width * 4 * sizeof(float));

    GaussBlurFromBottomToTop_SSE(Buffer, Width, Height, B0, B1, B2, B3);

    ConvertBGRAF2BGR8U_SSE(Buffer + 3 * Width * 4, Dest, Width, Height, Stride);

    _mm_free(Buffer);
}

  对于同样的3000*2000的彩色图像,SSE版本的代码耗时只有145ms的耗时,相对于一般性的C代码有约2.5倍左右的涨价,这也事出有因,因为我们在实践SSE的代码时时多处理了一个坦途的总结量的,可是同编译器自身的SSE优化220ms,唯有1.5倍的涨价了,这表明编译器的SSE优化能力或者十分强的。

   
 进一步的优化就是本身上边的注释掉的opemp相关代码,在ConvertBGR8U2BGRAF /
GaussBlurFromLeftToRight / GaussBlurFromRightToLeft / ConvertBGRAF2BGR8U
 4个函数中,直接助长简单的#pragma omp parallel
for就可以了,至于GaussBlurFromTopToBottom_SSE、
GaussBlurFromBottomToTop_SSE则是因为上下行之间数据的相关性,是心有余而力不足落实并行统计的,可是测试表示他们的耗时要比水平方向的少很多。

   
比如大家指定openmp使用2个线程,在上述机器上(四核的),纯C版本能优化到252ms,而纯SSE的只可以增强到100ms左右,编译器自身的SSE优化速度大约是150ms,基本上仍旧维持同一个级其余百分比。

   对于灰度图像,很惋惜,上述的水平方向上的优化措施就不管为力了,一种艺术就是把4行灰度像素混洗成一行,可是这么些操作不太方便用SSE实现,其它一种就是把水平方向的数量先转置,然后接纳垂直方向的SSE优化代码处理,完成在转置回去,最终对转置的数码再一次开展垂直方向SSE优化,当然转置的经过是可以依靠SSE的代码实现的,然而急需非常的一份内存,速度上恐怕和一般性的C相比就不会有那么多的升级换代了,这一个待有时光了再去测试。

   
 前前后后写这个大约也花了半个月的年月,写完了任何工艺流程,在倒过来看,真的是至极的简短,有的时候就是如此。

   
 本文并没有提供完整的可以从来实施的成套代码,需者自勤,提供一段C#的调用工程供有趣味的爱侣测试和比对(未利用OPENMP版本)。

   
 http://files.cnblogs.com/files/Imageshop/GaussBLur_Sample.rar

图片 20

     后记:后来测试发现,当半径参数较大时,无论是C版本依旧SSE版本都会现出局部语无伦次的通病,感觉像是溢出了,后来意识根本是当半径变大后,B0参数变得很小,以至于用float类型的多少来处理递归已经黔驴技穷确保丰富的精度了,解决的艺术是行使double类型,这对于C语言版本的话是秒秒的政工,而对于SSE版本则是较大的天灾人祸,double时换成AVX可能改动量不大,然而AVX的普及度以及…..,可是,一般景色下,半径不超过75时结果都仍然不错的,这对于绝大多数的采纳来说是十足的,半径75时,整个图像已经差不多没有其余的底细了,再大,区别也不大了。

   
 解决上述问题一个可行的方案就是拔取Deriche滤波器,我用该滤波器的float版本对大半径是不会产出问题的,并且也有相关SSE参考代码。

 图片 21

   后续小说:高斯模糊算法的完善优化过程分享(二)。

 

相关文章