SSE图像算法优化系列二十一：基于DCT变换图像去噪算法的进一步优化（100W像素30ms）。

时间：2022-10-03 02:58:39

　　在优化IPOL网站中基于DCT（离散余弦变换）的图像去噪算法（附源代码） 一文中，我们曾经优化过基于DCT变换的图像去噪算法，在那文所提供的Demo中，处理一副1000*1000左右的灰度噪音图像耗时约450ms，如果采用所谓的快速模式耗时约150ms，说实在的，这个速度确实还是有点慢，后续曾尝试用AVX优化，但是感觉AVX真的没有SSE用的方便，而且AVX里还有不少陷阱，本以为这个算法优化没有什么希望了，但前几日网友推荐了一片论文《Randomized Redundant DCT Efficent Denoising by Using Random Subsampling of DCT Patches》里提供了比较诱人的数据，于是我又对这个算法有了新的兴趣，那我们来看看这篇论文有何特点。

　　先来看看论文提供的图片和数据：

SSE图像算法优化系列二十一：基于DCT变换图像去噪算法的进一步优化（100W像素30ms）。

其中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 + )));

        Out[] = Data[] + Data[] + Data[] + Data[];                            //    这里是否还有更好的方式实现呢?

　　上述代码使用了SSE指令和FPU指令共用，不是不行，而是确实效率不是很好，以前做这个算法时对SIMD指令集的了解还不是很深刻，因此，现在看来这里有很好的解决方案，为表述方便，我们把原始的8x8的DCT变换C代码贴出如下：

//    普通C语言版本的8x1的一维DCT变换

__forceinline void IM_DCT1D_8x1_C(float *In, float *Out)

{

    for (int Y = , Index = ; Y < ; Y++)

    {

        float Sum = ;

        for (int X = ; X < ; X++)

        {

            Sum += In[X] * DCTbasis[Index++];

        }

        Out[Y] = Sum;

    }

}

　　其中的DCTbasis是预先定义好的64个元素的浮点数组。

　　Out的每个元素等于8个浮点数的和，但是是水平方向元素的和，我们知道SSE里有个指令叫_mm_hadd_ps，他就是执行的水平方向元素相加，他执行的类似下图的操作：

SSE图像算法优化系列二十一：基于DCT变换图像去噪算法的进一步优化（100W像素30ms）。

　　我们的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 + );

    __m128 Sum0_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + ));

    __m128 Sum0_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + ));

    __m128 Sum1_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + ));

    __m128 Sum1_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + ));

    __m128 Sum2_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + ));

    __m128 Sum2_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + ));

    __m128 Sum3_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + ));

    __m128 Sum3_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + ));

    __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 + , _mm_hadd_ps(Sum01, Sum23));

    __m128 Sum4_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + ));

    __m128 Sum4_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + ));

    __m128 Sum5_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + ));

    __m128 Sum5_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + ));

    __m128 Sum6_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + ));

    __m128 Sum6_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + ));

    __m128 Sum7_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + ));

    __m128 Sum7_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + ));

    __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 + , _mm_hadd_ps(Sum45, Sum67));

}

　　对于DCT逆变换，也做类似的处理，100万速度一下子就能提升到350ms，快速版本约为120ms，也是一个长足的进度。

　　 在DCT变换完成后，为了进行去噪，需要一个阈值的处理，过程，普通的C语言类似下述过程：

    for (int YY = ; YY < ; YY++)

        if (abs(DctOut[YY] )< Threshold) DctOut[YY]  = ;            //    阈值处理

　　这也完全可以转换成SSE处理：

//    DCT阈值处理

__forceinline void IM_DctThreshold8x8(float *DctOut, float Threshold)

