基于DCT变换图像去噪算法的进一步优化（100W像素30ms）

最新推荐文章于 2024-01-31 13:12:26 发布
最新推荐文章于 2024-01-31 13:12:26 发布
阅读量1.4k 收藏 7
点赞数 1
分类专栏： 机器视觉_opencv
原文链接：https://www.cnblogs.com/Imageshop/p/9584024.html
版权

本篇博文来自博主Imageshop，打赏或想要查阅更多内容可以移步至Imageshop。

转载自：https://www.cnblogs.com/Imageshop/p/9584024.html   侵删

　　在优化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的图像去噪时间只要几毫秒，不知道是真是假，反正挺受打击的。

算法Lab
关注
  • 1
    点赞
  • 踩
  • 7
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
专栏目录
【图像去噪】DCT图像去噪【含GUI Matlab源码 614期】
Matlab领域的博客
12-17 523
DCT图像去噪 完整的代码，方可运行；可提供运行操作视频！适合小白！
分块DCT.rar_DCT 去噪_DCT 图像 去噪_DCT去噪_denoising_分块DCT变换
07-14
实现了一个图像分块的DCT去噪算法，用于学习参考
【图像隐写】离散余弦变换DCT结合奇异值分解SVD灰色图像信息隐写（数字水印嵌入攻击提取 含NC PSNR）【含Matlab源码 3892期】
最新发布
Matlab领域的博客
01-31 767
离散余弦变换DCT结合奇异值分解SVD灰色图像信息隐写（数字水印嵌入攻击提取 含NC PSNR） 完整的代码，方可运行；可提供运行操作视频！适合小白！
MATLAB分析图像的离散余弦变换（DCT）
Chen Yuanshen的专栏
08-20 7万+
1. MATLAB的介绍以及所需函数的说明：  1.1 MATLAB  MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室）。是由美国mathworks 公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程
图像去噪及Matlab实现
热门推荐
Joey's Blog
01-06 10万+
图像去噪常用方法 图像去噪处理方法可分为空间域法和变换域法两大类。 基于离散余弦变换的图像去噪 一般而言，我们认为图像的噪声在离散余弦变换结果中处在其高频部分，而高频部分的幅值一般很小，利用这一性质，就可以实现去噪。然而，同时会失去图像的部分细节。 %读取图像 X=imread('wangshi.jpg'); X=rgb2gray(X); %读取图像尺寸 [m,n]=size(X
优化IPOL网站中基于DCT（离散余弦变换）的图像去噪算法（附源代码）
weixin_39450742的博客
07-01 460
本篇博文来自博主Imageshop，打赏或想要查阅更多内容可以移步至Imageshop。 转载自：https://www.cnblogs.com/Imageshop/p/4965192.html侵删 在您阅读本文前，先需要告诉你的是：即使是本文优化过的算法，DCT去噪的计算量依旧很大，请不要向这个算法提出实时运行的苛刻要求。 　　言归正传，在IPOL网站中有一篇基于DCT的图像去噪文章，具体的链接地址是：http://www.ipol.im/pub/art/2011/ys-dct/，IPOL...
DCT 图像去除噪声
01-18
linux下已经跑通，图像去噪声，参考readme.txt 进行编译。
【图像去噪】基于matlab GUI DCT图像去噪【含Matlab源码 614期】
订阅付费专栏（初级版），可赠送该付费专栏任意代码1份；
03-23 2259
DCT图像去噪 完整代码，直接运行，适合小白！可提供运行操作视频！
matlab DCT水印算法的四种改进
dzcera的博客
08-13 1571
DCT水印算法的四种改进
dct椒盐噪声去噪效果不好_自编码降噪效果到底怎么样？我用对比实验告诉你
weixin_29961459的博客
12-30 187
当你变强时，世界才会对你温柔以待。机器学习是目前人工智能最火的领域之一，信号的降噪是机器学习擅长的一个领域。今天讲的自编码就是机器学习的一个新研究方向，具体见前文《地震无监督自编码降噪怎样做出效果，我想清楚了几个关键问题》01 自编码是有监督学习首先要更正一个错误，自编码是监督学习，不是无监督学习。判断是否监督学习，重要的判断就是看是否使用标签。自编码学习不需要人工标注标签，它的标签就是从输入数据...
基于DCT 及高频屏蔽窗的小波图像去噪方法（清晰版）
08-26
摘要：为精确估计噪声方差，提出一种基于离散余弦变换及高频屏蔽窗的去噪方法，利用离散余弦变换（DCT）对 小波高频子带进行局部特征提取，在此基础上对由噪声引起的小波系数的幅值进行估计，从原小波系数中屏蔽噪声 部分的幅值，并用小波逆变换对剔除噪声分量后的小波系数进行恢复，从而得到去噪后的图像.该方法避开了噪声 方差的估计,并在去噪效果和保留细节方面要好于传统的阈值去噪方法.
图像去噪（matlab）.rar_DCT降噪matlab_PCA图像去噪_previous12j_图像 去噪_图像去噪
07-14
使用中值，均值，小波，DCT，PCA五种方法实现对图像的去噪处理。
DCT.zip_DCT去噪
09-23
离散余弦变换的图像去噪的编码及图像的保存和显示
图像滤波算法学习图像降噪、bilateral filtering、nonlocal means filtering、小波降噪、DCT降噪。
03-10
学习图像降噪、bilateral filtering、nonlocal means filtering、小波降噪、DCT降噪。 报告要求：对给出的图像进行上述处理，附上代码并且对处理结果给出展示。 报告形式，里面有附带的程序以及结果。
图像降噪、bilateral filtering、nonlocal means filtering、小波降噪、DCT降噪
03-10
学习图像降噪、bilateral filtering、nonlocal means filtering、小波降噪、DCT降噪。 报告要求：对给出的图像进行上述处理，附上代码并且对处理结果给出展示。 为程序
基于图像分块DCT变换的压缩重构算法matlab仿真+仿真操作录像
10-14
1.版本：matlab2021a，我录制了仿真操作录像，可以跟着操作出仿真结果 2.领域：压缩重构 3.内容：基于图像分块DCT变换的压缩重构算法matlab仿真 4.适合人群：本，硕等教研学习使用
基于DCT变换的数字水印算法
11-21
DCT离散余弦变换实现数字水印算法；在进行水印嵌入时，先对水印用了Arnold变换；包含了几种不同的水印攻击
【图像去噪】基于 DCT变换实现图像去噪matlab源码含 GUI.md
08-09
【图像去噪】基于 DCT变换实现图像去噪matlab源码含 GUI.md
基于DCT变换的图像压缩解压缩算法matlab仿真+仿真录像
10-14
1.版本：matlab2021a，我录制了仿真操作录像，可以跟着操作出仿真结果 2.领域：图像压缩解压缩 3.内容：基于DCT变换的图像压缩解压缩算法 4.适合人群：本，硕等教研学习使用
基于matlab的DCT变换细节增强算法
06-10
DCT（离散余弦变换）是一种常用的信号处理方法。基于DCT的细节增强算法可以通过提高图像的高频部分来增强图像的细节。以下是一个基于MATLAB的DCT变换细节增强算法的示例： 1. 读入图像并将其转换为灰度图像。 2. 对图像进行DCT变换。 ```matlab img_dct = dct2(img_gray); ``` 3. 对变换后的图像进行正则化。 ```matlab img_dct_norm = abs(img_dct) / max(abs(img_dct(:))); ``` 4. 设置一个阈值，将小于阈值的值置为0。 ```matlab threshold = 0.1; img_dct_norm(img_dct_norm < threshold) = 0; ``` 5. 对变换后的图像进行逆DCT变换，得到增强后的图像。 ```matlab img_enhanced = idct2(img_dct_norm); ``` 6. 显示原始图像和增强后的图像。 ```matlab figure; subplot(1, 2, 1), imshow(img_gray), title('Original Image'); subplot(1, 2, 2), imshow(img_enhanced), title('Enhanced Image'); ``` 该算法的基本思路是利用DCT变换将图像转换到频域，然后通过对频域信息进行处理来增强图像的细节，最后再进行逆变换得到增强后的图像。调整阈值可以控制增强程度和噪声的影响。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值