# SSE图像算法优化系列一：一段BGR2Y的SIMD代码解析。

一个同事在github上淘到一个基于SIMD的RGB转Y（彩色转灰度或者转明度）的代码，我抽了点时间看了下，顺便学习了一些SIMD指令，这里把学习过程中的一些理解和认识共享给大家。

我们首先说说普通的RGB2Y的代码：

```void RGB2Y(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride)
{
const int B_WT = int(0.114 * 256 + 0.5);
const int G_WT = int(0.587 * 256 + 0.5);
const int R_WT = 256 - B_WT - G_WT;            //     int(0.299 * 256 + 0.5);

for (int Y = 0; Y < Height; Y++)
{
unsigned char *LinePS = Src + Y * Stride;
unsigned char *LinePD = Dest + Y * Width;
for (int X = 0; X < Width; X++, LinePS += 3)
{
LinePD[X] = (B_WT * LinePS[0] + G_WT * LinePS[1] + R_WT * LinePS[2]) >> 8;
}
}
}```

上述代码的速度已经非常快了，在测试机上1920*1280的图像单次执行也只需要3.95ms左右，如果还需要优化，可以像下面这样模拟并行操作：

```void RGB2Y(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride)
{
const int B_WT = int(0.114 * 256 + 0.5);
const int G_WT = int(0.587 * 256 + 0.5);
const int R_WT = 256 - B_WT - G_WT;            //     int(0.299 * 256 + 0.5);

for (int Y = 0; Y < Height; Y++)
{
unsigned char *LinePS = Src + Y * Stride;
unsigned char *LinePD = Dest + Y * Width;
int X = 0;
for (; X < Width - 4; X += 4, LinePS += 12)
{
LinePD[X + 0] = (B_WT * LinePS[0] + G_WT * LinePS[1] + R_WT * LinePS[2]) >> 8;
LinePD[X + 1] = (B_WT * LinePS[3] + G_WT * LinePS[4] + R_WT * LinePS[5]) >> 8;
LinePD[X + 2] = (B_WT * LinePS[6] + G_WT * LinePS[7] + R_WT * LinePS[8]) >> 8;
LinePD[X + 3] = (B_WT * LinePS[9] + G_WT * LinePS[10] + R_WT * LinePS[11]) >> 8;
}
for (; X < Width; X++, LinePS += 3)
{
LinePD[X] = (B_WT * LinePS[0] + G_WT * LinePS[1] + R_WT * LinePS[2]) >> 8;
}
}
}```

我们来分析分析他提供的代码先，为了篇幅，这里仅仅贴出核心的我稍作修改的SIMD指令（SSE）：

