SSE图像算法优化系列九:灵活运用SIMD指令16倍提升Sobel边缘检测的速度(4000*3000的24位图像时间由480ms降低到30ms)。

  这半年多时间,基本都在折腾一些基本的优化,有很多都是十几年前的技术了,从随大流的角度来考虑,研究这些东西在很多人看来是浪费时间了,即不能赚钱,也对工作能力提升无啥帮助。可我觉得人类所谓的幸福,可以分为物质档次的享受,还有更为复杂的精神上的富有,哪怕这种富有只是存在于短暂的自我满足中也是值得的。 

     闲话少说, SIMD指令集,这个古老的东西,从第一代开始算起,也快有近20年的历史了,从最开始的MMX技术,到SSE,以及后来的SSE2、SSE3、SSE4、AVX以及11年以后的AVX2,逐渐的成熟和丰富,不过目前考虑通用性方面,AVX的辐射范围还是有限,大部分在优化时还是考虑使用128位的SSE指令集,我在之前的一系列文章中,也有不少文章涉及到了这个方面的优化了。

      今天我们来学习下Sobel算法的优化,首先,我们给出传统的C++实现的算法代码:

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

    unsigned char *RowCopy = (unsigned char *)malloc((Width + 2) * 3 * Channel);
    if (RowCopy == NULL)    return IM_STATUS_OUTOFMEMORY;

    unsigned char *First = RowCopy;
    unsigned char *Second = RowCopy + (Width + 2) * Channel;
    unsigned char *Third = RowCopy + (Width + 2) * 2 * Channel;

    memcpy(Second, Src, Channel);
    memcpy(Second + Channel, Src, Width * Channel);                                                    //    拷贝数据到中间位置
    memcpy(Second + (Width + 1) * Channel, Src + (Width - 1) * Channel, Channel);

    memcpy(First, Second, (Width + 2) * Channel);                                                    //    第一行和第二行一样

    memcpy(Third, Src + Stride, Channel);                                                    //    拷贝第二行数据
    memcpy(Third + Channel, Src + Stride, Width * Channel);
    memcpy(Third + (Width + 1) * Channel, Src + Stride + (Width - 1) * Channel, Channel);

    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePS = Src + Y * Stride;
        unsigned char *LinePD = Dest + Y * Stride;
        if (Y != 0)
        {
            unsigned char *Temp = First; First = Second; Second = Third; Third = Temp;
        }
        if (Y == Height - 1)
        {
            memcpy(Third, Second, (Width + 2) * Channel);
        }
        else
        {
            memcpy(Third, Src + (Y + 1) * Stride, Channel);
            memcpy(Third + Channel, Src + (Y + 1) * Stride, Width * Channel);                            //    由于备份了前面一行的数据,这里即使Src和Dest相同也是没有问题的
            memcpy(Third + (Width + 1) * Channel, Src + (Y + 1) * Stride + (Width - 1) * Channel, Channel);
        }
        if (Channel == 1)
        {
            for (int X = 0; X < Width; X++)
            {
                int GX = First[X] - First[X + 2] + (Second[X] - Second[X + 2]) * 2 + Third[X] - Third[X + 2];
                int GY = First[X] + First[X + 2] + (First[X + 1] - Third[X + 1]) * 2 - Third[X] - Third[X + 2];
                LinePD[X] = IM_ClampToByte(sqrtf(GX * GX + GY * GY + 0.0F));
            }
        }
        else
        {
            for (int X = 0; X < Width * 3; X++)
            {
                int GX = First[X] - First[X + 6] + (Second[X] - Second[X + 6]) * 2 + Third[X] - Third[X + 6];
                int GY = First[X] + First[X + 6] + (First[X + 3] - Third[X + 3]) * 2 - Third[X] - Third[X + 6];
                LinePD[X] = IM_ClampToByte(sqrtf(GX * GX + GY * GY + 0.0F));
            }
        }
    }
    free(RowCopy);
    return IM_STATUS_OK;
}

  代码很短,但是这段代码很经典,第一,这个代码支持In-Place操作,也就是Src和Dest可以是同一块内存,第二,这个代码本质就支持边缘。网络上很多参考代码都是只处理中间有效的区域。具体的实现细节我不愿意多述,自己看。

  那么Sobel的核心就是计算X方向的梯度GX和Y方向的梯度GY,最后有一个耗时的操作是求GX*GX+GY*GY的平方。

      上面这段代码,在不打开编译器的SSE优化和快速浮点计算的情况,直接使用FPU,对4000*3000的彩色图约需要480ms,当开启SSE后,大概时间为220ms ,因此系统编译器的SSE优化也很厉害,反编译后可以看到汇编里这样的部分:

