優化IPOL網站中基於DCT(離散余弦變換)的圖像去噪算法(附源代碼)。


     在您閱讀本文前,先需要告訴你的是:即使是本文優化過的算法,DCT去噪的計算量依舊很大,請不要向這個算法提出實時運行的苛刻要求。  

  言歸正傳,在IPOL網站中有一篇基於DCT的圖像去噪文章,具體的鏈接地址是:http://www.ipol.im/pub/art/2011/ys-dct/,IPOL網站的最大特點就是他的文章全部提供源代碼,而且可以基於網頁運行相關算法,得到結果。不過其里面的代碼本身是重實現論文的過程,基本不考慮速度的優化,因此,都相當的慢。

      這篇文章的原理也是非常簡單的,整個過程就是進行dct變換,然后在dct域進行硬閾值操作,再反變換回來。但是DCT變換不是針對整幅圖進行處理,而是基於每個像素點的領域(這里使用的8領域或者16領域),每次則移動一個像素。IPOL上提供的代碼函數也很少,但是一個關鍵的問題就是內存占用特別誇張,我們貼他的部分代碼:

 

// Transfer an image im of size width x height x channel to sliding patches of // size width_p x height_p xchannel. // The patches are stored in patches, where each ROW is a patch after being // reshaped to a vector.
void Image2Patches(vector<float>& im, vector< vector< vector< vector< float > > > >& patches, int width, int height, int channel, int width_p, int height_p) { int size1 = width * height; int counter_patch = 0; // Loop over the patch positions
    for (int j = 0; j < height - height_p + 1; j ++) for (int i = 0; i < width - width_p + 1; i ++) { int counter_pixel = 0; // loop over the pixels in the patch
            for (int kp = 0; kp < channel; kp++) for (int jp = 0; jp < height_p; jp ++) for (int ip = 0; ip < width_p; ip ++) { patches[counter_patch][kp][jp][ip] = im[kp*size1 + (j+jp)*width + i + ip]; counter_pixel ++; } counter_patch ++; } }

  看到這里的patches了,他的作用是保存每個點周邊的8*8領域的DCT變換的結果的,即使使用float類型的變量,也需要約Width * Height * 8 * 8 * sizeof(float)個字節的數組,假定寬度和高度都為1024的灰度圖,這個數據量為256MB,其實256MB的內存在現在機器來說毫無壓力,可這里要求是連續分布內存,則很有可能分配失敗,在外部表現的錯誤則是內存溢出。我們首先來解決這個問題。

  我們來看看原作者的代碼中patches的主要作用,見下面這部分代碼:

    // Loop over the patch positions
    for (int j = 0; j < height - height_p + 1; j ++) for (int i = 0; i < width - width_p + 1; i ++) { int counter_pixel = 0; // loop over the pixels in the patch
            for (int kp = 0; kp < channel; kp++) for (int jp = 0; jp < height_p; jp ++) for (int ip = 0; ip < width_p; ip ++) { im[kp*size1 + (j+jp)*width + i + ip] += patches[counter_patch][kp][jp][ip]; im_weight[kp*size1 + (j+jp)*width + i + ip] ++; counter_pixel ++; } counter_patch ++; } // Normalize by the weight
    for (int i = 0; i < size; i ++) im[i] = im[i] / im_weight[i];

  可見patches主要是為了保存每個點領域的DCT變換的數據,然后累加,上述四重循環外圍兩個是圖像的寬度和高度方向,內部兩重則是DCT的變換數據的行列方向,如果我們把DCT的行列方向的循環提到最外層,而把圖像的寬度和高度方向的循環放到內存,其一就是整個過程只需要一個8*8*sizeof(float)大小的重復利用的內存,第二就是正好符號了內部放大循環,外部放小循環的優化准則,在速度和內存占用上都有重大提高。

      我們來繼續優化,在獲取每個點的領域時,傳統的方式需要8*8個循環,那整個圖像就需要Width * Height * 8 * 8次了, 這個數據量太可觀了,在圖像處理中任意核卷積(matlab中conv2函數)的快速實現 一文中共享的代碼里提到了一種方式,大約只需要Width * Height * 8次讀取,並且其cache命中率還要高很多,具體的方式可參考本文附帶的代碼。

      繼續可以優化的地方就是8*8點的浮點DCT變換了。原文提供一維DCT變換的代碼如下:

