SSE图像算法优化系列六:OpenCv关于灰度积分图的SSE代码学习和改进。

  最近一直沉迷于SSE方面的优化,实在找不到想学习的参考资料了,就拿个笔记本放在腿上翻翻OpenCv的源代码,无意中看到了OpenCv中关于积分图的代码,仔细研习了一番,觉得OpenCv对SSE的灵活运用真的做的很好,这里记录下我对该段代码的品味并将其思路扩展到其他通道数的图像。

     该核心代码位于:Opencv 3.0\opencv\sources\modules\imgproc\src\sumpixels.cpp文件中。

     我们贴出最感兴趣的一部分代码以便分析:

    bool operator()(const uchar * src, size_t _srcstep,int * sum, size_t _sumstep,double * sqsum, size_t, int * tilted, size_t,Size size, int cn) const
    {
        if (sqsum || tilted || cn != 1 || !haveSSE2) return false;
        // the first iteration
        memset(sum, 0, (size.width + 1) * sizeof(int));
        __m128i v_zero = _mm_setzero_si128(), prev = v_zero;
        int j = 0;
        // the others
        for (int i = 0; i < size.height; ++i)
        {
            const uchar * src_row = src + _srcstep * i;
            int * prev_sum_row = (int *)((uchar *)sum + _sumstep * i) + 1;
            int * sum_row = (int *)((uchar *)sum + _sumstep * (i + 1)) + 1;
            sum_row[-1] = 0;
            prev = v_zero;
            j = 0;
            for ( ; j + 7 < size.width; j += 8)
            {
                __m128i vsuml = _mm_loadu_si128((const __m128i *)(prev_sum_row + j));
                __m128i vsumh = _mm_loadu_si128((const __m128i *)(prev_sum_row + j + 4));
                __m128i el8shr0 = _mm_loadl_epi64((const __m128i *)(src_row + j));
                __m128i el8shr1 = _mm_slli_si128(el8shr0, 1);
                __m128i el8shr2 = _mm_slli_si128(el8shr0, 2);
                __m128i el8shr3 = _mm_slli_si128(el8shr0, 3);
                vsuml = _mm_add_epi32(vsuml, prev);
                vsumh = _mm_add_epi32(vsumh, prev);
                __m128i el8shr12 = _mm_add_epi16(_mm_unpacklo_epi8(el8shr1, v_zero),
                                                 _mm_unpacklo_epi8(el8shr2, v_zero));
                __m128i el8shr03 = _mm_add_epi16(_mm_unpacklo_epi8(el8shr0, v_zero),
                                                 _mm_unpacklo_epi8(el8shr3, v_zero));
                __m128i el8 = _mm_add_epi16(el8shr12, el8shr03);
                __m128i el4h = _mm_add_epi16(_mm_unpackhi_epi16(el8, v_zero),
                                             _mm_unpacklo_epi16(el8, v_zero));
                vsuml = _mm_add_epi32(vsuml, _mm_unpacklo_epi16(el8, v_zero));
                vsumh = _mm_add_epi32(vsumh, el4h);
                _mm_storeu_si128((__m128i *)(sum_row + j), vsuml);
                _mm_storeu_si128((__m128i *)(sum_row + j + 4), vsumh);
                prev = _mm_add_epi32(prev, _mm_shuffle_epi32(el4h, _MM_SHUFFLE(3, 3, 3, 3)));
            }
            for (int v = sum_row[j - 1] - prev_sum_row[j - 1]; j < size.width; ++j)
                sum_row[j] = (v += src_row[j]) + prev_sum_row[j];
        }

     为了说明更方便,这里贴出我做的普通C语言的代码和重新优化后的SSE代码。

     普通C语言:

 void GetGrayIntegralImage(unsigned char *Src, int *Integral, int Width, int Height, int Stride)
 {
      memset(Integral, 0, (Width + 1) * sizeof(int));                    //    第一行都为0
      for (int Y = 0; Y < Height; Y++)
      {
          unsigned char *LinePS = Src + Y * Stride;
          int *LinePL = Integral + Y * (Width + 1) + 1;                 //    上一行位置            
          int *LinePD = Integral + (Y + 1) * (Width + 1) + 1;           //    当前位置,注意每行的第一列的值都为0
          LinePD[-1] = 0;                                               //    第一列的值为0
          for (int X = 0, Sum = 0; X < Width; X++)
          {
             Sum += LinePS[X];                                          //    行方向累加
             LinePD[X] = LinePL[X] + Sum;                               //    更新积分图
          }
     }
}

       优化后的SSE算法:

void GetGrayIntegralImage(unsigned char *Src, int *Integral, int Width, int Height, int Stride)
{
    memset(Integral, 0, (Width + 1) * sizeof(int));            //    第一行都为0
    int BlockSize = 8, Block = Width / BlockSize;
    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePS = Src + Y * Stride;
        int *LinePL = Integral + Y * (Width + 1) + 1;                //    上一行位置            
        int *LinePD = Integral + (Y + 1) * (Width + 1) + 1;          //    当前位置,注意每行的第一列的值都为0
        LinePD[-1] = 0;
        __m128i PreV = _mm_setzero_si128();
        __m128i Zero = _mm_setzero_si128();
        for (int X = 0; X < Block * BlockSize; X += BlockSize)
        {
            __m128i Src_Shift0 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(LinePS + X)), Zero);        //    A7 A6 A5 A4 A3 A2 A1 A0
            __m128i Src_Shift1 = _mm_slli_si128(Src_Shift0, 2);                                            //    A6 A5 A4 A3 A2 A1 A0 0     
            __m128i Src_Shift2 = _mm_slli_si128(Src_Shift1, 2);    //    移位改成基于Shift0,速度慢,Why?    //    A5 A4 A3 A2 A1 A0 0  0
            __m128i Src_Shift3 = _mm_slli_si128(Src_Shift2, 2);                                            //    A4 A3 A2 A1 A0 0  0  0
            __m128i Shift_Add12 = _mm_add_epi16(Src_Shift1, Src_Shift2);                                   //    A6+A5 A5+A4 A4+A3 A3+A2 A2+A1 A1+A0 A0+0  0+0
            __m128i Shift_Add03 = _mm_add_epi16(Src_Shift0, Src_Shift3);                                   //    A7+A4 A6+A3 A5+A2 A4+A1 A3+A0 A2+0  A1+0  A0+0    
            __m128i Low = _mm_add_epi16(Shift_Add12, Shift_Add03);                                         //    A7+A6+A5+A4 A6+A5+A4+A3 A5+A4+A3+A2 A4+A3+A2+A1 A3+A2+A1+A0 A2+A1+A0+0 A1+A0+0+0 A0+0+0+0
            __m128i High = _mm_add_epi32(_mm_unpackhi_epi16(Low, Zero), _mm_unpacklo_epi16(Low, Zero));    //    A7+A6+A5+A4+A3+A2+A1+A0  A6+A5+A4+A3+A2+A1+A0  A5+A4+A3+A2+A1+A0  A4+A3+A2+A1+A0
            __m128i SumL = _mm_loadu_si128((__m128i *)(LinePL + X + 0));
            __m128i SumH = _mm_loadu_si128((__m128i *)(LinePL + X + 4));
            SumL = _mm_add_epi32(SumL, PreV);
            SumL = _mm_add_epi32(SumL, _mm_unpacklo_epi16(Low, Zero));
            SumH = _mm_add_epi32(SumH, PreV);
            SumH = _mm_add_epi32(SumH, High);
            PreV = _mm_add_epi32(PreV, _mm_shuffle_epi32(High, _MM_SHUFFLE(3, 3, 3, 3)));
            _mm_storeu_si128((__m128i *)(LinePD + X + 0), SumL);
            _mm_storeu_si128((__m128i *)(LinePD + X + 4), SumH);
        }
        for (int X = Block * BlockSize, V = LinePD[X - 1] - LinePL[X - 1]; X < Width; X++)
        {
            V += LinePS[X];
            LinePD[X] = V + LinePL[X];
        }   }

  我们先来解释下这段代码的SSE优化过程吧。

 首先,用_mm_loadl_epi64一次性加载8个字节数据到XMM寄存器中,其中寄存器的高8位位0,此时寄存器的数据为:

      高位            0  0  0  0  0  0  0  0 A7 A6 A5 A4 A3 A2 A1 A0        低位   (8位)

     因为涉及到加法,并且最大为8个字节数据的加法,因此转换到16位数据类型,使用_mm_unpacklo_epi8结合zero即可实现。

     此时XMM寄存器内容变为:

           Src_Shift0    A7 A6 A5 A4 A3 A2 A1 A0    (16位)

     此后有3次移位分别得到:

            Src_Shift1    A6 A5 A4 A3 A2 A1 A0 0       (16位)
            Src_Shift2    A5 A4 A3 A2 A1 A0 0  0     (16位)
            Src_Shift3    A4 A3 A2 A1 A0 0  0  0         (16位)  通过_mm_add_epi16分别对4组16位数据进行8次相加:
            Shift_Add12   A6+A5 A5+A4 A4+A3 A3+A2 A2+A1 A1+A0 A0+0  0+0   (16位)
            Shift_Add03   A7+A4 A6+A3 A5+A2 A4+A1 A3+A0 A2+0  A1+0  A0+0   (16位)  
  再对他们进行相加:
        Low            A7+A6+A5+A4 A6+A5+A4+A3 A5+A4+A3+A2 A4+A3+A2+A1 A3+A2+A1+A0 A2+A1+A0+0 A1+A0+0+0 A0+0+0+0

     注意到低4位的16位数已经是连续相加的数据了,只要将他们转换为32位就可以直接使用。

