SSE圖像算法優化系列二十一：基於DCT變換圖像去噪算法的進一步優化（100W像素30ms）。

時間  2019-11-06
標籤 sse 圖像 算法 優化 系列 二十一 基於 dct 變換 進一步 100w 像素 30ms 简体版
原文   原文鏈接

　　在優化IPOL網站中基於DCT（離散餘弦變換）的圖像去噪算法（附源代碼） 一文中，咱們曾經優化過基於DCT變換的圖像去噪算法，在那文所提供的Demo中，處理一副1000*1000左右的灰度噪音圖像耗時約450ms，若是採用所謂的快速模式耗時約150ms，說實在的，這個速度確實仍是有點慢，後續曾嘗試用AVX優化，可是感受AVX真的沒有SSE用的方便，並且AVX裏還有很多陷阱，本覺得這個算法優化沒有什麼但願了，但前幾日網友推薦了一片論文《Randomized Redundant DCT Efficent Denoising by Using Random Subsampling of DCT Patches》裏提供了比較誘人的數據，因而我又對這個算法有了新的興趣，那咱們來看看這篇論文有何特色。html

 　　先來看看論文提供的圖片和數據：算法

      

      其中R-DCT即《DCT Image Denoising a Simple and Eective Image Denoising Algorithm》一文的經典算法，而RR-DCT則是論文提供的速度，23ms對768*512的彩色像來講應該說是比較快的了。縱觀全文的內容，其提出的加速的方式其實相似於我在博文中提出的取樣，另外他提出了使用多線程的方式來加速。編程

　　可是，其取樣的方式並非我在博文中提出的規則的隔行隔列的方式，而是一種更爲隨機但也不會太隨機的方式，稱爲Poisson-disk sampling (PDS) 取樣，這種取樣他能夠指定每一個取樣點最遠距離，而又不是徹底均勻分佈，這既能保證每一個像素都能被取樣到，又保證每一個像素不會被取樣過多。固然，能這樣作的主要緣由仍是基於DCT算法的去噪，好比8*8的塊，裏面存在較多的數據冗餘。數組

　　在這篇論文提供了相關代碼，日本人寫的，我覺代碼寫的我想吐，很難受，我如今還在吐，可是這個代碼也不是徹底一無可取，其中關於DCT的優化無心中給我繼續優化這個算法帶來了極大的但願。多線程

　　好，那我接着說個人優化思路，在以前的博文中，咱們在8x8的DCT優化中留下了一個疑問，如何獲得SSE變量裏的四個浮點數的累加值，即以下代碼：dom

        const float *Data = (const float *)∑
        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];                            //    這裏是否還有更好的方式實現呢?

　　上述代碼使用了SSE指令和FPU指令共用，不是不行，而是確實效率不是很好，之前作這個算法時對SIMD指令集的瞭解還不是很深入，所以，如今看來這裏有很好的解決方案，爲表述方便，咱們把原始的8x8的DCT變換C代碼貼出以下：ide

//    普通C語言版本的8x1的一維DCT變換
__forceinline void IM_DCT1D_8x1_C(float *In, float *Out)
{
    for (int Y = 0, Index = 0; Y < 8; Y++)
    {
        float Sum = 0;
        for (int X = 0; X < 8; X++)
        {
            Sum += In[X] * DCTbasis[Index++];
        }
        Out[Y] = Sum;
    }
}

　　其中的DCTbasis是預先定義好的64個元素的浮點數組。函數

　　Out的每一個元素等於8個浮點數的和，可是是水平方向元素的和，咱們知道SSE裏有個指令叫_mm_hadd_ps，他就是執行的水平方向元素相加，他執行的相似下圖的操做：post

　　咱們的Sum是8個浮點的相加，8個浮點2個SSE變量可保存下，執行四次水平加法就能夠獲得結果，同時咱們在考慮到結合後續的幾行數據，因而就有下面的優化代碼：測試

