优化IPOL网站中基于DCT（离散余弦变换）的图像去噪算法（附源代码）。

　　在您阅读本文前，先需要告诉你的是：即使是本文优化过的算法，DCT去噪的计算量依旧很大，请不要向这个算法提出实时运行的苛刻要求。

　　言归正传，在IPOL网站中有一篇基于DCT的图像去噪文章，具体的链接地址是：http://www.ipol.im/pub/art/2011/ys-dct/，IPOL网站的最大特点就是他的文章全部提供源代码，而且可以基于网页运行相关算法，得到结果。不过其里面的代码本身是重实现论文的过程，基本不考虑速度的优化，因此，都相当的慢。

这篇文章的原理也是非常简单的，整个过程就是进行dct变换，然后在dct域进行硬阈值操作，再反变换回来。但是DCT变换不是针对整幅图进行处理，而是基于每个像素点的领域（这里使用的8领域或者16领域），每次则移动一个像素。IPOL上提供的代码函数也很少，但是一个关键的问题就是内存占用特别夸张，我们贴他的部分代码：

// Transfer an image im of size width x height x channel to sliding patches of
// size width_p x height_p xchannel.
// The patches are stored in patches, where each ROW is a patch after being
// reshaped to a vector.
void Image2Patches(vector<float>& im,
                   vector< vector< vector< vector< float > > > >& patches,
                   int width, int height, int channel, int width_p,
                   int height_p)
{
    int size1 = width * height;

    int counter_patch = 0;

    // Loop over the patch positions
    for (int j = 0; j < height - height_p + 1; j ++)
        for (int i = 0; i < width - width_p + 1; i ++) {
            int counter_pixel = 0;
            // loop over the pixels in the patch
            for (int kp = 0; kp < channel; kp++)
                for (int jp = 0; jp < height_p; jp ++)
                    for (int ip = 0; ip < width_p; ip ++) {
                        patches[counter_patch][kp][jp][ip] =
                                         im[kp*size1 + (j+jp)*width + i + ip];
                        counter_pixel ++;
                    }
            counter_patch ++;
        }
}

　　看到这里的patches了，他的作用是保存每个点周边的8*8领域的DCT变换的结果的，即使使用float类型的变量，也需要约Width * Height * 8 * 8 * sizeof(float)个字节的数组，假定宽度和高度都为1024的灰度图，这个数据量为256MB，其实256MB的内存在现在机器来说毫无压力，可这里要求是连续分布内存，则很有可能分配失败，在外部表现的错误则是内存溢出。我们首先来解决这个问题。

　　我们来看看原作者的代码中patches的主要作用，见下面这部分代码：

    // Loop over the patch positions
    for (int j = 0; j < height - height_p + 1; j ++)
        for (int i = 0; i < width - width_p + 1; i ++) {
            int counter_pixel = 0;
            // loop over the pixels in the patch
            for (int kp = 0; kp < channel; kp++)
                for (int jp = 0; jp < height_p; jp ++)
                    for (int ip = 0; ip < width_p; ip ++) {
                        im[kp*size1 + (j+jp)*width + i + ip] +=
                                       patches[counter_patch][kp][jp][ip];
                        im_weight[kp*size1 + (j+jp)*width + i + ip] ++;
                        counter_pixel ++;
                    }
            counter_patch ++;
        }

    // Normalize by the weight
    for (int i = 0; i < size; i ++)
        im[i] = im[i] / im_weight[i];

　　可见patches主要是为了保存每个点领域的DCT变换的数据，然后累加，上述四重循环外围两个是图像的宽度和高度方向，内部两重则是DCT的变换数据的行列方向，如果我们把DCT的行列方向的循环提到最外层，而把图像的宽度和高度方向的循环放到内存，其一就是整个过程只需要一个8*8*sizeof(float)大小的重复利用的内存，第二就是正好符号了内部放大循环，外部放小循环的优化准则，在速度和内存占用上都有重大提高。

我们来继续优化，在获取每个点的领域时，传统的方式需要8*8个循环，那整个图像就需要Width * Height * 8 * 8次了， 这个数据量太可观了，在图像处理中任意核卷积(matlab中conv2函数)的快速实现 一文中共享的代码里提到了一种方式，大约只需要Width * Height * 8次读取，并且其cache命中率还要高很多，具体的方式可参考本文附带的代码。