``` 1     __m128i p1aL = _mm_mullo_epi16(_mm_cvtepu8_epi16(_mm_loadu_si128((__m128i *)(LinePS + 0))), _mm_setr_epi16(B_WT, G_WT, R_WT, B_WT, G_WT, R_WT, B_WT, G_WT));
2     __m128i p2aL = _mm_mullo_epi16(_mm_cvtepu8_epi16(_mm_loadu_si128((__m128i *)(LinePS + 1))), _mm_setr_epi16(G_WT, R_WT, B_WT, G_WT, R_WT, B_WT, G_WT, R_WT));
3     __m128i p3aL = _mm_mullo_epi16(_mm_cvtepu8_epi16(_mm_loadu_si128((__m128i *)(LinePS + 2))), _mm_setr_epi16(R_WT, B_WT, G_WT, R_WT, B_WT, G_WT, R_WT, B_WT));
4
5     __m128i p1aH = _mm_mullo_epi16(_mm_cvtepu8_epi16(_mm_loadu_si128((__m128i *)(LinePS + 8))), _mm_setr_epi16(R_WT, B_WT, G_WT, R_WT, B_WT, G_WT, R_WT, B_WT));
6     __m128i p2aH = _mm_mullo_epi16(_mm_cvtepu8_epi16(_mm_loadu_si128((__m128i *)(LinePS + 9))), _mm_setr_epi16(B_WT, G_WT, R_WT, B_WT, G_WT, R_WT, B_WT, G_WT));
7     __m128i p3aH = _mm_mullo_epi16(_mm_cvtepu8_epi16(_mm_loadu_si128((__m128i *)(LinePS + 10))), _mm_setr_epi16(G_WT, R_WT, B_WT, G_WT, R_WT, B_WT, G_WT, R_WT));
8
9     __m128i p1bL = _mm_mullo_epi16(_mm_cvtepu8_epi16(_mm_loadu_si128((__m128i *)(LinePS + 18))), _mm_setr_epi16(B_WT, G_WT, R_WT, B_WT, G_WT, R_WT, B_WT, G_WT));
10     __m128i p2bL = _mm_mullo_epi16(_mm_cvtepu8_epi16(_mm_loadu_si128((__m128i *)(LinePS + 19))), _mm_setr_epi16(G_WT, R_WT, B_WT, G_WT, R_WT, B_WT, G_WT, R_WT));
11     __m128i p3bL = _mm_mullo_epi16(_mm_cvtepu8_epi16(_mm_loadu_si128((__m128i *)(LinePS + 20))), _mm_setr_epi16(R_WT, B_WT, G_WT, R_WT, B_WT, G_WT, R_WT, B_WT));
12
13     __m128i p1bH = _mm_mullo_epi16(_mm_cvtepu8_epi16(_mm_loadu_si128((__m128i *)(LinePS + 26))), _mm_setr_epi16(R_WT, B_WT, G_WT, R_WT, B_WT, G_WT, R_WT, B_WT));
14     __m128i p2bH = _mm_mullo_epi16(_mm_cvtepu8_epi16(_mm_loadu_si128((__m128i *)(LinePS + 27))), _mm_setr_epi16(B_WT, G_WT, R_WT, B_WT, G_WT, R_WT, B_WT, G_WT));
15     __m128i p3bH = _mm_mullo_epi16(_mm_cvtepu8_epi16(_mm_loadu_si128((__m128i *)(LinePS + 28))), _mm_setr_epi16(G_WT, R_WT, B_WT, G_WT, R_WT, B_WT, G_WT, R_WT));
16
21     __m128i sclaL = _mm_srli_epi16(sumaL, 8);
22     __m128i sclaH = _mm_srli_epi16(sumaH, 8);
23     __m128i sclbL = _mm_srli_epi16(sumbL, 8);
24     __m128i sclbH = _mm_srli_epi16(sumbH, 8);
25     __m128i shftaL = _mm_shuffle_epi8(sclaL, _mm_setr_epi8(0, 6, 12, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1));
26     __m128i shftaH = _mm_shuffle_epi8(sclaH, _mm_setr_epi8(-1, -1, -1, 18, 24, 30, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1));
27     __m128i shftbL = _mm_shuffle_epi8(sclbL, _mm_setr_epi8(-1, -1, -1, -1, -1, -1, 0, 6, 12, -1, -1, -1, -1, -1, -1, -1));
28     __m128i shftbH = _mm_shuffle_epi8(sclbH, _mm_setr_epi8(-1, -1, -1, -1, -1, -1, -1, -1, -1, 18, 24, 30, -1, -1, -1, -1));
29     __m128i accumL = _mm_or_si128(shftaL, shftbL);
30     __m128i accumH = _mm_or_si128(shftaH, shftbH);
31     __m128i h3 = _mm_or_si128(accumL, accumH);
32     //__m128i h3 = _mm_blendv_epi8(accumL, accumH, _mm_setr_epi8(0, 0, 0, -1, -1, -1, 0, 0, 0, -1, -1, -1, 1, 1, 1, 1));
33     _mm_storeu_si128((__m128i *)(LinePD + X), h3);```

