文档章节

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

abcijkxyz
 abcijkxyz
发布于 2016/11/22 16:38
字数 2905
阅读 22
收藏 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滤波的效果上很多情况下也是相当不错的,想比NLM也毫不逊色,我们贴一些图片看下效果:

                         

    

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

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

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

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

 

  后记:  继续优化了下8*8点的DCT里SSE代码的处理方式,改变了累加的方向,速度提高30%;然后把64次阈值处理的代码也改成SSE优化,速度提高10%;在如果考虑进行隔行隔列取样然后在DCT进行累加,经过测试基本上看不出有什么效果上的区别,但是速度大约能提高3.5倍左右;如果考虑多线程的方式,比如开个双线程,速度约能提高0.8倍,如果CPU支撑AVX,则大概又能提高0.9倍,算来算去,我感觉可以实时了。

 

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

 

 

 

本文转载自:http://www.cnblogs.com/Imageshop/p/4965192.html

共有 人打赏支持
abcijkxyz
粉丝 63
博文 6196
码字总数 1876
作品 0
深圳
项目经理
私信 提问
SSE图像算法优化系列二十一:基于DCT变换图像去噪算法的进一步优化(100W像素30ms)。

  在优化IPOL网站中基于DCT(离散余弦变换)的图像去噪算法(附源代码) 一文中,我们曾经优化过基于DCT变换的图像去噪算法,在那文所提供的Demo中,处理一副1000*1000左右的灰度噪音图像耗...

Imageshop
2018/09/04
0
0
【DCT】OPENCV python 离散余弦变换问题

因为项目需要,目前在研究《基于离散余弦变换和区域生长的白粉虱图像分割算法》,使用的OPENCV python。 论文地址:%e5%9f%ba%e4%ba%8e%e7%a6%bb%e6%95%a3%e4%bd%99%e5%bc%a6%e5%8f%98%e6%8d...

majisong
2016/06/02
1K
0
离散余弦变换(DCT)的来龙去脉

1. 图像的二维离散变换   与一维的有限长离散非周期信号存在傅里叶变换(DFT)一样,图像作为一个二维离散信号同样存在着二维离散变换(注意这里是介绍一个通用的概念,二维离散变换,包括...

dugudaibo
2017/11/01
0
0
Java进阶(五十七)-基于感知哈希算法的pHash图像配准算法

Java进阶(五十七)-基于感知哈希算法的pHash图像配准算法   毕业论文提交之后,老师交给自己一项任务:图像配准,也就是给你两幅图像,通过系统来判定两幅图像是否为同一副图像。自己作为这...

sunhuaqiang1
2017/04/18
0
0
离散余弦变换(DCT)

DCT变换、DCT反变换、分块DCT变换 一、引言 DCT变换的全称是离散余弦变换(Discrete Cosine Transform),主要用于将数据或图像的压缩,能够将空域的信号转换到频域上,具有良好的去相关性的性...

li_wen01
2017/06/05
0
0

没有更多内容

加载失败,请刷新页面

加载更多

钢铁侠自曝:特斯拉“哨兵模式” 来了!

北京时间1月23日消息,特斯拉车主可能很快就可以查看和记录汽车在无人看管时受到的损坏。 一位特斯拉车主在Twitter上向特斯拉CEO埃隆·马斯克(Elon Musk)抱怨称,他的车门凹进去了一大块,希...

Linux就该这么学
31分钟前
3
0
Java基础教程,第十一讲,阻止继承与抽象类

上一节课我们学习了Java中的继承语法,以及如何继承中的方法重写和super()关键字,今天我们将学习一下如何阻止继承以及抽象类。 课程内容 阻止继承:final类和方法 有些时候我们希望阻止某个...

程序员补给栈
44分钟前
1
0
mount命令

mount命令可以将分区挂接到Linux的一个文件夹下,从而将分区和该目录联系起来,因此我们只要访问这个文件夹,就相当于访问该分区了。 mount [-t vfstype] [-o options] device dir loop:用来...

Danni3
50分钟前
1
0
推荐一款接口 API 设计神器!

今天栈长给大家推荐一款接口 API 设计神器,传说中的,牛逼哄洪的 Swagger,它到底是什么?今天为大家揭开谜底! Swagger是什么? 官网:https://swagger.io/ Swagger 如官网所示,它是最好的...

Java技术栈
今天
10
0
AMD直奔5nm!这一步棋下得妙

AMD今年将推出采用7nm工艺的第二代EPYC霄龙、第三代Ryzen锐龙处理器,其中后者已经在CES 2019上公开首秀,性能追评i9-9900K,功耗则低得多。 虽然被称为“女友”的GlobalFoundries临时决定放...

linuxCool
今天
3
0

没有更多内容

加载失败,请刷新页面

加载更多

返回顶部
顶部