最近一直沉迷於SSE方面的優化，實在找不到想學習的參考資料了，就拿個筆記本放在腿上翻翻OpenCv的源代碼，無意中看到了OpenCv中關於積分圖的代碼，仔細研習了一番，覺得OpenCv對SSE的靈活運用真的做的很好，這里記錄下我對該段代碼的品味並將其思路擴展到其他通道數的圖像。
該核心代碼位於：Opencv 3.0\opencv\sources\modules\imgproc\src\sumpixels.cpp文件中。
我們貼出最感興趣的一部分代碼以便分析:
bool operator()(const uchar * src, size_t _srcstep,int * sum, size_t _sumstep,double * sqsum, size_t, int * tilted, size_t,Size size, int cn) const { if (sqsum || tilted || cn != 1 || !haveSSE2) return false; // the first iteration memset(sum, 0, (size.width + 1) * sizeof(int)); __m128i v_zero = _mm_setzero_si128(), prev = v_zero; int j = 0; // the others for (int i = 0; i < size.height; ++i) { const uchar * src_row = src + _srcstep * i; int * prev_sum_row = (int *)((uchar *)sum + _sumstep * i) + 1; int * sum_row = (int *)((uchar *)sum + _sumstep * (i + 1)) + 1; sum_row[-1] = 0; prev = v_zero; j = 0; for ( ; j + 7 < size.width; j += 8) { __m128i vsuml = _mm_loadu_si128((const __m128i *)(prev_sum_row + j)); __m128i vsumh = _mm_loadu_si128((const __m128i *)(prev_sum_row + j + 4)); __m128i el8shr0 = _mm_loadl_epi64((const __m128i *)(src_row + j)); __m128i el8shr1 = _mm_slli_si128(el8shr0, 1); __m128i el8shr2 = _mm_slli_si128(el8shr0, 2); __m128i el8shr3 = _mm_slli_si128(el8shr0, 3); vsuml = _mm_add_epi32(vsuml, prev); vsumh = _mm_add_epi32(vsumh, prev); __m128i el8shr12 = _mm_add_epi16(_mm_unpacklo_epi8(el8shr1, v_zero), _mm_unpacklo_epi8(el8shr2, v_zero)); __m128i el8shr03 = _mm_add_epi16(_mm_unpacklo_epi8(el8shr0, v_zero), _mm_unpacklo_epi8(el8shr3, v_zero)); __m128i el8 = _mm_add_epi16(el8shr12, el8shr03); __m128i el4h = _mm_add_epi16(_mm_unpackhi_epi16(el8, v_zero), _mm_unpacklo_epi16(el8, v_zero)); vsuml = _mm_add_epi32(vsuml, _mm_unpacklo_epi16(el8, v_zero)); vsumh = _mm_add_epi32(vsumh, el4h); _mm_storeu_si128((__m128i *)(sum_row + j), vsuml); _mm_storeu_si128((__m128i *)(sum_row + j + 4), vsumh); prev = _mm_add_epi32(prev, _mm_shuffle_epi32(el4h, _MM_SHUFFLE(3, 3, 3, 3))); } for (int v = sum_row[j - 1] - prev_sum_row[j - 1]; j < size.width; ++j) sum_row[j] = (v += src_row[j]) + prev_sum_row[j]; }
為了說明更方便，這里貼出我做的普通C語言的代碼和重新優化后的SSE代碼。
普通C語言：
void GetGrayIntegralImage(unsigned char *Src, int *Integral, int Width, int Height, int Stride) { memset(Integral, 0, (Width + 1) * sizeof(int)); // 第一行都為0 for (int Y = 0; Y < Height; Y++) { unsigned char *LinePS = Src + Y * Stride; int *LinePL = Integral + Y * (Width + 1) + 1; // 上一行位置 int *LinePD = Integral + (Y + 1) * (Width + 1) + 1; // 當前位置，注意每行的第一列的值都為0 LinePD[-1] = 0; // 第一列的值為0 for (int X = 0, Sum = 0; X < Width; X++) { Sum += LinePS[X]; // 行方向累加 LinePD[X] = LinePL[X] + Sum; // 更新積分圖 } } }
優化后的SSE算法：
void GetGrayIntegralImage(unsigned char *Src, int *Integral, int Width, int Height, int Stride) { memset(Integral, 0, (Width + 1) * sizeof(int)); // 第一行都為0 int BlockSize = 8, Block = Width / BlockSize; for (int Y = 0; Y < Height; Y++) { unsigned char *LinePS = Src + Y * Stride; int *LinePL = Integral + Y * (Width + 1) + 1; // 上一行位置 int *LinePD = Integral + (Y + 1) * (Width + 1) + 1; // 當前位置，注意每行的第一列的值都為0 LinePD[-1] = 0; __m128i PreV = _mm_setzero_si128(); __m128i Zero = _mm_setzero_si128(); for (int X = 0; X < Block * BlockSize; X += BlockSize) { __m128i Src_Shift0 = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(LinePS + X)), Zero); // A7 A6 A5 A4 A3 A2 A1 A0 __m128i Src_Shift1 = _mm_slli_si128(Src_Shift0, 2); // A6 A5 A4 A3 A2 A1 A0 0 __m128i Src_Shift2 = _mm_slli_si128(Src_Shift1, 2); // 移位改成基於Shift0，速度慢，Why？ // A5 A4 A3 A2 A1 A0 0 0 __m128i Src_Shift3 = _mm_slli_si128(Src_Shift2, 2); // A4 A3 A2 A1 A0 0 0 0 __m128i Shift_Add12 = _mm_add_epi16(Src_Shift1, Src_Shift2); // A6+A5 A5+A4 A4+A3 A3+A2 A2+A1 A1+A0 A0+0 0+0 __m128i Shift_Add03 = _mm_add_epi16(Src_Shift0, Src_Shift3); // A7+A4 A6+A3 A5+A2 A4+A1 A3+A0 A2+0 A1+0 A0+0 __m128i Low = _mm_add_epi16(Shift_Add12, Shift_Add03); // A7+A6+A5+A4 A6+A5+A4+A3 A5+A4+A3+A2 A4+A3+A2+A1 A3+A2+A1+A0 A2+A1+A0+0 A1+A0+0+0 A0+0+0+0 __m128i High = _mm_add_epi32(_mm_unpackhi_epi16(Low, Zero), _mm_unpacklo_epi16(Low, Zero)); // A7+A6+A5+A4+A3+A2+A1+A0 A6+A5+A4+A3+A2+A1+A0 A5+A4+A3+A2+A1+A0 A4+A3+A2+A1+A0 __m128i SumL = _mm_loadu_si128((__m128i *)(LinePL + X + 0)); __m128i SumH = _mm_loadu_si128((__m128i *)(LinePL + X + 4)); SumL = _mm_add_epi32(SumL, PreV); SumL = _mm_add_epi32(SumL, _mm_unpacklo_epi16(Low, Zero)); SumH = _mm_add_epi32(SumH, PreV); SumH = _mm_add_epi32(SumH, High); PreV = _mm_add_epi32(PreV, _mm_shuffle_epi32(High, _MM_SHUFFLE(3, 3, 3, 3))); _mm_storeu_si128((__m128i *)(LinePD + X + 0), SumL); _mm_storeu_si128((__m128i *)(LinePD + X + 4), SumH); } for (int X = Block * BlockSize, V = LinePD[X - 1] - LinePL[X - 1]; X < Width; X++) { V += LinePS[X]; LinePD[X] = V + LinePL[X]; }
}
我們先來解釋下這段代碼的SSE優化過程吧。
首先，用_mm_loadl_epi64一次性加載8個字節數據到XMM寄存器中，其中寄存器的高8位位0，此時寄存器的數據為:
高位 0 0 0 0 0 0 0 0 A7 A6 A5 A4 A3 A2 A1 A0 低位 (8位）
因為涉及到加法，並且最大為8個字節數據的加法，因此轉換到16位數據類型，使用_mm_unpacklo_epi8結合zero即可實現。
此時XMM寄存器內容變為：
Src_Shift0 A7 A6 A5 A4 A3 A2 A1 A0 (16位）
此后有3次移位分別得到：
Src_Shift1 A6 A5 A4 A3 A2 A1 A0 0 （16位）
Src_Shift2 A5 A4 A3 A2 A1 A0 0 0 （16位）
Src_Shift3 A4 A3 A2 A1 A0 0 0 0 （16位）
通過_mm_add_epi16分別對4組16位數據進行8次相加：
Shift_Add12 A6+A5 A5+A4 A4+A3 A3+A2 A2+A1 A1+A0 A0+0 0+0 （16位）
Shift_Add03 A7+A4 A6+A3 A5+A2 A4+A1 A3+A0 A2+0 A1+0 A0+0 （16位）
再對他們進行相加：
Low A7+A6+A5+A4 A6+A5+A4+A3 A5+A4+A3+A2 A4+A3+A2+A1 A3+A2+A1+A0 A2+A1+A0+0 A1+A0+0+0 A0+0+0+0
注意到低4位的16位數已經是連續相加的數據了，只要將他們轉換為32位就可以直接使用。
而通過 __m128i High = _mm_add_epi32(_mm_unpackhi_epi16(Low, Zero), _mm_unpacklo_epi16(Low, Zero)); 這一句則可以把前面的高4位連續相加的值拼接起來得到：
High A7+A6+A5+A4+A3+A2+A1+A0 A6+A5+A4+A3+A2+A1+A0 A5+A4+A3+A2+A1+A0 A4+A3+A2+A1+A0
后面的操作則順理成章了。
注意到我核心的改動在於原始代碼中的el8shr12和el8shr03的計算中的_mm_unpacklo_epi8被消除了，而在el8shr0一句中增加了一個_mm_unpacklo_epi8，因此少了3次這個函數，很明顯這樣做是不會改變計算結果的。
另外源代碼中的部分_mm_add_epi16被我用_mm_add_epi32代替了，這主要是因為用_mm_add_epi32意義更明顯，而且由於高位數據為0，他們的執行結果不會有任何區別。
還有一點在測試時發現，如果Src_Shift2，Src_Shift3的移位是基於Src_Shift0，即使用如下代碼:
__m128i Src_Shift2 = _mm_slli_si128(Src_Shift0, 4); __m128i Src_Shift3 = _mm_slli_si128(Src_Shift0, 6);
速度會有較為明顯的下降，難道說移動的位數多少和CPU的耗時有關？
以上是灰度模式的算法，在我的筆記本電腦上，SSE優化后的語句雖然增加了很多，但是執行效率約能提升30%，不過在一些PC上，普通的C和SSE優化后卻沒有啥速度區別了，這也不知道是為什么了。
如果是針對24位或者32位圖像，基本的優化思想是一致的，不過有更多的細節需要自己注意。
24位或者32位圖像在任何機器配置上，速度都能有30%的提升的。
還是感覺這種算法用文字很難表述清楚，用代碼再加上自己的空間組合可能更能理解吧。
