SSE圖像算法優化系列十七：多個圖像處理中經常使用函數的SSE實現。

時間  2019-11-06
標籤 sse 圖像 算法 優化 系列 十七 多個 圖像處理 中經 常使 函數 實現 简体版
原文   原文鏈接

　　在作圖像處理的SSE優化時，也會常常遇到一些小的過程、數值優化等代碼，本文分享一些我的收藏或實現的代碼片斷給你們。c++

1、快速求對數運算git

　　對數運算在圖像處理中也是個常常會遇到的過程，特備是在一些數據壓縮和空間轉換時經常會用到，並且是個比較耗時的函數，標準的SSE庫裏並無提供該函數的實現，若是須要高精度的SSE版本，網絡上已經有了，參考：https://github.com/to-miz/sse_mathfun_extension/blob/master/sse_mathfun.h，這個的精度和標準庫的精度基本一致了，稍做整理後的代碼以下：github

//    對數函數的SSE實現，高精度版
inline __m128 _mm_log_ps(__m128 x)
{
    static const __declspec(align(16)) int _ps_min_norm_pos[4] = { 0x00800000, 0x00800000, 0x00800000, 0x00800000 };
    static const __declspec(align(16)) int _ps_inv_mant_mask[4] = { ~0x7f800000, ~0x7f800000, ~0x7f800000, ~0x7f800000 };
    static const __declspec(align(16)) int _pi32_0x7f[4] = { 0x7f, 0x7f, 0x7f, 0x7f };
    static const __declspec(align(16)) float _ps_1[4] = { 1.0f, 1.0f, 1.0f, 1.0f };
    static const __declspec(align(16)) float _ps_0p5[4] = { 0.5f, 0.5f, 0.5f, 0.5f };
    static const __declspec(align(16)) float _ps_sqrthf[4] = { 0.707106781186547524f, 0.707106781186547524f, 0.707106781186547524f, 0.707106781186547524f };
    static const __declspec(align(16)) float _ps_log_p0[4] = { 7.0376836292E-2f, 7.0376836292E-2f, 7.0376836292E-2f, 7.0376836292E-2f };
    static const __declspec(align(16)) float _ps_log_p1[4] = { -1.1514610310E-1f, -1.1514610310E-1f, -1.1514610310E-1f, -1.1514610310E-1f };
    static const __declspec(align(16)) float _ps_log_p2[4] = { 1.1676998740E-1f, 1.1676998740E-1f, 1.1676998740E-1f, 1.1676998740E-1f };
    static const __declspec(align(16)) float _ps_log_p3[4] = { -1.2420140846E-1f, -1.2420140846E-1f, -1.2420140846E-1f, -1.2420140846E-1f };
    static const __declspec(align(16)) float _ps_log_p4[4] = { 1.4249322787E-1f, 1.4249322787E-1f, 1.4249322787E-1f, 1.4249322787E-1f };
    static const __declspec(align(16)) float _ps_log_p5[4] = { -1.6668057665E-1f, -1.6668057665E-1f, -1.6668057665E-1f, -1.6668057665E-1f };
    static const __declspec(align(16)) float _ps_log_p6[4] = { 2.0000714765E-1f, 2.0000714765E-1f, 2.0000714765E-1f, 2.0000714765E-1f };
    static const __declspec(align(16)) float _ps_log_p7[4] = { -2.4999993993E-1f, -2.4999993993E-1f, -2.4999993993E-1f, -2.4999993993E-1f };
    static const __declspec(align(16)) float _ps_log_p8[4] = { 3.3333331174E-1f, 3.3333331174E-1f, 3.3333331174E-1f, 3.3333331174E-1f };
    static const __declspec(align(16)) float _ps_log_q1[4] = { -2.12194440e-4f, -2.12194440e-4f, -2.12194440e-4f, -2.12194440e-4f };
    static const __declspec(align(16)) float _ps_log_q2[4] = { 0.693359375f, 0.693359375f, 0.693359375f, 0.693359375f };

    __m128 one = *(__m128*)_ps_1;
    __m128 invalid_mask = _mm_cmple_ps(x, _mm_setzero_ps());
    /* cut off denormalized stuff */
    x = _mm_max_ps(x, *(__m128*)_ps_min_norm_pos);
    __m128i emm0 = _mm_srli_epi32(_mm_castps_si128(x), 23);

    /* keep only the fractional part */
    x = _mm_and_ps(x, *(__m128*)_ps_inv_mant_mask);
    x = _mm_or_ps(x, _mm_set1_ps(0.5f));

    emm0 = _mm_sub_epi32(emm0, *(__m128i *)_pi32_0x7f);
    __m128 e = _mm_cvtepi32_ps(emm0);
    e = _mm_add_ps(e, one);

    __m128 mask = _mm_cmplt_ps(x, *(__m128*)_ps_sqrthf);
    __m128 tmp = _mm_and_ps(x, mask);
    x = _mm_sub_ps(x, one);
    e = _mm_sub_ps(e, _mm_and_ps(one, mask));
    x = _mm_add_ps(x, tmp);

    __m128 z = _mm_mul_ps(x, x);
    __m128 y = *(__m128*)_ps_log_p0;
    y = _mm_mul_ps(y, x);
    y = _mm_add_ps(y, *(__m128*)_ps_log_p1);
    y = _mm_mul_ps(y, x);
    y = _mm_add_ps(y, *(__m128*)_ps_log_p2);
    y = _mm_mul_ps(y, x);
    y = _mm_add_ps(y, *(__m128*)_ps_log_p3);
    y = _mm_mul_ps(y, x);
    y = _mm_add_ps(y, *(__m128*)_ps_log_p4);
    y = _mm_mul_ps(y, x);
    y = _mm_add_ps(y, *(__m128*)_ps_log_p5);
    y = _mm_mul_ps(y, x);
    y = _mm_add_ps(y, *(__m128*)_ps_log_p6);
    y = _mm_mul_ps(y, x);
    y = _mm_add_ps(y, *(__m128*)_ps_log_p7);
    y = _mm_mul_ps(y, x);
    y = _mm_add_ps(y, *(__m128*)_ps_log_p8);
    y = _mm_mul_ps(y, x);

    y = _mm_mul_ps(y, z);
    tmp = _mm_mul_ps(e, *(__m128*)_ps_log_q1);
    y = _mm_add_ps(y, tmp);
    tmp = _mm_mul_ps(z, *(__m128*)_ps_0p5);
    y = _mm_sub_ps(y, tmp);
    tmp = _mm_mul_ps(e, *(__m128*)_ps_log_q2);
    x = _mm_add_ps(x, y);
    x = _mm_add_ps(x, tmp);
    x = _mm_or_ps(x, invalid_mask); // negative arg will be NAN

    return x;
}