继续可以优化的地方就是8*8点的浮点DCT变换了。原文提供一维DCT变换的代码如下：

for (int j = 0; j < 8; j ++)
{
    out[j] = 0;
    for (int i = 0; i < 8; i ++)
    {
        out[j] += in[i] * DCTbasis[j][i];
    }
}

就是两个循环，共64次乘法和64次加法，上面的DCTbasis为固定的DCT系数矩阵。

　　而一维的IDCT的变换代码为：

for (int j = 0; j < PATCHSIZE; j ++)
{
    out[j] = 0;
    for (int i = 0; i < PATCHSIZE; i ++)
    {
        out[j] += in[i] * DCTbasis[i][j];
    }
}

和DCT唯一的区别仅仅是DCTbasis的行列坐标相反。

这种代码一看就想到了有SSE进行优化，PATCHSIZE为8 正好是两个SSE浮点数m128的大小，乘法和加法都有对应的SSE函数，一次性可进行4个浮点加法和浮点乘法，效率当然会高很多，优化后的代码如下所示：

/// <summary>
/// 8*8的一维DCT变换及其逆变换。
/// </summary>
/// <param name="In">输入的数据。</param>
/// <param name="Out">输出的数据。</param>
/// <param name="Invert">是否进行逆变换。</param>
/// <remarks> 1、输入和输出不能相同，即不支持in-place操作。</remarks>
/// <remarks> 2、算法只直接翻译IPOL上的，利用了SSE加速。</remarks>
/// <remarks> 3、在JPEG压缩等程序中的8*8DCT变换里优化的算法乘法比较少，但不好利用SSE优化，我用那里的代码测试还比下面的慢。</remarks>
/// <remarks> 4、有关8*8 DCT优化的代码见：http://insurgent.blog.hexun.com/1797358_d.html。 </remarks>
void DCT8X81D(float *In, float *Out, bool Invert)
{
    __m128 Sum, A, B;
    const float *Data = (const float *)&Sum;
    A = _mm_load_ps(In);
    B = _mm_load_ps(In + 4);

    if (Invert == FALSE)
    {
        /*for (int Y = 0; Y < PATCHSIZE; Y ++)
        {
            Out[Y] = 0;
            for (int X = 0; X < PATCHSIZE; X ++)
            {
                Out[Y] += In[X] * DCTbasis[Y * PATCHSIZE + X];
            }
        }*/

        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];                            //    这里是否还有更好的方式实现呢?

        Sum = _mm_mul_ps(A, _mm_load_ps(DCTbasis + 8));
        Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(DCTbasis + 12)));        //    不用一个Sum变量，用多个是不是还能提高指令级别的并行啊
        Out[1] = Data[0] + Data[1] + Data[2] + Data[3];

        Sum = _mm_mul_ps(A, _mm_load_ps(DCTbasis + 16));
        Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(DCTbasis + 20)));
        Out[2] = Data[0] + Data[1] + Data[2] + Data[3];

        Sum = _mm_mul_ps(A, _mm_load_ps(DCTbasis + 24));
        Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(DCTbasis + 28)));
        Out[3] = Data[0] + Data[1] + Data[2] + Data[3];

        Sum = _mm_mul_ps(A, _mm_load_ps(DCTbasis + 32));
        Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(DCTbasis + 36)));
        Out[4] = Data[0] + Data[1] + Data[2] + Data[3];

        Sum = _mm_mul_ps(A, _mm_load_ps(DCTbasis + 40));
        Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(DCTbasis + 44)));
        Out[5] = Data[0] + Data[1] + Data[2] + Data[3];

        Sum = _mm_mul_ps(A, _mm_load_ps(DCTbasis + 48));
        Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(DCTbasis + 52)));
        Out[6] = Data[0] + Data[1] + Data[2] + Data[3];

        Sum = _mm_mul_ps(A, _mm_load_ps(DCTbasis + 56));
        Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(DCTbasis + 60)));
        Out[7] = Data[0] + Data[1] + Data[2] + Data[3];
    }
    else
    {
        /*for (int Y = 0; Y < PATCHSIZE; Y ++)
        {
            Out[Y] = 0;
            for (int X = 0; X < PATCHSIZE; X ++)
            {
                Out[Y] += In[X] * IDCTbasis[Y * PATCHSIZE + X];
            }
        }*/

        Sum = _mm_mul_ps(A, _mm_load_ps(IDCTbasis));
        Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(IDCTbasis + 4)));
        Out[0] = Data[0] + Data[1] + Data[2] + Data[3];

        Sum = _mm_mul_ps(A, _mm_load_ps(IDCTbasis + 8));
        Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(IDCTbasis + 12)));
        Out[1] = Data[0] + Data[1] + Data[2] + Data[3];

        Sum = _mm_mul_ps(A, _mm_load_ps(IDCTbasis + 16));
        Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(IDCTbasis + 20)));
        Out[2] = Data[0] + Data[1] + Data[2] + Data[3];

        Sum = _mm_mul_ps(A, _mm_load_ps(IDCTbasis + 24));
        Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(IDCTbasis + 28)));
        Out[3] = Data[0] + Data[1] + Data[2] + Data[3];

        Sum = _mm_mul_ps(A, _mm_load_ps(IDCTbasis + 32));
        Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(IDCTbasis + 36)));
        Out[4] = Data[0] + Data[1] + Data[2] + Data[3];

        Sum = _mm_mul_ps(A, _mm_load_ps(IDCTbasis + 40));
        Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(IDCTbasis + 44)));
        Out[5] = Data[0] + Data[1] + Data[2] + Data[3];

        Sum = _mm_mul_ps(A, _mm_load_ps(IDCTbasis + 48));
        Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(IDCTbasis + 52)));
        Out[6] = Data[0] + Data[1] + Data[2] + Data[3];

        Sum = _mm_mul_ps(A, _mm_load_ps(IDCTbasis + 56));
        Sum = _mm_add_ps(Sum, _mm_mul_ps(B, _mm_load_ps(IDCTbasis + 60)));
        Out[7] = Data[0] + Data[1] + Data[2] + Data[3];
    }
}