{

    __m128 Th = _mm_set_ps1(Threshold);

    __m128 DctOut0 = _mm_load_ps(DctOut + );

    __m128 DctOut1 = _mm_load_ps(DctOut + );

    _mm_store_ps(DctOut + , _mm_blend_ps(_mm_and_ps(DctOut0, _mm_cmpgt_ps(_mm_abs_ps(DctOut0), Th)), DctOut0, ));        //    // keep 00 coef.

    _mm_store_ps(DctOut + , _mm_and_ps(DctOut1, _mm_cmpgt_ps(_mm_abs_ps(DctOut1), Th)));

    __m128 DctOut2 = _mm_load_ps(DctOut + );

    __m128 DctOut3 = _mm_load_ps(DctOut + );

    _mm_store_ps(DctOut + , _mm_and_ps(DctOut2, _mm_cmpgt_ps(_mm_abs_ps(DctOut2), Th)));

    _mm_store_ps(DctOut + , _mm_and_ps(DctOut3, _mm_cmpgt_ps(_mm_abs_ps(DctOut3), Th)));

    __m128 DctOut4 = _mm_load_ps(DctOut + );

    __m128 DctOut5 = _mm_load_ps(DctOut + );

    _mm_store_ps(DctOut + , _mm_and_ps(DctOut4, _mm_cmpgt_ps(_mm_abs_ps(DctOut4), Th)));

    _mm_store_ps(DctOut + , _mm_and_ps(DctOut5, _mm_cmpgt_ps(_mm_abs_ps(DctOut5), Th)));

    __m128 DctOut6 = _mm_load_ps(DctOut + );

    __m128 DctOut7 = _mm_load_ps(DctOut + );

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

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

    __m128 DctOut8 = _mm_load_ps(DctOut + );

    __m128 DctOut9 = _mm_load_ps(DctOut + );

    _mm_store_ps(DctOut + , _mm_and_ps(DctOut8, _mm_cmpgt_ps(_mm_abs_ps(DctOut8), Th)));

    _mm_store_ps(DctOut + , _mm_and_ps(DctOut9, _mm_cmpgt_ps(_mm_abs_ps(DctOut9), Th)));

    __m128 DctOut10 = _mm_load_ps(DctOut + );

    __m128 DctOut11 = _mm_load_ps(DctOut + );

    _mm_store_ps(DctOut + , _mm_and_ps(DctOut10, _mm_cmpgt_ps(_mm_abs_ps(DctOut10), Th)));

    _mm_store_ps(DctOut + , _mm_and_ps(DctOut11, _mm_cmpgt_ps(_mm_abs_ps(DctOut11), Th)));

    __m128 DctOut12 = _mm_load_ps(DctOut + );

    __m128 DctOut13 = _mm_load_ps(DctOut + );

    _mm_store_ps(DctOut + , _mm_and_ps(DctOut12, _mm_cmpgt_ps(_mm_abs_ps(DctOut12), Th)));

    _mm_store_ps(DctOut + , _mm_and_ps(DctOut13, _mm_cmpgt_ps(_mm_abs_ps(DctOut13), Th)));

    __m128 DctOut14 = _mm_load_ps(DctOut + );

    __m128 DctOut15 = _mm_load_ps(DctOut + );

    _mm_store_ps(DctOut + , _mm_and_ps(DctOut14, _mm_cmpgt_ps(_mm_abs_ps(DctOut14), Th)));

    _mm_store_ps(DctOut + , _mm_and_ps(DctOut15, _mm_cmpgt_ps(_mm_abs_ps(DctOut15), Th)));

}

　　充分利用了_mm_cmpgt_ps比较函数的返回值和_mm_and_ps的位运算功能。

　　后续的权重累积以及DCT值的累加也一样可以用SSE优化，也是很简单的事情。

