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》里提供了比较诱人的数据，于是我又对这个算法有了新的兴趣，那我们来看看这篇论文有何特点。
先来看看论文提供的图片和数据：
其中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的图像去噪时间只要几毫秒，不知道是真是假，反正挺受打击的。