看上去篇幅好大，真的怀疑是否能比普通的C代码快，现给大家吃个定心丸吧，同样的机器和图片，以上述代码片段为核心的算法执行时间约为1.66ms, 比优化后的普通C代码还要快1倍多，好，先吃饭去了，早餐还没吃（12：06），咱们回来再继续分享。

好，吃了两个青菜，感觉不错，继续。

首先，代码一次性处理12个像素，我们用BGR序列表达出来如下:

B1  G1  R1  B2  G2  R2  B3  G3  R3  B4  G4  R4  B5  G5  R5  B6  G6  R6  B7  G7  R7  B8  G8  R8  B9  G9  R9  B10  G10  R10  B11  G11  R11  B12  G12  R12

SSE指令一次性能处理16个字节型数据，8个short类型的，或者4个int类型数据，考虑到计算过程中有乘法，综合效率和结果的准确性，我们选用short类型的作为主要的计算对象（这也就是说上述的权重定点化最大的放大范围为什么是256了，两个byte类型相乘不会超出short能表达的范围）。

那是如何利用SSE的批处理能力的呢，由于SSE的序列性，我们很难直接把图像数据相乘相加然后得到结果，上述代码的最大特点就是巧妙的从图像数据中构建了几个SSE的数据，然后在利用SSE指令进行一次性的处理，比如第一到第三行以及第17行就是实现了如下的功能：

(B1  G1  R1  B2  G2  R2  B3  G3)  x (B_WT  G_WT  R_WT  B_WT  G_WT  R_WT  B_WT  G_WT) +

(G1  R1  B2  G2  R2  B3  G3  R3)  x (G_WT  R_WT  B_WT  G_WT  R_WT  B_WT  G_WT  R_WT) +

(R1  B2  G2  R2  B3  G3  R3  B4)  x (R_WT  B_WT  G_WT  R_WT  B_WT  G_WT  R_WT  B_WT) =

(B1 x B_WT + G1 x G_WT + R1 x R_WT)    (G1 x G_WT + R1 x R_WT + B2 x B_WT)    (R1 x R_WT + B2 x B_WT + G2 x G_WT) +

(B2 x B_WT + G2 x G_WT + R2 x R_WT)    (G2 x G_WT + R2 x R_WT + B3 x B_WT)    (R2 x R_WT + B3 x B_WT + G3 x G_WT) +

(B3 x B_WT + G3 x G_WT + R3 x R_WT)    (G3 x G_WT + R3 x R_WT + B4 x B_WT)

注意到了上述式子中我用黑体字加粗的部分没有，这部分的结果就是我们需要获得的嘛，这是本SIMD代码的核心，好，现在一个函数一个函数的解释了吧。

_mm_setr_epi16这个实际上就是用已知数的8个16位数据来构造一个SSE整型数，仔细观察，这个代码的12个_mm_setr_epi16函数实际只有3个是不一样的，我曾经是这个把他们定义为一个全局的变量，在循环体内部直接使用，结果速度无区别。后面反汇编看看这些语句都便以为类似于这样的汇编码：

`5D48116C  movdqa      xmm7,xmmword ptr ds:[5D482110h]  `

ptr ds:[5D482110h] 这是个内存地址，也就是说他并不会在循环里每次都构造这个数据，而是直接从某个内存里读，这也就和我们在外部定义是一个意思了。

当然这主要是两个原因造成的，第一，我们这里的数据都是常数，因此每次循环内部不会变，编译器也是能够认识到这一点的，第二，由于这个SSE的变量比较多，已经大大的超过了SIMD内部的寄存器数量，因此，他需要用内存来缓存这些数据。

_mm_mullo_epi16 指令就是两个16位的乘法，注意不是用的_mm_mulhi_epi16，因为两个16位数相乘，一般要用32位数才能完整的保存结果，而_mm_mullo_epi16 是提取这个32位的低16位，我们这里前面已经明确了成绩的结果是不会超出short类型的，因此，所以只取低16位就已经完全保留了所有的信息。

