# 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);
_mm_unpacklo_epi8(el8shr2, v_zero));
_mm_unpacklo_epi8(el8shr3, v_zero));
_mm_unpacklo_epi16(el8, v_zero));
_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 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));
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优化过程吧。

高位            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位就可以直接使用。

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次这个函数，很明显这样做是不会改变计算结果的。

还有一点在测试时发现，如果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 条评论

## 相关文章

935

3236

1132

### CG实验5 简单光照明模型

(1) 示范代码为立方体在一束平行光照射下的漫反射光照效果。结合示范代码，学习掌握简单光照明模型的基本原理与实现； (2) 修改示范代码，给出不同光照参数...

1323

1872

841

30110

2508

2433

572