數學本質上優化能比代碼層次的優化帶來更大的效益，真心感謝曆史上那些偉大的數學家，讓DCT去噪的速度能進一步提高三四倍。
在優化IPOL網站中基于DCT（離散餘弦變換）的圖像去噪算法（附源代碼） 一文中，我們曾經優化過基于DCT變換的圖像去噪算法，在那文所提供的Demo中，處理一副1000*1000左右的灰階噪音圖像耗時約450ms，如果采用所謂的快速模式耗時約150ms，說實在的，這個速度确實還是有點慢，後續曾嘗試用AVX優化，但是感覺AVX真的沒有SSE用的友善，而且AVX裡還有不少陷阱，本以為這個算法優化沒有什麼希望了，但前幾日網友推薦了一片論文《Randomized Redundant DCT Efficent Denoising by Using Random Subsampling of DCT Patches》裡提供了比較誘人的資料，于是我又對這個算法有了新的興趣，那我們來看看這篇論文有何特點。
先來看看論文提供的圖檔和資料：
其中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變量裡的四個浮點數的累加值，即如下代碼：
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代碼貼出如下：
// 普通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，他就是執行的水準方向元素相加，他執行的類似下圖的操作：
我們的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的圖像去噪時間隻要幾毫秒，不知道是真是假，反正挺受打擊的。