for (int j = 0; j < 8; j ++) { out[j] = 0; for (int i = 0; i < 8; i ++) { out[j] += in[i] * DCTbasis[j][i]; } }

     就是兩個循環,共64次乘法和64次加法,上面的DCTbasis為固定的DCT系數矩陣。

  而一維的IDCT的變換代碼為:

for (int j = 0; j < PATCHSIZE; j ++) { out[j] = 0; for (int i = 0; i < PATCHSIZE; i ++) { out[j] += in[i] * DCTbasis[i][j]; } }

      和DCT唯一的區別僅僅是DCTbasis的行列坐標相反。

      這種代碼一看就想到了有SSE進行優化,PATCHSIZE為8 正好是兩個SSE浮點數m128的大小,乘法和加法都有對應的SSE函數,一次性可進行4個浮點加法和浮點乘法,效率當然會高很多,優化后的代碼如下所示:

/// <summary>
/// 8*8的一維DCT變換及其逆變換。 /// </summary>
/// <param name="In">輸入的數據。</param>
/// <param name="Out">輸出的數據。</param>
/// <param name="Invert">是否進行逆變換。</param>
/// <remarks> 1、輸入和輸出不能相同,即不支持in-place操作。</remarks>
/// <remarks> 2、算法只直接翻譯IPOL上的,利用了SSE加速。</remarks>
/// <remarks> 3、在JPEG壓縮等程序中的8*8DCT變換里優化的算法乘法比較少,但不好利用SSE優化,我用那里的代碼測試還比下面的慢。</remarks>
/// <remarks> 4、有關8*8 DCT優化的代碼見:http://insurgent.blog.hexun.com/1797358_d.html</remarks>
void DCT8X81D(float *In, float *Out, bool Invert) { __m128 Sum, A, B; const float *Data = (const float *)&Sum; A = _mm_load_ps(In); B = _mm_load_ps(In + 4); if (Invert == FALSE) { /*for (int Y = 0; Y < PATCHSIZE; Y ++) { Out[Y] = 0; for (int X = 0; X < PATCHSIZE; X ++) { Out[Y] += In[X] * DCTbasis[Y * PATCHSIZE + X]; } }*/ Sum = _mm_mul_ps(A, _mm_load_ps(DCTbasis)); Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(DCTbasis + 4))); Out[0] = Data[0] + Data[1] + Data[2] + Data[3];                            // 這里是否還有更好的方式實現呢?
 Sum = _mm_mul_ps(A, _mm_load_ps(DCTbasis + 8)); Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(DCTbasis + 12)));        // 不用一個Sum變量,用多個是不是還能提高指令級別的並行啊
        Out[1] = Data[0] + Data[1] + Data[2] + Data[3]; Sum = _mm_mul_ps(A, _mm_load_ps(DCTbasis + 16)); Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(DCTbasis + 20))); Out[2] = Data[0] + Data[1] + Data[2] + Data[3]; Sum = _mm_mul_ps(A, _mm_load_ps(DCTbasis + 24)); Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(DCTbasis + 28))); Out[3] = Data[0] + Data[1] + Data[2] + Data[3]; Sum = _mm_mul_ps(A, _mm_load_ps(DCTbasis + 32)); Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(DCTbasis + 36))); Out[4] = Data[0] + Data[1] + Data[2] + Data[3]; Sum = _mm_mul_ps(A, _mm_load_ps(DCTbasis + 40)); Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(DCTbasis + 44))); Out[5] = Data[0] + Data[1] + Data[2] + Data[3]; Sum = _mm_mul_ps(A, _mm_load_ps(DCTbasis + 48)); Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(DCTbasis + 52))); Out[6] = Data[0] + Data[1] + Data[2] + Data[3]; Sum = _mm_mul_ps(A, _mm_load_ps(DCTbasis + 56)); Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(DCTbasis + 60))); Out[7] = Data[0] + Data[1] + Data[2] + Data[3]; } else { /*for (int Y = 0; Y < PATCHSIZE; Y ++) { Out[Y] = 0; for (int X = 0; X < PATCHSIZE; X ++) { Out[Y] += In[X] * IDCTbasis[Y * PATCHSIZE + X]; } }*/ Sum = _mm_mul_ps(A, _mm_load_ps(IDCTbasis)); Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(IDCTbasis + 4))); Out[0] = Data[0] + Data[1] + Data[2] + Data[3]; Sum = _mm_mul_ps(A, _mm_load_ps(IDCTbasis + 8)); Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(IDCTbasis + 12))); Out[1] = Data[0] + Data[1] + Data[2] + Data[3]; Sum = _mm_mul_ps(A, _mm_load_ps(IDCTbasis + 16)); Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(IDCTbasis + 20))); Out[2] = Data[0] + Data[1] + Data[2] + Data[3]; Sum = _mm_mul_ps(A, _mm_load_ps(IDCTbasis + 24)); Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(IDCTbasis + 28))); Out[3] = Data[0] + Data[1] + Data[2] + Data[3]; Sum = _mm_mul_ps(A, _mm_load_ps(IDCTbasis + 32)); Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(IDCTbasis + 36))); Out[4] = Data[0] + Data[1] + Data[2] + Data[3]; Sum = _mm_mul_ps(A, _mm_load_ps(IDCTbasis + 40)); Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(IDCTbasis + 44))); Out[5] = Data[0] + Data[1] + Data[2] + Data[3]; Sum = _mm_mul_ps(A, _mm_load_ps(IDCTbasis + 48)); Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(IDCTbasis + 52))); Out[6] = Data[0] + Data[1] + Data[2] + Data[3]; Sum = _mm_mul_ps(A, _mm_load_ps(IDCTbasis + 56)); Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(IDCTbasis + 60))); Out[7] = Data[0] + Data[1] + Data[2] + Data[3]; } }

   當然,簡單的循環並不是效率最高的算法,在標准的JPG壓縮中就有着8*8的DCT變換,哪里的計算量有着更少的乘法和加法,在 http://insurgent.blog.hexun.com/1797358_d.html 中提出代碼里,只有32次乘法和更少的加法,但是由於這個代碼的向量性很差,是很難用SSE進行優化的,我實測的結果時他的代碼比我用SSE優化后的速度慢。

     當進行2維的DCT的時候,其步驟為對每行先進行行方向的一維DCT,然后對結果轉置,在對轉置后的數據進行行方向一維DCT,得到的結果再次轉置則為2維DCT。8*8的轉置雖然直接實現基本不存在cache miss的問題,不過還是用有關的SSE來實現未嘗不是個好主意,在intrinsic中提供了一個4*4浮點轉置的宏_MM_TRANSPOSE4_PS,我們對這個宏稍作修改,修改后的代碼如下:

