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

最新推荐文章于 2024-01-05 10:22:29 发布

原创

最新推荐文章于 2024-01-05 10:22:29 发布 · 344 阅读

2 ·

CC 4.0 BY-SA版权

文章标签：

#算法 #人工智能 #机器学习

本文介绍了一种基于SSE指令集优化的DCT图像去噪算法，通过对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;
        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(DctOut