课程设计自适应的快速非局部医学图像去噪算法.docx
《课程设计自适应的快速非局部医学图像去噪算法.docx》由会员分享，可在线阅读，更多相关《课程设计自适应的快速非局部医学图像去噪算法.docx（23页珍藏版）》请在冰点文库上搜索。
课程设计自适应的快速非局部医学图像去噪算法
中北大学
课程设计说明书
学生姓名：
陈杰学号：
1005084122
学院：
信息与通信工程学院
专业：
生物医学工程
题目：
自适应的快速非局部医学图像去噪算法
指导教师：
张权职称:
讲师
2013年6月28日
中北大学
课程设计任务书
12/13学年第二学期
学院：
信息与通信工程学院
专业：
生物医学工程
学生姓名：
陈杰学号：
1005084122
课程设计题目：
自适应的快速非局部医学图像去噪算法
起迄日期：
2013年6月7日～2013年6月28日
课程设计地点：
201、503
指导教师：
张权
系主任：
王浩全
下达任务书日期:
2013年6月7日
课程设计任务书
1．设计目的：
综合运用医学图像处理的相关知识，基于理论推导，设计一种自适应的快速非局部医学图像去噪算法，利用MATLAB作为编程工具进行计算机实现，从而加深对所学知识的理解。
2．设计内容和要求（包括原始数据、技术参数、条件、设计要求等）：
设计一种自适应的快速非局部医学图像去噪算法，以改善图像的视觉效果、利于临床诊断。
3．设计工作任务及工作量的要求〔包括课程设计计算说明书（论文）、图纸、实物样品等〕：
（1）查阅文献资料，了解相关知识；
（2）确定设计方法；
（3）编写MATLAB程序，并上机调试；
（4）撰写课程设计说明书。
课程设计任务书
4．主要参考文献：
●要求按国标GB7714—87《文后参考文献著录规则》书写，例：
1傅承义，陈运泰，祁贵中.地球物理学基础.北京：
科学出版社，1985
（5篇以上）
5．设计成果形式及要求：
提供课程设计说明书一份；
MATLAB程序；
6．工作计划及进度：
6月7日~6月12日：
查资料；
6月13日~6月19日：
确定设计方案；
6月20日~6月25日：
编写程序，上机调试；
6月26日~6月27日：
完成课程设计说明书；
6月28日：
答辩。
系主任审查意见：
签字：
年月日
目录
摘要·················································1
一课题题目及背景·······························2
1.1课题题目·······································2
1.2课题背景·······································2
二课题目的及意义···································2
三算法简介·········································2
3.1非局部图像去噪算法·····························3
3.2噪声与滤波参数h的关系·························4
3.3最佳滤波参数h的确定···························5
3.4高斯噪声标准差的估计···························6
3.5计算复杂度的降低·······························8
四matlab程序简介···································9
4.1程序实现概要···································9
4.2具体参数的确取值·······························9
五图像处理结果及分析·······························9
5.1图像处理结果···································9
5.2结果分析·······································11
六课题总结及心得体会·······························12
七参考文献·········································12
八附录·············································12
摘要
由于图像获取、传输等过程中常常会产生噪声,造成图像质量下降,因此图像去噪作为图像处理的基本技术之一。
本课题既是基于matlab的自适应快速非局部的医学图像去噪算法。
本课题主要在非局部均值去噪（NL-Means）基础上对算法做适当改进与简化处理，提出一种自适应的非局部均值去噪方法，以更加方便的应用到临床医学图像处理中。
本文中对算法做了简介，以及介绍了参数h和a的简单估算方法。
通过实验说明本算法达到了预期效果，完成了课题要求。
在本文最后附录中我们给出了具体的matlab实现程序，以方便验证。
一、课题题目及背景
1.1、课题题目：
自适应的快速非局部医学图像去噪算法
1.2、课题背景
由于图像的采集获取及传输过程中图像质量会受到各种因素的影响,常常会产生噪声造成图像质量下降,因此图像去噪技术对作为图像处理的基本技术,对于获得高质量的图像至关重要.现实中的数字图像在数字化和传输过程中受到成像设备与外部环境噪声干扰等影响，称为含噪图像或噪声图像。
减少数字图像中噪声的过程称为图像去噪。
实际中常见的噪声主要有加性噪声和乘性噪声。
常用去噪方法：
均值滤波器，自适应维纳滤波器，中值滤波器，形态学噪声滤除器，小波去噪。
在临床医学应用中,医学图像在疾病诊断的应用及其广泛,现在已成为临床诊断及治疗的过程中必不可少的部分,并且对其依赖性极高.因此,图像处理在医学中的应用非常重要,图像的质量有可能会影响到疾病的诊断甚至是误诊,对医学图像的高效的可靠去噪处理在医学应用中也极为关键.本文既是基于此原因对原有的医学图像去噪算法作出改进,以利于在实际中的应用。
二、课题目的及意义
本课程目的主要是综合运用医学图像处理的相关知识，基于理论推导，设计一种自适应的快速非局部医学图像去噪算法，利用MATLAB作为编程工具进行计算机实现，从而加深对所学知识的理解。
本文中我们对非局部均值（NL-Means）图像去噪算法进行了改进,提出一种定量计算法滤波参数最优值的方法,由噪声图像估计噪声方差,进而由噪声方差与图像方差估计滤波参数h。
其中根据局部区域加权欧氏距离的对称性,将算法中复杂度最高的两像素间距离计算由两次降为一次,从而在不损失性能的条件下使计算复杂度降低到原来的一半左右,在不影响图像处理质量的情况下节省了处理时间,方便了临床的应用。
三、算法简介
在算法自适应方面,有些算法利用了噪声方差、图像方差等整体信息,有些算法利用了局部灰度相似性、局部直方图等图像局部信息。
噪声水平自适应的非局部图像去噪算法,根据噪声方差和图像方差估计最优h参数,利用了图像的整体信息,去噪效果接近NL-Means所能达到的最优性能。
NL-Means算法的另一个不足之处是计算量较大,原因是在计算每两个像素之间距离时,考虑了其周围一定区域内的多个像素,造成了处理时间较长的后果。
我们根据欧氏距离的对称性,对基于局部图像块计算的任意两点之间的距离只进行一次,算法整体时间复杂度约降为原算法的一半。
3.1、非局部图像去噪算法
由于NL-Means算法在对每个像素的加权平滑过程中考虑了局部结构的相似性,取得了很高的滤噪效果。
其基本思想是:
当前像素值由图像中所有与它结构相似的像素加权平均得到。
对于每个像素的权值,采用以它为中心图像子块（一般取7x7）与以当前像素为中心子块之间的高斯加权欧氏距离来计算,权值设为此距离的负指数函数值。
这样做的好处是在估计当前像素值时,局部结构上与它相似的像素权重较大,而结构相似像素上叠加的噪声是随机的,因而通过加权均值可有效去除噪声。
其处理过程如下:
设噪声图像为v={v（i）|iI},去除噪声后图像为NL[v],对每个像素i的值通过下式加权得到
w（i,j）为依赖于像素i与像素j相似程度（高斯加权欧氏距离）的权值:
式中,||v（Ni）-v（Nj）||
a表示像素i与像素j之间以它们中心的子块内像素之间高斯（标准差为a）加权距离的平方,v（Ni）表示i周围的局部子块像素集合,h为滤波参数。
权值介于0和1之间,且满足
实际计算过程中,为避免计算量太大,参与加权的并非图像中所有像素,而是利用它周围一定大小的区域。
图1
NL-means算法执行过程示意图
对于算法的实际执行，还需要设置两个窗口尺寸，一个是像素邻域窗尺寸KxK，一个是像素邻域窗搜索范围的窗口大小LxL，即在LxL大小的区域里面选择像素的邻域大小为KxK执行NL-means算法，KxK的窗口在LxL大小的区域里滑动，根据区域的相似性确定区域中心像素灰度值的贡献权值。
如右图所示。
依据算法原理，像素邻域及其搜索范围应是整幅图像，之所以选择此两窗口大小，是为了提高计算效率，如果在整个图像区域内计算，那么执行效率太低，工程实用性较差.对于常见的噪声较大的图像，一般取K=7，取L=21就足够了，对于低噪声的图像K=3，L=7基本就能满足降噪要求。
K和L可以根据图像噪声的水平而做适当调整。
在算法中有两个待选参数:
α和h，参数α由选定像素邻域的窗口大小决定，加高斯卷积核意味着距离中心像素越远，权值就会越小，反之则越大。
一般α取邻域窗口大小的1-10倍。
而滤波参数h由噪声的标准差决定，一般取噪声标准差σ的10-15倍。
Buades通过理论分析和实验结果表明,NL-means算法在主客观性能上都优于常见的图像去噪算法,如高斯滤波、各向异性滤波、总误差最小化、邻域滤波等等。
但算法也存在两个不足之处:
关键滤波参数h如何设定不明确,算法计算复杂度较高。
为此我们提出了ANL-means自适应快速非局部算法。
它根据噪声方差和图像方差估计最优h参数,利用了图像的整体信息,去噪效果接近NL-means所能达到的最优性能。
同时根据局部区域加权欧式距离的对称性，将算法中复杂度最高的两像素间距离计算由两次降为一次从而在不损失性能的条件下算法降低一半。
3.2、噪声与滤波参数h的关系
为了定量地分析噪声水平与去噪算法滤波参数的关系,以及这些参数是如何与图像内容相关的。
选择几个典型的图像进行了较大范围的滤波参数变化实验,结果如图1所示。
图中不同曲线对应不同的图像,噪声标准差为20,横坐标为h参数的变化过程,纵坐标为去噪结果与原始图像的均方误差（MSE）。
图2去噪性能随滤波参数h的变化
从图中可以看出:
（1）h参数确实存在一个最优参数,太大或大小都会造成性能下降;
（2）尽管图像内容不同,但在同样噪声水平下最优h参数值基本一致。
3.3、最佳滤波参数h的确定
Buades在NL-Means算法中采用与噪声标准差线性正比关系确定滤波参数h,即h=10
。
通过理论和实验分析得出h应与噪声方差
有近似线性正比关系,并受到噪声图像方差的影响。
从式
（2）可以看出,为了不损失噪声之外的信息,当两个像素的加权距离大于一定值时,权值应接近于零,由于式
（2）中的权值函数为典型的高斯函数（均值为0,标准差为h/
）,可能通过改变h值来调节权值分布。
根据高斯分布的特点,距离大于标准偏差某一倍数的所有权值之和可以估计出来,如距离大于3h/
的权值之和小于1%。
如果设定了可以忽略的权值和,便可以得到这一倍数值。
现假设这一倍数值为
则不被忽略的距离应满足:
另一方面,在假设噪声为高斯分布时,与当前像素点及其邻域像素灰度相同的点在叠加了的高斯噪声（
）后会与当前像素值有一些差异,这些差异决定了式（5）左边的距离值。
同理,根据高斯分布的特点,应该有一定比例的像素落在与当前考虑像素值的加权绝对差值不超过噪声标准差
倍的范围之内。
当设定可忽略权值和与计算加权距离时可忽略噪声像素百分比（由于这里考虑是加权距离,并非精确的像素百分比）相同时,即有
由此可以得出滤波参数h与噪声水平的近似关系:
此时考虑了所有与当前像素加权绝对差值不超过
的像素点,并且被忽略权值和小于某一个由
确定的比率。
根据高斯分布,
取值为2和3时,对应比例分别为95.4%和99.7%。
以上只是从噪声水平进行了分析,但是
的取值与图像本身的平滑程度也有一定关系。
当
值较大时,包含了较高比例的噪声点,但同时
也可能会将图像中本应与当前像素取值不同的点包含进来,因此
值并非越大越好。
实验中发现取统计的权值或噪声像素比例约为98%比较合适,此时
进一步分析可知,当图像内容比较平缓（标准差较小）时,像素之间差值较小,
应该取较小的值；相反,图像内容比较剧烈（标准差较大）时,像素之间差值较大,
可以取较大的值。
因此,可以根据图像本身标准差对这一估计进行进一步的调整。
由于图像内容较为复杂,标准差所包含的信息有限,准确的定量关系很难确定。
但可以结合实验分析给出一个比较好的估计。
实验发现,在图像本身标准差小于一定值时,对h取值影响较小。
对于灰度为256级的图像,小于70的标准差对h值影响较小；另外,当噪声水平较高时,图像标准差的增大对取值影响也较大。
因此,对于灰度级数为256的图像采用下式估计h值:
式中,
为噪声图像的标准差。
当采用这一参数对噪声图像进行滤波时,可以达到近似最优的效果。
h参数值过小会使许多噪声没有被考虑;反之,h参数值过大,会平滑掉一些超出噪声范围的图像差异,从而造成图像过于平滑、丢失边缘等细节信息。
3.4、高斯噪声标准差的估计
噪声水平与去噪算法参数有密切的关系。
要达到较好的去噪效果需要估计一幅图像中的噪声水平,对于高斯白噪声,均值为0,只需估计其标准差参数。
比较简单的噪声估计方法包括基于图像块的方法、基于滤波器的方法,以及二者相结合的方法。
基于图像块的方法将图像分为许多小块,计算其各自方差,以若干最小方差的均值作为估计结果；基于滤波器的方法先对图像进行一次平滑,再计算原图与平滑后图像的差别,由此差值图像估计噪声水平;二者相结合的方法先通过基于图像块方法选择一些平缓块,再利用滤波器方法根据这些块估计噪声水平。
由于基于滤波器的方法首先要选择适当的滤波器和参数,而这些参数本身又与噪声水平相关,使得噪声的估计较为困难。
另外,基于滤波器的方法计算复杂度也高于基于图像块的噪声估计方法。
因此,采用基于图像块的噪声估计方法。
将图像划分为M
M的无重叠小块,对每一块计算其方差。
在图像本身灰度基本一致的块中,其方差可近似为噪声方差,为避免个别噪声干扰,采用总块数1%（不少于两个子块）的最小方差均值作为估计结果。
实验中采用固定25
25的子块计算区域方差。
另外采用基于Immerker提出的噪声方差估计方法。
其滤波器设计理念是：
边缘类的图象结构具有很强的二阶微分成份，噪声估计器应当对拉普拉斯算子或者相应的微分算子比较敏感，固此，我们就可以利用两个模板L1和L2之间的差别来估计出图象的噪声方差。
如果将这两个模板设计为近似图象的拉普拉斯算子，则L1和L2为
。
噪声估计模板N可有下式得出
，
模板N的均值为零。
如果假设每一个象素的标准偏差为
的话，模板的方差可以表示为：
我们用
表示在图象I中（x，y）位置应用模N后的值。
因此图象，中的噪声方差可以由下式计算（W，H分别表示图宽和图高）：
那么，所求的图象噪声绝对偏差可表示为：
由此可得到较为理想的噪声方差估计，计算也较为简化，另外由于图像的不同可根据情况设计较为理想的模板。
本文由于主要研究项目为设计去噪算法，故噪声方差的估计只给出估计方法并不设计具体处理程序。
3.5、计算复杂度的降低
NL-Means算法中,复杂度最高的是两像素之间加权高斯距离的计算,因为每个像素都要计算其周围区域内所有像素与它的距离值,并按此距离计算出加权的权值。
Mahmoudi对NL-Means算法的计算复杂度进行了改进,计算周围像素的权值时,只有满足以下条件的像素才考虑:
（1）以两个像素为中心子块的灰度均值之比介于一定范围内；
（2）两个像素点的梯度均小于某一阈值,且方向差异小于某一阈值。
但这种方法对计算复杂度的降低是靠减少计算加权像素值时的邻域像素个数来达到的,如果算法中几个相关阈值确定不当,会造成去噪性能上的损失。
这里仅从计算方式的角度来降低计算复杂度,不会造成任何性能上的损失。
考虑到距离d（i,j）=∀v（Ni）-v（Nj）∀22,a是对称的,即d（i,j）=d（j,i）,可以对任意两个像素的距离只计算一次,同时累加到两个像素点各自的累积加权像素值和、及累积权值和之中。
设i点为当前像素,j点为计算距离及权值的像素,计算得到的权值为W（i,j）,则在i点的累积加权像素值和中加入v（j）*W（ij）,同时在j点的累积加权像素值和中加入v（i）*W（ij）,并在i和j的累加权值和中分别加入W（ij）。
图像边缘像素要做一些特殊处理。
这样,在累加过程中,总的计算量可降低为原来的一半左右,而估计噪声水平的计算复杂度与去噪算法本身相比是非常小的。
四、matlab程序简介
4.1、程序实现概要
在去噪处理的matlab程序中我们先定义了函数functionIout=NLMeans（Iin,a,h,K,L），然后选取典型的图像读入，进行加噪声处理，并显示原图和去噪图像。
然后通过对高斯核a,邻域窗口K,搜索窗口L，以及过滤参数h的定义实现对滤波去噪的控制。
程序的主体通过四层for循环的嵌套对图像进行处理滤波。
最后输出具体结果。
具体matlab实现程序见附录。
4.2、具体参数的取值
高斯核的标准差a的确定可以根据大量实验选取处理效果较好的值，在其他值不变的情况下a值越大处理图像相对越平滑。
邻域窗口K,搜索窗口L的确定具有相关性。
通过实验可知两者值如果太小或者太接近处理结果较为模糊，太大则对图像处理质量贡献不大，但是处理时间较长。
滤波参数h通过下式得到：
其中
为噪声标准差，
为噪声图像的标准差。
各参数的取值相互影响，实际应用中要根据具体情况而定具体值。
五、图像处理结果及分析
5.1、图像处理结果
原图像：
加高斯噪声后的图像：
加方差为0.005的高斯噪声
加方差为0.01的高斯噪声
对方差为0.01噪声图像去噪处理结果：
窗口K=5，L=19，a=1的处理结果
窗口K=7，L=19，a=1的处理结果
窗口K=7，L=21，a=1.5的处理结果
窗口K=7，L=21，a=1的处理结果
对方差为0.005的噪声图像去噪处理的结果：
窗口K=7，L=19，a=1的处理结果
窗口K=5，L=19，a=1的处理结果
窗口K=7，L=21，a=1.5处理结果
窗口K=7，L=21，a=1的处理结果
5.2、结果分析
由上处理结果对比可知：
图像的处理结果受到多个因素影响窗口的选择，高斯核函数的取值都会影响处理质量，图像中所加噪声的大小也会对结果造成较大影响。
这里列举了添加标准差分别为0.01和0.005的高斯噪声的处理，通过与原图像以及两种情况下处理结果比较可知噪声标准差越大经相同处理后细节可能丧失越多。
对于窗口的选择，两窗口较为接近时处理较为模糊，窗口较小结果也会不尽理想。
窗口较大则对质量贡献不大，但是会大大增加处理时间，这里邻域窗口选择5和7，搜索窗口选取19和21时不论时间还是处理质量都在可接受范围内。
对于高斯核函数a的选择,相同条件下a值较大处理也较平滑，但是太大也会是图像失真严重，在这里选择1和1.5较为合适。
六、课题总结及心得体会
经过几天的努力，这次的课程设计终于圆满完成了。
由于缺乏医学图像去噪处理的知识，我们通过查找资料和询问老师，终于克服各种难关，顺利完成了这次的课程设计。
这次的课程设计，使我们进一步掌握了医学图像的相关知识，更多的是学到了很多课本里学不到的知识，达到了学以致用的目的，避免了华而不实的局面。
通过这次的课程设计，使我们深刻了解到知识储备的重要性，也更加坚定了我们努力学习，不断充实自我的决心和信心。
七、参考文献
1.秦襄培、郑贤中.MATLAB图像处理宝典.北京：
电子工业出版社，2011
2.陈显毅.图像配准技术及其MATLAB编程实现.北京：
电子工业出版社，2009
3.张强、王正林.精通MATLAB图像处理.北京:
电子工业出版社.2009
4.赵书兰.MATLAB数字图像处理与分析实例教程.北京.北京工业出版社.2009
5.王志明、张丽.自适应的快速非局部图像去噪算法.中国图像图形学报.第14卷.第4期,2009
6.张权、桂志国、刘祎.医学图像的自适应非局部均值去噪算法.计算机工程.第38卷.第7期.2012
7.许俊玲、罗晓军.基于邻域相似性的快速非局部滤波算法.商丘职业技术学院学报.第8卷.第44期.2009
八、附录
functionIout=NLMeans（Iin,a,h,K,L）%定义功能函数
I=imread（'图像名.类型'）;
sa=0.005;
Iin=imnoise（I,'gaussian',0,sa）;
Figure,imshow（I）;%显示原图片
title（'原图像'）;
figure,imshow（Iin）;%显示高斯加噪图片
title（'加噪图像'）;
K=9;%邻域窗尺寸，为奇数
L=27;%邻域窗搜索范围的窗口，为奇数
a=1.5;%高斯内核标准差
q=std2（Iin）;%图像标准差
n=max（0,（q-70））/10;
h=（3.30+n）*sa;%过滤的程度
[nx,ny]=size（Iin）;
Iout=zeros（nx,ny）;
fori=1:
nx;
forj=1:
ny
%初始化邻域窗口
KxStart=-floor（K/2）;
KxEnd=floor（K/2）;
KyStart=-floor（K/2）;
KyEnd=floor（K/2）;
ifj==254
end
%修改窗口附近防止出边界
KxStart=max（1,i+KxStart）-i;
KxEnd=min（nx,i+KxEnd）-i;
KyStart=max（1,j+KyStart）-j;
KyEnd=min（ny,j+KyEnd）-j;
%初始化搜索窗口
LxStart=-floor（L/2）;
LxEnd=floor（L/2）;
LyStart=-floor（L/2）;
LyEnd=floor（L/2）;
升级会员 