而通过 __m128i High = _mm_add_epi32(_mm_unpackhi_epi16(Low, Zero), _mm_unpacklo_epi16(Low, Zero)); 这一句则可以把前面的高4位连续相加的值拼接起来得到:

       High                  A7+A6+A5+A4+A3+A2+A1+A0  A6+A5+A4+A3+A2+A1+A0  A5+A4+A3+A2+A1+A0  A4+A3+A2+A1+A0

  后面的操作则顺理成章了。

       注意到我核心的改动在于原始代码中的el8shr12和el8shr03的计算中的_mm_unpacklo_epi8被消除了,而在el8shr0一句中增加了一个_mm_unpacklo_epi8,因此少了3次这个函数,很明显这样做是不会改变计算结果的。

       另外源代码中的部分_mm_add_epi16被我用_mm_add_epi32代替了,这主要是因为用_mm_add_epi32意义更明显,而且由于高位数据为0,他们的执行结果不会有任何区别。

  还有一点在测试时发现,如果Src_Shift2,Src_Shift3的移位是基于Src_Shift0,即使用如下代码:

__m128i Src_Shift2 = _mm_slli_si128(Src_Shift0, 4);    
__m128i Src_Shift3 = _mm_slli_si128(Src_Shift0, 6);

  速度会有较为明显的下降,难道说移动的位数多少和CPU的耗时有关?

      以上是灰度模式的算法,在我的笔记本电脑上,SSE优化后的语句虽然增加了很多,但是执行效率约能提升30%,不过在一些PC上,普通的C和SSE优化后却没有啥速度区别了,这也不知道是为什么了。

      如果是针对24位或者32位图像,基本的优化思想是一致的,不过有更多的细节需要自己注意。

      24位或者32位图像在任何机器配置上,速度都能有30%的提升的。

      还是感觉这种算法用文字很难表述清楚,用代码再加上自己的空间组合可能更能理解吧。

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏Android小菜鸡

