HyzhangAllrightsreserved.基于迭代(自动阈值)算法的医学图像增强方法1算法原理介绍1.1基于灰度阈值的图像分割原理采用灰度阈值对图像进行分割是图像分割的基本方法之一。通过设定灰度阈值,把图像的象素点按灰度等级进行分类,把图像分割成若干子图像。在实际应用中,最常用的一种分割方法是将图像分割成高灰度区和低灰度区两部分,组成一幅二值图像。该分割方法按式1对图像进行操作。(1)其中,T为预先设定的灰度阈值,f为输入图像,g为输出图像。基于灰度阈值对图像进行分割的基本条件是式(1)中阈值T的选择。虽然采用人工方法往往可以获得较理想的阈值,但是在很多情况下,需要计算机自动完成阈值的选择,这样就要求有合适的算法对图像的灰度直方图进行分析,选择合适的阈值。在实际中常采用迭代算法或Ostu法对阈值进行自动计算[1],本文仅介绍迭代算法。迭代算法的基本思想是:首先设定一个阈值的估计值;采用一定的算法反复对该估计值进行修正,保证每次修正后的结果都优于前一次;当进行一定次数的修正之后,结果趋于收敛,即相邻两次的结果的差异较小,当该差异小到可接受范围时,表明一个理想的阈值已经求出。最后利用该阈值按式(1)对图像进行操作,即完成了图像的自动阈值分割。但是,一些图片由于照度不均、阴影、对比度差异等,使得如果采用同一阈值对整幅图片进行处理(即全局阈值)时会出现不兼容图像各处的情况,使得分割效果变差。这时,可以考虑将图片分割成若干子图片,将每个子图片按自动阈值算法进行处理,然后再将各个子图片的处理结果合并成整体结果输出,该方法称局部阈值或动态阈值法。1.2基于图像自动阈值分割的边缘检测采用全局阈值或局部阈值获得的二值图像可以方便地应用在图像的边缘检测当中。由于图像已经转换为二值图像,所以对图像中边缘信息的提取较为方便,仅需判断某像素是区域内部点还是边缘点即可。实际操作中,对二值图像中黑色(或白色)点的四邻域进行判断,若该点的四邻域均为黑色(或白色),即可判断该点为区域内部点而不是边缘点。依次对图像中每个象素的四邻域进行判断,即可得到图像的边缘信息。HyzhangAllrightsreserved.2程序设计本次实验中,程序采用Matlab语言编写。采用基于迭代的自动阈值算法分别实现了对图像的全局阈值分割、局部阈值分割,并利用全局阈值分割的结果,实现图片的边缘检测,增强了图片中的边缘信息。2.1迭代算法的实现[2]根据1.1节中介绍的迭代算法思想,采用如下语句实现利用迭代算法找出用于图像分割的阈值。t=mean(gray_image(:));%设置估计值is_done=false;%迭代完成标志位count=0;%迭代计数归零while~is_done%迭代循环r1=find(gray_image=t);%按前次结果t对图像二分r2=find(gray_imaget);temp1=mean(block(r1));%求r1区域均值ifisnan(temp1);%保证结果非NaNtemp1=0;endtemp2=mean(block(r2));%求r2区域均值ifisnan(temp2)%保证结果非NaNtemp2=0;endt_new=(temp1+temp2)/2;is_done=abs(t_new-t)1;%检测两次迭代计算结果的差值t=t_new;%更新结果count=count+1;%迭代计数增1ifcount=1000%判断是否迭代次数过多Error='Error:Cannotfindtheidealthreshold.'%出错提示Return%退出程序endend该程序中,首先将估计值设为灰度图像的灰度均值。该估计值对最终处理结果不产生影响,但该值越接近最终结果,则迭代计算的次数越少,增加程序效率。程序中采用is_done标志来决定是否继续进行迭代计算,is_done标志由判断两次迭代结果的差值决定,本程序中设定若两次迭代结果相差小于1则视为迭代结束。在某些极端情况下,该算法无法找出理想阈值,迭代计算将陷入死循环,所以在程序的迭代循环中,除了进行迭代结果的计算外,还增加了迭代次数的检测,当迭代运算进行了1000次仍没有找出理想阈值时,视为迭代失败,退出程序。同时,还对区域求均值得操作结果作了检验,如果结果为NaN,就将结果更正为0,这样增强了整个程序的健壮性,防止程序陷入死循环。2.2全局阈值图像分割的实现对图像进行全局阈值分割按照图1的流程进行。HyzhangAllrightsreserved.读取图像并转换为灰度图像输出原图及其灰度直方图利用求出的阈值对图像进行二值化处理求出用于图像分割的阈值输出结果及其灰度直方图图1全局阈值图像分割程序流程其中,自动求阈值采用2.1节中的程序实现。对图像进行二值化处理采用如下程序进行。[m,n]=size(gray_image);%获得原图尺寸result=zeros(m,n);%申请内存空间用来存放结果result(r1)=255;%将图像二值化处理resule(r2)=0;经过上述程序,得到的二值化图像就存放在了数组result中。全局阈值图像分割程序的整体代码见附录I。2.3局部阈值图像分割的实现局部阈值图像分割的原理与全局阈值分割相似,只是在使用2.1节的程序前,将图像分割成若干子图像进行处理,在处理后,在将各自的结果拼接起来。该部分程序流程图见图2。读取图像并转换为灰度图像输出原图读取图像块输入分块大小对图像块进行自动灰度阈值分割将该结果块合并入结果图像是否为最后一个图像块否是输出结果图2局部阈值图像分割程序流程HyzhangAllrightsreserved.该程序中,对每个图像块进行自动阈值分割的程序与2.1节中的程序相同,不再赘述。采如下程序对图像进行分块与合并结果。block_size=input('Blocksize=');%输入分块大小fori=1:block_size:m%依次对图像分块并处理forj=1:block_size:nif((i+block_size)m)&&((j+block_size)n)block=gray_image(i:end,j:end);%右下角区块elseif((i+block_size)m)&&((j+block_size)=n)block=gray_image(i:end,j:j+block_size-1);%最右列区块elseif((i+block_size)=m)&&((j+block_size)n)block=gray_image(i:i+block_size-1,j:end);%最下行区块elseblock=gray_image(i:i+block_size-1,j:j+block_size-1);%普通区块end…………%该部分为基于迭代的灰度阈值分割程序,同2.1节,略if((i+block_size)m)&&((j+block_size)n)result(i:end,j:end)=block;%右下角区块elseif((i+block_size)m)&&((j+block_size)=n)result(i:end,j:j+block_size-1)=block;%最右列区块elseif((i+block_size)=m)&&((j+block_size)n)result(i:i+block_size-1,j:end)=block;%最下行区块elseresult(i:i+block_size-1,j:j+block_size-1)=block;%普通区块endendend在程序中,由于无法保证图片的尺寸能被输入的区块大小整除,即最右列和最下行区块可能为非完整区块,所以在分块时要对这些区块特殊处理。程序中采用if语句来判定区块类型,并选择合适的操作。该部分完整代码见附录II。2.4利用二值图像进行图像边缘的检测利用迭代算法得到了图像的二值图像后,可利用该二值图像进行图像边缘的检测。从二值图像中提取图像的边缘较为方便,既对图像中黑色或白色区域中每一个像素的四邻域进行判断即可判断该像素是区域内点还是边缘点。用于边缘检测的代码如下:edge=zeros(m,n);fork=2:1:m-1forj=2:1:n-1ifresult(k,j)==255%白色区域if((result(k,j)==255)&&(result(k+1,j)==255)&&(result(k-1,j)==255)&&(result(k,j+1)==255)&&(result(k,j-1)==255))%判断四邻域edge(k,j)=255;%区域内点elseHyzhangAllrightsreserved.edge(k,j)=0;%边缘点endelseedge(k,j)=255;%其他区域endendend经过上述代码,将在edge数组中得到图像的边缘。将得到边缘合并到原图,即可完成图像的边缘增强。采用如下代码完成图像边缘与原图像的合并。mix=gray_image;mix=uint8(mix);fork=1:1:mforj=1:1:nifedge(k,j)==0mix(k,j)=0;endendend经上述程序,最终的结果存放在mix数组中。该部分完整程序见附录III。3结果分析本实验中,选用头部MRI照片(图3)用于评估上述程序的性能。图3脑瘤患者头部MRI照片HyzhangAllrightsreserved.图4图3的灰度直方图3.1全局阈值分割结果采用基于迭代算法的全局阈值分割对图3进行操作,经2次迭代运算,得到结果图5,自动求得的分割阈值为53,即图4中的谷点。图5全局阈值分割得到的二值图像3.2局部阈值分割结果采用局部阈值对图像进行分割,分割的结果与选择的子图像区块大小有直接关系。分别以10×10、50×50和100×100为区块大小进行处理,结果分别见图6(a)、6(b)和6(c)。HyzhangAllrightsreserved.图6局部阈值分割结果从图6的结果中可以看出,区块的大小对局部阈值分割的结果有很大的影响。当选择区块较小时,由于很多区块仅处于背景区域中,所以将背景中的噪声被放大了。3.3边缘检测及增强的结果对利用结果图5,进行边缘检测,得到图像的边缘见图7。图7边缘检测的结果将图7与原图合并,得到最终结果图8。HyzhangAllrightsreserved.图8边缘增强结果及该图局部放大4结论基于迭代算法能够叫高效地自动求出灰度阈值,对于大多数图像,迭代算法的结果都是收敛的,即能在有限迭代次数内求理想阈值。该迭代算法可以用于全局阈值分割和局部阈值分割处理中。对于一般的图像,采用全局阈值分割能获得较好的分割效果,并且利用分割后得到的二值图像可以方便地提取出图像的边缘信息,用于图像边缘提取与增强,得到较清晰的图像轮廓。但是,对于由于曝光不均等原因导致的不同位置灰度差异过大的图片,采用全局阈值分割得到的结果往往不理想。对于这样的图片,可以采用局部阈值分割进行处理。采用局部阈值分割时,分割的结果受对图片进行分块的大小的影响。分块时,图片的分块不宜过小,因为如果分块过小,每个区块内部的像素过少,统计规律不明显,分割效果不好,且会对背景中的噪声敏感。参考文献[1]陈冬岚,刘京南,余玲玲.几种图像分割阈值选取方法的比较与研究[J].电气技术与自动化,2003,(1):77-80.[2]姚敏.数字图像处理(M).北京,机械工业出版社,2006.HyzhangAllrightsreserved.附录I:全局阈值分割程序original_image=imread('w:\head.jpg');gray_image=rgb2gray(original_image);subplot(1,2,1);imshow(original_image);subplot(1,2,2);imhist(gray_image);gray_image=double(gray_image);t=mean(gray_image(: