码迷，mamicode.com
首页 > 编程语言 > 详细

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

时间：2015-11-14 23:29:40      阅读：460      评论：0      收藏：0      [点我收藏+]

标签：

　　在您阅读本文前，先需要告诉你的是：即使是本文优化过的算法，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 转载请保留本行信息**********************

 

 

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

标签：

原文地址：http://www.cnblogs.com/Imageshop/p/4965192.html
踩
(0)
赞
(0)
   
举报
评论 一句话评论（0）
登录后才能评论！
分享档案
更多>
2021年07月29日 (22)
2021年07月28日 (40)
2021年07月27日 (32)
2021年07月26日 (79)
2021年07月23日 (29)
2021年07月22日 (30)
2021年07月21日 (42)
2021年07月20日 (16)
2021年07月19日 (90)
2021年07月16日 (35)
周排行
mamicode.com排行更多图片
更多
友情链接
兰亭集智   国之画   百度统计   站长统计   阿里云   chrome插件   新版天听网
关于我们 - 联系我们 - 留言反馈
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码！