第6次课 医学图像恢复技术

1、医医 学学 图图 像像 处处 理理Medical Image Processing陈家益陈家益()广东医学院 信息工程学院 计算机科学教研室主要内容：主要内容：5.2 医学图像中的噪声及滤除医学图像中的噪声及滤除5.3 系统退化模型系统退化模型5.1 图像恢复的基本概念图像恢复的基本概念第五章 医学图像恢复技术（3）噪声图像的恢复3）顺序统计滤波器修正后阿尔法均值滤波器：中值滤波器： g(x,y)=median(f(x,y)最大值滤波器： g(x,y)=max(f(x,y)最小值滤波器： g(x,y)=min(f(x,y)中点滤波器：g(x,y)=(1/2)*(min(f(x,y)+max(f

2、(x,y)( , )1g( , )( , )xys tSx yf s tmnd （3）噪声图像的恢复3）顺序统计滤波器修正后阿尔法均值滤波器：适用于混合噪声中值滤波器： 适用于脉冲噪声最大值滤波器： 适用于胡椒噪声，滤除暗点最小值滤波器： 适用于盐噪声，滤除亮点中点滤波器：适用于高斯噪声和均匀随机分布的噪声（3）噪声图像的恢复修正后阿尔法均值滤波 设g(x,y)表示带噪图像，Sxy表示图像中点(x,y)的m*n邻域。对邻域Sxy的所有像素g(s,t)排序，去掉d/2个最低灰度值和d/2个最高灰度值。用gr(s,t)表示剩余的m*n-d个像素。由这些像素的平均值形成的滤波器称为修正后阿尔法均值滤

3、波器。( , )1( , )( , )xyrs tSf x yg s tmnd （3）噪声图像的恢复修正后阿尔法均值滤波 其中0dmn-1，当d=0,修正后阿尔法均值滤波退化为算术均值滤波；当d=mn-1，滤波退化为中值滤波。当d取中间的其他值时，修正后阿尔法均值滤波适用于处理包括多种噪声的混合情况。( , )1( , )( , )xyrs tSf x yg s tmnd 修正后阿尔法均值滤波function f=alphatrim(g,m,n,d)g=im2double(g);f=imfilter(g,ones(m,n),symmetric);for k=1:d/2 f=imsubtract

4、(f,ordfilt2(g,round(k),ones(m,n),symmetric);endfor k=(m*n-(d/2)+1):m*n f=imsubtract(f,ordfilt2(g,round(k),ones(m,n),symmetric);endf=f/(m*n-d); （3）噪声图像的恢复修正后阿尔法均值滤波 若对修正后阿尔法均值滤波稍加改变，就可以得到K邻近均值（中值）滤波。它是在Sxy中选择K个与待处理像素的灰度差为最小的像素，再用这K个像素的均值（中值）替换原来的像素，因此有mn-d=K.( , )1g( , )( , )xys tSx yf s tmnd （3）噪声图像

5、的恢复修正后阿尔法均值滤波 为了减少平滑后图像的模糊，在对当前像素进行平滑时，应尽量保留边界点和细节点，而选择非边界点非细节点进行平滑。修正后阿尔法均值滤波和K邻近均值（中值）滤波都是对与当前像素的灰度差较小的邻域像素进行均值或中值处理，具有一定的边界点和细节点保持效果。（3）噪声图像的恢复中值滤波器g=medfilt2(f,m n, symmetric);最大值滤波g=ordfilt2(f,m*n,ones(m,n),symmetric);最小值滤波g=ordfilt2(f,1,ones(m,n),symmetric);（3）噪声图像的恢复中点滤波g1=ordfilt2(f,1,ones(m

6、,n),symmetric);g2=ordfilt2(f,m*n,ones(m,n),symmetric);g=imlincomb(0.5,g1,0.5,g2);（3）噪声图像的恢复4）自适应滤波器 自适应中值滤波算法(Adaptive Median Filtering Algorithm, AMF)是一种中值滤波，使用一个邻域区域作为滤波窗口，在滤波过程中会根据一定的设定条件改变滤波窗口的大小。（3）噪声图像的恢复4）自适应滤波器 另外，自适应中值滤波算法主要通过判断当前滤波窗口的中心像素点的性质来决定滤波的输出：当判断窗口中心的像素是脉冲噪声时，将中心像素用滤波窗口中的像素排序后的中值代替

7、；否则不改变当前的像素值，直接输出。（3）噪声图像的恢复4）自适应滤波器 自适应中值滤波算法分为两步进行：噪声检测和噪声滤波，基本步骤如下： Wij表示正被处理的、中心在(i,j)处的滤波窗口； S0为默认的初始窗口大小，Smax为允许的最大窗口大小。 fij为(i,j)处的灰度值，fmin为Wij中的最小灰度值，fmax为Wij中的最大灰度值，fmed为Wij中的灰度中值。 （3）噪声图像的恢复4）自适应滤波器 自适应中值滤波算法工作在两个层面，表示为LevelA和LevelB。LevelA: 如果fminfmedfmax，转向LevelB，否则，增大窗口尺寸；如果窗口WijWmax，重复L

8、evelA，否则输出fmed。LevelB: 如果fminfijfmax，输出fij，否则输出fmed。（3）噪声图像的恢复function f=adpmedian(g,Smax)if(Smax1.)endf=g; f(:)=0; alreadyProcessed=false(size(g);for k=3:2:Smax zmin=ordfilt2(g,1,ones(k,k),symmetric); zmax=ordfilt2(g,k*k,ones(k,k),symmetric); zmed=medfilt2(g,k k,symmetric); processUsingLevelB=(zmed

9、zmin) & (zmaxzmed) & alreadyProcessed; zB=(gzmin) & (zmaxg);（3）噪声图像的恢复 outputZxy=processUsingLevelB & zB; outputZmed=processUsingLevelB & zB; f(outputZxy)=g(outputZxy); f(outputZmed)=zmed(outputZmed); alreadyProcessed=alreadyProcessed| processUsingLevelB; if all(alreadyProcessed(:) break; endendf(al

10、readyProcessed)=zmed(alreadyProcessed);（1）常见的系统退化函数 实际获得的图像不仅包含噪声，还包含各种因素引起的图像退化，如相机聚焦不佳、系统散射等。对于线性移不变成像系统，退化因素可由退化函数h(x,y)反映。（1）常见的系统退化函数 成像系统的系统函数或退化函数h(x,y)又称点扩散函数 (PSF)。输出的退化图像可以表示为原图像与点扩散函数的卷积，再加上系统的噪声。( , )( , )*( , )( , )g x yh x yf x yn x y（1）常见的系统退化函数 图像恢复的效果取决于对图像退化过程的先验知识所掌握的准备程度，知道关于退化函数

11、h和噪声n的信息越多，对原图像的估计就越准确。 对于退化图像的恢复，关键是确定或得到系统的点扩散函数h。（1）常见的系统退化函数 成像系统的退化函数一般可以通过对成像系统退化过程的先验知识得到。 常用的估计方法有：模型估计法图像观察估计法实验估计法（1）常见的系统退化函数图像观察估计法 利用退化图像自身携带的信息对退化模型进行估计，它通过选择图像中具有简单结构或较强特征的局部区域（如点源的模糊图像），并构造一个具有相同大小和特征、但没有退化的估计图像（点源）通过解卷积来得到退化函数的近似估计。（1）常见的系统退化函数图像观察估计法 如果能确定原图像的某部位有一个清晰的点，假定噪声可忽略，则使得

12、该点退化的模糊图像就是h(x,y)。 从原图像中找到一块信号内容较强的区域，构建出该区域的估计图像，进而大致推出退化函数h(x,y):( , )( , )( , )SSSgx yh x yfx y （1）常见的系统退化函数图像观察估计法图像观察估计法特点和应用场合： 没有任何关于退化函数 h 的知识，只能从图像本身来收集信息。 选择强信号区域（如高对比度区域）进行处理，是为了降低噪声的影响。 处理过程复杂，仅在特殊环境下使用，如恢复具有历史价值的老照片。（1）常见的系统退化函数实验估计法 实验估计法通过使用或设计一个与原图像退化过程相似的成像系统装置或过程，并利用这个系统或过程对点源成像，从而

13、得到原系统点扩散函数的较准确估计。（1）常见的系统退化函数实验估计法实验估计法特点和应用场合： 可以使用与获取退化图像设备一致或相似的装置时使用。 用上述系统对一个冲激成像，冲激可由一个亮点来模拟，该点应尽可能亮，以便将噪声的影响降低到可以忽略的程度。 对冲激的成像结果即为PSF函数。（1）常见的系统退化函数模型估计法 模型估计法基于对不同成像过程基本原理和物理特性的了解和建模。模型估计法特点和应用场合： 一般需要从基本物理原理推到数学模型。 模型建立后应用范围广。（1）常见的系统退化函数 可通过调用函数fspecial创建常见的退化滤波器。h=fspecial(type，parameter)

14、;常用的滤波器：type parametergaussian hsize，sigmasobelprewitt（1）常见的系统退化函数常用的滤波器：type parameterlaplacian alphalog hsize，sigmaaverage hsizedisk radiusMotion len, theta（1）常见的系统退化函数用退化滤波器退化图像：h=fspecial(motion,10,60);g=imfilter(f,h); h=fspecial(disk,10);g=imfilter(f,h);（2）逆滤波 对于线性移不变系统，有退化模型 为了方便讨论和求解，用顺序堆叠的方式

15、将二维图形表示为一维，上述退化模型可转化为：原图像的恢复就转化为求解上面的线性方程组。( , )( , )*( , )( , )g x yh x yf x yn x ygHfn（2）逆滤波 即 在无噪声的理想情况下，有： 原图像可以通过退化图像与系统矩阵的逆估计出来，这就是逆滤波。1()fHgn 1fHg （2）逆滤波 实验证明，当图像的噪声为零或非常小时，逆滤波的效果很好。 但在H很小或为零的情况下，逆滤波的效果会受到很大影响，因为当H很小或为零时，噪声会被放大。这就意味着，退化图像中小的噪声干扰在H接近零时会被迅速放大，甚至淹没原图像信号。（2）逆滤波 解决上述问题的方法是，避开H的零点或

16、极小值点。 可有两种途径：一、在H等于零或接近于零值的地方，人为设置H的值（如一个固定的小值）。二、使图像的恢复在有限区域内进行，以避免H的零点或极小值点的影响。（3）维纳滤波 根据逆滤波的理论，从统计学的观点出发，提出了维纳滤波器。 维纳滤波恢复的目标是寻找一个滤波器，使得恢复后的图像与原图像的均方差最小，维纳滤波器又称最小均方误差滤波器。（3）维纳滤波 在图像信号接近平稳随机过程，噪声和图像不相关且均值为零的情况下，复原图像可通过下式估计： H表示系统退化矩阵，T表示矩阵的转置，Rf和Rn表示原图像和噪声图像的自相关函数。 即为维纳滤波器。11()TTfnfH HR RH g11()TTf

17、nH HR RH （3）维纳滤波 维纳滤波器 能够自动抑制噪声的放大。 当H为零时，由于括号内第二项的存在，逆不为零，而此时滤波器的值为零。 当图像信噪比很高或不存在噪声时，括号内第二项很小或为零（噪声的自相关函数远小于原图像或为零），此时维纳滤波器成为逆滤波器。 当噪声的自相关函数远远大于原图像时，滤波器的值接近零，从而避免了对噪声的放大。 11()TTfnH HR RH （3）维纳滤波 噪声和原图像的自相关函数不容易获得，滤波器常用一个K值代替。 称为有参维纳滤波。在滤波过程中，通过调节K值，得到较好的恢复。1()TTfH HKH g （3）维纳滤波 通过调用函数 deconvwnr 可实

18、现维纳滤波。 g= deconvwnr(f,PSF,NCORR,FCORR) 或者 g= deconvwnr(f,PSF,NSR) PSF为点扩散函数，NCORR, FCORR为可选参数，分别为噪声和原图像的自相关函数。NSR为可选参数，表示信噪比，即前面所说的K值，缺省值为0,此时维纳滤波变成逆滤波。（3）维纳滤波图像的功率谱：Sf=abs(fft2(f).2;图像的平均功率谱：fA=sum(Sf(:)/prod(size(f);图像的信噪比：NSR=nA/fA;图像的自相关函数: FCORR=fftshift(real(ifft2(Sf);（3）维纳滤波逆滤波实例：f=phantom(2048);psf=fspecial(motion,10,45);motionf=imfilter(f,psf);gaus=imnoise(motionf