第17和20行直接就是把每个元素相乘后的结果在相加，21和24行明显就是归一化的过程，进行移位，明显移位后的8个16位数只有低8位具有有效数字，高八位必然为0。

第25行到第32行其实把结果重新拼接到一个完整的SSE变量的过程，我们以第25句为例， shftaL 变量就正好记录了我们上面举例的那个结果，我们看到黑体加粗的结果是我们需要的，并且他应该位于真正结果的前3个字节的位置，因此，我们用_mm_shuffle_epi8指令把他们提到前面去，注意这个指令的前三个数字是0， 6， 12。为-1的部分都会变为0.

h3变量原代码是用_mm_blendv_epi8指令实现的，我觉得考虑这里的特殊性，用_mm_or_si128实现也是无妨的。

_mm_storeu_si128把处理的结果写入到目标内存中，注意，这里会多写了4个字节的内存数据（128 - 12 * 8），但是我们后面又会把他们重新覆盖掉，但是有一点要注意，就是如果是最后一行数据，在某些情况下超出的这个几个字节就已经不属于你这个进程该管理的范围了， 这个时候就会出现OOM错误，因此一种简单的方式就是在宽度方向上的循环终止条件设置为： X < Width - 12；这样剩余的像素用普通的算法处理即可避免这种问题出现。

上述代码中，一条SSE指令能同时执行8个short类型的计算，那为什么最后的提速只有1倍多一点呢，这其实很好解释，我们看到前面的计算中，计算出的8个累加值里只有3个是有效的，而其他的结果对我们来说毫无意义，并且在计算完之后，还有其他的一些合并操作，因此，最终只能获得这样的收益。

从个人理解来看，他这里一次性只处理了12个像素，其考虑的主要因素是最后一个_mm_blendv_epi8指令的方便，实际上稍作修改同时处理15个像素是可以的（小于16的最大的3的倍数）。

代码中还有一些其他的技巧，有些数字可能读者自己去看看那个指令的意义后会更加清晰，这些不太是语言能够完美表达清楚的。当然，这段代码在书写艺术上还是有很大的改良空间的，有兴趣的读者可以自行研究下。

我把这些代码稍做整理提供给大家使用和测试，我也相信，上述指令肯定不是最佳的SIMD实现方式，比如改变指令顺序让指令级又可以并行等手段也许也是有效的提速方式之一。

最后一点就是我有个疑问，在我提供的代码执行后，如果先使用SSE测试，后使用AVX测试，SSE的速度和上述报告数据差不多，但是一旦点了AVX测试后，在点SSE测试，SSE的速度就骤然下降很多，甚至比普通C的都要慢，我水平很有限，实在是不知道这是什么道理，烦请有知道的告知。

本笔记创建于2016年1月8日即将离开南京之际，特此纪念。

74 篇文章36 人订阅

0 条评论

## 相关文章

4116

### 基础野：细说浮点数

Brief　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　   本来只打算理解JS中0.1 + 0.2 == 0.30000000000000004...

1768

2916

2973

### 最短路（Floyd算法的动态规划本质）- HDU 2544

Floyd–Warshall（简称Floyd算法）是一种著名的解决任意两点间的最短路径（All Paris Shortest Paths，APS...

631

1162

3174

26410

### PhiloGL学习（2）——骚年，让我们荡起双桨

前言 上一篇文章中简单介绍了PhiloGL框架如何上手、GLSL语言以及简单的绘制一个方块（见PhiloGL学习（1）——场景创建及二维方块加载）。本文很简单...

3387

### Go语言 －浮点数

Go提供了两种size的浮点数，float32和float64。它们的算术规范是由IEEE754国际标准定义，现代CPU都实现了这个规范。 浮点数能够表示的范...

2894