　 当然，简单的循环并不是效率最高的算法，在标准的JPG压缩中就有着8*8的DCT变换，哪里的计算量有着更少的乘法和加法，在 http://insurgent.blog.hexun.com/1797358_d.html 中提出代码里，只有32次乘法和更少的加法，但是由于这个代码的向量性很差，是很难用SSE进行优化的，我实测的结果时他的代码比我用SSE优化后的速度慢。

当进行2维的DCT的时候，其步骤为对每行先进行行方向的一维DCT，然后对结果转置，在对转置后的数据进行行方向一维DCT,得到的结果再次转置则为2维DCT。8*8的转置虽然直接实现基本不存在cache miss的问题，不过还是用有关的SSE来实现未尝不是个好主意，在intrinsic中提供了一个4*4浮点转置的宏_MM_TRANSPOSE4_PS，我们对这个宏稍作修改，修改后的代码如下：

//    http://stackoverflow.com/questions/5200338/a-cache-efficient-matrix-transpose-program
//    http://stackoverflow.com/questions/16737298/what-is-the-fastest-way-to-transpose-a-matrix-In-c
//    https://msdn.microsoft.com/en-us/library/4d3eabky(v=vs.90).aspx
//    高效的矩阵转置算法，速度约为直接编写的4倍, 整理时间 2015.10.12
inline void TransposeSSE4x4(float *Src, float *Dest)
{
    __m128 Row0 = _mm_load_ps(Src);
    __m128 Row1 = _mm_load_ps(Src + 8);
    __m128 Row2 = _mm_load_ps(Src + 16);
    __m128 Row3 = _mm_load_ps(Src + 24);

    //        _MM_TRANSPOSE4_PS(Row0, Row1, Row2, Row3);                            //    或者使用这个SSE的宏

    __m128 Temp0    = _mm_shuffle_ps(Row0, Row1, _MM_SHUFFLE(1, 0, 1, 0));      //    Row0[0] Row0[1] Row1[0] Row1[1]
    __m128 Temp2    = _mm_shuffle_ps(Row0, Row1, _MM_SHUFFLE(3, 2, 3, 2));      //    Row0[2] Row0[3] Row1[2] Row1[3]
    __m128 Temp1    = _mm_shuffle_ps(Row2, Row3, _MM_SHUFFLE(1, 0, 1, 0));      //    Row2[0] Row2[1] Row3[0] Row3[1]
    __m128 Temp3    = _mm_shuffle_ps(Row2, Row3, _MM_SHUFFLE(3, 2, 3, 2));        //    Row2[2] Row2[3] Row3[2] Row3[3]          

    Row0 = _mm_shuffle_ps(Temp0, Temp1, _MM_SHUFFLE(2, 0, 2, 0));                //    Row0[0] Row1[0] Row2[0] Row3[0]
    Row1 = _mm_shuffle_ps(Temp0, Temp1, _MM_SHUFFLE(3, 1, 3, 1));                //    Row0[1] Row1[1] Row2[1] Row3[1]
    Row2 = _mm_shuffle_ps(Temp2, Temp3, _MM_SHUFFLE(2, 0, 2, 0));                //    Row0[2] Row1[2] Row2[2] Row3[2]
    Row3 = _mm_shuffle_ps(Temp2, Temp3, _MM_SHUFFLE(3, 1, 3, 1));                //    Row0[3] Row1[3] Row2[3] Row3[3]             

    _mm_store_ps(Dest, Row0);
    _mm_store_ps(Dest + 8, Row1);
    _mm_store_ps(Dest + 16, Row2);
    _mm_store_ps(Dest + 24, Row3);
}

