前言
ImageShop博主是我的Idol。我会逐渐出一些SSE优化一些图像领域的算法的文章，并且将代码开源到Github，希望可以为我多多star。https://github.com/BBuf/Image-processing-algorithm-Speed
算法原理
RGB转灰度图没啥说的了，这里的重点在于使用SSE优化，但为了和普通以及多线程优化的程序做一个速度对比这里也提供了一下对应的程序实现，我在实现这个代码的时候主要参考了ImageShop博主的这篇文章，https://www.cnblogs.com/Imageshop/p/6261719.html ，他只提供了一些代码片段，我依靠自己的逻辑和理解完整实现了整个代码并做了速度测试，我认为指令集加速无论是在图像领域还是在CV领域都是非常重要的，所以我也愿意将代码分享出来，接下来我们就开始今天的主题吧。
朴素的RGB转灰度算法
//1920 * 1080 3.71ms
void RGB2Y_1(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;
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;
}
}
}
这个我不用解释了，只说速度测试，在我的PC(core-i7-3770处理器的速度为3.95ms)。
四路循环展开算法
//1920 * 1080 2.0ms
void RGB2Y_2(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;
}
}
}
这个代码的速度为3.4ms处理一张1920*1080的图片。
使用SSE优化的算法
void RGB2Y_3(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 - 12; X += 12, LinePS += 36) {
__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)); //1
__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)); //2
__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)); //3
__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));
__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));
__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));
__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));
__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));
__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));
__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));
__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));
__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));
__m128i sumaL = _mm_add_epi16(p3aL, _mm_add_epi16(p1aL, p2aL));
__m128i sumaH = _mm_add_epi16(p3aH, _mm_add_epi16(p1aH, p2aH));
__m128i sumbL = _mm_add_epi16(p3bL, _mm_add_epi16(p1bL, p2bL));
__m128i sumbH = _mm_add_epi16(p3bH, _mm_add_epi16(p1bH, p2bH));
__m128i sclaL = _mm_srli_epi16(sumaL, 8);
__m128i sclaH = _mm_srli_epi16(sumaH, 8);
__m128i sclbL = _mm_srli_epi16(sumbL, 8);
__m128i sclbH = _mm_srli_epi16(sumbH, 8);
__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));
__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));
__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));
__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));
__m128i accumL = _mm_or_si128(shftaL, shftbL);
__m128i accumH = _mm_or_si128(shftaH, shftbH);
__m128i h3 = _mm_or_si128(accumL, accumH);
//__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));
_mm_storeu_si128((__m128i *)(LinePD + X), h3);
}
for (; X < Width; X++, LinePS += 3) {
LinePD[X] = (B_WT * LinePS[0] + G_WT * LinePS[1] + R_WT * LinePS[2]) >> 8;
}
}
}
这个算法速度是1.4ms，是普通算法(OpenCV自带CvtColor)的2-3倍速度。指令和原理讲解请看ImageShop大牛的博客：https://www.cnblogs.com/Imageshop/p/6261719.html 。我这里也尝试自己来描述一下：
算法原理
SSE优化的代码一次可以处理12个像素，每个像素有3个值(B, G, R)，所以一次相当于可以处理36个数。我们把这个用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能表达的最大数值，所以这里使用short作为计算对象。然后我们使用SSE指令读取数据然后进行一些相乘操作，例如这几行代码：
__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));
__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));
__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));
以及这行
__m128i sumaL = _mm_add_epi16(p3aL, _mm_add_epi16(p1aL, p2aL));
就实现了这样一个过程：
(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)
上面得到了一个sse的包含8个short类型数的向量，其中还黑色加粗部分就是本次我们要的第一个答案（一共处理12个像素，现在处理了3个像素，得到了第一个答案）。
指令集解释
_mm_loadu_si128就是把之后16个字节的数据读入到一个SSE寄存器中,注意由于任意位置的图像数据内存地址肯定不可能都满足SIMD16字节对齐的规定，因此这里不是用的_mm_load_si128指令。_mm_cvtepu8_epi16指令则把这16个字节的低64位的8个字节数扩展为8个16位数据，这样做主要是为了上述的乘法进行准备的。_mm_setr_epi16这个实际上就是用已知数的8个16位数据来构造一个SSE整型数。_mm_mullo_epi16指令就是两个16位的乘法，注意不是用的_mm_mulhi_epi16，因为两个16位数相乘，一般要用32位数才能完整的保存结果，而_mm_mullo_epi16是提取这个32位的低16位，我们这里前面已经明确了成绩的结果是不会超出short类型的，因此，所以只取低16位就已经完全保留了所有的信息。_mm_srli_epi16这个就是移位操作，相当于普通算法中的>>8。_mm_shuffle_epi8把保存数字的位提到前面去，这样方便我们后面进行合并为一个完整的sse变量。举个例子，这里第一个变量的位置为什么是0,6,12呢，因为最后计算得到的变量高位是没有信息的，我们只使用了低8位，而sse中的内存排布是大概是这样子，可以把黑色的看成地位，所以这里下表是0,6,12，而不是想当然的0,3,6。
- _mm_storeu_si128把处理的结果写入到目标内存中，注意，这里会多写了4个字节的内存数据（128 - 12 * 8），但是我们后面又会把他们重新覆盖掉，但是有一点要注意，就是如果是最后一行数据，在某些情况下超出的这个几个字节就已经不属于你这个进程该管理的范围了， 这个时候就会出现OOM错误，因此一种简单的方式就是在宽度方向上的循环终止条件设置为： X < Width - 12；这样剩余的像素用普通的算法处理即可避免这种问题出现。
剩下的可以自行查看SSE指令进行学习，也没什么好说的了。
京公网安备 11010502036488号