使用SSE優化兩種顔色空間的轉換，測試1080P的圖像能達到455fps左右，可顯著提高視訊算法的實時性。
在顔色空間系列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); Return Value
Adds the signed 32-bit integer results pairwise and packs the 4 signed 32-bit integer results.
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
到年底了，沒錢回家過年，望各位路過的大神施舍點。