本质上说，这些优化都是小打小闹，提高不了多少速度，当然还有一些可以优化的地方，比如权重和累加和的更新，最后的累加和和权重的相除得到结果等等都有有关的SSE函数可以使用。

还有一个可以优化的地方就是，在高度方向上前后两个像素8*8领域 在进行2D的DCT变换时，其第一次行方向上的DCT变换有7行的结果是可以重复利用的，如果利用这一点，则可以获得约15%的速度提示。

　   综合以上各种优化手段，在I5的机器上一幅512*512 的灰度图像大约处理用时为100ms左右 ，比起IPOL的demo运行速度快了很多倍了。

如果进行隔行隔列取样然后在DCT进行累加，经过测试基本上看不出有什么效果上的区别，但是速度大约能提高3倍左右，这部分代码大家可以自己实现下。

DCT滤波的效果上很多情况下也是相当不错的，想比NLM也毫不逊色，我们贴一些图片看下效果：

    

    

　　　　　　　　　　　　       噪音图像　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　去噪后效果（Sigma = 10)

为显示方便，上面的图像都是设置了一定程度的缩放。

共享我改写的这个代码供打击共同学习提高，如果哪位能进一步提高算法的速度 ，也希望不吝赐教。

　　代码下载链接：http://files.cnblogs.com/files/Imageshop/DCT_Denosing.rar

****************************作者： laviewpbt   时间： 2015.11.14    联系QQ:  33184777 转载请保留本行信息**********************

时间： 2024-10-18 01:16:59

优化IPOL网站中基于DCT（离散余弦变换）的图像去噪算法（附源代码）。的相关文章

JPEG压缩原理与DCT离散余弦变换——有实际的数据演示

http://blog.csdn.net/newchenxf/article/details/51719597 1 前言 JPEG是joint Photographic Experts Group(联合图像专家组)的缩写,文件后辍名为"．jpg"或"．jpeg". jpg图片可以说是最常见的图片格式了,基本上你的自拍照,要么是png的,要么就是jpeg的了.(有关jpeg和png的区别,请参考我的另一博文[jpeg 与 png 图片格式的区别]) 但它是一种有损压缩

隐马尔可夫模型中基于比例因子的前向算法(java实现)