Android语音录制,语音发送

这是一个录音的例子,可用于IM的语音发送,OA的语音留言等。 首先我们需要引入权限:

1572
来自专栏木头编程 - moTzxx

PHP算法 [杨辉三角的求解]

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/u011415782/article/de...

922
来自专栏深度学习之tensorflow实战篇

R语言函数的含义与用法,实现过程解读

R的源起 R是S语言的一种实现。S语言是由 AT&T贝尔实验室开发的一种用来进行数据探索、统计分析、作图的解释型语言。最初S语言的实现版本主要是S-PLUS。S...

36312
来自专栏数据科学学习手札

(数据科学学习手札42)folium进阶内容介绍

  在上一篇(数据科学学习手札41)中我们了解了folium的基础内容,实际上folium在地理信息可视化上的真正过人之处在于其绘制图像的高度可定制化上,本文就...

2354
来自专栏一心无二用,本人只专注于基础图像算法的实现与优化。

一个简单的统计图像主颜色的算法(C#源代码)

      前段日子有朋友咨询了下分析图像主颜色的算法,我对这一块也没有什么深入的研究,参考了一些小代码,然后自己写了一个很简单的小工具,现共享给大家。    ...

2505
来自专栏天天

background、转换、过渡

743
来自专栏程序你好

CSharp代码示例每日一讲: 在GDI+中使用画笔和画刷

701
来自专栏Java成神之路

JavaUtil_02_二维码的生成与解析

 zxing-core-1.7.jar  :   http://viralpatel.net/blogs/download/jar/zxing-core-1.7...

672
来自专栏Android常用基础

自定义View(四)-动画- Interpolator与Evaluator

Interpolator插值器之前我们已经接触过了,而Evaluator好像我们还没有将,这是属性动画中俩个比较中的两个知识点,弄清楚它们有助于我们更好的使用与...

1062
来自专栏漫漫深度学习路

pytorch学习笔记(十二):详解 Module 类

Module 是 pytorch 提供的一个基类,每次我们要 搭建 自己的神经网络的时候都要继承这个类,继承这个类会使得我们 搭建网络的过程变得异常简单。 本文...

3266

扫码关注云+社区