C++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独指令,代码简洁了还要速度为重新快了。

  C++尝试这样修改,整个的精打细算过程时间减少至了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,有趣味之对象可看看。

     
 C++ 1

  欢迎点赞和打赏。

   
  C++ 2