直接上干货哈,其他子算法,后续补上. 1 package jxutcm.edu.cn.hmm.model;  2   3 import jxutcm.edu.cn.hmm.bean.HMMHelper;  4 import jxutcm.edu.cn.util.TCMMath;  5   6 /**  7  * [改进后的前向算法]  8  * [带比例因子修正的前向算法 :计算观察序列的概率 ]  9  * [注意] 改进后,就没必要使用后向算法来求观测序列概率了,直接利用中间比例因子scal

深圳网站优化：网站上线前要做哪些准备？

随着搜索引擎的普及,越来越多的人通过搜索引擎进行学习,学习网站搭建也已经不算什么难事,但是有很多朋友网站还没完善好就匆匆将网站上线,导致搜索引擎对网站的评分大大降低,甚至因为网站上线后经常改动,导致网站被降权,犯这种错误的朋友并不少见,所以今天我们就来聊一聊,网站上线前要做哪些准备? 网站上线的底层需求 1.一个合适的域名:域名就相当于网站在搜索引擎中的门牌号,一个没有域名的网站是不能展现到互联网上的,小编推荐在大型的域名供应商购买主流的域名,比如.com.cn.net,不要图便宜去购买中文域名

基于离散余弦变换的同图复制的检测技术

本文是一则学习小结 图像篡改的背景及意义 ??目前图像已经成为人类社会中必不可少的一部分,人们的日常生活中到处都可以见到图像.特别在医学.商业.军事.情报.学术研究.法律和新闻领域中,图像作为原始事件或现象的真实记录,具有信息载体和数字证明的作用.然而人们渐渐发现图像的真实性已经不再可靠了.篡改图像出现的频率越来越大,要分辨出它们也越来越困难. 已知检测技术 数字水印 易损水印,半易损水印和鲁棒水印 水印在图像篡改之前存在 要求所有设备都带有水印装置是不可能的 盲检测技术 不依赖数字水印,是基于

离散余弦变换(C++实现)

理论部分转载自这篇blog: http://blog.csdn.net/luoweifu/article/details/8214959  该blog给出的是java代码,我用c++将其实现了. 理论: 图像处理中常用的正交变换除了傅里叶变换外,还有其他一些有用的正交变换,其中离散余弦就是一种.离散余弦变换表示为DCT( Discrete Cosine Transformation),常用于图像处理和图像识别等. 一维离散余弦变换 正变换                            

让那些为Webkit优化的网站也能适配IE10（转载）

转载地址:http://www.w3cplus.com/css3/adapting-your-webkit-optimized-site-for-internet-explorer-10.html 特别声明:此篇文章由David根据Charles Morris的英文文章原名<Adapting your WebKit-optimized site for Internet Explorer 10>进行翻译,整个译文带有我们自己的理解与思想,如果译得不好或不对之处还请同行朋友指点.如需转载此译文,

C++实现离散余弦变换（参数为二维指针）

http://www.cnblogs.com/scut-linmaojiang/p/5013590.html 写在前面 到目前为止已经阅读了相当一部分的网格水印等方面的论文了,但是论文的实现进度还没有更上,这个月准备挑选一些较为经典的论文,将其中的算法实现.在实现论文的过程中,发现论文中有用到一些空域转频率域的算法.因此也就想到了实现一下离散余弦变换.虽然本文的代码和网上很多已有的代码很类似,思路都没有太多的差别,但是本文有一个比较重要的改进.具体的说,网上现有DCT算法输入的是一个固定的二维数

根据自己的需要，把别人开发好的东西搬过来，优化and重构，在优化的过程中，甚至也会弄出一套全新的东西（转）

赵海平在今年三月份来到阿里,听毕玄(他现任主管)说去年五六月份就跟赵海平聊上了.有人问:为啥 BAT 三大巨头,你看中了阿里巴巴?在今天现场达一千多人的分享中赵海平给出了回复:“因为百度和腾讯没找我呗~”,他笑道,“百度以搜索为核心,优化了很多年了,估计也没啥可以优化的了:而腾讯除了 QQ 和微信,也没什么大型应用(别跟人家说哦)”.这不是原话哈,赵海平还是相当谦虚并且能言的,思维很开阔,两个小时的分享内容丰富,时不时还插两个故事,起初进场的手机和电脑都很自觉的收起来了~旁边的同事侃道:“高 P

怎样搭建一个完整的便于优化的网站

在网站搭建初期就完成所有便于优化的框架和内容能够缩短整个网站优化时间,一个完整的便于优化的网站包含:满足用户需求的程序和内容.便于优化的路径.优先层次的布局.附加价值的内容.长尾关键词布局.合理的内链设置. 学习要点 一.符合用户需求的程序和内容 1.程序.在做网站之前要考察什么程序是符合这个网站的,比如养羊,搜索养羊出现在百度首页的是养羊吧,那么说明搜索养羊的用户有讨论需求,那么论坛类型更加适合网站.但是如果搜索婚纱摄影,出现在百度第一名的是百度图片,就说明图片类的程序更加适合这种类型的网站.