　　看上去有一大堆代碼，不過實測這個的速度越是標準庫（本文是指啓動加強指令集選項設置爲：未設置，設計上編譯器在此種狀況下會自動設置爲SSE2加強，這能夠從反編譯logf函數看到，所以，這裏的速度比較還不是和純Fpu實現的比較）的2倍，若是稍微下降點精度，好比_ps_log_p5到_ps_log_p8之間的代碼，還能提升點速度。算法

　　另外，在不少場合咱們還可使用另一種低精度的log函數，其C代碼以下所示：編程

//https://stackoverflow.com/questions/9411823/fast-log2float-x-implementation-c
inline float IM_Flog(float val)
{
    union
    {
        float val;
        int x;
    } u = { val };
    float log_2 = (float)(((u.x >> 23) & 255) - 128);
    u.x &= ~(255 << 23);
    u.x += (127 << 23);
    log_2 += ((-0.34484843f) * u.val + 2.02466578f) * u.val - 0.67487759f;
    return log_2 * 0.69314718f;
}

　　這個函數大概有小數點後2位精度。網絡

　　上述代碼大約也是標準函數的2倍速度左右。可是上述函數是能夠向量化的，咱們來嘗試實現。函數

　　咱們首先來看聯合體，其實這個東西就是兩個東西佔同一個內存空間，而後外部用不一樣的規則去讀取他，在SSE裏，有着豐富的cast函數，他也是幹這個事情的，好比這裏的聯合體就能夠用_mm_castps_si128來轉換，而實際上這個Intrinsic並不會產生任何的彙編語句。性能

　　那麼後面的那些移位、或運算、非運算、加減乘除之類的就是直接翻譯了，毫無難處，完整的代碼以下所示：學習

