在SSE圖像算法優化系列五：超高速指數模糊算法的實現和優化（10000*10000在100ms左右實現） 一文中，我曾經說過優化后的ExpBlur比BoxBlur還要快，那個時候我比較的BoxBlur算法是通過積分圖+SSE實現的，我在09年另外一個博客賬號上曾經提供過一篇這個文章彩色圖像高速模糊之懶惰算法，里面也介紹了一種快速的圖像模糊算法，這個算法的執行時間基本也是和半徑無關的。在今年的SSE優化學習之路上我曾經也考慮過將該算法使用SSE實現，但當時覺得這個算法逐像素同時逐行都是前后依賴的（單純的逐像素依賴算法我在指數模糊里有提到如何用SSE優化），是無法用SSE處理的，一直沒考慮，直到最近有朋友提出某個基於局部局方差的算法希望能提速，我又再次觸發靈感，終於將這種算法也實現的指令集實現，並且測試速度比積分圖快二倍，比解析opencv中Box Filter的實現並提出進一步加速的方案（源碼共享）這篇文章的速度快3倍，比opencv的cvSmooth函數快5倍，在一台老舊的I3筆記本上處理3000*2000的灰度圖達到了6ms的速度，本文分享該優化過程並提供灰度版本的優化代碼供大家學習和討論。
在彩色圖像高速模糊之懶惰算法一文附帶的代碼中（VB6的代碼）是針對24位的圖像，為了討論方便，我們先貼出8位灰度的C++的代碼：
1 /// <summary> 2 /// 按照Tile方式進行數據的擴展，得到擴展后在原尺寸中的位置，以0為下標。 3 /// </summary> 4 int IM_GetMirrorPos(int Length, int Pos) 5 { 6 if (Pos < 0) 7 return -Pos; 8 else if (Pos >= Length) 9 return Length + Length - Pos - 2; 10 else 11 return Pos; 12 } 13 14 void FillLeftAndRight_Mirror_C(int *Array, int Length, int Radius) 15 { 16 for (int X = 0; X < Radius; X++) 17 { 18 Array[X] = Array[Radius + Radius - X]; 19 Array[Radius + Length + X] = Array[Radius + Length - X - 2]; 20 } 21 } 22 23 int SumofArray_C(int *Array, int Length) 24 { 25 int Sum = 0; 26 for (int X = 0; X < Length; X++) 27 { 28 Sum += Array[X]; 29 } 30 return Sum; 31 } 32 33 /// <summary> 34 /// 將整數限幅到字節數據類型。 35 /// </summary> 36 inline unsigned char IM_ClampToByte(int Value) // 現代PC還是這樣直接寫快些 37 { 38 if (Value < 0) 39 return 0; 40 else if (Value > 255) 41 return 255; 42 else 43 return (unsigned char)Value; 44 //return ((Value | ((signed int)(255 - Value) >> 31)) & ~((signed int)Value >> 31)); 45 } 46 47 int IM_BoxBlur_C(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Radius) 48 { 49 int Channel = Stride / Width; 50 if ((Src == NULL) || (Dest == NULL)) return IM_STATUS_NULLREFRENCE; 51 if ((Width <= 0) || (Height <= 0) || (Radius <= 0)) return IM_STATUS_INVALIDPARAMETER; 52 if (Radius < 1) return IM_STATUS_INVALIDPARAMETER; 53 if ((Channel != 1) && (Channel != 3) && (Channel != 4)) return IM_STATUS_NOTSUPPORTED; 54 55 Radius = IM_Min(IM_Min(Radius, Width - 1), Height - 1); // 由於鏡像的需求，要求半徑不能大於寬度或高度-1的數據 56 57 int SampleAmount = (2 * Radius + 1) * (2 * Radius + 1); 58 float Inv = 1.0 / SampleAmount; 59 60 int *ColValue = (int *)malloc((Width + Radius + Radius) * (Channel == 1 ? Channel : 4) * sizeof(int)); 61 int *ColOffset = (int *)malloc((Height + Radius + Radius) * sizeof(int)); 62 if ((ColValue == NULL) || (ColOffset == NULL)) 63 { 64 if (ColValue != NULL) free(ColValue); 65 if (ColOffset != NULL) free(ColOffset); 66 return IM_STATUS_OUTOFMEMORY; 67 } 68 for (int Y = 0; Y < Height + Radius + Radius; Y++) 69 ColOffset[Y] = IM_GetMirrorPos(Height, Y - Radius); 70 71 if (Channel == 1) 72 { 73 for (int Y = 0; Y < Height; Y++) 74 { 75 unsigned char *LinePD = Dest + Y * Stride; 76 if (Y == 0) 77 { 78 memset(ColValue + Radius, 0, Width * sizeof(int)); 79 for (int Z = -Radius; Z <= Radius; Z++) 80 { 81 unsigned char *LinePS = Src + ColOffset[Z + Radius] * Stride; 82 for (int X = 0; X < Width; X++) 83 { 84 ColValue[X + Radius] += LinePS[X]; // 更新列數據 85 } 86 } 87 } 88 else 89 { 90 unsigned char *RowMoveOut = Src + ColOffset[Y - 1] * Stride; // 即將減去的那一行的首地址 91 unsigned char *RowMoveIn = Src + ColOffset[Y + Radius + Radius] * Stride; // 即將加上的那一行的首地址 92 for (int X = 0; X < Width; X++) 93 { 94 ColValue[X + Radius] -= RowMoveOut[X] - RowMoveIn[X]; // 更新列數據 95 } 96 } 97 FillLeftAndRight_Mirror_C(ColValue, Width, Radius); // 鏡像填充左右數據 98 int LastSum = SumofArray_C(ColValue, Radius * 2 + 1); // 處理每行第一個數據 99 LinePD[0] = IM_ClampToByte(LastSum * Inv); 100 for (int X = 1; X < Width; X++) 101 { 102 int NewSum = LastSum - ColValue[X - 1] + ColValue[X + Radius + Radius]; 103 LinePD[X] = IM_ClampToByte(NewSum * Inv); 104 LastSum = NewSum; 105 } 106 } 107 } 108 else if (Channel == 3) 109 { 110 111 } 112 else if (Channel == 4) 113 { 114 115 } 116 free(ColValue); 117 free(ColOffset); 118 return IM_STATUS_OK; 119 }
以前沒意識到，就這么簡單的代碼用C寫后能產生速度也是很誘人的，3000*2000的圖能做到39ms，如果在編譯選項里勾選編譯器的“啟用增強指令集：流式處理 SIMD 擴展 2 (/arch:SSE2)”， 則系統會對上述具有浮點計算的部分使用相關的SIMD指令優化，如下圖所示：
這個時候3000*2000的圖能做到25ms,,基本上接近我改進的OPENCV的代碼的速度了。
簡單的描述下各函數的作用先。
IM_GetMirrorPos函數主要是得到某一個位置Pos按照鏡像的方式處理時在Length方向的坐標，這里Pos可以為負值，這個主要是用來獲得后期的坐標偏移的。
FillLeftAndRight_Mirror_C主要是用來獲取兩邊鏡像數據的（直接獲取，不調用IM_GetMirrorPos函數），比如比如Array原始數據為 ***abcdefgh*** (Length = 8, Radius = 3)，則結果為dcbabcdefghgfe。
SumofArray_C主要是計算一個數組的所有的元素的和。
IM_BoxBlur_C函數內部即為模糊的函數體，采用的優化思路整體和任意半徑中值濾波（擴展至百分比濾波器）O(1)時間復雜度算法的原理、實現及效果是一致的。當半徑為R時，方框模糊的值等於以某點為中心，左右上下各擴展R像素的的范圍內所有像素的綜合，像素總個數為（2*R+1）*（2*R+1）個，當然我們也可以把他分成（2*R+1）列，每列有（2*R+1）個元素本例的優化方式我們就是把累加數據分成一列一列的，充分利用重復信息來達到速度提升。
我們定義一個Radius + Width + Radius的內存數據ColValue，用來保存列方向的累加數據，對於第一行數據，需要做特殊處理，也就是用原始的方式計算一行像素所有元素在列方向上的和，
詳見78至於86行代碼，當然這里只計算了中間Width范圍內的數據，當不是第一行時，如下圖左邊所示，新的累加值只需減去移出的哪一行像素值同時加上移入的一行像素值，詳見90到96
行代碼。
上面只計算了中間Width范圍內的累加值，為了方便后續代碼的編寫以及使用SSE優化，下面的FillLeftAndRight_Mirror_C主要作用就是填充左邊和右邊分別填充數據，而且是按照鏡像的方式進行填充。
在更新了上述累加值后，我們開始處理計算均值了，對於每行的第一個點，SumofArray_C計算了前2*R + 1個列的累加值，這個累加值就是該點周邊（2*R+1）*（2*R+1）個元素的累積和，對於一行的其他像素，其實就類似於行方向列累加值的更新，減去移出的加入新進的列，如下圖右側圖所示，102到104行即實現了該過程。
原理基本上就是這樣，這個算法占用的額外內存很明顯很少，但是不支持In-Place操作。
為了追求速度，我們把整個過程能用SSE優化的地方都用SSE優化。
首先是第79至86行的數據，這個很容易優化：
for (int Z = -Radius; Z <= Radius; Z++) { unsigned char *LinePS = Src + ColOffset[Z + Radius] * Stride; int BlockSize = 8, Block = Width / BlockSize; for (int X = 0; X < Block * BlockSize; X += BlockSize) { int *DestP = ColValue + X + Radius; __m128i Sample = _mm_cvtepu8_epi16(_mm_loadl_epi64((__m128i *)(LinePS + X))); _mm_storeu_si128((__m128i *)DestP, _mm_add_epi32(_mm_loadu_si128((__m128i *)DestP), _mm_cvtepi16_epi32(Sample))); _mm_storeu_si128((__m128i *)(DestP + 4), _mm_add_epi32(_mm_loadu_si128((__m128i *)(DestP + 4)), _mm_unpackhi_epi16(Sample, _mm_setzero_si128()))); } // 處理剩余不能被SSE優化的數據 }
用_mm_loadl_epi64加載8個字節數據到內存，並用_mm_cvtepu8_epi16將其轉換為16位的__m128i變量，后面再把低4位和高4位的數據分別轉換成32位的，然后用_mm_add_epi32累加，注意后面我轉換函數用了兩種不同的方式，因為這里的數據絕對都是正數，因此也是可以使用_mm_cvtepi16_epi32和_mm_unpackhi_epi16組合Zero實現。
再來看看92到95行代碼，這個也很簡單。
int BlockSize = 8, Block = Width / BlockSize; __m128i Zero = _mm_setzero_si128(); for (int X = 0; X < Block * BlockSize; X += BlockSize) { int *DestP = ColValue + X + Radius; __m128i MoveOut = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(RowMoveOut + X)), Zero); __m128i MoveIn = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i *)(RowMoveIn + X)), Zero); __m128i Diff = _mm_sub_epi16(MoveIn, MoveOut); // 注意這個有負數也有正數的，有負數時轉換為32位是不能用_mm_unpackxx_epi16體系的函數 _mm_storeu_si128((__m128i *)DestP, _mm_add_epi32(_mm_loadu_si128((__m128i *)DestP), _mm_cvtepi16_epi32(Diff))); _mm_storeu_si128((__m128i *)(DestP + 4), _mm_add_epi32(_mm_loadu_si128((__m128i *)(DestP + 4)), _mm_cvtepi16_epi32(_mm_srli_si128(Diff, 8)))); } // 處理剩余不能被SSE優化的數據
這里也是一次性處理8個像素，這里我使用了另外一種轉換技巧來把8位轉換為16位，但是后面的Diff因為有正有負，要轉換為32位就必須使用_mm_cvtepi16_epi32來轉換，不能用unpack系列組合函數來實現，因為unpack不會移動符號位，我在這里吃了好幾次虧了。
接着是FillLeftAndRight_Mirror_C函數的優化，改寫如下：
void FillLeftAndRight_Mirror_SSE(int *Array, int Length, int Radius) { int BlockSize = 4, Block = Radius / BlockSize; for (int X = 0; X < Block * BlockSize; X += BlockSize) { __m128i SrcV1 = _mm_loadu_si128((__m128i *)(Array + Radius + Radius - X - 3)); __m128i SrcV2 = _mm_loadu_si128((__m128i *)(Array + Radius + Length - X - 5)); _mm_storeu_si128((__m128i *)(Array + X), _mm_shuffle_epi32(SrcV1, _MM_SHUFFLE(0, 1, 2, 3))); _mm_storeu_si128((__m128i *)(Array + Radius + Length + X), _mm_shuffle_epi32(SrcV2, _MM_SHUFFLE(0, 1, 2, 3))); } // 處理剩余不能被SSE優化的數據 }
鏡像就是倒序的過程，直接使用SSE的shuffle函數很方便實現。
計算數組的累加和優化也方便。
int SumofArray_SSE(int *Array, int Length) { int BlockSize = 8, Block = Length / BlockSize; __m128i Sum1 = _mm_setzero_si128(); __m128i Sum2 = _mm_setzero_si128(); for (int X = 0; X < Block * BlockSize; X += BlockSize) { Sum1 = _mm_add_epi32(Sum1, _mm_loadu_si128((__m128i *)(Array + X + 0))); Sum2 = _mm_add_epi32(Sum2, _mm_loadu_si128((__m128i *)(Array + X + 4))); } int Sum = SumI32(_mm_add_epi32(Sum1, Sum2)); // 處理剩余不能被SSE優化的數據 return Sum; }
使用兩個__m128i 變量主要是為了充分利用XMM寄存器，其中SumI32函數如下，主要是為了計算__m128i 內四個整數的累加值。
// Convert 4 32-bit integer values to 4 unsigned char values.
void _mm_storesi128_4char(__m128i Src, unsigned char *Dest) { __m128i T = _mm_packs_epi32(Src, Src); T = _mm_packus_epi16(T, T); *((int*)Dest) = _mm_cvtsi128_si32(T); }
代碼不解釋。
最后就是100到104行的代碼的轉換。
int BlockSize = 4, Block = (Width - 1) / BlockSize; __m128i OldSum = _mm_set1_epi32(LastSum); __m128 Inv128 = _mm_set1_ps(Inv); for (int X = 1; X < Block * BlockSize + 1; X += BlockSize) { __m128i ColValueOut = _mm_loadu_si128((__m128i *)(ColValue + X - 1)); __m128i ColValueIn = _mm_loadu_si128((__m128i *)(ColValue + X + Radius + Radius)); __m128i ColValueDiff = _mm_sub_epi32(ColValueIn, ColValueOut); // P3 P2 P1 P0 __m128i Value_Temp = _mm_add_epi32(ColValueDiff, _mm_slli_si128(ColValueDiff, 4)); // P3+P2 P2+P1 P1+P0 P0 __m128i Value = _mm_add_epi32(Value_Temp, _mm_slli_si128(Value_Temp, 8)); // P3+P2+P1+P0 P2+P1+P0 P1+P0 P0 __m128i NewSum = _mm_add_epi32(OldSum, Value); OldSum = _mm_shuffle_epi32(NewSum, _MM_SHUFFLE(3, 3, 3, 3)); // 重新賦值為最新值 __m128 Mean = _mm_mul_ps(_mm_cvtepi32_ps(NewSum), Inv128); _mm_storesi128_4char(_mm_cvtps_epi32(Mean), LinePD + X); }
以前認為這個算法難以使用SSE優化，主要難處就在這里，但是在學習了Opencv的積分圖技術時，這個問題就迎刃而解了，進一步分析還發現Opencv的代碼寫的並不完美，還有改進的空間，見上述代碼，使用兩次對同一數據移位就可以獲取累加，由P3 P2 P1 P0變為P3+P2+P1+P0 P2+P1+P0 P1+P0 P0。我們稍微分析一下。
__m128i ColValueDiff = _mm_sub_epi32(ColValueIn, ColValueOut); 這句代碼得到了移入和移出序列的差值，我們計為； P3 P2 P1 P0 （高位到低位）由於算法的累加特性，如果說OldSum的原始值為A3 A3 A3 A3，那么新的NewSum就應該是：
A3+P3+P2+P1+P0 A3+P2+P1+P0 A3+P1+P0 A3+P0；
__m128i Value_Temp = _mm_add_epi32(ColValueDiff, _mm_slli_si128(ColValueDiff, 4)); 這句的_mm_slli_si128得到中間結果 P2 P1 P0 0, 再和P3 P2 P1 P0相加得到
Value_Temp = P3+P2 P2+P1 P1+P0 P0
__m128i Value = _mm_add_epi32(Value_Temp, _mm_slli_si128(Value_Temp, 8));這句的_mm_slli_si128得到中間結果P1+P0 P0 0 0,再和P3+P2 P2+P1 P1+P0 P0相加得到：
Value = P3+P2+P1+P0 P2+P1+P0 P1+P0 P0
好簡單的過程。
OldSum = _mm_shuffle_epi32(NewSum, _MM_SHUFFLE(3, 3, 3, 3)); 這一句為什么要這樣寫，恐怕還是讀者自己體會比較好，似乎不好用文字表達。
最后一個_mm_storesi128_4char是我自己定義的一個將1個__m128i里面的4個32位整數轉換為字節數並保存到內存中的函數，詳見附件代碼。
至於24位圖像的優化，是個比較尷尬的處境，因為SSE一次性處理4個32位數，而24位每個像素只有3個分量，這種情況一般還是把他擴展到4個分量，比如說ColValue就可以改成4通道的，然后累積和也需要處理成4通道的，這當然需要一定的奇巧淫技，這里我不想多談，留點東西給自己。當然也可以考慮先把24位分解到3個灰度內存，然后利用灰度的算法進行計算，后面在合成到24位，這個分解和合成的過程也是可以使用SSE加速的，如果你結合多線程，還可以把3個灰度過程的處理並行化，這樣除了多占用內存外，速度比至二級處理24位要快（因為直接處理算法無法並行的，前后依賴的緣故）。另外就是最后在計算列累積求平均值的過程也變得更加自然了，不會出現灰度那種要在__mm128i內部進行累加的過程，而是直接得兩個SSE變量累加。
還說一點，現在大部分的CPU都支持AVX256了，還可以使用AVX進一步加速，似乎代碼該起來也不是很難，有興趣的朋友可以自己試試。
可以說，目前這個速度已經基本上達到了CPU的極限了，但是測試過IPP的速度，似乎比這個還要快點，不排除他使用了AVX，也不排除他使用多核的資源。
這個的優化對於BoxBlur來說是重要的一步，但是更重要的是他可以運用在很多場合，比如圖像的局部均方差計算，也可以使用類似的技術進行加速，兩幅圖像大的局部平方差也是可以這樣優化的，后續我會簡單的談下他在這方面加速的應用。
源代碼下載：https://files.cnblogs.com/files/Imageshop/FastBlur.rar
彩色圖工程：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar
不好意思，圖太小，速度為0ms了。
