SSE圖像算法優化系列十三：超高速BoxBlur算法的實現和優化（Opencv的速度的五倍）

時間 2019-11-06

標籤 sse 圖像算法優化系列十三超高速 boxblur 實現 opencv 速度五倍简体版

原文原文鏈接

在SSE圖像算法優化系列五：超高速指數模糊算法的實現和優化（10000*10000在100ms左右實現）一文中，我曾經說過優化後的ExpBlur比BoxBlur還要快，那個時候我比較的BoxBlur算法是經過積分圖+SSE實現的，我在09年另一個博客帳號上曾經提供過一篇這個文章彩色圖像高速模糊之懶惰算法，裏面也介紹了一種快速的圖像模糊算法，這個算法的執行時間基本也是和半徑無關的。在今年的SSE優化學習之路上我曾經也考慮過將該算法使用SSE實現，但當時以爲這個算法逐像素同時逐行都是先後依賴的（單純的逐像素依賴算法我在指數模糊裏有提到如何用SSE優化），是沒法用SSE處理的，一直沒考慮，直到最近有朋友提出某個基於局部局方差的算法但願能提速，我又再次觸發靈感，終於將這種算法也實現的指令集實現，而且測試速度比積分圖快二倍，比解析opencv中Box Filter的實現並提出進一步加速的方案（源碼共享）這篇文章的速度快3倍，比opencv的cvSmooth函數快5倍，在一臺老舊的I3筆記本上處理3000*2000的灰度圖達到了6ms的速度，本文分享該優化過程並提供灰度版本的優化代碼供你們學習和討論。html

　　在彩色圖像高速模糊之懶惰算法一文附帶的代碼中（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的代碼的速度了。多線程

　　簡單的描述下各函數的做用先。ide

　　IM_GetMirrorPos函數主要是獲得某一個位置Pos按照鏡像的方式處理時在Length方向的座標，這裏Pos能夠爲負值，這個主要是用來得到後期的座標偏移的。　　　　　　函數

　　FillLeftAndRight_Mirror_C主要是用來獲取兩邊鏡像數據的（直接獲取，不調用IM_GetMirrorPos函數），好比好比Array原始數據爲 ***abcdefgh*** (Length = 8, Radius = 3)，則結果爲dcbabcdefghgfe。post

　　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了。