59AD12F8  movd        xmm0,ecx  
59AD12FC  cvtdq2ps    xmm0,xmm0  
59AD12FF  sqrtss      xmm0,xmm0  
59AD1303  cvttss2si   ecx,xmm0  

 可见他是调用的单浮点数的sqrt优化。

    由于该Sobel的过程最后是把数据用图像的方式记录下来,因此,IM_ClampToByte(sqrtf(GX * GX + GY * GY + 0.0F))可以用查表的方式来实现。简单改成如下的版本, 避免了浮点计算。

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

    unsigned char *RowCopy = (unsigned char *)malloc((Width + 2) * 3 * Channel);
    if (RowCopy == NULL)    return IM_STATUS_OUTOFMEMORY;

    unsigned char *First = RowCopy;
    unsigned char *Second = RowCopy + (Width + 2) * Channel;
    unsigned char *Third = RowCopy + (Width + 2) * 2 * Channel;

    memcpy(Second, Src, Channel);
    memcpy(Second + Channel, Src, Width * Channel);                                                    //    拷贝数据到中间位置
    memcpy(Second + (Width + 1) * Channel, Src + (Width - 1) * Channel, Channel);

    memcpy(First, Second, (Width + 2) * Channel);                                                    //    第一行和第二行一样

    memcpy(Third, Src + Stride, Channel);                                                    //    拷贝第二行数据
    memcpy(Third + Channel, Src + Stride, Width * Channel);
    memcpy(Third + (Width + 1) * Channel, Src + Stride + (Width - 1) * Channel, Channel);

    int BlockSize = 16, Block = (Width * Channel) / BlockSize;

    unsigned char Table[65026];
    for (int Y = 0; Y < 65026; Y++)        Table[Y] = (sqrtf(Y + 0.0f) + 0.5f);

    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePS = Src + Y * Stride;
        unsigned char *LinePD = Dest + Y * Stride;
        if (Y != 0)
        {
            unsigned char *Temp = First; First = Second; Second = Third; Third = Temp;
        }
        if (Y == Height - 1)
        {
            memcpy(Third, Second, (Width + 2) * Channel);
        }
        else
        {
            memcpy(Third, Src + (Y + 1) * Stride, Channel);
            memcpy(Third + Channel, Src + (Y + 1) * Stride, Width * Channel);                            //    由于备份了前面一行的数据,这里即使Src和Dest相同也是没有问题的
            memcpy(Third + (Width + 1) * Channel, Src + (Y + 1) * Stride + (Width - 1) * Channel, Channel);
        }
        if (Channel == 1)
        {
            for (int X = 0; X < Width; X++)
            {
                int GX = First[X] - First[X + 2] + (Second[X] - Second[X + 2]) * 2 + Third[X] - Third[X + 2];
                int GY = First[X] + First[X + 2] + (First[X + 1] - Third[X + 1]) * 2 - Third[X] - Third[X + 2];
                LinePD[X] = Table[IM_Min(GX * GX + GY * GY, 65025)];
            }
        }
        else
        {
            for (int X = 0; X < Width * 3; X++)
            {
                int GX = First[X] - First[X + 6] + (Second[X] - Second[X + 6]) * 2 + Third[X] - Third[X + 6];
                int GY = First[X] + First[X + 6] + (First[X + 3] - Third[X + 3]) * 2 - Third[X] - Third[X + 6];
                LinePD[X] = Table[IM_Min(GX * GX + GY * GY, 65025)];
            }
        }
    }
    free(RowCopy);
    return IM_STATUS_OK;
}

对4000*3000的彩色图约需要180ms,比系统的SSE优化快了40ms,而这个过程完全无浮点计算,因此,可以知道计算GX和GY的耗时在本例中也占据了相当大的部分。

这样的过程最适合于SSE处理了。

我们分析之。

