在顏色空間系列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
到年底了，沒錢回家過年，望各位路過的大神施舍點。