__forceinline void IM_DCT1D_8x1_SSE(float *In, float *Out)
{
    __m128 In1 = _mm_loadu_ps(In);
    __m128 In2 = _mm_loadu_ps(In + 4);

    __m128 Sum0_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + 0));
    __m128 Sum0_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + 4));
    __m128 Sum1_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + 8));
    __m128 Sum1_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + 12));
    __m128 Sum2_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + 16));
    __m128 Sum2_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + 20));
    __m128 Sum3_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + 24));
    __m128 Sum3_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + 28));

    __m128 Sum01 = _mm_hadd_ps(_mm_hadd_ps(Sum0_0, Sum0_1), _mm_hadd_ps(Sum1_0, Sum1_1));            //    使用_mm_hadd_ps提升速度
    __m128 Sum23 = _mm_hadd_ps(_mm_hadd_ps(Sum2_0, Sum2_1), _mm_hadd_ps(Sum3_0, Sum3_1));

    _mm_storeu_ps(Out + 0, _mm_hadd_ps(Sum01, Sum23));

    __m128 Sum4_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + 32));
    __m128 Sum4_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + 36));
    __m128 Sum5_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + 40));
    __m128 Sum5_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + 44));
    __m128 Sum6_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + 48));
    __m128 Sum6_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + 52));
    __m128 Sum7_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + 56));
    __m128 Sum7_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + 60));

    __m128 Sum45 = _mm_hadd_ps(_mm_hadd_ps(Sum4_0, Sum4_1), _mm_hadd_ps(Sum5_0, Sum5_1));
    __m128 Sum67 = _mm_hadd_ps(_mm_hadd_ps(Sum6_0, Sum6_1), _mm_hadd_ps(Sum7_0, Sum7_1));

    _mm_storeu_ps(Out + 4, _mm_hadd_ps(Sum45, Sum67));
}

 　　對於DCT逆變換，也作相似的處理，100萬速度一會兒就能提高到350ms，快速版本約爲120ms，也是一個長足的進度。

　　 在DCT變換完成後，爲了進行去噪，須要一個閾值的處理，過程，普通的C語言相似下述過程：

    for (int YY = 0; YY < 64; YY++)    
        if (abs(DctOut[YY] )< Threshold) DctOut[YY]  = 0;            //    閾值處理

　　這也徹底能夠轉換成SSE處理：

//    DCT閾值處理
__forceinline void IM_DctThreshold8x8(float *DctOut, float Threshold)
{
    __m128 Th = _mm_set_ps1(Threshold);

    __m128 DctOut0 = _mm_load_ps(DctOut + 0);
    __m128 DctOut1 = _mm_load_ps(DctOut + 4);
    _mm_store_ps(DctOut + 0, _mm_blend_ps(_mm_and_ps(DctOut0, _mm_cmpgt_ps(_mm_abs_ps(DctOut0), Th)), DctOut0, 1));        //    // keep 00 coef.
    _mm_store_ps(DctOut + 4, _mm_and_ps(DctOut1, _mm_cmpgt_ps(_mm_abs_ps(DctOut1), Th)));

    __m128 DctOut2 = _mm_load_ps(DctOut + 8);
    __m128 DctOut3 = _mm_load_ps(DctOut + 12);
    _mm_store_ps(DctOut + 8, _mm_and_ps(DctOut2, _mm_cmpgt_ps(_mm_abs_ps(DctOut2), Th)));
    _mm_store_ps(DctOut + 12, _mm_and_ps(DctOut3, _mm_cmpgt_ps(_mm_abs_ps(DctOut3), Th)));

    __m128 DctOut4 = _mm_load_ps(DctOut + 16);
    __m128 DctOut5 = _mm_load_ps(DctOut + 20);
    _mm_store_ps(DctOut + 16, _mm_and_ps(DctOut4, _mm_cmpgt_ps(_mm_abs_ps(DctOut4), Th)));
    _mm_store_ps(DctOut + 20, _mm_and_ps(DctOut5, _mm_cmpgt_ps(_mm_abs_ps(DctOut5), Th)));

    __m128 DctOut6 = _mm_load_ps(DctOut + 24);
    __m128 DctOut7 = _mm_load_ps(DctOut + 28);

    _mm_store_ps(DctOut + 24, _mm_and_ps(DctOut6, _mm_cmpgt_ps(_mm_abs_ps(DctOut6), Th)));
    _mm_store_ps(DctOut + 28, _mm_and_ps(DctOut7, _mm_cmpgt_ps(_mm_abs_ps(DctOut7), Th)));

    __m128 DctOut8 = _mm_load_ps(DctOut + 32);
    __m128 DctOut9 = _mm_load_ps(DctOut + 36);
    _mm_store_ps(DctOut + 32, _mm_and_ps(DctOut8, _mm_cmpgt_ps(_mm_abs_ps(DctOut8), Th)));
    _mm_store_ps(DctOut + 36, _mm_and_ps(DctOut9, _mm_cmpgt_ps(_mm_abs_ps(DctOut9), Th)));

    __m128 DctOut10 = _mm_load_ps(DctOut + 40);
    __m128 DctOut11 = _mm_load_ps(DctOut + 44);
    _mm_store_ps(DctOut + 40, _mm_and_ps(DctOut10, _mm_cmpgt_ps(_mm_abs_ps(DctOut10), Th)));
    _mm_store_ps(DctOut + 44, _mm_and_ps(DctOut11, _mm_cmpgt_ps(_mm_abs_ps(DctOut11), Th)));

    __m128 DctOut12 = _mm_load_ps(DctOut + 48);
    __m128 DctOut13 = _mm_load_ps(DctOut + 52);
    _mm_store_ps(DctOut + 48, _mm_and_ps(DctOut12, _mm_cmpgt_ps(_mm_abs_ps(DctOut12), Th)));
    _mm_store_ps(DctOut + 52, _mm_and_ps(DctOut13, _mm_cmpgt_ps(_mm_abs_ps(DctOut13), Th)));

    __m128 DctOut14 = _mm_load_ps(DctOut + 56);
    __m128 DctOut15 = _mm_load_ps(DctOut + 60);
    _mm_store_ps(DctOut + 56, _mm_and_ps(DctOut14, _mm_cmpgt_ps(_mm_abs_ps(DctOut14), Th)));
    _mm_store_ps(DctOut + 60, _mm_and_ps(DctOut15, _mm_cmpgt_ps(_mm_abs_ps(DctOut15), Th)));
}