第一来看一看这两句:

  int GX = First[X] - First[X + 2] + (Second[X] - Second[X + 2]) * 2 + Third[X] - Third[X + 2];
  int GY = First[X] + First[X + 2] + (First[X + 1] - Third[X + 1]) * 2 - Third[X] - Third[X + 2];

 里面涉及到了8个不同的像素,考虑计算的特性和数据的范围,在内部计算时这个int可以用short代替,也就是要把加载的字节图像数据转换为short类型先,这样SSE优化方式为用8个SSE变量分别记录8个连续的像素风量的颜色值,每个颜色值用16位数据表达。

  这可以使用_mm_unpacklo_epi8配合_mm_loadl_epi64实现:

    __m128i FirstP0 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(First + X)), Zero);
    __m128i FirstP1 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(First + X + 3)), Zero);
    __m128i FirstP2 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(First + X + 6)), Zero);

    __m128i SecondP0 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(Second + X)), Zero);
    __m128i SecondP2 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(Second + X + 6)), Zero);

    __m128i ThirdP0 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(Third + X)), Zero);
    __m128i ThirdP1 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(Third + X + 3)), Zero);
    __m128i ThirdP2 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(Third + X + 6)), Zero);

  接着就是搬积木了,用SSE的指令代替普通的C的函数指令实现GX和GY的计算。

    __m128i GX16 = _mm_abs_epi16(_mm_add_epi16(_mm_add_epi16(_mm_sub_epi16(FirstP0, FirstP2), _mm_slli_epi16(_mm_sub_epi16(SecondP0, SecondP2), 1)), _mm_sub_epi16(ThirdP0, ThirdP2)));
    __m128i GY16 = _mm_abs_epi16(_mm_sub_epi16(_mm_add_epi16(_mm_add_epi16(FirstP0, FirstP2), _mm_slli_epi16(_mm_sub_epi16(FirstP1, ThirdP1), 1)), _mm_add_epi16(ThirdP0, ThirdP2)));

找个时候的GX16和GY16里保存的是8个16位的中间结果,由于SSE只提供了浮点数的sqrt操作,我们必须将它们转换为浮点数,那么这个转换的第一步就必须是先将它们转换为int的整形数,这样,就必须一个拆成2个,即:

    __m128i GX32L = _mm_unpacklo_epi16(GX16, Zero);
    __m128i GX32H = _mm_unpackhi_epi16(GX16, Zero);
    __m128i GY32L = _mm_unpacklo_epi16(GY16, Zero);
    __m128i GY32H = _mm_unpackhi_epi16(GY16, Zero);  

  接着又是搬积木了:

    __m128i ResultL = _mm_cvtps_epi32(_mm_sqrt_ps(_mm_cvtepi32_ps(_mm_add_epi32(_mm_mullo_epi32(GX32L, GX32L), _mm_mullo_epi32(GY32L, GY32L)))));
    __m128i ResultH = _mm_cvtps_epi32(_mm_sqrt_ps(_mm_cvtepi32_ps(_mm_add_epi32(_mm_mullo_epi32(GX32H, GX32H), _mm_mullo_epi32(GY32H, GY32H)))));

  整形转换为浮点数,最后计算完之后又要将浮点数转换为整形数。

  最后一步,得到8个int型的结果,这个结果有要转换为字节类型的,并且这些数据有可能会超出字节所能表达的范围,所以就需要用到SSE自带的抗饱和向下打包函数了,如下所示:

_mm_storel_epi64((__m128i *)(LinePD + X), _mm_packus_epi16(_mm_packus_epi32(ResultL, ResultH), Zero));

 Ok, 一切搞定了,还有一些细节处理自己慢慢补充吧。

运行,哇,只要37ms了,速度快了N倍,可结果似乎和其他方式实现的不一样啊,怎么回事。

我也是找了半天都没有找到问题所在,后来一步一步的测试,最终问题定位在16位转换为32位整形那里去了。

通常,我们都是对像素的字节数据进行向上扩展,他们都是正数,所以用unpack之类的配合zero把高8位或高16位的数据填充为0就可以了,但是在本例中,GX16或者GY16很有可能是负数,而负数的最高位是符号位,如果都填充为0,则变为正数了,明显改变原始的数据了,所以得到了错误的结果。

那如何解决问题呢,对于本例,很简单,因为后面只有一个平方操作,因此,对GX先取绝对值是不会改变计算的结果的,这样就不会出现负的数据了,修改之后,果然结果正确。

还可以继续优化,我们来看最后的计算GX*GX+GY*GY的过程,我们知道,SSE3提供了一个_mm_madd_epi16指令,其作用为:

