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，有兴趣的朋友可以看看。

　　欢迎点赞和打赏。

时间： 2024-10-22 04:35:20

SSE再学习：灵活运用SIMD指令16倍提升Sobel边缘检测的速度（4000*3000的24位图像时间由480ms降低到30ms）。的相关文章

UNIX再学习 -- 文件I/O

在 UNIX/Linux 系统中,一切皆文件,这句话想必都有听过.对于文件的操作几乎适用于所有的设备,这也就看出了文件操作的重要性了.在C语言再学习部分有讲过标准I/O文件操作,参看:C语言再学习 -- 文件 下面我们来讲解下系统文件I/O的. 一.文件描述符 1.文件描述符简介 首先从文件描述符开始讲起.因为,对于内核而言,所有打开的文件都是通过文件描述符引用的.那么文件描述符到底是什么? 文件描述符(file descriptor)通常是一个小的非负整数,内核用以标识一个特定进程正在访问的文

基于8086CPU微处理器的汇编学习之MOV指令

汇编指令:MOV的作用是往某个寄存器中存入数值. 格式:mov  寄存器名,数值                数值-->寄存器 mov  寄存器A,存器寄B          B-->A PS:必须前后位数匹配,如: mov   ah,bx     ;error   ah is 8 bit,bx is 16 bit mov   ah, bh    ;right    ah and bh all is 8  bit mov   cx,dx     ;right     cx and dx al

学习AngularJs:Directive指令用法（完整版）

这篇文章主要学习AngularJs:Directive指令用法,内容很全面,感兴趣的小伙伴们可以参考一下 本教程使用AngularJs版本:1.5.3 AngularJs GitHub: https://github.com/angular/angular.js/ AngularJs下载地址:https://angularjs.org/ 摘要:Directive(指令)笔者认为是AngularJ非常强大而有有用的功能之一.它就相当于为我们写了公共的自定义DOM元素或LASS属性或ATTR属性,并

JDBC再学习

JDBC是规范,地球人都知道. 啥是规范呢?反正我说不好,真要让我说的话,就是SUN制订了一大堆接口,然后你要是想实现一些功能就要去实现这些接口,他要是也想要实现这些功能也得去老实儿的实现这些接口. JDBC就是这些接口们,java.sql包下面有好多个接口文件,这些接口文件就是所谓的规范,标准. 无论Oracle,MySql,还是DB2,SqlServer都实现了这些接口.这样一来我们只需要针对着jdk中的接口编程就可以了. 记得上学的时候,最讨厌的就是JDBC,因为就这里需要记一大串东西,第

Java反射再学习

在最初学习Java的时候觉得反射真的好难,并不是技术负责,而是思想复杂,无法接受.随着工作经验的增多,今日偶然间又看见某智的一个视频,感觉茅塞顿开.顺便在此系统整理一下反射的知识. 一言以蔽之:反射就是将Java类的各个组成部分转换为对应的Java对象. 我们知道,一切皆对象,那么这个“一切”必然也包含了Java类啊,Java类也是一种事物,那么他是什么的对象呢?毫无疑问,Java类是Class类的对象.(PS:那么Class类又是谁的对象呢?求大神指教?这问题貌似无穷无尽啊 %>_<% )

PHP再学习5——RESTFul框架 远程控制LED

0.前言 去年(2013年)2月第一次接触yeelink平台,当时该平台已经运行了一些时间也吸引了不少极客.试想自己也将投身IoT(物联网)行业,就花了些时间研究了它.陆陆续续使用和研究了一年,大致围绕两个问题展开——1.yeelink平台如何使用,2.如何构造一个功能简单些的yeelink平台.    [PHP学习笔记——索引博文] 本文将讨论如何构造一个简单restful架构平台(该平台有点像yeelink,不过功能比yeelink少的多),并结合树莓派实现LED的远程控制(网络控制).构建

UNIX再学习 -- 环境变量

之前讲gcc编译的时候,参看:C语言再学习 -- GCC编译过程 提到过静态库和共享库,那时只是简单的讲了下它们相关的编译链接,接下来就该详细介绍它们了.不过再讲解之前还需了解一下编程相关的环境变量. 一.环境变量 参看:百度百科--环境变量 环境变量一般是指在操作系统中用来指定操作系统运行环境的一些参数,如:临时文件夹位置和系统文件夹位置等. 环境变量时在操作系统中一个具有特定名字的对象,它包含了一个或者多个应用程序所将使用的信息. 1.Windows下的环境变量 (1)环境变量配置 右击我的

vue2.0学习(一)-内部指令

vue2.0学习(一)-内部指令 GitHub源码 https://github.com/sunnyguyan/VueDemo 1.v-if v-else v-show v-if:是vue的一个内部指令,指令用在我们的html中 v-if用来判断是否加载html的DOM,比如我们模拟一个用户登录状态,在用户登录后现实用户名称 <div v-if='isLogin'>你好,guyan</div> 完整的html代码 <!DOCTYPE html> <html lan

深度学习 GPU环境 Ubuntu 16.04 + Nvidia GTX 1080 + Python 3.6 + CUDA 9.0 + cuDNN 7.1 + TensorFlow 1.6 环境配置

本节详细说明一下深度学习环境配置,Ubuntu 16.04 + Nvidia GTX 1080 + Python 3.6 + CUDA 9.0 + cuDNN 7.1 + TensorFlow 1.6. Python 3.6 首先安装 Python 3.6,这里使用 Anaconda 3 来安装,下载地址:https://www.anaconda.com/download/#linux,点击 Download 按钮下载即可,这里下载的是 Anaconda 3-5.1 版本,如果下载速度过慢可以选