圖像處理之卷積---任意卷積核的快速實現

時間  2019-11-20
標籤 圖像處理 任意 快速 實現 简体版
原文   原文鏈接

　　 卷積實際上是圖像處理中最基本的操做，咱們常見的一些算法好比：均值模糊、高斯模糊、銳化、Sobel、拉普拉斯、prewitt邊緣檢測等等一些和領域相關的算法，均可以經過卷積算法實現。只不過因爲這些算法的卷積矩陣的特殊性，通常不會直接實現它，而是經過一些優化的手段讓計算量變小。可是有些狀況下卷積矩陣的元素值無甚規律或者有特殊要求，沒法經過常規手段優化，這個時候只能經過原始的方式實現。所以，如何快速的實現圖像的任意卷積矩陣操做也有必要作適當的研究。html

      目前，經過友人共享或本身搜索找到的一片關於任意覈算法優化的文章有： Reshuffling: A Fast Algorithm for Filtering with Arbitrary Kernels，改文章稱可以提升原始程序速度的40%左右，可是原始的程序是如何寫的也還不明白。算法

      在matlab中有幾個函數都與圖像卷積有關，好比imfilter就能夠實現卷積，或者 conv2也行，他們的速度都是至關快的，好比3000*3000的灰度圖，卷積矩陣大小爲15*15，在I5的CPU上運行時間只要170ms左右，至關的給力。編程

      在Celery的博客中，也提到了他的優化後的conv2和matlab至關甚至快於matlab，詳見http://blog.csdn.net/celerychen2009/article/details/38852105。ide

      因爲matlab的代碼中使用到了IPL庫進行加速，目前我寫的Conv2函數還沒法作到和其至關，對於任何核速度約爲matlab的一半。函數

      簡單的記錄下我在作卷積過程當中用到的優化吧。佈局

      原始的卷積的實現須要四重循環，簡單的表達以下：測試

複製代碼 
for (Y = 0; Y < Height; Y++)
{
    for (X = 0; X < Width; X++)
    {
        Index = .....;
        Sum = 0;
        for (XX = 0; XX < ConvW; XX++)
        {
            for (YY = 0; YY < ConvH; YY++)
            {
                Index1 = ..... ;
                Index2 = ..... ;
                Sum += Conv[Index1] * Pixel[Index2];
            }
        }
        Dest[Index] = Sum / Weight;
    }
}
複製代碼

　　當卷積矩陣較大時，計算量將會很大，並且因爲程序中的內存訪問很頻繁，cache miss現象比較嚴重，所以效率極爲低下。優化

     個人優化方法主要包括如下幾個方面： 編碼

     一：使用SSE進行乘法計算，因爲SSE能夠一次性進行4個單精度浮點數的計算，所以能夠有明顯的速度提高。spa

     二：經過適當的處理方式，對每一個取樣點周邊的卷積矩陣內的元素進行集中，使得每移動一個像素點不會須要從內存中進行大量的搜索工做。

     具體來講實現過程以下：

           一、爲了使用SSE的優點，首先將卷積矩陣進行調整，調整卷積矩陣一行的元素個數，使其爲不小於原始值的4的整數倍，而且讓新的卷積矩陣的內存佈局符合SSE相關函數的16字節對齊的要求。

　　實現代碼以下：

複製代碼 
float *Conv16 = (float *)_mm_malloc(PadConvLine * ConvH * sizeof(float), 16);                        //    保存16字節對齊的卷積矩陣，以方便使用SSE
                                                
for(Y = 0; Y < ConvH; Y++) 
{
    memcpy (Conv16 + Y * PadConvLine, Conv->Data.F + Y * ConvW , ConvW * sizeof(float));            //    複製卷積矩陣的數據
    memset(Conv16 + Y * PadConvLine + ConvW, 0, (PadConvLine - ConvW) * sizeof(float));                //    把冗餘部分的卷積數據設置爲0
}
複製代碼

    其中PadConvLine = Pad4(ConvW) 以及Pad4的原型爲： #define Pad4(bits) (((bits) + 3) / 4 * 4)；

    注意_mm_malloc函數分配的內存中的值是隨機值，對於擴展的部分必定要填充0，不然就會破壞卷積的結果。

    那麼若是咱們也同時得到了須要被卷積的部分數據的話（卷積核確定和卷積矩陣同樣大小，且也應該是16字節對齊的），能夠用以下的SSE的代碼進行乘法計算：