　　充分利用了_mm_cmpgt_ps比較函數的返回值和_mm_and_ps的位運算功能。

　　後續的權重累積以及DCT值的累加也同樣能夠用SSE優化，也是很簡單的事情。

       最後的DCT累計值和權重相除便可獲得最終的結果值，即以下代碼：

        for (int Y = 0, Index = 0; Y < Height; Y++)
        {
            unsigned char *LinePD = Dest + Y * Stride;
            for (int X = 0; X < Width; X++)
            {
                LinePD[X] = ClampToByte(Sum[Index] / Weight[Index]);        
                Index++;
            }
        }

　　也是很是容易能產生高效的並行化的代碼的。

//    完整的累加值和權重相除
__forceinline void IM_SumDivWeight2uchar8x8(float *Sum, float *Weight, unsigned char *Dest, int Width, int Height, int Stride)
{
    int BlockSize = 8, Block = Width / BlockSize;
    for (int Y = 0, Index = 0; Y < Height; Y++)
    {
        unsigned char *LinePD = Dest + Y * Stride;
        for (int X = 0; X < Block * BlockSize; X += BlockSize, Index += BlockSize)
        {
            __m128 Div1 = _mm_div_ps(_mm_loadu_ps(Sum + Index + 0), _mm_loadu_ps(Weight + Index + 0));        //    即便Sum或Weight是16字節對齊的，但若是Width是奇數，因爲最後幾個像素的做用，使得並非每次加載數據時都是16字節對齊的內存，因此只能用loadu而不能用load
            __m128 Div2 = _mm_div_ps(_mm_loadu_ps(Sum + Index + 4), _mm_loadu_ps(Weight + Index + 4));
            __m128i Result = _mm_packus_epi32(_mm_cvttps_epi32(Div1), _mm_cvttps_epi32(Div2));
            Result = _mm_packus_epi16(Result, Result);
            _mm_storel_epi64((__m128i *)(LinePD + X), Result);
        }
        for (int X = Block * BlockSize; X < Width; X++, Index++)
        {
            LinePD[X] = IM_ClampToByte((int)(Sum[Index] / Weight[Index] + 0.4999999f));
        }
    }
}

　　使用_mm_packus_epi16等飽和運算的相關指令能有效的避免跳轉，從而適當的加快速度。

       那麼另一點，在以前博文中提到了須要進行DCT變換原始數據的獲取，是經過相似滑塊的方式處理的，其實咱們感受不是須要，咱們能夠藉助SSE直接將8個字節數組轉換爲浮點，而後保存。

//    將8個字節數據轉換位浮點數
__forceinline void IM_Convert8ucharTo8float(unsigned char *Src, float *Dest)
{
    __m128i SrcV = _mm_loadl_epi64((__m128i *)Src);
    _mm_storeu_ps(Dest + 0, _mm_cvtepi32_ps(_mm_shuffle_epi8(SrcV, _mm_setr_epi8(0, -1, -1, -1, 1, -1, -1, -1, 2, -1, -1, -1, 3, -1, -1, -1))));
    _mm_storeu_ps(Dest + 4, _mm_cvtepi32_ps(_mm_shuffle_epi8(SrcV, _mm_setr_epi8(4, -1, -1, -1, 5, -1, -1, -1, 6, -1, -1, -1, 7, -1, -1, -1))));
}

　　更進一步，咱們能夠把他和DCT1D的代碼結合起來，直接實現由字節數據計算出對應的DCT1D結果，這樣減小中間的一次保存一次加載時間，以下所示：