//    http://stackoverflow.com/questions/5200338/a-cache-efficient-matrix-transpose-program
//    http://stackoverflow.com/questions/16737298/what-is-the-fastest-way-to-transpose-a-matrix-In-c
//    https://msdn.microsoft.com/en-us/library/4d3eabky(v=vs.90).aspx
//    高效的矩陣轉置算法,速度約為直接編寫的4倍, 整理時間 2015.10.12
inline void TransposeSSE4x4(float *Src, float *Dest) 
{
    __m128 Row0 = _mm_load_ps(Src);
    __m128 Row1 = _mm_load_ps(Src + 8);
    __m128 Row2 = _mm_load_ps(Src + 16);
    __m128 Row3 = _mm_load_ps(Src + 24);

    //        _MM_TRANSPOSE4_PS(Row0, Row1, Row2, Row3);                            //    或者使用這個SSE的宏

    __m128 Temp0    = _mm_shuffle_ps(Row0, Row1, _MM_SHUFFLE(1, 0, 1, 0));      //    Row0[0] Row0[1] Row1[0] Row1[1]    
    __m128 Temp2    = _mm_shuffle_ps(Row0, Row1, _MM_SHUFFLE(3, 2, 3, 2));      //    Row0[2] Row0[3] Row1[2] Row1[3]        
    __m128 Temp1    = _mm_shuffle_ps(Row2, Row3, _MM_SHUFFLE(1, 0, 1, 0));      //    Row2[0] Row2[1] Row3[0] Row3[1]        
    __m128 Temp3    = _mm_shuffle_ps(Row2, Row3, _MM_SHUFFLE(3, 2, 3, 2));        //    Row2[2] Row2[3] Row3[2] Row3[3]          

    Row0 = _mm_shuffle_ps(Temp0, Temp1, _MM_SHUFFLE(2, 0, 2, 0));                //    Row0[0] Row1[0] Row2[0] Row3[0]             
    Row1 = _mm_shuffle_ps(Temp0, Temp1, _MM_SHUFFLE(3, 1, 3, 1));                //    Row0[1] Row1[1] Row2[1] Row3[1]                     
    Row2 = _mm_shuffle_ps(Temp2, Temp3, _MM_SHUFFLE(2, 0, 2, 0));                //    Row0[2] Row1[2] Row2[2] Row3[2]                
    Row3 = _mm_shuffle_ps(Temp2, Temp3, _MM_SHUFFLE(3, 1, 3, 1));                //    Row0[3] Row1[3] Row2[3] Row3[3]             

    _mm_store_ps(Dest, Row0);
    _mm_store_ps(Dest + 8, Row1);
    _mm_store_ps(Dest + 16, Row2);
    _mm_store_ps(Dest + 24, Row3);
}

     本質上說,這些優化都是小打小鬧,提高不了多少速度,當然還有一些可以優化的地方,比如權重和累加和的更新,最后的累加和和權重的相除得到結果等等都有有關的SSE函數可以使用。

     還有一個可以優化的地方就是,在高度方向上前后兩個像素8*8領域 在進行2D的DCT變換時,其第一次行方向上的DCT變換有7行的結果是可以重復利用的,如果利用這一點,則可以獲得約15%的速度提示。

   綜合以上各種優化手段,在I5的機器上一幅512*512 的灰度圖像大約處理用時為100ms左右 ,比起IPOL的demo運行速度快了很多倍了。

     DCT濾波的效果上很多情況下也是相當不錯的,想比NLM也毫不遜色,我們貼一些圖片看下效果:

                         

    

                        噪音圖像                                                                                            去噪后效果(Sigma = 10)

     為顯示方便,上面的圖像都是設置了一定程度的縮放。

     共享我改寫的這個代碼供大家共同學習提高,如果哪位能進一步提高算法的速度 ,也希望不吝賜教。

  代碼下載鏈接:http://files.cnblogs.com/files/Imageshop/DCT_Denosing.rar

 

  后記:  繼續優化了下8*8點的DCT里SSE代碼的處理方式,改變了累加的方向,速度提高30%;然后把64次閾值處理的代碼也改成SSE優化,速度提高10%;在如果考慮進行隔行隔列取樣然后在DCT進行累加,經過測試基本上看不出有什么效果上的區別,但是速度大約能提高3.5倍左右;如果考慮多線程的方式,比如開個雙線程,速度約能提高0.8倍,如果CPU支撐AVX,則大概又能提高0.9倍,算來算去,我感覺可以實時了。

 ****************************作者: laviewpbt   時間: 2015.11.14    聯系QQ:  33184777 轉載請保留本行信息**********************

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM