對圖像算法進行SSE優化時，有很多常用的過程，本文列舉了十幾個例子，有些是很高效的，有些是很常用的，對研究圖像的朋友有一定的幫助。
在做圖像處理的SSE優化時，也會經常遇到一些小的過程、數值優化等代碼，本文分享一些個人收藏或實作的代碼片段給大家。
一、快速求對數運算
對數運算在圖像進行中也是個經常會遇到的過程，特備是在一些資料壓縮和空間轉換時常常會用到，而且是個比較耗時的函數，标準的SSE庫裡并沒有提供該函數的實作，如果需要高精度的SSE版本，網絡上已經有了，參考：https://github.com/to-miz/sse_mathfun_extension/blob/master/sse_mathfun.h，這個的精度和标準庫的精度基本一緻了，稍作整理後的代碼如下：
// 對數函數的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倍的優勢。
二、快速求幂運算
一般圖像程式設計中有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代碼你會了嗎？
三、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()));
} 六、将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));
} 七、讀取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快點，但是要注意循環的結束為止。
八、整除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);
} 九、求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類型也是使用類似的過程的。
10、求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);
} 十一、其他一些優化技巧
在http://www.alfredklomp.com/programming/sse-intrinsics/ 以及 http://www.itkeyword.com/doc/0326039046115117x827/c++-sse2-intrinsics-comparing-unsigned-integers等網站上還有很多參考的資料，希望大家自己去學習下。
![Log抓取和分析-BugReport[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)