//    直接由字節數據執行一維8X8的DCT變換
__forceinline void IM_DCT1D_8x1_SSE(unsigned char *In, float *Out)
{
    __m128i SrcV = _mm_loadl_epi64((__m128i *)In);
    __m128 In1 = _mm_cvtepi32_ps(_mm_shuffle_epi8(SrcV, _mm_setr_epi8(0, -1, -1, -1, 1, -1, -1, -1, 2, -1, -1, -1, 3, -1, -1, -1)));
    __m128 In2 = _mm_cvtepi32_ps(_mm_shuffle_epi8(SrcV, _mm_setr_epi8(4, -1, -1, -1, 5, -1, -1, -1, 6, -1, -1, -1, 7, -1, -1, -1)));

    __m128 Sum0_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + 0));
    __m128 Sum0_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + 4));
    __m128 Sum1_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + 8));
    __m128 Sum1_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + 12));
    __m128 Sum2_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + 16));
    __m128 Sum2_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + 20));
    __m128 Sum3_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + 24));
    __m128 Sum3_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + 28));

    __m128 Sum01 = _mm_hadd_ps(_mm_hadd_ps(Sum0_0, Sum0_1), _mm_hadd_ps(Sum1_0, Sum1_1));
    __m128 Sum23 = _mm_hadd_ps(_mm_hadd_ps(Sum2_0, Sum2_1), _mm_hadd_ps(Sum3_0, Sum3_1));

    _mm_storeu_ps(Out + 0, _mm_hadd_ps(Sum01, Sum23));

    __m128 Sum4_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + 32));
    __m128 Sum4_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + 36));
    __m128 Sum5_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + 40));
    __m128 Sum5_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + 44));
    __m128 Sum6_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + 48));
    __m128 Sum6_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + 52));
    __m128 Sum7_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + 56));
    __m128 Sum7_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + 60));

    __m128 Sum45 = _mm_hadd_ps(_mm_hadd_ps(Sum4_0, Sum4_1), _mm_hadd_ps(Sum5_0, Sum5_1));
    __m128 Sum67 = _mm_hadd_ps(_mm_hadd_ps(Sum6_0, Sum6_1), _mm_hadd_ps(Sum7_0, Sum7_1));

    _mm_storeu_ps(Out + 4, _mm_hadd_ps(Sum45, Sum67));
}

 　　結合在上一篇博文中咱們談到的充分利用兩次DCT2D變換中第一次DCT1D有7行結果是重複的現象，咱們新的一個版本的優化代碼出爐：

int IM_DCT_Denoising_02(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, float Sigma, bool FastMode)
{
    int Channel = Stride / Width;
    if ((Src == NULL) || (Dest == NULL))                                            return IM_STATUS_NULLREFRENCE;
    if ((Width <= 7) || (Height <= 7) || (Sigma <= 0))                                return IM_STATUS_INVALIDPARAMETER;
    if ((Channel != 1) && (Channel != 3))                                            return IM_STATUS_NOTSUPPORTED;

    int Status = IM_STATUS_OK;

    if (Channel == 3)
    {
                ////
    }
    else
    {
        int Step = (FastMode == true ? 2 : 1);
        float Threshold = 3 * Sigma;
        float *Dct1D = (float *)_mm_malloc(8 * Height * sizeof(float), 32);
        float *DctOut = (float *)_mm_malloc(64 * sizeof(float), 32);
        float *Sum = (float *)_mm_malloc(Width * Height * sizeof(float), 32);
        float *Weight = (float *)_mm_malloc(Width * Height * sizeof(float), 32);

        if ((Dct1D == NULL) || (DctOut == NULL) || (Sum == NULL) || (Weight == NULL))
        {
            Status = IM_STATUS_OUTOFMEMORY;
            goto FreeMemory;
        }
        memset(Sum, 0, Width * Height * sizeof(float));
        for (int X = 0; X < Width - 7; X += Step)
        {
            for (int Y = 0; Y < Height; Y++)
            {
                IM_DCT1D_8x1_SSE(Src + Y * Stride + X, Dct1D + Y * 8);
            }
            for (int Y = 0; Y < Height - 7; Y += Step)
            {
                IM_DCT2D_8x8_With1DRowDCT_SSE(Dct1D + Y * 8, DctOut);                                    //    DCT變換
                IM_DctThreshold8x8(DctOut, Threshold);                                                    //    閾值處理
                IM_IDCT2D_8x8_SSE(DctOut, DctOut);                                                        //    在反變換回來
                IM_UpdateSumAndWeight8x8(Sum + Y * Width + X, Weight + Y * Width + X, DctOut, Width);    //    更新權重和閾值
            }
        }
        IM_SumDivWeight2uchar8x8(Sum, Weight, Dest, Width, Height, Stride);
    FreeMemory:
        if (Dct1D != NULL)        _mm_free(Dct1D);
        if (DctOut != NULL)        _mm_free(DctOut);
        if (Sum != NULL)        _mm_free(Sum);
        if (Weight != NULL)        _mm_free(Weight);
        return Status;
    }
}