最后的DCT累计值和权重相除即可得到最终的结果值，即如下代码：

        for (int Y = , Index = ; Y < Height; Y++)

        {

            unsigned char *LinePD = Dest + Y * Stride;

            for (int X = ; 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 = , Block = Width / BlockSize;

    for (int Y = , Index = ; Y < Height; Y++)

    {

        unsigned char *LinePD = Dest + Y * Stride;

        for (int X = ; X < Block * BlockSize; X += BlockSize, Index += BlockSize)

        {

            __m128 Div1 = _mm_div_ps(_mm_loadu_ps(Sum + Index + ), _mm_loadu_ps(Weight + Index + ));        //    即使Sum或Weight是16字节对齐的，但如果Width是奇数，由于最后几个像素的作用，使得并不是每次加载数据时都是16字节对齐的内存，所以只能用loadu而不能用load

            __m128 Div2 = _mm_div_ps(_mm_loadu_ps(Sum + Index + ), _mm_loadu_ps(Weight + Index + ));

            __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 + , _mm_cvtepi32_ps(_mm_shuffle_epi8(SrcV, _mm_setr_epi8(, -, -, -, , -, -, -, , -, -, -, , -, -, -))));

    _mm_storeu_ps(Dest + , _mm_cvtepi32_ps(_mm_shuffle_epi8(SrcV, _mm_setr_epi8(, -, -, -, , -, -, -, , -, -, -, , -, -, -))));

}

　　更进一步，我们可以把他和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(, -, -, -, , -, -, -, , -, -, -, , -, -, -)));

    __m128 In2 = _mm_cvtepi32_ps(_mm_shuffle_epi8(SrcV, _mm_setr_epi8(, -, -, -, , -, -, -, , -, -, -, , -, -, -)));

    __m128 Sum0_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + ));

    __m128 Sum0_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + ));

    __m128 Sum1_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + ));

    __m128 Sum1_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + ));

    __m128 Sum2_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + ));

    __m128 Sum2_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + ));

    __m128 Sum3_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + ));

    __m128 Sum3_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + ));

    __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 + , _mm_hadd_ps(Sum01, Sum23));

    __m128 Sum4_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + ));

    __m128 Sum4_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + ));

    __m128 Sum5_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + ));

    __m128 Sum5_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + ));

    __m128 Sum6_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + ));

    __m128 Sum6_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + ));

    __m128 Sum7_0 = _mm_mul_ps(In1, _mm_load_ps(DCTbasis + ));

    __m128 Sum7_1 = _mm_mul_ps(In2, _mm_load_ps(DCTbasis + ));

    __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 + , _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 <= ) || (Height <= ) || (Sigma <= ))                                return IM_STATUS_INVALIDPARAMETER;

    if ((Channel != ) && (Channel != ))                                            return IM_STATUS_NOTSUPPORTED;

    int Status = IM_STATUS_OK;

    if (Channel == )

    {

                ////

    }

    else

    {

        int Step = (FastMode == true ?  : );

        float Threshold =  * Sigma;

        float *Dct1D = (float *)_mm_malloc( * Height * sizeof(float), );

        float *DctOut = (float *)_mm_malloc( * sizeof(float), );

        float *Sum = (float *)_mm_malloc(Width * Height * sizeof(float), );

        float *Weight = (float *)_mm_malloc(Width * Height * sizeof(float), );

        if ((Dct1D == NULL) || (DctOut == NULL) || (Sum == NULL) || (Weight == NULL))

        {

            Status = IM_STATUS_OUTOFMEMORY;

            goto FreeMemory;

        }

        memset(Sum, , Width * Height * sizeof(float));

        for (int X = ; X < Width - ; X += Step)

        {

            for (int Y = ; Y < Height; Y++)

            {

                IM_DCT1D_8x1_SSE(Src + Y * Stride + X, Dct1D + Y * );

            }

            for (int Y = ; Y < Height - ; Y += Step)

            {

                IM_DCT2D_8x8_With1DRowDCT_SSE(Dct1D + Y * , 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 = ; i < ; i++)

    {

        const float mx00 = src[] + src[];

        const float mx01 = src[] + src[];

        const float mx02 = src[] + src[];

        const float mx03 = src[] + src[];

        const float mx04 = src[] - src[];

        const float mx05 = src[] - src[];

        const float mx06 = src[] - src[];

        const float mx07 = src[] - src[];

        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.353553390593274f * (mx08 + mx09);

        dst[] = 0.353553390593274f * (mx0c + mx0d);

        dst[] = 0.461939766255643f*mx0a + 0.191341716182545f*mx0b;

        dst[] = 0.707106781186547f * (mx10 - mx11);

        dst[] = 0.353553390593274f * (mx08 - mx09);

        dst[] = 0.707106781186547f * (mx10 + mx11);

        dst[] = 0.191341716182545f*mx0a - 0.461939766255643f*mx0b;

        dst[] = 0.353553390593274f * (mx0e + mx0f);

        dst += ;

        src += ;

    }

}

　　我靠，怎么这么复杂，比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 = ; Index < ; Index++)

    {

        __m128 ms0 = _mm_load_ps(In);

        __m128 ms1 = _mm_load_ps(In + );

        __m128 ms2 = _mm_load_ps(In + );

        __m128 ms3 = _mm_load_ps(In + );

        __m128 ms4 = _mm_load_ps(In + );

        __m128 ms5 = _mm_load_ps(In + );

        __m128 ms6 = _mm_load_ps(In + );

        __m128 ms7 = _mm_load_ps(In + );

        __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 + , _mm_mul_ps(c0353, _mm_add_ps(mx08, mx09)));

        _mm_store_ps(Out + , _mm_mul_ps(c0353, _mm_add_ps(mx0c, mx0d)));

        _mm_store_ps(Out + , _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 + , _mm_mul_ps(c0707, _mm_sub_ps(mx10, mx11)));

        _mm_store_ps(Out + , _mm_mul_ps(c0353, _mm_sub_ps(mx08, mx09)));

        _mm_store_ps(Out + , _mm_mul_ps(c0707, _mm_add_ps(mx10, mx11)));

        _mm_store_ps(Out + , _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 + , _mm_mul_ps(c0353, _mm_add_ps(mx0e, mx0f)));

        In += ;

        Out += ;

    }

}

