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

基于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)

為顯示友善，上面的圖像都是設定了一定程度的縮放。

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

　　代碼下載下傳連結：https://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 轉載請保留本行資訊**********************

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

繼續閱讀

JPEG圖像密寫研究（一） JPEG圖像檔案結構JPEG圖像密寫研究（一） JPEG圖像檔案結構

DCT 離散餘弦變換原理與應用講解

DCT離散餘弦變換程式設計

蚊式噪聲

HEVC變換編碼實作代碼詳解

基于DCT和SVD的圖像哈希水印算法（論文閱讀）論文背景論文貢獻

二維離散變換DCT的一點簡單了解

python求音頻的梅爾倒譜系數寫在前面一、必備的關于音頻的知識二、音頻的梅爾倒譜系數的含義和求解步驟三、求梅爾倒譜系數詳細過程四、代碼總和

程式人生 - DCT、AT、CVT 到底哪個好？

一維2N浮點FDCT變換以及反變換

DCT變換公式一維DCT一維IDCT

JPEG編碼中的DCT與量化

DCT變換、DCT反變換、分塊DCT變換

MATLAB實作圖像的DCT、DFT變換與SSIM品質評估

DCT變換的基函數與基圖像1. 圖像變換的基函數2. 圖像變換的基圖像3.如何繪制基圖像