　　這樣速度也不會比以前的慢，測試100W像素大約320ms，快速模式90ms，並且佔用內存相比以前更小。

　　進一步考慮，其實在整個過程當中Weight的值是能夠預測，也就是說不一樣位置的Weight值實際上是固定的，咱們沒有必要去計算，在圖像中心，這個值必定等於64，若是是快速模式，則等於16，惟一須要注意處理的是邊緣部分的像素，他們的權重是變化的。

　　這樣處理的好處是即減小了Weight的權重部分佔用的內存，也減小了計算量，並且部分計算能夠用乘法代替除法，速度也會有顯著提高。

       如此改進後測試100W像素大約270ms，快速模式78ms，內存佔用更小。

　　到此爲止，彷佛已經沒有什麼可以進一步獲取速度的提高了，我也差點絕望了，由於最爲耗時的DCT部分的進一步提高從C的實現上看，對應的SSE語句彷佛已經到了極致。

　　翻閱《Randomized Redundant DCT Efficent Denoising by Using Random Subsampling of DCT Patches》一文的官方網站，下載了其對應的代碼，在代碼裏找到了dct_simd.cpp模塊，裏面的代碼也很是之亂，但觀察其有16x16,8x8,4x4的DCT的多種實現，會不會比IPOL裏這種直接乘的更快呢，找到了我關心的8x8的DCT部分代碼，發現了其C算法的1D實現：

static void fdct81d_GT(const float *src, float *dst)
{
    for (int i = 0; i < 2; i++)
    {
        const float mx00 = src[0] + src[7];
        const float mx01 = src[1] + src[6];
        const float mx02 = src[2] + src[5];
        const float mx03 = src[3] + src[4];
        const float mx04 = src[0] - src[7];
        const float mx05 = src[1] - src[6];
        const float mx06 = src[2] - src[5];
        const float mx07 = src[3] - src[4];
        const float mx08 = mx00 + mx03;
        const float mx09 = mx01 + mx02;
        const float mx0a = mx00 - mx03;
        const float mx0b = mx01 - mx02;
        const float mx0c = 1.38703984532215f*mx04 + 0.275899379282943f*mx07;
        const float mx0d = 1.17587560241936f*mx05 + 0.785694958387102f*mx06;
        const float mx0e = -0.785694958387102f*mx05 + 1.17587560241936f*mx06;
        const float mx0f = 0.275899379282943f*mx04 - 1.38703984532215f*mx07;
        const float mx10 = 0.353553390593274f * (mx0c - mx0d);
        const float mx11 = 0.353553390593274f * (mx0e - mx0f);
        dst[0] = 0.353553390593274f * (mx08 + mx09);
        dst[1] = 0.353553390593274f * (mx0c + mx0d);
        dst[2] = 0.461939766255643f*mx0a + 0.191341716182545f*mx0b;
        dst[3] = 0.707106781186547f * (mx10 - mx11);
        dst[4] = 0.353553390593274f * (mx08 - mx09);
        dst[5] = 0.707106781186547f * (mx10 + mx11);
        dst[6] = 0.191341716182545f*mx0a - 0.461939766255643f*mx0b;
        dst[7] = 0.353553390593274f * (mx0e + mx0f);
        dst += 4;
        src += 4;
    }
}

　　我靠，怎麼這麼複雜，比IPOL的複雜多了，代碼也多了不少，這不可能比IPOL快，但仔細的數一數，這個算法只用了24個加減法及20次乘法，而IPOL裏的則使用了56次加法及64次乘法，相比之下IPOL的計算量明顯更大。

　　上述代碼其實有錯誤，很明顯，若是執行這個for循環，則指針明顯超過了邊界，因此實際是沒有這個for的，這裏有，主要是後面的SSE優化確實須要這個for循環，做者沒有刪除他而已。

　　這樣一個算法說實在的不是很適合使用SSE指令的，不是說作不到，也能夠作到，用shuffle指令分別把單個數據分配到XMM寄存器中，而後執行操做，只是這樣作可能拔苗助長，反而速度更慢。