r0 := (a0 * b0) + (a1 * b1)
r1 := (a2 * b2) + (a3 * b3)
r2 := (a4 * b4) + (a5 * b5)
r3 := (a6 * b6) + (a7 * b7)

  如果我们能把GX和GY的数据拼接成另外两个数据:

         GXYL   =   GX0  GY0  GX1  GY1  GX2  GY2  GX3  GY3

         GXYH   =   GX4  GY4  GX5  GY5  GX6  GY6  GX7  GY7

  那么调用_mm_madd_epi16(GXYL ,GXYL )和_mm_madd_epi16(GXYH ,GXYH )不就是能得到和之前一样的结果了吗,而这个拼接SSE已经有现成的函数了即:

__m128i GXYL = _mm_unpacklo_epi16(GX16, GY16);
__m128i GXYH = _mm_unpackhi_epi16(GX16, GY16);

  这样就把原来需要的10个指令变为了4个指令,代码简洁了并且速度也更快了。

尝试如此修改,整个的计算过程时间减少到了32ms左右。

另外,还有一个可以优化的地方就是借用  _mm_maddubs_epi16  函数实现像素之间的加减乘除和扩展。

这个函数的作用如下:

r0 := SATURATE_16((a0 * b0) + (a1 * b1))
r1 := SATURATE_16((a2 * b2) + (a3 * b3))
...
r7 := SATURATE_16((a14 * b14) + (a15 * b15))

  他的第一个参数是16个无符号的字节数据,第二个参数是16个有符号的char数据。

配合unpack使用类似上面的技术就可以一次性处理16个字节的像素简加减了,这样做整个过程大概能再加速2ms,达到最终的30ms。

源代码地址:http://files.cnblogs.com/files/Imageshop/Sobel.rar (其中的SSE代码请按照本文的思路自行整理。)

http://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar,这里是一个我全部用SSE优化的图像处理的Demo,有兴趣的朋友可以看看。

欢迎点赞和打赏。

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏林冠宏的技术文章

经典面试问题: Top K 之 ---- 海量数据找出现次数最多或,不重复的。

作者:林冠宏 / 指尖下的幽灵 掘金:https://juejin.im/user/587f0dfe128fe100570ce2d8 博客:htt...

49915
来自专栏iKcamp

翻译连载 | 附录 A:Transducing(下)-《JavaScript轻量级函数式编程》 |《你不知道的JS》姊妹篇

原文地址:Functional-Light-JS 原文作者:Kyle Simpson-《You-Dont-Know-JS》作者 JavaScript 轻量级函数...

1808
来自专栏Android机器圈

数据结构与算法----数学应用之一元多项式

PS:上一篇说了线性表的顺序表和链式表表达,该片就写一下应用到现实数学中去,一元多项式的加减。

382
来自专栏陈树义

BloomFilter算法

Bloom filter 是由 Howard Bloom 在 1970 年提出的二进制向量数据结构,它具有很好的空间和时间效率,被用来检测一个元素是不是集合中的...

3588
来自专栏C/C++基础

C++IO流简介

输入输出(IO)是指计算机同任何外部设备之间的数据传递。常见的输入输出设备有文件、键盘、打印机、屏幕等。数据可以按记录(或称数据块)的方式传递,也可以 流的方式...

913
来自专栏lonelydawn的前端猿区

lambda+reduce的一句艰深代码

一句话一脸懵逼 某天晚上看到一句lambda+reduce 组合的代码,看的头都炸了,愣是没看懂,不过也可能因为稀疏的脑神经经过一天的摧残已经运转不动了,这两天...

1778
来自专栏抠抠空间

python常见模块之time模块

一、模块简介 在Python中,通常有这几种方式来表示时间: 时间戳(timestamp):通常来说,时间戳表示的是从1970年1月1日00:00:00开始按...

3227
来自专栏用户2442861的专栏

数据库系统——B+树索引

http://blog.csdn.net/cjfeii/article/details/10858721

751
来自专栏swag code

抽象类与抽象方法

在我们抽象实例对象的时候,有这样一种情况,往上层抽象时就会发现很难描述对象的属性和行为,比如“形状” ,其方法计算面积怎么计算?正方形知道怎么计算,长方形也知道...

643
来自专栏信数据得永生

JavaScript 编程精解 中文第三版 五、高阶函数

27810

扫码关注云+社区