複製代碼 
float MultiplySSE(float *Kernel, float *Conv, int Length)
{
    int Block;  
    const float *Data;                        // 將SSE變量上的多個數值合併時所用指針.
    float Sum = 0;
    if (Length > 16)                        //    能夠進行四次SSE計算，測試代表，這個仍是快些的    
    {
        const int BlockWidth = 4 * 4;        // 塊寬. SSE寄存器能一次處理4個float，而後循環展開4次.
        Block = Length / BlockWidth;        // 塊數.    
        float *KernelP = Kernel, *ConvP = Conv;                // SSE批量處理時所用的指針.
        
        __m128 Sum0 = _mm_setzero_ps();         // 求和變量。SSE賦初值0
        __m128 Sum1 = _mm_setzero_ps();
        __m128 Sum2 = _mm_setzero_ps();
        __m128 Sum3 = _mm_setzero_ps();

        for(int I = 0; I < Block; I++)
        {
            Sum0 = _mm_add_ps(Sum0, _mm_mul_ps(_mm_load_ps(KernelP), _mm_load_ps(ConvP)));                    // SSE單精浮點緊縮加法
            Sum1 = _mm_add_ps(Sum1, _mm_mul_ps(_mm_load_ps(KernelP + 4), _mm_load_ps(ConvP + 4)));
            Sum2 = _mm_add_ps(Sum2, _mm_mul_ps(_mm_load_ps(KernelP + 8), _mm_load_ps(ConvP + 8)));
            Sum3 = _mm_add_ps(Sum3, _mm_mul_ps(_mm_load_ps(KernelP + 12), _mm_load_ps(ConvP + 12)));
            KernelP += BlockWidth;
            ConvP += BlockWidth;
        }
        
        Sum0 = _mm_add_ps(Sum0, Sum1);    // 兩兩合併(0~1).
        Sum2 = _mm_add_ps(Sum2, Sum3);    // 兩兩合併(2~3).
        Sum0 = _mm_add_ps(Sum0, Sum2);    // 兩兩合併(0~2).

        Data = (const float *)&Sum0;
        Sum = Data[0] + Data[1] + Data[2] + Data[3];

        Length = Length - Block * BlockWidth;            // 剩餘數量.
    }
    if (Length != 0)
    {
        const int BlockWidth = 4;                        //    程序已經保證了數量必然是4的倍數 
        Block = Length / BlockWidth;        
        float *KernelP = Kernel, *ConvP = Conv;                
        __m128 Sum0 = _mm_setzero_ps();        

        for(int I = 0; I < Block; I++)
        {
            Sum0 = _mm_add_ps(Sum0, _mm_mul_ps(_mm_load_ps(KernelP), _mm_load_ps(ConvP)));        
            KernelP += BlockWidth;
            ConvP += BlockWidth;
        }

        Data = (const float *)&Sum0;
        Sum += Data[0] + Data[1] + Data[2] + Data[3];
    }
    return Sum;
}
複製代碼

　　當卷積矩陣（擴充後）的元素數量大於16時，咱們採用了4路並行的SSE乘法實現，我在I3的CPU上測試時，2路SSE和4路SSE已經沒有啥大的區別了，而在I5的CPU上則4路仍是有較爲明顯的提升，所以採用4路SSE同時運行。固然1路SSE確定仍是比2路慢。另外，若是元素的數量少於16或者大於16但不能被16整除，那麼餘下的部分因爲先前的擴充，剩餘元素數量也確定是4的倍數，所以能夠用單路的SSE實現。 這也是編碼上的技巧。

     二、前面提到了須要被卷積的部分數據，這部分如何快速的獲取呢。觀察最原始的4重循環，其內部的2重即爲獲取須要被卷積的部分，可是這裏其實有不少問題。第一：因爲卷積取樣時必然有部分取樣點的座標在原始圖像的有效範圍外，所以必須進行判斷，耗時。第二：一樣爲了使用SSE,也必須把取樣的數據放在和擴充的卷積矩陣同樣大小的內存中。這裏我先貼出個人代碼在進行解釋具體的實現：