　　可是咱們知道這樣的算法，每一行之間的數據是互不干涉的，若是咱們考慮8行數據一塊兒算，而後先把他們轉置，這樣就能夠充分利用SSE的特性了，而且這個時候基本上就是把上述的C裏面的運算符換成對等的SIMD運算符，以下所示：

//    8*8轉置浮點數據的一維DCT變換，即對每行分別執行一維的DCT變換。
__forceinline void IM_DCT1D_8x8_T_GT_SSE(float *In, float *Out)
{
    __m128 c0353 = _mm_set1_ps(0.353553390593274f);
    __m128 c0707 = _mm_set1_ps(0.707106781186547f);
    for (int Index = 0; Index < 2; Index++)
    {
        __m128 ms0 = _mm_load_ps(In);
        __m128 ms1 = _mm_load_ps(In + 8);
        __m128 ms2 = _mm_load_ps(In + 16);
        __m128 ms3 = _mm_load_ps(In + 24);
        __m128 ms4 = _mm_load_ps(In + 32);
        __m128 ms5 = _mm_load_ps(In + 40);
        __m128 ms6 = _mm_load_ps(In + 48);
        __m128 ms7 = _mm_load_ps(In + 56);

        __m128 mx00 = _mm_add_ps(ms0, ms7);
        __m128 mx01 = _mm_add_ps(ms1, ms6);
        __m128 mx02 = _mm_add_ps(ms2, ms5);
        __m128 mx03 = _mm_add_ps(ms3, ms4);
        __m128 mx04 = _mm_sub_ps(ms0, ms7);
        __m128 mx05 = _mm_sub_ps(ms1, ms6);
        __m128 mx06 = _mm_sub_ps(ms2, ms5);
        __m128 mx07 = _mm_sub_ps(ms3, ms4);
        __m128 mx08 = _mm_add_ps(mx00, mx03);
        __m128 mx09 = _mm_add_ps(mx01, mx02);
        __m128 mx0a = _mm_sub_ps(mx00, mx03);
        __m128 mx0b = _mm_sub_ps(mx01, mx02);

        __m128 mx0c = _mm_add_ps(_mm_mul_ps(_mm_set1_ps(1.38703984532215f), mx04), _mm_mul_ps(_mm_set1_ps(0.275899379282943f), mx07));
        __m128 mx0d = _mm_add_ps(_mm_mul_ps(_mm_set1_ps(1.17587560241936f), mx05), _mm_mul_ps(_mm_set1_ps(+0.785694958387102f), mx06));
        __m128 mx0e = _mm_add_ps(_mm_mul_ps(_mm_set1_ps(-0.785694958387102f), mx05), _mm_mul_ps(_mm_set1_ps(+1.17587560241936f), mx06));
        __m128 mx0f = _mm_add_ps(_mm_mul_ps(_mm_set1_ps(0.275899379282943f), mx04), _mm_mul_ps(_mm_set1_ps(-1.38703984532215f), mx07));
        __m128 mx10 = _mm_mul_ps(c0353, _mm_sub_ps(mx0c, mx0d));
        __m128 mx11 = _mm_mul_ps(c0353, _mm_sub_ps(mx0e, mx0f));

        _mm_store_ps(Out + 0, _mm_mul_ps(c0353, _mm_add_ps(mx08, mx09)));
        _mm_store_ps(Out + 8, _mm_mul_ps(c0353, _mm_add_ps(mx0c, mx0d)));
        _mm_store_ps(Out + 16, _mm_add_ps(_mm_mul_ps(_mm_set1_ps(0.461939766255643f), mx0a), _mm_mul_ps(_mm_set1_ps(0.191341716182545f), mx0b)));
        _mm_store_ps(Out + 24, _mm_mul_ps(c0707, _mm_sub_ps(mx10, mx11)));
        _mm_store_ps(Out + 32, _mm_mul_ps(c0353, _mm_sub_ps(mx08, mx09)));
        _mm_store_ps(Out + 40, _mm_mul_ps(c0707, _mm_add_ps(mx10, mx11)));
        _mm_store_ps(Out + 48, _mm_add_ps(_mm_mul_ps(_mm_set1_ps(0.191341716182545f), mx0a), _mm_mul_ps(_mm_set1_ps(-0.461939766255643f), mx0b)));
        _mm_store_ps(Out + 56, _mm_mul_ps(c0353, _mm_add_ps(mx0e, mx0f)));
        
        In += 4;
        Out += 4;
    }
}

　　注意這裏就有個for循環了，因爲SSE一次性只能處理4個像素，一次須要2個才能處理一行。