inline __m128 _mm_flog_ps(__m128 x)
{
    __m128i I = _mm_castps_si128(x);
    __m128 log_2 = _mm_cvtepi32_ps(_mm_sub_epi32(_mm_and_si128(_mm_srli_epi32(I, 23), _mm_set1_epi32(255)), _mm_set1_epi32(128)));
    I = _mm_and_si128(I, _mm_set1_epi32(-2139095041));        //    255 << 23
    I = _mm_add_epi32(I, _mm_set1_epi32(1065353216));        //    127 << 23
    __m128 F = _mm_castsi128_ps(I);
    __m128 T = _mm_add_ps(_mm_mul_ps(_mm_set1_ps(-0.34484843f), F), _mm_set1_ps(2.02466578f));
    T = _mm_sub_ps(_mm_mul_ps(T, F), _mm_set1_ps(0.67487759f));
    return _mm_mul_ps(_mm_add_ps(log_2, T), _mm_set1_ps(0.69314718f));
}

　　通過實測，這個速度能夠達到標準庫的7到8倍的優點。測試

2、快速求冪運算

　　通常圖像編程中有log出現的地方就會有exp出現，所以exp的優化也尤其重要，一樣在sse_mathfun.h中也有exp的優化（還有sin,cos的SSE優化語句呢），我這裏就不貼那個的代碼了，咱們一樣關注下用聯合體實現的近似快速算法，其C代碼以下所示：

inline float IM_Fexp(float Y)            
{
    union
    {
        double Value;
        int X[2];
    } V;
    V.X[1] = (int)(Y * 1512775 + 1072632447 + 0.5F);
    V.X[0] = 0;
    return (float)V.Value;
}

　　測試這個和標準的exp庫函數速度竟然差很少，不曉得爲啥，但咱們來試下他的SSE優化版本了。

　V.X[1] = (int)(Y * 1512775 + 1072632447 + 0.5F);這句話沒啥難度，直接翻譯就能夠了，注意幾個強制類型轉化就能夠了，以下所示：
__m128i T = _mm_cvtps_epi32(_mm_add_ps(_mm_mul_ps(Y, _mm_set1_ps(1512775)), _mm_set1_ps(1072632447)));

　　因爲咱們想一次性處理4個float類型的數據，所以也就須要4個union的空間，這樣就須要2個__m128i變量來保存數據，每一個XMM寄存器的數據應該分別爲：

　　T1    0    T0    0         +     T3   0    T2    0 　　　　（高位----》低位）

　　這個可使用unpack來實現，具體以下：

    __m128i TL = _mm_unpacklo_epi32(_mm_setzero_si128(), T);
    __m128i TH = _mm_unpackhi_epi32(_mm_setzero_si128(), T);

　　最後咱們認爲__m128i裏的數據是double數據，直接一個cast就能夠了，而後由於咱們只須要單精度的數據，再使用_mm_cvtpd_ps將double轉換爲float類型，注意這個時候還須要將他們鏈接再一塊兒造成一個完整的__m128變量，最終的代碼以下：

inline __m128 _mm_fexp_ps(__m128 Y)
{
    __m128i T = _mm_cvtps_epi32(_mm_add_ps(_mm_mul_ps(Y, _mm_set1_ps(1512775)), _mm_set1_ps(1072632447)));
    __m128i TL = _mm_unpacklo_epi32(_mm_setzero_si128(), T);
    __m128i TH = _mm_unpackhi_epi32(_mm_setzero_si128(), T);
    return _mm_movelh_ps(_mm_cvtpd_ps(_mm_castsi128_pd(TL)), _mm_cvtpd_ps(_mm_castsi128_pd(TH)));
}

　　實測這個的提速大概有10倍。

　　若是要求double的exp，其SSE代碼你會了嗎？

3、pow函數的優化。

　　一種經常使用的近似算法以下所示：

inline float IM_Fpow(float a, float b) 
{
    union
    {
        double Value;
        int X[2];
    } V;
    V.X[1] = (int)(b * (V.X[1] - 1072632447) + 1072632447);
    V.X[0] = 0;
    return (float)V.Value;
}

　　和exp很相似，留給有興趣的人本身實現。

四：兩個求倒數函數的優化誤區

　　SSE提供了連個快速求倒數的函數，_mm_rcp_ps，_mm_rsqrt_ps，他們都是近似值，只有12bit的精度，若是想經過他們獲得精確的倒數值，須要牛頓 - 拉弗森方法，好比利用_mm_rcp_ps求精確倒數的代碼以下：

__forceinline __m128 _mm_prcp_ps(__m128 a)
{
    __m128 rcp = _mm_rcp_ps(a);            //    此函數只有12bit的精度.                    
    return _mm_sub_ps(_mm_add_ps(rcp, rcp), _mm_mul_ps(a, _mm_mul_ps(rcp, rcp)));    //    x1 = x0 * (2 - d * x0) = 2 * x0 - d * x0 * x0，使用牛頓 - 拉弗森方法這種方法能夠提升精度到23bit
}

