SSE圖像算法優化系列一：一段BGR2Y的SIMD代碼解析。

     一個同事在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日即將離開南京之際，特此記念。