　　注意，這裏基本上就是純C的語法直接翻譯到SIMD指令，沒有任何其餘的技巧。

　　至於DCT的逆變換，也是同樣的道理，能夠參考相關代碼。

　　注意上述計算是得到了轉置後的1D計算結果，若是須要得到真正的8x1的DCT結果，必須還要進行一次轉置，即：

__forceinline void IM_DCT1D_8x8_GT_SSE(float *In, float *Out)
{
    __declspec(align(16)) float Temp1[8 * 8];
    __declspec(align(16)) float Temp2[8 * 8];
    IM_Transpose8x8(In, Temp1);
    IM_DCT1D_8x8_T_GT_SSE(Temp1, Temp2);
    IM_Transpose8x8(Temp2, Out);
}

　　那麼這樣按照2D轉置的寫法，此算法的2D轉置過程爲：

__forceinline void IM_DCT2D_8x8_GT_SSE(float *In, float *Out)
{
    __declspec(align(16)) float Temp1[8 * 8];
    __declspec(align(16)) float Temp2[8 * 8];
    IM_DCT1D_8x8_GT_SSE(In, Temp1);
    IM_Transpose8x8(Temp1, Temp2);
    IM_DCT1D_8x8_GT_SSE(Temp2, Temp1);
    IM_Transpose8x8(Temp1, Out);
}

　　把其中的IM_DCT1D_8x8_GT_SSE展開後合併有關內存獲得：

 1 __forceinline void IM_DCT2D_8x8_GT_SSE(float *In, float *Out)
 2 {
 3     __declspec(align(16)) float Temp1[8 * 8];
 4     __declspec(align(16)) float Temp2[8 * 8];
 5 
 6     IM_Transpose8x8(In, Temp1);
 7     IM_DCT1D_8x8_T_GT_SSE(Temp1, Temp2);
 8     IM_Transpose8x8(Temp2, Temp1);
 9     
10     IM_Transpose8x8(Temp1, Temp2);
11 
12     IM_Transpose8x8(Temp2, Temp1);
13     IM_DCT1D_8x8_T_GT_SSE(Temp1, Temp2);
14     IM_Transpose8x8(Temp2, Temp1);
15 
16     IM_Transpose8x8(Temp1, Out);
17 }

　　咱們看到第10和第12行，明顯是個相互抵消的過程，徹底能夠山粗，第14行和第16也有是抵消的過程，所以，最終的2D的8x8 的DCT能夠表達以下：

__forceinline void IM_DCT2D_8x8_GT_SSE(float *In, float *Out)
{
    __declspec(align(16)) float Temp1[8 * 8];
    __declspec(align(16)) float Temp2[8 * 8];

    IM_Transpose8x8(In, Temp1);
    IM_DCT1D_8x8_T_GT_SSE(Temp1, Temp2);
    IM_Transpose8x8(Temp2, Temp1);

    IM_DCT1D_8x8_T_GT_SSE(Temp1, Out);
}

 　　一樣的，咱們能夠考慮利用列方向相鄰DCT的的重複信息，這樣能夠到的又一個版本的優化代碼：

int IM_DCT_Denoising_8x8(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, float Sigma, bool FastMode)
{
    int Channel = Stride / Width;
    if ((Src == NULL) || (Dest == NULL))                                            return IM_STATUS_NULLREFRENCE;
    if ((Width <= 7) || (Height <= 7) || (Sigma <= 0))                                return IM_STATUS_INVALIDPARAMETER;
    if ((Channel != 1) && (Channel != 3))                                            return IM_STATUS_NOTSUPPORTED;

    int Status = IM_STATUS_OK;
    if (Channel == 3)
    {

    }
    else
    {
        int Step = (FastMode == true ? 2 : 1);
        float Threshold = 3 * Sigma;
        float *DctIn = (float *)_mm_malloc(8 * Height * sizeof(float), 32);
        float *DctOut = (float *)_mm_malloc(64 * sizeof(float), 32);
        float *Sum = (float *)_mm_malloc(Width * Height * sizeof(float), 32);
        if ((DctIn == NULL) || (DctOut == NULL) || (Sum == NULL))
        {
            Status = IM_STATUS_OUTOFMEMORY;
            goto FreeMemory;
        }
        memset(Sum, 0, Width * Height * sizeof(float));
        for (int X = 0; X < Width - 7; X += Step)
        {
            for (int Y = 0; Y < Height; Y++)
            {
                IM_Convert8ucharTo8float(Src + Y * Stride + X, DctIn + Y * 8);                        //    把一整列的字節數據轉換爲浮點數
            }
            int BlockSize = 8, Block = Height / BlockSize;
            for (int Y = 0; Y < Block * BlockSize; Y += BlockSize)                                    //    利用他快速計算
            {
                IM_Transpose8x8(DctIn + Y * 8, DctOut);
                IM_DCT1D_8x8_T_GT_SSE(DctOut, DctOut);
                IM_Transpose8x8(DctOut, DctIn + Y * 8);
            }
            for (int Y = Block * BlockSize; Y < Height; Y++)
            {
                IM_DCT1D_8x1_GT_C(DctIn + Y * 8, DctIn + Y * 8);
            }
            for (int Y = 0; Y < Height - 7; Y += Step)
            {
                IM_DCT1D_8x8_T_GT_SSE(DctIn + Y * 8, DctOut);                                            //    DCT變換
                IM_DctThreshold8x8(DctOut, Threshold);                                                    //    閾值處理
                IM_IDCT2D_8x8_GT_SSE(DctOut, DctOut);                                                    //    在反變換回來
                IM_UpdateSum8x8(Sum + Y * Width + X, DctOut, Width);                                    //    更新權重和閾值
            }
        }
        IM_SumDivWeight2uchar8x8(Sum, Dest, Width, Height, Stride, FastMode);
    FreeMemory:
        if (DctIn != NULL)        _mm_free(DctIn);
        if (DctOut != NULL)        _mm_free(DctOut);
        if (Sum != NULL)        _mm_free(Sum);
        return Status;
    }
}

