一個同事在github上淘到一個基於SIMD的RGB轉Y（彩色轉灰度或者轉明度）的代碼，我抽了點時間看了下，順便學習了一些SIMD指令，這里把學習過程中的一些理解和認識共享給大家。
github上相關代碼見鏈接：https://github.com/komrad36/RGB2Y，這哥們還有其他一些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; } } }
很簡單，就是R/G/B分量分別乘以各自的系數得到亮度值，注意這個系數是歸一化的，為了速度考慮，我們采用了定點化措施，但是注意R_WT最后的量化方法，為了保證系數定點化后的歸一性，最佳方式就是用他們的理論之和減去其他的系數。
上述代碼的速度已經非常快了，在測試機上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; } } }
即采用4路並行，這樣同樣的圖大概需要3.40ms，稍有提高。
基本上這樣的速度能夠滿足所有場合的需求，但是也有一些極端的條件，比如一些用於高清視頻的美容方面等，每個過程能提高1ms都是很有必要的，因此，我一直在關注這方面的優化算法，正好komrad36提供了這方面的SIMD代碼。
我們來分析分析他提供的代碼先，為了篇幅，這里僅僅貼出核心的我稍作修改的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 17 __m128i sumaL = _mm_add_epi16(p3aL, _mm_add_epi16(p1aL, p2aL)); 18 __m128i sumaH = _mm_add_epi16(p3aH, _mm_add_epi16(p1aH, p2aH)); 19 __m128i sumbL = _mm_add_epi16(p3bL, _mm_add_epi16(p1bL, p2bL)); 20 __m128i sumbH = _mm_add_epi16(p3bH, _mm_add_epi16(p1bH, p2bH)); 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代碼的核心，好，現在一個函數一個函數的解釋了吧。
第一到第15行都是一樣的過程，其核心就是把字節數據讀入並和相應的權重相乘。_mm_loadu_si128就是把之后16個字節的數據讀入到一個SSE寄存器中,注意由於任意位置的圖像數據內存地址肯定不可能都滿足SIMD16字節對齊的規定，因此這里不是用的_mm_load_si128指令。而_mm_cvtepu8_epi16指令則把這16個字節的低64位的8個字節數擴展為8個16位數據，這樣做主要是為了上述的乘法進行准備的。那么這個其實也有其他的方式實現，比如使用_mm_unpacklo_epi8和 _mm_unpackhi_epi8配合_mm_setzero_si128也可以實現這樣的效果，而且可以節省后面的某一句_mm_loadu_si128指令，不過實測速度無區別。
_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的倍數）。
代碼中還有一些其他的技巧，有些數字可能讀者自己去看看那個指令的意義后會更加清晰，這些不太是語言能夠完美表達清楚的。當然，這段代碼在書寫藝術上還是有很大的改良空間的，有興趣的讀者可以自行研究下。
同時，komrad36還提供了基於AVX的指令算法，執行耗時大概是1.15ms，和普通的算法比提速比約為3:1，思路和SSE基本是一樣的。
我把這些代碼稍做整理提供給大家使用和測試，我也相信，上述指令肯定不是最佳的SIMD實現方式，比如改變指令順序讓指令級又可以並行等手段也許也是有效的提速方式之一。
最后一點就是我有個疑問，在我提供的代碼執行后，如果先使用SSE測試，后使用AVX測試，SSE的速度和上述報告數據差不多，但是一旦點了AVX測試后，在點SSE測試，SSE的速度就驟然下降很多，甚至比普通C的都要慢，我水平很有限，實在是不知道這是什么道理，煩請有知道的告知。
源代碼下載地址：http://files.cnblogs.com/files/Imageshop/FastRGB2Y.rar
本筆記創建於2016年1月8日即將離開南京之際，特此紀念。