　　可是實測這個還不如直接用_mm_div_ps的速度，即便是下面的函數：

__forceinline __m128 _mm_fdiv_ps(__m128 a, __m128 b)
{
    return _mm_mul_ps(a, _mm_rcp_ps(b));
}

　　彷佛速度也不夠好，並且精度還低了。

　　特別低，若是使用_mm_rcp_ps和_mm_rsqrt_ps聯合求近似sqrt，即以下代碼，速度好像還慢了，真搞不明白爲何。

__forceinline __m128 _mm_fsqrt_ps(__m128 a)
{
    return _mm_rcp_ps(_mm_rsqrt_ps(a));
}

五：避免除數爲0時的得到無效效果

　　在SSE指令中，沒有提供整數的除法指令，不知道這是爲何，因此整數除法通常只能借用浮點版本的指令，同時，除法存在的一個問題就是若是除數爲0，可能會觸發異常，不過SSE在這種狀況下不會拋出異常，可是咱們應該避免，避免的方式有不少，好比判斷若是除數爲0，就作特殊處理，或者若是除數爲0就除以一個很小的數，不過大部分的需求是，除數爲0，則返回0，此時就可使用下面的SSE指令代替_mm_div_ps:

//    四個浮點數的除法a/b，若是b中某個份量爲0，則對應位置返回0值
inline __m128 _mm_divz_ps(__m128 a, __m128 b)
{
    __m128 Mask = _mm_cmpeq_ps(b, _mm_setzero_ps());
    return _mm_blendv_ps(_mm_div_ps(a, b), _mm_setzero_ps(), Mask);
}

　　即先把除數和0進行比較，而後在把_mm_div_ps的返回值中，除數爲0的部分用0代替，固然，這會帶來必定的性能降低。

　　實際上，利用位運算，上述代碼還能夠稍做優化以下：

inline __m128 _mm_divz_ps(__m128 a, __m128 b)
{
    return _mm_and_ps(_mm_div_ps(a, b), _mm_cmpneq_ps(b, _mm_setzero_ps()));
}

6、將4個32位整數轉換爲字節數並保存

　　不少狀況下，咱們的算法計算須要將字節類型擴展到32位，計算完成後再保存的字節數據中，這個時候使用SSE的話是沒有直接的指令的，不過SSE4提供了一條_mm_cvtsi128_si32指令，能夠將XMM寄存器的4個int32數據轉換爲4個字節數據並保存到一個普通的寄存器中，所以，咱們只要調用幾個合適的pack語句就能夠實現這個功能了，以下所示：

//    將4個32位整形變量數據打包到4個字節數據中
inline void _mm_storesi128_4char(unsigned char *Dest, __m128i P)
{
    __m128i T = _mm_packs_epi32(P, P);
    *((int *)Dest) = _mm_cvtsi128_si32(_mm_packus_epi16(T, T));
}

7、讀取12個字節數到一個XMM寄存器中

　　XMM寄存器是16個字節大小的，並且SSE的不少計算是以4的整數倍字節位單位進行的，可是在圖像處理中，70%狀況下處理的是彩色的24位圖像，即一個像素佔用3個字節，若是直接使用load指令載入數據，一次性可載入5加1/3個像素，這對算法的處理是很不方便的，通常情況下都是加載4個像素，即12個字節，而後擴展成16個字節（給每一個像素增長一個Alpha值），咱們固然能夠直接使用load加載16個字節，而後每次跳過12個字節在進行load加載，可是其實也可使用下面的加載12個字節的函數：

//    從指針p處加載12個字節數據到XMM寄存器中，寄存器最高32位清0
inline __m128i _mm_loadu_epi96(const __m128i * p)
{
    return _mm_unpacklo_epi64(_mm_loadl_epi64(p), _mm_cvtsi32_si128(((int *)p)[2]));
}

　　保存固然也能夠只保存XMM的低12位：

//     將寄存器Q的低位12個字節數據寫入到指針P中。
inline void _mm_storeu_epi96(__m128i *P, __m128i Q)
{
    _mm_storel_epi64(P, Q);
    ((int *)P)[2] = _mm_cvtsi128_si32(_mm_srli_si128(Q, 8));
}

　　不過實際測試，可能仍是直接使用_mm_loadu_si128和_mm_storeu_si128快點，可是要注意循環的結束爲止。

8、整除255

