BM3D算法实现图像降噪.docx
《BM3D算法实现图像降噪.docx》由会员分享，可在线阅读，更多相关《BM3D算法实现图像降噪.docx（24页珍藏版）》请在冰豆网上搜索。
BM3D算法实现图像降噪
广东工业大学
《数字图像、数字信号处理及应用》
实验报告
题目图像处理综合实验
院、系（部）自动化学院
专业及班级
学号
姓名
日期
1实验目的
理解图像平滑和锐化的基本方法；
了解图像复原的基本方法；
综合使用多种方法改善图像质量。
2实验要求
（1）利用Photoshop之类的图像处理工具软件，尝试对Moon.bmp尽量改善图像质量。
（2）综合采用各种合理的方法，编写程序（C/C++/OpenCV、MATLAB、Python……均可）对Moon.bmp进行图像质量改善，实现以下目标的权衡折中：
a．降低噪声；
b.增大对比度；
c.锐化增强。
3实验设备
安装了VC6/VS2010、PSC6、MATLAB的PC机
4实验原理
4.1利用拉普拉斯算子实现图像锐化
锐化处理的目的是突出图像中的细节或者增强被模糊了的细节。
锐化处理可以用空间微分来完成。
微分算子的响应强度与图像在该点的突变程度有关，图像微分增强了边缘和其他突变（如噪声）而消弱了灰度变化缓慢的区域。
由于我们处理的是数字量，最小灰度级的变化是有限的，而且变化发生的最短距离是在两个相邻像素之间。
因此，
用一阶微分的差值定义一元函数f（x）的二阶微分:
二元图像函数f（x,y）的拉普拉斯变换定义为：
离散方式：
X方向：
Y方向：
故二维拉普拉斯数字实现由以上两个分量相加：
因此拉普拉斯算子用于图像增强的基本方法如下:
拉普拉斯算子处理后的图像如下：
图1（拉普拉斯处理图像对比图。
处理前左，处理后右）
4.2利用分段线性函数实现对比度扩展
分段线性变换函数的对比度拉伸相对于直方图均衡（直方图均衡只能按照统计特性进行变换）可以更加灵活地控制输出灰度直方图的分布，可以有选择地拉伸某段灰度区间，以改善输出图像。
如果一幅图像灰度集中在较暗的区域而导致图像偏暗，可以用灰度拉伸功能来扩展（斜率>1）物体的灰度区间以改善图像；如果图像灰度集中在较亮的区域而导致图像偏亮，也可以用灰度拉伸功能来压缩（斜率<1）物体灰度区间以改善图像质量。
如图2所示，线性函数分为3段，转折点在（c，a）和（d，b）。
从（0，0）到（c，a）段的斜率为
；从（c，a）到（d，b）段的斜率为
；从（d，b）到（Mf,Mg）段的斜率为
。
所以分段函数的表达式为：
图2（分段线性函数示意图）
4.3余弦变换（DCT）
离散余弦变换（DiscreteCosineTransform）是与傅里叶变换相关的一种变换，它类似于离散傅里叶变换（DiscreteFourierTransform）,但是只使用实数。
离散余弦变换相当于一个长度大概是它两倍的离散傅里叶变换，是对实信号定义的一种变换，变换后在频域中得到的也是一个实信号。
相比DFT，DCT可以减少一半以上的计算。
DCT还有一个很重要的性质（能量集中特性）：
大多书自然信号（声音、图像）的能量都集中在离散余弦变换后的低频部分，因而DCT在（声音、图像）数据压缩、图像处理等方面得到了广泛的使用。
●二维余弦变换为：
其中f（x,y）是空间域二维向量之元素，x,y=0,1,2,......N-1；F（u,v）是变换系数阵列之元素。
式中表示的阵列为N×N
●二维余弦逆变换为：
式中的符号意义同正变换式一样
4.4BM3D降噪算法（BlockMatching3DFilterAlgorithm）
一些传统的图像视频去噪算法,会在滤除噪声的同时引入人工噪声或对图像有很大的模糊效果,存在很大的局限性。
而基于块匹配和三维变换域滤波（BM3D）的图像视频去噪算法则采用了不同的去噪策略。
通过搜索相似块并在变换域进行滤波，得到块评估值，最后对图像中每个点进行加权得到最终去噪效果。
BM3D算法不仅有一个较高的信噪比，而且视觉效果也很好。
因此研究者提出了很多基于BM3D的去噪方法，例如:
基于余弦变换的BM3D去噪，基于Anscombe变换域BM3D滤波等等。
BM3D算法的大致算法流程为：
第1步.初始估计
1）逐块估计。
对含噪图像中的每一块
i）分组。
找到它的相似块然后把它们聚集到一个三维数组。
ii）联合硬阈值。
对形成的三维数组进行三维变换，通过对变换域的系数进
行硬阈值处理减弱噪声，然后逆变换得到组中所有图像块的估计值，然后
把这些估计值返回到他们的原始位置。
2）聚集。
对得到的有重叠的块估计，通过对他们进行加权平均得到真实图像的
基础估计。
第2步.最终估计
1）逐块估计。
对基础估计图像中的每一块
i）分组。
通过块匹配找到与它相似的相似块在基础估计图像中的位置，通
过这些位置得到两个三维数组，一个是从含噪图像中得到的，一个是从基
础估计图像中得到的。
ii）联合维纳滤波。
对形成的两个三维数组均进行三维变换，以基础估计图
像中的能量谱作为能量谱对含噪三维数组进行维纳滤波，然后逆变换得到
组中所有图像块的估计，然后把这些估计值返回到他们的原始位置。
2）聚集。
对得到的有重叠的局部块估计，通过对他们进行加权平均得到真实图
像的最终估计。
算法大致流程图如下：
图3（BM3D算法大致流程图）
5软件设计
5.1总体设计
本次实验主要用到对图像进行了去噪、对比度扩展、锐化增强的处理。
其中图像去噪最终使用的是BM3D去噪算法；对比度增强使用的是分段线性对比增强算法；锐化增强是使用拉普拉斯算子实现的。
主程序（main函数）的流程图如下图。
主函数开始时，先读取如工程目录下面的BMP格式文件的数据。
接着按照用户给定的去噪残差进行BM3D算法去噪。
BM3D去噪之后的数据再进行分段对比度增强。
最后，再对数据进行对比度增强处理。
图4（主程序流程图）
BM3D算法、对比度增强、锐化增强函数在下一小节（详细设计）里有更加详细的分析，所以这里就不介绍了。
5.2详细设计
1、锐化增强函数
对于锐化增强算子，如果使用普通锐化算子，即
由于算子h为常量，所以并不能控制锐化的程度和效果。
为了解决这一点，采用MATLAB提供的拉普拉斯锐化算子。
在MATLAB命令行中键入：
h=fspecial（‘laplacian’,0.2）;
即可得到算子的值。
在MATLAB命令行键入：
editfspecial
即可查看MATLAB中拉普拉斯算子的详细计算方式。
详细的拉普拉斯算子计算方式如下：
其中，
；
；
；
为用户控制锐化效果的变量，取值在0~1之间（不包含0和1）。
详细的代码和注释如下：
voidsharp（constBITMAPINFOHEADER*pInfoHead,BYTE*pdata,constDWORDh_len,constdoublegamma）{
inti,j;
doubleG;
/*先根据gamma值算出h1,h2,h3*/
constdoubleh1=gamma/（gamma+1）,
h2=1+4/（gamma+1）,
h3=（1-gamma）/（gamma+1）;
/*创建临时存储空间，用于存储锐化结果*/
BYTE*ptemp=（BYTE*）malloc（sizeof（BYTE）*pInfoHead->biHeight*pInfoHead->biWidth）;
for（i=1;ibiHeight-1;i++）{
for（j=1;jbiWidth-1;j++）{
/*做卷积*/
G=pdata[i*h_len+j]*h2-（（pdata[（i-1）*h_len+j-1]+pdata[（i-1）*h_len+j+1]
+pdata[（i+1）*h_len+j-1]+pdata[（i+1）*h_len+j+1]）*h1）
-（（pdata[（i-1）*h_len+j]+pdata[i*h_len+j-1]+pdata[（i+1）*h_len+j]
+pdata[i*h_len+j+1]）*h3）;
/*防止溢出*/
if（G>255）
ptemp[i*h_len+j]=255;
elseif（G<0）
ptemp[i*h_len+j]=0;
else
ptemp[i*h_len+j]=（BYTE）G;
}
}
/*结果复制到输出*/
for（i=1;ibiHeight-1;i++）{
for（j=1;jbiWidth-1;j++）{
pdata[i*h_len+j]=ptemp[i*h_len+j];
}
}
free（ptemp）;
}
2、中值、均值滤波函数
虽然最后的图像去噪算法采用的是BM3D算法，但是因为后面要把BM3D去噪算法与中值滤波算法作比较，所以在这里还是给出详细的设计。
中值滤波详细代码及注释如下：
voidmedfilt（constBITMAPINFOHEADER*pInfoHead,BYTE*pdata,constDWORDh_len）{
inti,j,k,z,index;
BYTEtemp,nPixel[9];
/*创建临时存储空间，用于中值滤波结果*/
BYTE*ptemp=（BYTE*）malloc（sizeof（BYTE）*pInfoHead->biHeight*pInfoHead->biWidth）;
for（i=1;ibiHeight-1;i++）{
for（j=1;jbiWidth-1;j++）{
/*先排列9个像素*/
nPixel[0]=pdata[（i-1）*h_len+（j-1）];
nPixel[1]=pdata[（i-1）*h_len+j];
nPixel[2]=pdata[（i-1）*h_len+（j+1）];
nPixel[3]=pdata[i*h_len+（j-1）];
nPixel[4]=pdata[i*h_len+j];
nPixel[5]=pdata[i*h_len+（j+1）];
nPixel[6]=pdata[（i+1）*h_len+（j-1）];
nPixel[7]=pdata[（i+1）*h_len+j];
nPixel[8]=pdata[（i+1）*h_len+（j+1）];
/*找中值，这里只需找5次*/
for（k=0;k<5;k++）{
index=0;
for（z=0;z<9-k;z++）
if（nPixel[index]index=z;
temp=nPixel[8-k];
nPixel[8-k]=nPixel[index];
nPixel[index]=temp;
}
ptemp[i*h_len+j]=nPixel[4];
}
}
/*复制结果到输出*/
for（i=1;ibiHeight-1;i++）
for（j=1;jbiWidth-1;j++）
pdata[i*h_len+j]=ptemp[i*h_len+j];
free（ptemp）;
}
均值滤波算法的详细设计及代码：
voidjunfilt（constBITMAPINFOHEADER*pInfoHead,BYTE*pdata,constDWORDh_len）{
inti,j,temp;
/*创建临时存储空间，用于均值滤波结果*/
BYTE*ptemp=（BYTE*）malloc（sizeof（BYTE）*pInfoHead->biHeight*pInfoHead->biWidth）;
for（i=1;ibiHeight-1;i++）{
for（j=1;jbiWidth-1;j++）{
temp=0;
temp+=pdata[（i-1）*h_len+（j-1）];
temp+=pdata[（i-1）*h_len+j];
temp+=pdata[（i-1）*h_len+（j+1）];
temp+=pdata[i*h_len+（j-1）];
temp+=pdata[i*h_len+j];
temp+=pdata[i*h_len+（j+1）];
temp+=pdata[（i+1）*h_len+（j-1）];
temp+=pdata[（i+1）*h_len+j];
temp+=pdata[（i+1）*h_len+（j+1）];
ptemp[i*h_len+j]=（BYTE）（temp/9）;
}
}
/*复制结果到输出*/
for（i=1;ibiHeight-1;i++）{
for（j=1;jbiWidth-1;j++）{
pdata[i*h_len+j]=ptemp[i*h_len+j];
}
}
free（ptemp）;
}
3、分段线性对比度增强函数
因为分段线性对比度增强的详细介绍已经在实验原理部分给出，所以只给出详细代码和注释。
/*指定对比度增强区间的结构体*/
typedefstructadjust_range{
doublein_low;
doublein_high;
doubleout_low;
doubleout_high;
}ADJ_RANGE;
voidadjust（constBITMAPINFOHEADER*pInfoHead,BYTE*pdata,constDWORDh_len,constADJ_RANGE*ar）{
/*
总共分成三段，三段的变换范围为：
第一段：
（0-in_low_c）->（0-out_low_c）
第二段：
（in_low_c-in_high_c）->（out_low_c-out_high_c）
第三段：
（in_high_c-255）->（out_high_c-255）
三段的斜率分别为：
第一段：
h1
第二段：
h2
第三段：
h3
*/
constBYTEin_low_c=（BYTE）ar->in_low*255;
constBYTEin_high_c=（BYTE）ar->in_high*255;
constBYTEout_low_c=（BYTE）ar->out_low*255;
constBYTEout_high_c=（BYTE）ar->out_high*255;
constdoubleh1=ar->out_low/ar->in_low;
constdoubleh2=（ar->out_high-ar->out_low）/（ar->in_high-ar->in_low）;
constdoubleh3=（1-ar->out_high）/（1-ar->in_high）;
inti,j;
doubletemp;
for（i=1;ibiHeight-1;i++）{
for（j=1;jbiWidth-1;j++）{
/*分段线性增强*/
if（pdata[i*h_len+j]temp=h1*pdata[i*h_len+j];
elseif（pdata[i*h_len+j]>in_high_c）
temp=out_high_c+（pdata[i*h_len+j]-in_high_c）*h3;
else
temp=out_low_c+（pdata[i*h_len+j]-in_low_c）*h2;
/*防止溢出*/
if（temp>255）
pdata[i*h_len+j]=255;
elseif（temp<0）
pdata[i*h_len+j]=0;
else
pdata[i*h_len+j]=（BYTE）temp;
}
}
}
4、BM3D去噪函数
对于这次图像复原实验来讲，去噪绝对是复原的首要部分。
因为去噪效果不好的话，接下来的对比度扩展和图像锐化处理会把之前经过去噪处理的图片的噪声又放大了出来；或者要控制对比度扩展和图像锐化的程度，使得处理的效果不太理想。
本人在尝试了很久的空域滤波之后，发现空域滤波的方法是没有办法避免对图像模糊化的，所以决定另外寻找比较好的去噪算法。
经过决定使用BM3D去噪算法。
BM3D去噪算法是这次实验的重头戏，不管是从去噪的效果，还是从去噪后对图像的模糊程度来看，BM3D算法要比中值滤波和均值滤波的效果好多了（BM3D算法与中值滤波、均值滤波的实验效果对比图在“第6章测试与分析”给出）。
不过，BM3D去噪算法卓越的去噪优点牺牲的是比较多的计算时间。
我现在用过的最快的BM3D算法也要花去大概22秒，这还是使用芬兰大学研究所公布的MATLAB动态链接库；而我自己的C语言BM3D去噪算法对一个普通大小的彩色图像处理在Release模式下要花去大概3-4分钟。
BM3D算法这个速度对于DSP图像处理系统来说是难以接受的（不过已经有研究者开始把该算法移植到FPGA上面，并且获得了成功）。
下面是本次使用的C语言版本的BM3D算法的流程框图：
图5（本次实验BM3D算法的流程框图）
如图5所示，BM3D算法的流程主要分为两步：
第1步.初始估计
1）分块。
图像进行分块（块大小为8×8最佳）。
2）余弦变换。
对每一块进行二维余弦变换。
3）分组。
找到它的相似块然后把它们聚集到一个三维数组。
4）联合硬阈值。
对形成的三维数组进行三维变换，通过对变换域的系数进
行硬阈值处理减弱噪声，然后逆变换得到组中所有图像块的估计值，然后
把这些估计值返回到他们的原始位置。
5）聚集。
对得到的有重叠的块估计，通过对他们进行加权平均得到真实图像的
基础估计图像。
第2步.最终估计
1）分块。
对基础估计图像中的每一块进行分块。
2）分组。
通过块匹配找到基础估计图像块与原始图像相似的相似块在基础估
计图像中的位置。
通过这些位置得到两个三维数组，一个是从原始含噪图
像中得到的，一个是从基础估计图像中得到的。
3）联合维纳滤波。
对形成的两个三维数组均进行三维变换，以基础估计图
像中的能量谱作为能量谱对含噪三维数组进行维纳滤波，然后逆变换得到
组中所有图像块的估计，然后把这些估计值返回到他们的原始位置。
2）聚集。
对得到的有重叠的局部块估计进行凯撒滤波之后，通过对他们进行加
权平均得到真实图像的最终估计。
由于C语言开发的BM3D算法的代码长度有大概600多行，为了不影响实验报告的易读性，把该算法的详细代码放在附录部分。
6测试与分析
6.1测试步骤
1、进入工程目录，用VS2010打开工程（.sln）。
2、把解决方案改为Release模式（否则会造成运行时间多出大概一倍）。
3、依次点击“调试”->“开始执行（不调试）”或者“Ctrl+F5”开始图像处理。
4、等待4-5分钟，命令行显示“totalfinish!
”。
如图6所示。
5、到工程目录下的“BM3D”夹下打开“out.bmp”查看处理结果。
图6
由于要进行彩印，对“Moon.bmp”图像处理的结果存放于附录中。
6.2比较中值、均值、BM3D滤波信噪比
使用MATLAB测试中值、均值、BM3D滤波三个算法的信噪比提高量。
信噪比提高量的公式如下：
其中，c（x,y）为退化图像；x（x,y）为原始图像；a（x,y）为还原图像。
计算信噪比提高量的代码如下：
%@imgx：
原始图像
%@imgc：
退化图像
%@imga：
还原图像
functionsnri=snri_cal（imgx,imgc,imga）
c_x=imgc-imgx;
a_x=imga-imgx;
S=sum（sum（（c_x-mean（mean（c_x）））.^2））;
N=sum（sum（（a_x-mean（mean（a_x）））.^2））;
snri=10*log10（S/N）;
end
接下来先使用MATLAB准备一幅加噪声的图片，编写的MATLAB代码如下：
Img=imread（‘cameraman.tif’）;
Imgn=imnoise（cameraman,'gaussian',0,0.002）;
imwrite（Imgn,‘cameramen.bmp’）;
接着在VS2010环境下运行程序，分别对噪声图像“cameram.bmp”进行中值、均值、BM3D滤波处理。
处理完成之后，编写代码计算信噪比提高量，代码如下：
imgx=imread（'cameraman.tif'）;%原始图像
imgc=imnoise（imgx,'gaussian'）;%退化图像
imga=imread（'中值滤波.bmp'）;%还原图像
snri_cal（imgx,imgc,imga）%中值滤波的信噪比提高量
imga=imread（'均值滤波.bmp'）;%还原图像
snri_cal（imgx,imgc,imga）%均值滤波的信噪比提高量
imga=imread（'BM3D滤波.bmp'）;%还原图像
snri_cal（imgx,imgc,imga）%BM3D滤波的信噪比提高量
以上代码输出的结果如下：
中值滤波
均值滤波
BM3D滤波
信噪比提高量（db）
4.903
升级会员 