相關文章：
優化IPOL網站中基於DCT（離散餘弦變換）的圖像去噪算法（附源代碼）。
SSE圖像算法優化系列二十一：基於DCT變換圖像去噪算法的進一步優化（100W像素30ms）。
這個算法2015年優化過一版，2018年又優化過一版，2016年初又來回訪一下，感覺還有優化空間，繼續又折騰了將近一週，速度又有進一步的提升，下面是這個三個版本在同一台電腦上的速度比較：
對於小圖，其實2026版和2018年速度差異不是很大，到了大圖，新版約有25%左右的提升。
那麼2026版優化的渠道主要有兩處，第一個還是內部DCT變換本身的優化，以8*8的DCT為例，2015版DCT1D的代碼如下所示：
inline void IM_DCT1D_8x1_GT_C(float *In, float *Out) { const float mx00 = In[0] + In[7]; const float mx01 = In[1] + In[6]; const float mx02 = In[2] + In[5]; const float mx03 = In[3] + In[4]; const float mx04 = In[0] - In[7]; const float mx05 = In[1] - In[6]; const float mx06 = In[2] - In[5]; const float mx07 = In[3] - In[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); Out[0] = 0.353553390593274f * (mx08 + mx09); Out[1] = 0.353553390593274f * (mx0c + mx0d); Out[2] = 0.461939766255643f * mx0a + 0.191341716182545f * mx0b; Out[3] = 0.707106781186547f * (mx10 - mx11); Out[4] = 0.353553390593274f * (mx08 - mx09); Out[5] = 0.707106781186547f * (mx10 + mx11); Out[6] = 0.191341716182545f * mx0a - 0.461939766255643f * mx0b; Out[7] = 0.353553390593274f * (mx0e + mx0f); }
可以看到Out[0]到Out[7]這個8個數據前面都還有很多係數，共有26次加法及22次乘法，實際上可以把他們和前面的係數融合到一起，修改後數據如下所示：
inline void IM_DCT1D_1x8_C(float* In, float* Out) { float V00 = In[0] + In[7]; float V01 = In[0] - In[7]; float V02 = In[1] + In[6]; float V03 = In[1] - In[6]; float V04 = In[2] + In[5]; float V05 = In[2] - In[5]; float V06 = In[3] + In[4]; float V07 = In[3] - In[4]; float V08 = 0.353553390593274f * (V00 + V06); float V09 = 0.353553390593274f * (V02 + V04); float V10 = V00 - V06; float V11 = V02 - V04; float V12 = 0.490392640201616f * V01 + 0.097545161008064f * V07; float V13 = 0.415734806151273f * V03 + 0.277785116509801f * V05; float V14 = 0.415734806151273f * V05 - 0.277785116509801f * V03; float V15 = 0.097545161008064f * V01 - 0.490392640201616f * V07; float V16 = 0.707106781186547f * (V12 - V13); float V17 = 0.707106781186547f * (V14 - V15); Out[0] = V08 + V09; Out[1] = V12 + V13; Out[2] = 0.461939766255643f * V10 + 0.191341716182545f * V11; Out[3] = V16 - V17; Out[4] = V08 - V09; Out[5] = V16 + V17; Out[6] = 0.191341716182545f * V10 - 0.461939766255643f * V11; Out[7] = V14 + V15; }
修改後的代碼只有26次加法及16次乘法，減少了6次乘法。
加速的另外一個部分就是，原先的方法是沿着列方向更新權重和累加值，新的算法更改為行方向更新權重和累加值，這裏自然的就變為了緩存友好的方式了。當然其中，會增加一部分內存佔用。
原先的處理方式如下：
for (int X = 0; X < Width - 7; X += Step) { for (int Y = 0; Y < Height; Y++) { IM_Convert8ucharTo8float_PureC(Src + Y * Stride + X, DctIn + Y * 8); // 把一整列的字節數據轉換為浮點數 } for (int Y = 0; Y < Height; Y++) { IM_DCT1D_8x1_GT_C(DctIn + Y * 8, DctIn + Y * 8); } for (int Y = 0; Y < Height - 7; Y += Step) { IM_DCT2D_8x8_With1DRowDCT_GS_PureC(DctIn + Y * 8, DctOut); // DCT變換 IM_DctThreshold8x8_PureC(DctOut, Threshold); // 閾值處理 IM_IDCT2D_8x8_GS_PureC(DctOut, DctOut); // 在反變換回來 IM_UpdateSum8x8_PureC(Sum + Y * Width + X, DctOut, Width); // 更新權重和閾值 } } IM_SumDivWeight2uchar8x8_PureC(Sum, Dest, Width, Height, Stride, FastMode);
更改後的方案為：
// 先把所有的字節數據轉換為浮點數 for (int Y = 0; Y < Height; Y++) { IM_ConvertU8To32F_PureC(Src + Y * Stride, SrcF + Y * Width, Width); } // 因為後面進行了以y為起點，連續8個元素的行方向的DCT變換，所以最後7個點也是有取樣的 for (int Y = 0; Y < Height - 7; Y += Step) { float* SumL = Sum + Y * Width; // 8行數據的一維列DCT變換，保存到DctIn中 IM_DCT1D_8xW_C_PureC(SrcF + Y * Width, DctIn, Width); // 把列DCT變換的數據轉置下 IM_TransposeF_PureC(DctIn, DctInT, Width, 8); // 進行了以X為起點，X方向連續7個列方向的DCT列變換，所以列方向最後7個元素也能取樣到 for (int X = 0; X < Width - 7; X += Step) { // 可以重複利用行方向的轉置數據 IM_DCT1D_8x8_C_PureC(DctInT + X * 8, DctOut); // 轉換後的數據還要轉置 IM_Transpose8x8F_W8_PureC(DctOut, DctOut); // 閾值處理 IM_DctThreshold8x8_PureC(DctOut, Threshold); // 在反變換回來 IM_IDCT2D_8x8_PureC(DctOut, DctOut); // 更新權重和閾值 IM_UpdateSum8x8_PureC(SumL + X, DctOut, Width); } } IM_SumDivWeight2uchar8x8_PureC(Sum, Dest, Width, Height, Stride, FastMode);
具體的細節還需要大家自己慢慢悟。
對於16*16的DCT，本身的計算量就特別誇張，原先論文配套的代碼也有提供，但是同樣的道理那裏的代碼存在大量的可節省計算，很多乘法部分可以合併的。優化後的DCT16*16共需要46次乘法72次加法，16*16的去噪程度更強，同樣的Sigma參數結果會更加模糊，為了提高16*16的速度，也可以仿照8*8的方式，加大取樣間隔，實際測試間隔調整為2個像素，對結果基本沒影響，如果把間隔調整為4個像素，在加大的Sigma時，可以較為明顯的看到局部的方格，這個應該是不能忍受的， 實際測試如果把間距調整為3個像素，其實結果和原始的還是能承受的，而且加速也很客觀。
論文裏提到的棋格取樣，雖然有一定的理論正確性，但是實際不太可操作，首先是獲取滿足網格佈置的座標本身就是一個比較慢的過程，第二個即使獲取了網格座標，因為這些座標分佈基本不是按照先X方向增加，再Y方向增加的佈局方式佈置的，所以為了利用後續的加速技巧，還要先對他們進行排序，這個增加的耗時也是非常客觀的，所以我個人覺得那個只具有理論意義。
另外一個就是關於SSE和AVX加速的部分，這裏AVX起到的加速作用不是特別明顯，核心原因估計還是內存加載和保存佔比在整個過程中比較多，而且在計算部分，也曾嘗試用FMA相關指令替代mul+add，也不見有明顯的收益。 另外，對於算法中經常使用的8*8轉置，AVX版本如果學SSE那種直接處理，加速就更有限了。但是AVX裏有不同的port，如果充分利用這個則還不錯。
多線程：
因為這個DCT算法其實也是領域算法，所以不太好直接使用多線程進行處理，但是還是有辦法，對於灰度圖像，一種方法是按照高度或者寬度方向分塊，但是分塊的時候需要注意塊和塊之間必須有重疊，以保證邊緣的不分被充分處理，雖然重疊會增加那麼一丟丟的耗時，但是這不算啥，而對於彩色圖像，因為要把彩色圖分單通道後再調用單通道的處理算法，所以一個自然的想法就是分解後的單通道之間開三個線程並行就可以了。一個簡單的代碼如下：
int IM_DCT_Denoising_8x8_MultiThread_SSE(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 (Width < 256) { return IM_DCT_Denoising_8x8_SSE(Src, Dest, Width, Height, Stride, Sigma, FastMode); } else { if (Channel == 3) { int StatusA = IM_STATUS_OK, StatusB = IM_STATUS_OK, StatusC = IM_STATUS_OK; unsigned char* Blue = (unsigned char*)malloc(Width * Height * sizeof(unsigned char)); unsigned char* Green = (unsigned char*)malloc(Width * Height * sizeof(unsigned char)); unsigned char* Red = (unsigned char*)malloc(Width * Height * sizeof(unsigned char)); if ((Blue == NULL) || (Green == NULL) || (Red == NULL)) { Status = IM_STATUS_OUTOFMEMORY; goto FreeResource; } IM_SplitRGB(Src, Blue, Green, Red, Width, Height, Stride, 0); #pragma omp parallel sections num_threads(3) { #pragma omp section { StatusA = IM_DCT_Denoising_8x8_SSE(Blue, Blue, Width, Height, Width, Sigma, FastMode); if (StatusA != IM_STATUS_OK) Status = StatusA; } #pragma omp section { StatusB = IM_DCT_Denoising_8x8_SSE(Green, Green, Width, Height, Width, Sigma, FastMode); if (StatusB != IM_STATUS_OK) Status = StatusB; } #pragma omp section { StatusC = IM_DCT_Denoising_8x8_SSE(Red, Red, Width, Height, Width, Sigma, FastMode); if (StatusC != IM_STATUS_OK) Status = StatusC; } } if (Status != IM_STATUS_OK) goto FreeResource; IM_CombineRGB(Blue, Green, Red, Dest, Width, Height, Stride, 0); FreeResource: if (Blue != NULL) free(Blue); if (Green != NULL) free(Green); if (Red != NULL) free(Red); return Status; } else if (Channel == 1) { // 大圖分三個線程處理，因為在10年的電腦上，可能很多還是4核，要留一個核給操作系統比較好 int BlockSize = Width / 3; int BlockA_W = BlockSize + 7; // 1/3寬度在加上右側多7個像素，能完整的計算出左側BlockSize大小的信息 int BlockB_W = BlockSize + 14; // 1/3寬度在加上左右兩側各7個像素，能完整的計算出中間的BlockSize大小的信息 int BlockC_W = Width - BlockSize - BlockSize + 7; // 1/3寬度在加上左側多7個像素，能完整的計算出右側BlockSize大小的信息（注意寬度不一定能整除，所以要用總寬度減去2倍的1/3大小 int StatusA = IM_STATUS_OK, StatusB = IM_STATUS_OK, StatusC = IM_STATUS_OK; unsigned char* BlockA = (unsigned char*)malloc(BlockA_W * Height * sizeof(unsigned char)); unsigned char* BlockB = (unsigned char*)malloc(BlockB_W * Height * sizeof(unsigned char)); unsigned char* BlockC = (unsigned char*)malloc(BlockC_W * Height * sizeof(unsigned char)); if ((BlockA == NULL) || (BlockB == NULL) || (BlockC == NULL)) { Status = IM_STATUS_OUTOFMEMORY; goto FreeMemory; } #pragma omp parallel sections num_threads(3) { #pragma omp section { // 拷貝左側數據到臨時內存 for (int Y = 0; Y < Height; Y++) { memcpy(BlockA + Y * BlockA_W, Src + Y * Stride, BlockA_W * sizeof(unsigned char)); } StatusA = IM_DCT_Denoising_8x8_SSE(BlockA, BlockA, BlockA_W, Height, BlockA_W, Sigma, FastMode); if (StatusA != IM_STATUS_OK) { Status = StatusA; } else { // 把數據複製回去 for (int Y = 0; Y < Height; Y++) { memcpy(Dest + Y * Stride, BlockA + Y * BlockA_W, BlockSize * sizeof(unsigned char)); } } } #pragma omp section { for (int Y = 0; Y < Height; Y++) { memcpy(BlockB + Y * BlockB_W, Src + Y * Stride + BlockSize - 7, BlockB_W * sizeof(unsigned char)); } StatusB = IM_DCT_Denoising_8x8_SSE(BlockB, BlockB, BlockB_W, Height, BlockB_W, Sigma, FastMode); if (StatusB != IM_STATUS_OK) { Status = StatusB; } else { for (int Y = 0; Y < Height; Y++) { memcpy(Dest + Y * Stride + BlockSize, BlockB + Y * BlockB_W + 7, BlockSize * sizeof(unsigned char)); } } } #pragma omp section { for (int Y = 0; Y < Height; Y++) { memcpy(BlockC + Y * BlockC_W, Src + Y * Stride + 2 * BlockSize - 7, BlockC_W * sizeof(unsigned char)); } StatusC = IM_DCT_Denoising_8x8_SSE(BlockC, BlockC, BlockC_W, Height, BlockC_W, Sigma, FastMode); if (StatusC != IM_STATUS_OK) { Status = StatusC; } else { for (int Y = 0; Y < Height; Y++) { memcpy(Dest + Y * Stride + 2 * BlockSize, BlockC + Y * BlockC_W + 7, (Width - 2 * BlockSize) * sizeof(unsigned char)); } } } } FreeMemory: if (BlockA != NULL) free(BlockA); if (BlockB != NULL) free(BlockB); if (BlockC != NULL) free(BlockC); return Status; } } return Status; }
以上代碼灰度圖也是使用了3個線程，注意在中間塊部分呢，兩邊都要進行適當的擴展，否則會形成一條比較明顯的分界線，多線程加速也不是線性，開三個線程也是無法得到3倍加速的，大約為2倍左右的速度提升。
就這麼個算法，我把所有的可能旋向都做了實現，輕輕鬆鬆就弄了進7000行代碼，也是會折騰。
現在年齡大了，感覺有的時候確實是折騰不動了，新的技術在不斷髮展，而我們依舊趴在老舊的技術上啃食。 不過呢也無所謂，20年前的CS1.6，佔地1942，星際爭霸現在玩起來還是有那種感覺。
該算法在最新的Demo中已經更新，相關界面如上圖所示，下載地址：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar
