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

时间:2024-01-21 13:11:17

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