複製代碼 
IS_RET __stdcall Conv2(TImage *Src, TMatrix *Conv, TImage *Dest, EdgeMode Edge)
{
    if (Src == NULL || Dest == NULL || Conv == NULL) return IS_RET_ERR_PARA;
    if (Src->Width != Dest->Width || Src->Height != Dest->Height || Src->BitCount != Dest->BitCount || Src->Stride != Dest->Stride) return IS_RET_ERR_PARA;
    if (Src->Scan0 == NULL || Dest->Scan0 == NULL || Conv->Data.F == NULL) return IS_RET_ERR_MEM;
    if (Conv->Width < 1 || Conv->Height < 1) return IS_RET_ERR_PARA;
    
    int Width = Src->Width, Height = Src->Height, Stride = Src->Stride;
    int ConvW = Conv->Width, ConvH = Conv->Height;
    unsigned char *PtSrc = Src->Scan0, *PtDest = Dest->Scan0;


    if (Src->BitCount == 24)
    {

    }
    else
    {
        int Left = ConvW / 2, Top = ConvH / 2, Right = ConvW - Left - 1, Bottom = ConvH - Top - 1, ExpHeight = Height + ConvH - 1;        //    注意核中心那個元素不用擴展，好比核的寬度爲3，則只要左右各擴展一個像素就能夠了
        int PadConvLine = Pad4(ConvW), Length = PadConvLine * ConvH;
        int X, Y, IndexD, IndexE, IndexK, ExpStride;
        float *CurKer, Inv, Sum = 0;
        unsigned char *PtExp, *PtDest;

        TImage *Expand;
        IS_RET Ret = GetPadImage(Src, &Expand, Left, Right, Top, Bottom, Edge);                                //    獲得擴展後的數據，能夠提速和方便編程，可是多佔用一分內存
        if (Ret != IS_RET_OK) return Ret;
        
        PtExp = Expand->Scan0; PtDest = Dest->Scan0; ExpStride = Expand->Stride;
        
        for (X = 0; X < ConvH * ConvW; X ++) Sum += Conv->Data.F[X];
        Inv = (Sum == 0 ? 1: 1 / Sum);                                                                        //    若是卷積舉證的和爲0，則設置爲1

        float *Conv16 = (float *)_mm_malloc(PadConvLine * ConvH * sizeof(float), 16);                        //    保存16字節對齊的卷積矩陣，以方便使用SSE
        float *Kernel = (float *)_mm_malloc(PadConvLine * ExpHeight * sizeof(float), 16);                    //    保存16字節對齊的卷積核矩陣，以方便使用SSE
                                        
        for(Y = 0; Y < ConvH; Y++) 
        {
            memcpy (Conv16 + Y * PadConvLine, Conv->Data.F + Y * ConvW , ConvW * sizeof(float));            //    複製卷積矩陣的數據
            memset(Conv16 + Y * PadConvLine + ConvW, 0, (PadConvLine - ConvW) * sizeof(float));                //    把冗餘部分的卷積數據設置爲0
        }
        
        for (Y = 0; Y < ExpHeight; Y++)
        {
            IndexE = Y * ExpStride;
            CurKer = Kernel + Y * PadConvLine;                        //    計算第一列全部像素將要取樣的卷積核數據
            for (X = 0; X < ConvW; X++)
            {
                CurKer[X] = PtExp[IndexE++];
            }
        }
        for (X = 0 ; X < Width ; X ++)
        {
            if (X != 0)                                                //    若是不是第一列，須要更新卷積核的數據
            {
                memcpy(Kernel, Kernel + 1, (PadConvLine * ExpHeight - 1) * sizeof(float));    //    往前移動一個數據
                IndexK = ConvW - 1 ;
                IndexE = IndexK + X;
                for (Y = 0; Y < ExpHeight; Y++)
                {
                    Kernel[IndexK] = PtExp[IndexE];        //    只要刷新下一個元素
                    IndexK += PadConvLine;
                    IndexE += ExpStride;
                }
            }
            
            CurKer = Kernel;    IndexD = X;
            for (Y = 0; Y < Height; Y ++)                            //    沿列的方向進行更新
            {
                PtDest[IndexD] = Clamp((int)( MultiplySSE(Conv16, CurKer, Length) * Inv + 0.5));        //    直接把函數放在這裏也沒有啥提速的，注意改函數不會被內聯的
                CurKer += PadConvLine;
                IndexD += Stride;
            }
        }
        _mm_free(Conv16);
        _mm_free(Kernel);
        FreeImage(Expand);
        return IS_RET_OK;
    }
}
複製代碼

     對於第一個問題，解決的方式很簡答，即用空間換時間，新建一副（Width + ConvW - 1, Height + ConvH -1)大小的圖像，而後四周的ConvW及ConvH的像素用邊緣的值或者邊緣鏡像的值填充，正中間的則用原來的圖複製過來，這樣操做後進行取樣時再也不原圖取樣，而在這福擴展的圖中取樣，就避免了座標判斷等if語句的跳轉耗時了，上GetPadImage即實現了改功能。

     第二個問題則須要有必定的實現技巧，咱們分配一塊PadConvLine * (Height + ConvH - 1) 大小的內存，而後計算原圖第一列像素串聯起來的須要卷積的部分的數據，這一部分代碼如上述44-52行所示。有了這樣的數據，若是須要計算第一列的卷積結果，則很簡單了，每跳過一列則把被卷積的數據起點增長PadConvLine個元素，在調用上述MultiplySSE函數得到卷積結果。接着則計算第二列像素的卷積值，此時須要總體更新這一列像素串聯起來的須要被卷積的數據，更新也很簡單，就是把原來的數據總體向左移動一個像素，這個能夠用memcpy快速實現，而後在填充入新進來的那個元素，就ok了，接着就是再次調用MultiplySSE函數，如此重複下去。

     通過編碼測試，對於3000*3000的灰度圖，15*15的核在I5的CPU上的測試平均結果爲360ms,比matlab的慢了一半。

     最後說明一點，不少人都說用FFT能夠快速的實現卷積，而且是O(1)的，我比較贊成後半句，可是前面半句是絕對的有問題的，至少在覈小於50*50時，FFT實現的卷積不會比直接實現塊。要知道FFT的計算量實際上是很大的。

     http://www.cnblogs.com/Imageshop/p/4126753.html

http://blog.csdn.net/celerychen2009/article/details/38852105

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