　　整除255在圖像處理是很是很是常見的操做，前面說了，SSE裏沒有整數除法的指令，若是轉換到浮點在除那就慢的死了，通常狀況下若是要求精度不高可使用右移8位實現，若是非要精確值可使用以下的C代碼：

//    計算整數整除255的四捨五入結果。
inline int IM_Div255(int V)
{
    return (((V >> 8) + V + 1) >> 8);        //    彷佛V能夠是負數
}

　　翻譯爲SSE爲：

//    返回16位無符號整形數據整除255後四捨五入的結果： x = ((x + 1) + (x >> 8)) >> 8:
inline __m128i _mm_div255_epu16(__m128i x)
{
    return _mm_srli_epi16(_mm_adds_epu16(_mm_adds_epu16(x, _mm_set1_epi16(1)), _mm_srli_epi16(x, 8)), 8);
}

9、求XMM寄存器內全部元素的累加值

　　這也是個常見的需求，咱們可能把某個結果重複的結果保存在寄存器中，最後結束時在把寄存器中的每一個元素想加，你固然能夠經過訪問__m128i變量的內部的元素實現，可是聽說這樣會下降循環內的優化，一種方式是直接用SSE指令實現，好比對8個有符號的short類型的相加代碼以下所示：

//    8個有符號的16位的數據相加的和。
//    https://stackoverflow.com/questions/31382209/computing-the-inner-product-of-vectors-with-allowed-scalar-values-0-1-and-2-usi/31382878#31382878
inline int _mm_hsum_epi16(__m128i V)                            //    V7 V6 V5 V4 V3 V2 V1 V0
{
    //    V = _mm_unpacklo_epi16(_mm_hadd_epi16(V, _mm_setzero_si128()), _mm_setzero_si128());    也能夠用這句，_mm_hadd_epi16彷佛對計算結果超出32768能得到正確結果
    __m128i T = _mm_madd_epi16(V, _mm_set1_epi16(1));   //    V7+V6                        V5+V4            V3+V2    V1+V0
    T = _mm_add_epi32(T, _mm_srli_si128(T, 8));            //    V7+V6+V3+V2                    V5+V4+V1+V0        0        0        
    T = _mm_add_epi32(T, _mm_srli_si128(T, 4));            //    V7+V6+V3+V2+V5+V4+V1+V0        V5+V4+V1+V0        0        0    
    return _mm_cvtsi128_si32(T);                        //    提取低位    
}

　　對於epi32或者ps類型也是使用相似的過程的。

十、求16個字節的最小值

　　好比咱們要求一個字節序列的最小值，咱們確定會使用_mm_min_epi8這樣的函數保存每隔16個字節的最小值，這樣最終咱們獲得16個字節的一個XMM寄存器，整個序列的最小值確定在這個16個字節裏面，這個時候咱們能夠巧妙的借用下面的SSE語句獲得這16個字節的最小值：

//    求16個字節數據的最小值, 只能針對字節數據。
inline int _mm_hmin_epu8(__m128i a)
{
    __m128i L = _mm_unpacklo_epi8(a, _mm_setzero_si128());
    __m128i H = _mm_unpackhi_epi8(a, _mm_setzero_si128());
    return _mm_extract_epi16(_mm_min_epu16(_mm_minpos_epu16(L), _mm_minpos_epu16(H)), 0);
}

　　SSE3提供了_mm_minpos_epu16函數，他能獲取8個無符號數的的最小值及其最小值的索引，放置在寄存器的低16和低32位，咱們把字節數據擴展到16位，而後在經過兩次比較就能夠得到相應的最小值了。

　　那若是是求最大值呢，惋惜SSE沒有提供_mm_maxpos_epu16函數，可是也無妨，稍微修改下上面的代碼就能夠了，以下所示：

//    求16個字節數據的最大值, 只能針對字節數據。
inline int _mm_hmax_epu8(__m128i a)
{
    __m128i b = _mm_subs_epu8(_mm_set1_epi8(255), a);
    __m128i L = _mm_unpacklo_epi8(b, _mm_setzero_si128());
    __m128i H = _mm_unpackhi_epi8(b, _mm_setzero_si128());
    return 255 - _mm_extract_epi16(_mm_min_epu16(_mm_minpos_epu16(L), _mm_minpos_epu16(H)), 0);
}

11、其餘一些優化技巧

　　在http://www.alfredklomp.com/programming/sse-intrinsics/ 以及 http://www.itkeyword.com/doc/0326039046115117x827/c++-sse2-intrinsics-comparing-unsigned-integers等網站上還有不少參考的資料，但願你們本身去學習下。

相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程