　　果真不出所料，執行速度大爲提高，100W圖執行時間120ms，快速模式時間38ms左右。

　　看樣子數學自己的優化比你用其餘方式來搞確實更有效果啊，真心佩服這些數學高手。

　　那麼接着看，在《Randomized Redundant DCT Efficent Denoising by Using Random Subsampling of DCT Patches》所提供的代碼後面，還提供了另一種DCT的計算方法，其參考論文的名字是 Practical fast 1-D DCT algorithms with 11 multiplications，沒看錯，11次乘法，比前面的20次乘法的GT算法又要少9次，其使用的加法數據是26次，比GT算法稍微多點。趕快實踐，測試結果代表速度速度爲：

　　100W圖執行時間105ms，快速模式時間30ms左右。

　　再次感謝萬能的數學啊。

　　最後我也嘗試了進行4X4的DCT變換，看看效果和速度如何，實踐代表，4X4的DCT去噪在一樣的Sigma參數下，效果不是很好，可是速度呢要比8X8的快速模式要快那麼一丁點。

　　在《Randomized Redundant DCT Efficent Denoising by Using Random Subsampling of DCT Patches》論文中提到的多線程並行處理，其中計算DCT和閾值確實能夠並行（這個時候就不易考慮1D的DCT的重複計算事項了），可是權重和累加值的計算不是很好弄，論文的建議是將圖像然高度方向分塊，而後作邊界擴展，我認爲這樣是可行的，但編碼稍顯繁瑣，未有去進行嘗試，對於彩色圖像，論文裏提到用Color Decorrelation，說實在的沒看懂，我就直接分解成三個獨立通道，而後每一個通道獨立計算，最後再合成了。固然這裏能夠方便的用openmp的section來進行並行化處理，惟一的缺點就是佔用內存會比串行的大3倍多，這就看應用場合了。

　   在說明一點，通過測試8*8的快速DCT去噪效果和原始的從視覺上基本看不出什麼差別，但速度大約相差3到4倍，這是由於這種基於DCT的計算存在不少的冗餘量，快速算法只取其中的1/4數據，亦可基本知足去噪需求，而基於4X4的去噪，因爲其涉及的領域範圍太小難以有效的去除噪音。所以，8x8快速算法更具備實際應用價值。

　　說道最後，嘗試了一下AVX的代碼，這個算法使用AVX編程也很簡單，特別是DCT過程對一個的AVX函數只須要把那個for循環去除，而後裏面的相關SIMD指令更改成_mm256就能夠了，Transpose8x8則用AVX實現更爲直接，閾值和累加值的計算過程修改也較爲方便，編譯時在編譯器裏的代碼生成->啓用加強指令集裏選擇：高級矢量擴展 (/arch:AVX)，編碼AVX和SSE切換時的週期耗時，結果代表速度大概還能提高一點點，但僅僅是一點點，100W會對的快速模式大概能到27ms，所以對AVX的表現仍是很失望的。

　　提供一個測試效果圖：

　　     

　   Demo下載地址：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar

      不過令我沮喪的是，聽說在手機上使用高通的芯片的Dwt去噪針對20M的圖像去噪時間只要幾毫秒，不知道是真是假，反正挺受打擊的。

      

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