在優化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的圖像去噪時間只要幾毫秒，不知道是真是假，反正挺受打擊的。
