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

時間 2019-11-06

標籤優化 ipol 網站基於 dct 離散餘弦變換圖像算法源代碼欄目網站開發简体版

原文原文鏈接

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

　　言歸正傳，在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)大小的重複利用的內存，第二就是正好符號了內部放大循環，外部放小循環的優化準則，在速度和內存佔用上都有重大提升。oop

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

繼續能夠優化的地方就是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> 一、輸入和輸出不能相同，即不支持in-place操做。</remarks>
/// <remarks> 二、算法只直接翻譯IPOL上的，利用了SSE加速。</remarks>
/// <remarks> 三、在JPEG壓縮等程序中的8*8DCT變換裏優化的算法乘法比較少，但很差利用SSE優化，我用那裏的代碼測試還比下面的慢。</remarks>
/// <remarks> 四、有關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 轉載請保留本行信息**********************