SSE图像算法优化系列十五：YUV/XYZ和RGB空间相互转化的极速实现(此后老板不用再担心算法转到其他空间通道的耗时了)。
在颜色空间系列1: RGB和CIEXYZ颜色空间的转换及相关优化和颜色空间系列3: RGB和YUV颜色空间的转换及优化算法两篇文章中我们给出了两种不同的颜色空间的相互转换之间的快速算法的实现代码，但是那个是C#版本的，为了比较方便，我们这里提供C版本的代码，以RGB转到YUV空间的代码为例：
void RGBToYUV(unsigned char *RGB, unsigned char *Y, unsigned char *U, unsigned char *V, int Width, int Height, int Stride) { const int Shift = 15; const int HalfV = 1 << (Shift - 1); const int Y_B_WT = 0.114f * (1 << Shift), Y_G_WT = 0.587f * (1 << Shift), Y_R_WT = (1 << Shift) - Y_B_WT - Y_G_WT; const int U_B_WT = 0.436f * (1 << Shift), U_G_WT = -0.28886f * (1 << Shift), U_R_WT = -(U_B_WT + U_G_WT); const int V_B_WT = -0.10001 * (1 << Shift), V_G_WT = -0.51499f * (1 << Shift), V_R_WT = -(V_B_WT + V_G_WT); for (int YY = 0; YY < Height; YY++) { unsigned char *LinePS = RGB + YY * Stride; unsigned char *LinePY = Y + YY * Width; unsigned char *LinePU = U + YY * Width; unsigned char *LinePV = V + YY * Width; for (int XX = 0; XX < Width; XX++, LinePS += 3) { int Blue = LinePS[0], Green = LinePS[1], Red = LinePS[2]; LinePY[XX] = (Y_B_WT * Blue + Y_G_WT * Green + Y_R_WT * Red + HalfV) >> Shift; LinePU[XX] = ((U_B_WT * Blue + U_G_WT * Green + U_R_WT * Red + HalfV) >> Shift) + 128; LinePV[XX] = ((V_B_WT * Blue + V_G_WT * Green + V_R_WT * Red + HalfV) >> Shift) + 128; } } }
上述代码和颜色空间系列3: RGB和YUV颜色空间的转换及优化算法中的有所不同，但是应该说更加合理，注意Y_R_WT/U_R_WT/V_R_WT 的书写方式有所不同，这主要是为了保证定点化后的系数不会放大误差，同时注意计算时每一项还增加了HalfV值，这主要是为了保证对计算结果的四舍五入。
现代的CPU都很强劲，找了一台普通的I5电脑测试，1080P的图平均转换时间为10ms，大概能得到100fps的转换速率，但是从实际应用来讲，在很多场合这个耗时还是多了点，如果需要处理1080P这样的高清视频，考虑到综合因素，单帧的总处理时间不宜超过30ms,所以还有必要进一步提高这个算法的速度。
凭直觉，上述代码用GPU实现应该有很高的并行性，不知道最后速度会如何。
俺不会GPU，只会用SSE进行优化，先试一试把，一种最直接最朴实的写法如下：
void RGBToYUVSSE(unsigned char *RGB, unsigned char *Y, unsigned char *U, unsigned char *V, int Width, int Height, int Stride) { const int Shift = 15; const int HalfV = 1 << (Shift - 1); const int Y_B_WT = 0.114f * (1 << Shift), Y_G_WT = 0.587f * (1 << Shift), Y_R_WT = (1 << Shift) - Y_B_WT - Y_G_WT; const int U_B_WT = 0.436f * (1 << Shift), U_G_WT = -0.28886f * (1 << Shift), U_R_WT = -(U_B_WT + U_G_WT); const int V_B_WT = -0.10001 * (1 << Shift), V_G_WT = -0.51499f * (1 << Shift), V_R_WT = -(V_B_WT + V_G_WT); __m128i Weight_YB = _mm_set1_epi32(Y_B_WT), Weight_YG = _mm_set1_epi32(Y_G_WT), Weight_YR = _mm_set1_epi32(Y_R_WT); __m128i Weight_UB = _mm_set1_epi32(U_B_WT), Weight_UG = _mm_set1_epi32(U_G_WT), Weight_UR = _mm_set1_epi32(U_R_WT); __m128i Weight_VB = _mm_set1_epi32(V_B_WT), Weight_VG = _mm_set1_epi32(V_G_WT), Weight_VR = _mm_set1_epi32(V_R_WT); __m128i C128 = _mm_set1_epi32(128); __m128i Half = _mm_set1_epi32(HalfV); __m128i Zero = _mm_setzero_si128(); const int BlockSize = 16, Block = Width / BlockSize; for (int YY = 0; YY < Height; YY++) { unsigned char *LinePS = RGB + YY * Stride; unsigned char *LinePY = Y + YY * Width; unsigned char *LinePU = U + YY * Width; unsigned char *LinePV = V + YY * Width; for (int XX = 0; XX < Block * BlockSize; XX += BlockSize, LinePS += BlockSize * 3) { __m128i Src1, Src2, Src3, Blue, Green, Red; Src1 = _mm_loadu_si128((__m128i *)(LinePS + 0)); Src2 = _mm_loadu_si128((__m128i *)(LinePS + 16)); Src3 = _mm_loadu_si128((__m128i *)(LinePS + 32)); Blue = _mm_shuffle_epi8(Src1, _mm_setr_epi8(0, 3, 6, 9, 12, 15, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1)); Blue = _mm_or_si128(Blue, _mm_shuffle_epi8(Src2, _mm_setr_epi8(-1, -1, -1, -1, -1, -1, 2, 5, 8, 11, 14, -1, -1, -1, -1, -1))); Blue = _mm_or_si128(Blue, _mm_shuffle_epi8(Src3, _mm_setr_epi8(-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 1, 4, 7, 10, 13))); Green = _mm_shuffle_epi8(Src1, _mm_setr_epi8(1, 4, 7, 10, 13, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1)); Green = _mm_or_si128(Green, _mm_shuffle_epi8(Src2, _mm_setr_epi8(-1, -1, -1, -1, -1, 0, 3, 6, 9, 12, 15, -1, -1, -1, -1, -1))); Green = _mm_or_si128(Green, _mm_shuffle_epi8(Src3, _mm_setr_epi8(-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 2, 5, 8, 11, 14))); Red = _mm_shuffle_epi8(Src1, _mm_setr_epi8(2, 5, 8, 11, 14, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1)); Red = _mm_or_si128(Red, _mm_shuffle_epi8(Src2, _mm_setr_epi8(-1, -1, -1, -1, -1, 1, 4, 7, 10, 13, -1, -1, -1, -1, -1, -1))); Red = _mm_or_si128(Red, _mm_shuffle_epi8(Src3, _mm_setr_epi8(-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 0, 3, 6, 9, 12, 15))); // 以上操作把16个连续像素的像素顺序由 B G R B G R B G R B G R B G R B G R B G R B G R B G R B G R B G R B G R B G R B G R B G R B G R // 更改为适合于SIMD指令处理的连续序列 B B B B B B B B B B B B B B B B G G G G G G G G G G G G G G G G R R R R R R R R R R R R R R R R // 并分别保存于三个SSE变量中，并且没任何的冗余数据。 __m128i Blue16L = _mm_unpacklo_epi8(Blue, Zero); __m128i Blue16H = _mm_unpackhi_epi8(Blue, Zero); __m128i Blue32LL = _mm_unpacklo_epi16(Blue16L, Zero); __m128i Blue32LH = _mm_unpackhi_epi16(Blue16L, Zero); __m128i Blue32HL = _mm_unpacklo_epi16(Blue16H, Zero); __m128i Blue32HH = _mm_unpackhi_epi16(Blue16H, Zero); __m128i Green16L = _mm_unpacklo_epi8(Green, Zero); __m128i Green16H = _mm_unpackhi_epi8(Green, Zero); __m128i Green32LL = _mm_unpacklo_epi16(Green16L, Zero); __m128i Green32LH = _mm_unpackhi_epi16(Green16L, Zero); __m128i Green32HL = _mm_unpacklo_epi16(Green16H, Zero); __m128i Green32HH = _mm_unpackhi_epi16(Green16H, Zero); __m128i Red16L = _mm_unpacklo_epi8(Red, Zero); __m128i Red16H = _mm_unpackhi_epi8(Red, Zero); __m128i Red32LL = _mm_unpacklo_epi16(Red16L, Zero); __m128i Red32LH = _mm_unpackhi_epi16(Red16L, Zero); __m128i Red32HL = _mm_unpacklo_epi16(Red16H, Zero); __m128i Red32HH = _mm_unpackhi_epi16(Red16H, Zero); // 以上操作把三个SSE变量里的字节数据分别提取到12个包含4个int类型的数据的SSE变量里，以便后续的乘积操作不溢出 __m128i LL_Y = _mm_srai_epi32(_mm_add_epi32(_mm_add_epi32(_mm_mullo_epi32(Blue32LL, Weight_YB), _mm_add_epi32(_mm_mullo_epi32(Green32LL, Weight_YG), _mm_mullo_epi32(Red32LL, Weight_YR))),Half), Shift); __m128i LH_Y = _mm_srai_epi32(_mm_add_epi32(_mm_add_epi32(_mm_mullo_epi32(Blue32LH, Weight_YB), _mm_add_epi32(_mm_mullo_epi32(Green32LH, Weight_YG), _mm_mullo_epi32(Red32LH, Weight_YR))), Half), Shift); __m128i HL_Y = _mm_srai_epi32(_mm_add_epi32(_mm_add_epi32(_mm_mullo_epi32(Blue32HL, Weight_YB), _mm_add_epi32(_mm_mullo_epi32(Green32HL, Weight_YG), _mm_mullo_epi32(Red32HL, Weight_YR))), Half), Shift); __m128i HH_Y = _mm_srai_epi32(_mm_add_epi32(_mm_add_epi32(_mm_mullo_epi32(Blue32HH, Weight_YB), _mm_add_epi32(_mm_mullo_epi32(Green32HH, Weight_YG), _mm_mullo_epi32(Red32HH, Weight_YR))), Half), Shift); _mm_storeu_si128((__m128i*)(LinePY + XX), _mm_packus_epi16(_mm_packus_epi32(LL_Y, LH_Y), _mm_packus_epi32(HL_Y, HH_Y))); // 4个包含4个int类型的SSE变量重新打包为1个包含16个字节数据的SSE变量 // 以上操作完成：Y[0 - 15] = (Y_B_WT * Blue[0 - 15]+ Y_G_WT * Green[0 - 15] + Y_R_WT * Red[0 - 15] + HalfV) >> Shift; __m128i LL_U = _mm_add_epi32(_mm_srai_epi32(_mm_add_epi32(_mm_add_epi32(_mm_mullo_epi32(Blue32LL, Weight_UB), _mm_add_epi32(_mm_mullo_epi32(Green32LL, Weight_UG), _mm_mullo_epi32(Red32LL, Weight_UR))), Half), Shift), C128); __m128i LH_U = _mm_add_epi32(_mm_srai_epi32(_mm_add_epi32(_mm_add_epi32(_mm_mullo_epi32(Blue32LH, Weight_UB), _mm_add_epi32(_mm_mullo_epi32(Green32LH, Weight_UG), _mm_mullo_epi32(Red32LH, Weight_UR))), Half), Shift), C128); __m128i HL_U = _mm_add_epi32(_mm_srai_epi32(_mm_add_epi32(_mm_add_epi32(_mm_mullo_epi32(Blue32HL, Weight_UB), _mm_add_epi32(_mm_mullo_epi32(Green32HL, Weight_UG), _mm_mullo_epi32(Red32HL, Weight_UR))), Half), Shift), C128); __m128i HH_U = _mm_add_epi32(_mm_srai_epi32(_mm_add_epi32(_mm_add_epi32(_mm_mullo_epi32(Blue32HH, Weight_UB), _mm_add_epi32(_mm_mullo_epi32(Green32HH, Weight_UG), _mm_mullo_epi32(Red32HH, Weight_UR))), Half), Shift), C128); _mm_storeu_si128((__m128i*)(LinePU + XX), _mm_packus_epi16(_mm_packus_epi32(LL_U, LH_U), _mm_packus_epi32(HL_U, HH_U))); // 以上操作完成：U[0 - 15] = ((U_B_WT * Blue[0 - 15]+ U_G_WT * Green[0 - 15] + U_R_WT * Red[0 - 15] + HalfV) >> Shift) + 128; __m128i LL_V = _mm_add_epi32(_mm_srai_epi32(_mm_add_epi32(_mm_add_epi32(_mm_mullo_epi32(Blue32LL, Weight_VB), _mm_add_epi32(_mm_mullo_epi32(Green32LL, Weight_VG), _mm_mullo_epi32(Red32LL, Weight_VR))), Half), Shift), C128); __m128i LH_V = _mm_add_epi32(_mm_srai_epi32(_mm_add_epi32(_mm_add_epi32(_mm_mullo_epi32(Blue32LH, Weight_VB), _mm_add_epi32(_mm_mullo_epi32(Green32LH, Weight_VG), _mm_mullo_epi32(Red32LH, Weight_VR))), Half), Shift), C128); __m128i HL_V = _mm_add_epi32(_mm_srai_epi32(_mm_add_epi32(_mm_add_epi32(_mm_mullo_epi32(Blue32HL, Weight_VB), _mm_add_epi32(_mm_mullo_epi32(Green32HL, Weight_VG), _mm_mullo_epi32(Red32HL, Weight_VR))), Half), Shift), C128); __m128i HH_V = _mm_add_epi32(_mm_srai_epi32(_mm_add_epi32(_mm_add_epi32(_mm_mullo_epi32(Blue32HH, Weight_VB), _mm_add_epi32(_mm_mullo_epi32(Green32HH, Weight_VG), _mm_mullo_epi32(Red32HH, Weight_VR))), Half), Shift), C128); _mm_storeu_si128((__m128i*)(LinePV + XX), _mm_packus_epi16(_mm_packus_epi32(LL_V, LH_V), _mm_packus_epi32(HL_V, HH_V))); // 以上操作完成：V[0 - 15] = ((V_B_WT * Blue[0 - 15]+ V_G_WT * Green[0 - 15] + V_R_WT * Red[0 - 15] + HalfV) >> Shift) + 128; } for (int XX = Block * BlockSize; XX < Width; XX++, LinePS += 3) // 每行剩余的像素按照正常的方式处理 { int Blue = LinePS[0], Green = LinePS[1], Red = LinePS[2]; LinePY[XX] = (Y_B_WT * Blue + Y_G_WT * Green + Y_R_WT * Red + HalfV) >> Shift; LinePU[XX] = ((U_B_WT * Blue + U_G_WT * Green + U_R_WT * Red + HalfV) >> Shift) + 128; LinePV[XX] = ((V_B_WT * Blue + V_G_WT * Green + V_R_WT * Red + HalfV) >> Shift) + 128; } } }
（显示器不是宽屏的朋友请下载代码在VS里浏览，不然会看晕的）。
基本上就是按照普通的C的代码翻译到SSE版本，按F5编译执行，测试速度1080P的平均每帧约4.5ms，大概是普通C语言的2.2倍，为什么没到四倍，可能是我写的不到位，但是我认为主要的原因还是这里面有蛮多的数据拆解和数据类型的转换占用了一定的CPU周期。但总的来说还是有相当大的进步的。
稍微分析下程序以帮助不太懂SSE的朋友。
前面的几句加载和shffule语句主要是把BGR格式的图像数据拆解为单独的通道数据，其详细的原理见（真心描述的清楚)：
SSE图像算法优化系列八：自然饱和度(Vibrance)算法的模拟实现及其SSE优化（附源码，可作为SSE图像入门，Vibrance算法也可用于简单的肤色调整）。
中间一部分unpack代码主要的主要作用是把字节数据扩展到32位的int类型数，因为后面的乘法结果是只能用32位才能表达完整的。
那么后面的_mm_srai_epi32、_mm_add_epi32、_mm_mullo_epi32等就是直接的C代码的翻译了，只是把普通的晕算法用函数名代替了而已。
从理论上说，上述代码中的 const int Shift = 15;这个常量最大可以取到23，这个与最后的计算精度有一定的关系，但是也是影响不大。
2.2倍的提速，就这样放弃了吗，不，这不科学，也不是我Imageshop的风格，继续加油。
每当这个时候，我总是习惯性的再翻一遍Intel的SSE帮助文档，也许那个里面还藏着其他的葵花宝典。
当我们定位到_mm_madd_epi16这个函数时，虽然他只比普通的add函数多了一个m，但是当看到他下面对该函数功能的解释时，那真的眼前一亮，看下MSDN的说明：
Multiplies the 8 signed 16-bit integers from a by the 8 signed 16-bit integers from b.
__m128i _mm_madd_epi16 (__m128i a, __m128i b);
r0 := (a0 * b0) + (a1 * b1)
r1 := (a2 * b2) + (a3 * b3)
r2 := (a4 * b4) + (a5 * b5)
r3 := (a6 * b6) + (a7 * b7)
注意到这个函数内部实际执行了8次乘法和4次加法，唯一区别的就是参数要求是signed short类型，那能应用到我们的例子中吗？
首先我们来看Y变量的计算式：
LinePY[XX] = (Y_B_WT * Blue + Y_G_WT * Green + Y_R_WT * Red + HalfV) >> Shift;
LinePU[XX] = ((U_B_WT * Blue + U_G_WT * Green + U_R_WT * Red + HalfV) >> Shift) + 128;
LinePV[XX] = ((V_B_WT * Blue + V_G_WT * Green + V_R_WT * Red + HalfV) >> Shift) + 128;
注意到其中的HalfV = 1 << (Shift - 1)，稍微改写一下，我们得到：
LinePY[XX] = (Y_B_WT * Blue + Y_G_WT * Green + Y_R_WT * Red + 1 * HalfV) >> Shift;
LinePU[XX] = (U_B_WT * Blue + U_G_WT * Green + U_R_WT * Red + 257 * HalfV) >> Shift;
LinePV[XX] = (V_B_WT * Blue + V_G_WT * Green + V_R_WT * Red + 257 * HalfV) >> Shift;
至于算式中的257是怎么得到的，如果你不明白，回家问下你正在上小学的儿子吧。
这样改写目的很明显，我们正好把括号内的计算配对成了2组可利用_mm_madd_epi16函数进行计算的形式，那么他在数据类型方面能符合_mm_madd_epi16的需求吗，或者说我们能改造他使得我们的数据满足要求不。
_mm_madd_epi16的文档中有提到其作用是：Multiplies the 8 signed 16-bit integers from a by the 8 signed 16-bit integers from b，即参与计算的必须是带符号的16位数据，首先图像数范围[0,255]，肯定满足要求，同时注意到上述系数绝对值都小于1，因此只要(1 << Shift)不大于32768，即Shift最大取值15时，时能够满足需求的，而Shift=15时的精度也足以满足实际的运用需求。
注意我们改写的最后一项，比如257 * HalfV，这里我们可以认为HalfV是该项的权重，257是系数。
还有一点，也是最重要的一点，主要到没有r0 := (a0 * b0) + (a1 * b1)这里的系数是交叉的，比如我们认为变量b里面保存的是交叉的B和G分分量的权重，那么a变量里就应该交叉保存了Blue和Green的值，这时有两个途径：
1、我们上述代码里已经获得了Blue和Green分量的连续排列变量，这个时候只需要使用unpacklo和unpackhi就能分别获取低8位和高8位的交叉结果。
2、注意到获取Blue和Green分量的连续排列变量时是用的shuffle指令，我们也可以采用不同的shuffle系数直接获取交叉后的结果啊。
毫无疑问，第二种方式要快多了。
好，我们贴出这样改进后的结果（似乎不用宽屏也能显示全了）：
void RGBToYUVSSE(unsigned char *RGB, unsigned char *Y, unsigned char *U, unsigned char *V, int Width, int Height, int Stride) {const int Shift = 15; // 这里没有绝对值大于1的系数，最大可取2^15次方的放大倍数。 const int HalfV = 1 << (Shift - 1); const int Y_B_WT = 0.114f * (1 << Shift), Y_G_WT = 0.587f * (1 << Shift), Y_R_WT = (1 << Shift) - Y_B_WT - Y_G_WT, Y_C_WT = 1; const int U_B_WT = 0.436f * (1 << Shift), U_G_WT = -0.28886f * (1 << Shift), U_R_WT = -(U_B_WT + U_G_WT), U_C_WT = 257; const int V_B_WT = -0.10001 * (1 << Shift), V_G_WT = -0.51499f * (1 << Shift), V_R_WT = -(V_B_WT + V_G_WT), V_C_WT = 257; __m128i Weight_YBG = _mm_setr_epi16(Y_B_WT, Y_G_WT, Y_B_WT, Y_G_WT, Y_B_WT, Y_G_WT, Y_B_WT, Y_G_WT); __m128i Weight_YRC = _mm_setr_epi16(Y_R_WT, Y_C_WT, Y_R_WT, Y_C_WT, Y_R_WT, Y_C_WT, Y_R_WT, Y_C_WT); __m128i Weight_UBG = _mm_setr_epi16(U_B_WT, U_G_WT, U_B_WT, U_G_WT, U_B_WT, U_G_WT, U_B_WT, U_G_WT); __m128i Weight_URC = _mm_setr_epi16(U_R_WT, U_C_WT, U_R_WT, U_C_WT, U_R_WT, U_C_WT, U_R_WT, U_C_WT); __m128i Weight_VBG = _mm_setr_epi16(V_B_WT, V_G_WT, V_B_WT, V_G_WT, V_B_WT, V_G_WT, V_B_WT, V_G_WT); __m128i Weight_VRC = _mm_setr_epi16(V_R_WT, V_C_WT, V_R_WT, V_C_WT, V_R_WT, V_C_WT, V_R_WT, V_C_WT); __m128i Half = _mm_setr_epi16(0, HalfV, 0, HalfV, 0, HalfV, 0, HalfV); __m128i Zero = _mm_setzero_si128(); int BlockSize = 16, Block = Width / BlockSize; for (int YY = 0; YY < Height; YY++) { unsigned char *LinePS = RGB + YY * Stride; unsigned char *LinePY = Y + YY * Width; unsigned char *LinePU = U + YY * Width; unsigned char *LinePV = V + YY * Width; for (int XX = 0; XX < Block * BlockSize; XX += BlockSize, LinePS += BlockSize * 3) { __m128i Src1 = _mm_loadu_si128((__m128i *)(LinePS + 0)); __m128i Src2 = _mm_loadu_si128((__m128i *)(LinePS + 16)); __m128i Src3 = _mm_loadu_si128((__m128i *)(LinePS + 32)); __m128i BGL = _mm_shuffle_epi8(Src1, _mm_setr_epi8(0, 1, 3, 4, 6, 7, 9, 10, 12, 13, 15, -1, -1, -1, -1, -1)); BGL = _mm_or_si128(BGL, _mm_shuffle_epi8(Src2, _mm_setr_epi8(-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 0, 2, 3, 5, 6))); __m128i BGH = _mm_shuffle_epi8(Src2, _mm_setr_epi8(8, 9, 11, 12, 14, 15, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1)); BGH = _mm_or_si128(BGH, _mm_shuffle_epi8(Src3, _mm_setr_epi8(-1, -1, -1, -1, -1, -1, 1, 2, 4, 5, 7, 8, 10, 11, 13, 14))); __m128i RCL = _mm_shuffle_epi8(Src1, _mm_setr_epi8(2, -1, 5, -1, 8, -1, 11, -1, 14, -1, -1, -1, -1, -1, -1, -1)); RCL = _mm_or_si128(RCL, _mm_shuffle_epi8(Src2, _mm_setr_epi8(-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 1, -1, 4, -1, 7, -1))); __m128i RCH = _mm_shuffle_epi8(Src2, _mm_setr_epi8(10, -1, 13, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1)); RCH = _mm_or_si128(RCH, _mm_shuffle_epi8(Src3, _mm_setr_epi8(-1, -1, -1, -1, 0, -1, 3, -1, 6, -1, 9, -1, 12, -1, 15, -1))); __m128i BGLL = _mm_unpacklo_epi8(BGL, Zero); __m128i BGLH = _mm_unpackhi_epi8(BGL, Zero); __m128i RCLL = _mm_or_si128(_mm_unpacklo_epi8(RCL, Zero), Half); // HalfV是大于255的，只能在此处进行合并 __m128i RCLH = _mm_or_si128(_mm_unpackhi_epi8(RCL, Zero), Half); __m128i BGHL = _mm_unpacklo_epi8(BGH, Zero); __m128i BGHH = _mm_unpackhi_epi8(BGH, Zero); __m128i RCHL = _mm_or_si128(_mm_unpacklo_epi8(RCH, Zero), Half); __m128i RCHH = _mm_or_si128(_mm_unpackhi_epi8(RCH, Zero), Half); __m128i Y_LL = _mm_srai_epi32(_mm_add_epi32(_mm_madd_epi16(BGLL, Weight_YBG), _mm_madd_epi16(RCLL, Weight_YRC)), Shift); __m128i Y_LH = _mm_srai_epi32(_mm_add_epi32(_mm_madd_epi16(BGLH, Weight_YBG), _mm_madd_epi16(RCLH, Weight_YRC)), Shift); __m128i Y_HL = _mm_srai_epi32(_mm_add_epi32(_mm_madd_epi16(BGHL, Weight_YBG), _mm_madd_epi16(RCHL, Weight_YRC)), Shift); __m128i Y_HH = _mm_srai_epi32(_mm_add_epi32(_mm_madd_epi16(BGHH, Weight_YBG), _mm_madd_epi16(RCHH, Weight_YRC)), Shift); _mm_storeu_si128((__m128i*)(LinePY + XX), _mm_packus_epi16(_mm_packus_epi32(Y_LL, Y_LH), _mm_packus_epi32(Y_HL, Y_HH))); __m128i U_LL = _mm_srai_epi32(_mm_add_epi32(_mm_madd_epi16(BGLL, Weight_UBG), _mm_madd_epi16(RCLL, Weight_URC)), Shift); __m128i U_LH = _mm_srai_epi32(_mm_add_epi32(_mm_madd_epi16(BGLH, Weight_UBG), _mm_madd_epi16(RCLH, Weight_URC)), Shift); __m128i U_HL = _mm_srai_epi32(_mm_add_epi32(_mm_madd_epi16(BGHL, Weight_UBG), _mm_madd_epi16(RCHL, Weight_URC)), Shift); __m128i U_HH = _mm_srai_epi32(_mm_add_epi32(_mm_madd_epi16(BGHH, Weight_UBG), _mm_madd_epi16(RCHH, Weight_URC)), Shift); _mm_storeu_si128((__m128i*)(LinePU + XX), _mm_packus_epi16(_mm_packus_epi32(U_LL, U_LH), _mm_packus_epi32(U_HL, U_HH))); __m128i V_LL = _mm_srai_epi32(_mm_add_epi32(_mm_madd_epi16(BGLL, Weight_VBG), _mm_madd_epi16(RCLL, Weight_VRC)), Shift); __m128i V_LH = _mm_srai_epi32(_mm_add_epi32(_mm_madd_epi16(BGLH, Weight_VBG), _mm_madd_epi16(RCLH, Weight_VRC)), Shift); __m128i V_HL = _mm_srai_epi32(_mm_add_epi32(_mm_madd_epi16(BGHL, Weight_VBG), _mm_madd_epi16(RCHL, Weight_VRC)), Shift); __m128i V_HH = _mm_srai_epi32(_mm_add_epi32(_mm_madd_epi16(BGHH, Weight_VBG), _mm_madd_epi16(RCHH, Weight_VRC)), Shift); _mm_storeu_si128((__m128i*)(LinePV + XX), _mm_packus_epi16(_mm_packus_epi32(V_LL, V_LH), _mm_packus_epi32(V_HL, V_HH))); } for (int XX = Block * BlockSize; XX < Width; XX++, LinePS += 3) // 每行剩余的像素按照正常的方式处理 { int Blue = LinePS[0], Green = LinePS[1], Red = LinePS[2]; LinePY[XX] = (Y_B_WT * Blue + Y_G_WT * Green + Y_R_WT * Red + Y_C_WT * HalfV) >> Shift; LinePU[XX] = (U_B_WT * Blue + U_G_WT * Green + U_R_WT * Red + U_C_WT * HalfV) >> Shift; LinePV[XX] = (V_B_WT * Blue + V_G_WT * Green + V_R_WT * Red + V_C_WT * HalfV) >> Shift; } } }
速度如何呢，结果是2.2ms，相比纯C语言约有4.5倍的提升。
接下来我们简单的谈下YUV到RGB空间的优化，普通的C语言如下：
void YUVToRGB(unsigned char *Y, unsigned char *U, unsigned char *V, unsigned char *RGB, int Width, int Height, int Stride) { const int Shift = 15; const int HalfV = 1 << (Shift - 1); const int B_Y_WT = 1 << Shift, B_U_WT = 2.03211f * (1 << Shift), B_V_WT = 0; const int G_Y_WT = 1 << Shift, G_U_WT = -0.39465f * (1 << Shift), G_V_WT = -0.58060f * (1 << Shift); const int R_Y_WT = 1 << Shift, R_U_WT = 0, R_V_WT = 1.13983 * (1 << Shift); for (int YY = 0; YY < Height; YY++) { unsigned char *LinePD = RGB + YY * Stride; unsigned char *LinePY = Y + YY * Width; unsigned char *LinePU = U + YY * Width; unsigned char *LinePV = V + YY * Width; for (int XX = 0; XX < Width; XX++, LinePD += 3) { int YV = LinePY[XX], UV = LinePU[XX] - 128, VV = LinePV[XX] - 128; LinePD[0] = ClampToByte(YV + ((B_U_WT * UV + HalfV) >> Shift)); LinePD[1] = ClampToByte(YV + ((G_U_WT * UV + G_V_WT * VV + HalfV) >> Shift)); LinePD[2] = ClampToByte(YV + ((R_V_WT * VV + HalfV) >> Shift)); } } }
虽然内部的加减乘除运算减少了，但是由于需要进行范围的Clamp，因此纯C版本的实际耗时要比RGB转到YUV多，大约在11.5ms左右。
同样的道理，我们还准备使用_mm_madd_epi16来优化，但是这里的情况就稍微复杂了一点。
第一，上面的C代码为了速度考虑，对YV分量没有乘以系数，因为他乘以系数和后面的移位抵消了，但是我们还是考虑把他们展开，以LinePD[0]为例：
LinePD[0] = ClampToByte(YV + ((B_U_WT * UV + HalfV) >> Shift));
展开：
LinePD[0] = ClampToByte(YV * (1 << Shift) + B_U_WT * UV + (1 << (Shift - 1))) >> Shift))
= ClampToByte((YV + 0.5) * (1 << Shift) + B_U_WT * UV) >> Shift))
= ClampToByte((YV * 2 + 1) * ((1 << Shift) >> 1) + B_U_WT * UV) >> Shift))
为什么这样做呢，其目的就是想只用一次_mm_madd_epi16就可以达到结果，此时相当于我们把B_Y_WT = 1 << Shift改为了B_Y_WT = (1 << Shift) >> 1，同时像素值进行了乘2和加1，同样的道理LinePD[1]也可以进行同样的处理。
LinePD[1]展开后你们觉得还要不要说，是RGB到YUV里是相同的道理。
至于ClampToByte，SSE自带了抗饱和处理，packs或者packus直接实现了这个功能，这也是个额外提速的地方。
但是注意，由于YUV到RGB的转换系数里有个数据2.03211f ，大于2了，因此考虑_mm_madd_epi16的应用范围这里的Shift最大也只能取为13了，但是似乎也足够了。编码实测SSE优化后的代码速度也在2.2ms左右，提速比达到了5.2倍，没有比RGB到SSE快，是因为中间增加了一些加法计算。
这样计算下来，两个颜色空间相互转换的时间针对1080P的图，大约需要4.4ms，通常情况下这种转换是为了让算法只在Y通道做处理了，处理完后新的Y通道值和原有图的UV通道组合起来，在转换回RGB图像，这样就把本来需要在RGB三通道都处理的算法简化到一个通道处理了，在对Y通道做的变化不太大的时候，这种处理方式的处理结果和直接在RGB空间做处理时效果差不多，但是如果不考虑空间转换的耗时，那么速度确能提高三倍，这里的空间转化前后耗时4.4ms，对整体帧率的影响就非常小了，比图我是用的基于均方差的磨皮就可以采用这种方式优化从而达到实时的效果。
在如上的实际应用中，考虑到内存占用，我们也没有必要把UV通道的数据计算出来，而只计算Y通道的数据，这样把上述计算Y通道的代码单独提取出来，然后再写个void ModifyYComponent(unsigned char *Src, unsigned char *NewY, unsigned char *Dest, int Width, int Height, int Stride)这样的函数把Y通道替换掉，这样一个连续的写法其中的-128和+128这项又可以免去处理，当然具体怎么写内部还是有很多学问值得研究的，只是感觉上没有太多人去使用SSE了，这里就不浪费自己的宝贵时间了。
至于XYZ和RGB空间的互换，由于他们的3个系数都没有为0的，只要注意相关shift值得取值合理，优化方式都例同RGB到XYZ的过程，这里不予以赘述。
其实还有很多细节的东西再最初写的时候也经过了各种尝试，只有通过自己编写才能体会中这中间的快乐和痛苦，看别人的代码你看到的只是最后的成品。
回顾下今年一年的文章，基本上就是在学习SIMD优化，确实SIMD的优化确实在某些场合有着很大的提速价值，在手机端也有NEON指令集，基本上和SSE2都有对应的关系，不过我一心只关注PC的优化，这个市场现在是越来越小了，虽然在航空、医疗等工业场景还是有一定的应用场景。手机编程太痛苦了，都一个老人家了，懒得再去折腾了。
快过年了，考虑下2018年，作为一个编程自由人，我初步考虑下2018可能研究下16位图像的增强、HDR技术等等，至于深度学习，还是算了吧，知道有那么个东西就OK了。
提前祝大家2018年新年快乐。
本文相关代码下载：https://files.cnblogs.com/files/Imageshop/YUV_RGB.rar
基于本文的优化应用于磨皮的效果见：彩色图工程：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar
到年底了，没钱回家过年，望各位路过的大神施舍点。