　　注意这里就有个for循环了，由于SSE一次性只能处理4个像素，一次需要2个才能处理一行。

　　注意，这里基本上就是纯C的语法直接翻译到SIMD指令，没有任何其他的技巧。

　　至于DCT的逆变换，也是一样的道理，可以参考相关代码。

　　注意上述计算是获得了转置后的1D计算结果，如果需要获得真正的8x1的DCT结果，必须还要进行一次转置，即：

__forceinline void IM_DCT1D_8x8_GT_SSE(float *In, float *Out)

{

    __declspec(align()) float Temp1[ * ];

    __declspec(align()) float Temp2[ * ];

    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()) float Temp1[ * ];

    __declspec(align()) float Temp2[ * ];

    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展开后合并有关内存得到：

 __forceinline void IM_DCT2D_8x8_GT_SSE(float *In, float *Out)

 {

     __declspec(align()) float Temp1[ * ];

     __declspec(align()) float Temp2[ * ];

     IM_Transpose8x8(In, Temp1);

     IM_DCT1D_8x8_T_GT_SSE(Temp1, Temp2);

     IM_Transpose8x8(Temp2, Temp1);

     IM_Transpose8x8(Temp1, Temp2);

     IM_Transpose8x8(Temp2, Temp1);

     IM_DCT1D_8x8_T_GT_SSE(Temp1, Temp2);

     IM_Transpose8x8(Temp2, Temp1);

     IM_Transpose8x8(Temp1, Out);

 }

　　我们看到第10和第12行，明显是个相互抵消的过程，完全可以山粗，第14行和第16也有是抵消的过程，因此，最终的2D的8x8 的DCT可以表达如下：

__forceinline void IM_DCT2D_8x8_GT_SSE(float *In, float *Out)

{

    __declspec(align()) float Temp1[ * ];

    __declspec(align()) float Temp2[ * ];

    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 <= ) || (Height <= ) || (Sigma <= ))                                return IM_STATUS_INVALIDPARAMETER;

    if ((Channel != ) && (Channel != ))                                            return IM_STATUS_NOTSUPPORTED;

    int Status = IM_STATUS_OK;

    if (Channel == )

    {

    }

    else

    {

        int Step = (FastMode == true ?  : );

        float Threshold =  * Sigma;

        float *DctIn = (float *)_mm_malloc( * Height * sizeof(float), );

        float *DctOut = (float *)_mm_malloc( * sizeof(float), );

        float *Sum = (float *)_mm_malloc(Width * Height * sizeof(float), );

        if ((DctIn == NULL) || (DctOut == NULL) || (Sum == NULL))

        {

            Status = IM_STATUS_OUTOFMEMORY;

            goto FreeMemory;

        }

        memset(Sum, , Width * Height * sizeof(float));

        for (int X = ; X < Width - ; X += Step)

        {

            for (int Y = ; Y < Height; Y++)

            {

                IM_Convert8ucharTo8float(Src + Y * Stride + X, DctIn + Y * );                        //    把一整列的字节数据转换为浮点数

            }

            int BlockSize = , Block = Height / BlockSize;

            for (int Y = ; Y < Block * BlockSize; Y += BlockSize)                                    //    利用他快速计算

            {

                IM_Transpose8x8(DctIn + Y * , DctOut);

                IM_DCT1D_8x8_T_GT_SSE(DctOut, DctOut);

                IM_Transpose8x8(DctOut, DctIn + Y * );

            }

            for (int Y = Block * BlockSize; Y < Height; Y++)

            {

                IM_DCT1D_8x1_GT_C(DctIn + Y * , DctIn + Y * );

            }

            for (int Y = ; Y < Height - ; Y += Step)

            {

                IM_DCT1D_8x8_T_GT_SSE(DctIn + Y * , 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的表现还是很失望的。

　　提供一个测试效果图：

　　SSE图像算法优化系列二十一：基于DCT变换图像去噪算法的进一步优化（100W像素30ms）。     SSE图像算法优化系列二十一：基于DCT变换图像去噪算法的进一步优化（100W像素30ms）。

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

不过令我沮丧的是，据说在手机上使用高通的芯片的Dwt去噪针对20M的图像去噪时间只要几毫秒，不知道是真是假，反正挺受打击的。

SSE图像算法优化系列二十一：基于DCT变换图像去噪算法的进一步优化（100W像素30ms）。

相关文章