ch8- 图象复原2.ppt

1、逆滤波比较简单，但没有清楚地说明如何处理噪声，而维纳滤波综合了退化函数和噪声统计特性两个方面进行复原处理.,5 维纳滤波器,逆滤波方法不能完全恢复原始信号f(x,y)，而只能求出f(x,y)的一个估计值,希望找到一种方法，在有噪声条件下，从退化图像 g(x,y)复原出f(x,y)的估计值，该估计值符合一定的准则。,维纳滤波 (Wiener filtering)=最小均方差滤波,维纳滤波是最常用的图像恢复方法,基于维纳滤波的图像恢复方法是1967年提出的,C.W. Helstrom, “Image restoration by the method of lest sqaures,” Journ

2、al of the Optical Scoiety of America, vol.57, no.3, pp.297-303, 1967.,C.W.Helstrom, This weeks citation classic, 1982 1967-1982年SCI引用超过125次.,N.Wiener, “The extrapolation, interpolation and smoothing of stationary time series”, New York: Wiely, 1949.,在均方误差值最小的准则下得到的 称为对f(x,y)的最小二乘方估计。,按照该准则得到的滤波器叫维纳滤

3、波器。,目标：使得复原后图像 与原始图像 的均方 误差最小：,因此维纳滤波器又称为最小均方差滤波器,线性滤波：寻找点扩散函数hw(x,y)，使得,则有,由Andrews和Hunt推导满足这一要求的传递函数为：,这里， 是成像系统传递函数H(u,v)的复共轭； Sn(u,v) 是噪声功率谱； Sf (u,v)是输入图像的功率谱。,或者：,计算退化图像g(x,y)的二维Fourier变换G(u,v),Wiener滤波的过程:,计算理想图像的频谱估计,计算点扩展函数h(x,y)的二维Fourier变换H(u,v),计算退化 图像和噪声的功率谱Sf(u,v),Sn(u,v),计算滤波器HW(u,v),

4、求反Fourier变换,这一方法有如下特点： （1）当H（u,v）0或幅值很小时，分母不为零，不会造成严重的运算误差。 （2）在信噪比高的图像中，即Sn(u,v)Sf(u,v),维纳滤波复原法特点：,如果没有噪声，就成为逆滤波,(3) 当理想图像功率谱Sf (u,v)0)时 ,表明我们不可能从全是噪声的图像中恢复出任何有意义的信号,（4）往往未退化图像的功率谱Sf (u,v)难以知道，用下式近似表示：,维纳滤波函数：,测试维纳滤波效果：,逆滤波和维纳滤波的比较,维纳滤波的结果非常接近原始图像,比逆滤波要好,逆滤波和维纳滤波的比较,(a)运动模糊及均值为0方差为650的加性高斯噪声污染的图像 (

5、b) 逆滤波的结果 (c) 维纳滤波的结果 (d)-(f) 噪声幅度的方差比(a)小1个数量级 (g)-(i) 噪声幅度的方差比(a)小5个数量级,逆滤波(inverse filtering),线性代数恢复,图像恢复方法,频域图像恢复,维纳滤波(Wiener filtering),代数逆滤波,奇异值伪逆滤波。,线性代数恢复在1970s提出,B.R.Hunt, “A matrix theory proof the discrete convolution theorem,” IEEE Trans. Audio Electroacoust , vol.19, no.4, pp.285-288, 1

6、971.,B.R.Hunt, “The application of constrained least squares estimation to image restoration by digital computer,” IEEE Trans. Computers, vol.22, no.9, pp.805-812, 1973.,约束最小二乘滤波,Matrix-vector形式,投影法。,6 代数逆滤波,无约束代数逆滤波 (unconstrained restoration),已知退化图像的向量形式g和退化矩阵H,则无约束逆滤波恢复的的图像为:,结果,设恢复的图像为 ，如果不考虑噪声，

7、用 表示恢复误差。则代数逆滤波的目的是最小化目标函数,即,推导:,使导数为零,无约束代数逆滤波 (unconstrained restoration),约束最小二乘滤波 (constrained least squares restoration),令Q为f的线性矩阵算子，约束最小二乘恢复问题就是在满足约束条件 下，使 最小化的问题。,例如选择Q为有限差分矩阵，使得二阶差分的能量最小,约束最小二乘滤波 (constrained least squares restoration),f (x, y)在(x, y)处的二阶微分,用拉格朗日法求,微分，,可以用来调节以满足约束条件。,约束最小二乘滤波

8、 (constrained least squares restoration),推导,clear F=checkerboard(8); subplot(2,2,1),imshow(F,); title(a) 原始图像,FontSize,12) PSF=fspecial(motion,7,45); MF=imfilter(F,PSF,circular); noise=imnoise(zeros(size(F),gaussian,0,0.001); MEN=double(MF)+noise; subplot(2,2,2),imshow(MEN,); title(b) 运动模糊+高斯噪声后的图像,

9、FontSize,10) NSR=sum(noise(:).2)/sum(MEN(:).2); subplot(2,2,3),imshow(deconvwnr(MEN,PSF),); %matlab自带维纳滤波函数，若噪声功率谱设为0，则为逆滤波 title(c) 逆滤波后的图像,FontSize,10) subplot(2,2,4),imshow(deconvwnr(MEN,PSF,NSR),);%维纳滤波 title(d) 维纳滤波后的图像,FontSize,10),Matlab约束滤波函数：deconvreg,奇异值（SVD）伪逆滤波 ( SVD=Singular Value Decom

10、position),把退化矩阵H进行SVD分解，并得到H的广义逆矩阵。从而可以把逆滤波表示成迭代形式，通过控制迭代次数，可以在恢复效果较好时停止迭代，提供了一种人机交互机会。,基本思想,其中U的列向量是AAT的特征值， V的列向量是ATA的特征值； 是对角矩阵，对角线上的值i称为H的奇异值。,退化矩阵H的SVD分解为,于是恢复图像为,则退化矩阵H的广义逆矩阵为,写成迭代形式,7 图象灰度校正和几何畸变校正,一、灰度校正 假设理想图象为f（x,y），由于灰度失真因 子D（x,y）的影响，实际得到的图象为g（x,y）， g（x,y）=D（x,y）f（x,y） 灰度校正是要从畸变的图象g（x,y）中

11、复原原 始图象f（x,y）。一种最直接的方法是用光密度计 测量出被拍摄景物中的某一部分区域S内真实的灰 （亮）度数值C，而对应的图象灰度为gc（x,y）， 则 代入图象复原方程则有,三、几何畸变校正 在图象获取或显示过程中可能产生图象的几何失真， 如下图所示,从上图可以看出，几何畸变是将无失真坐标系中函 数f（x,y）变换到另外一个坐标上，例如，原先在（x,y）点上的象素（灰度）变化到（u,v），在图象上 反映有些位置被挤压，而另一些位置被扩张。我们希望 找到这两个坐标系之间的关系。,几何基准图像的坐标系统用(x, y)来表示 需要校正的图像的坐标系统用(x, y)表示,设两个图像坐标系统之间

12、的关系用解析式表示,通常h1(x,y)和h2(x,y)用多项式来表示：,通常用线性畸变来近似较小的几何畸变,更精确一些可以用二次型来近似,若基准图像为f(x,y)，畸变图像为g(x,y)，对于景物上的同一个点，假定其灰度不变，则,若已知h1(x,y)和h2(x,y),通常用线性畸变来近似较小的几何畸变,零级内插,双线性插值,则：,h1(x,y)和h2(x,y)未知,通常用已知的多对对应点来确定系数a, b,线性畸变,可由基准图找出三个点(r1, s1), (r2, s2), (r3, s3)与畸变图像上三个点(u1, v1), (u2, v2), (u3, v3)一一对应。,将对应点代入，有：

13、,解联立方程组，得出6个系数。,此时,二次型畸变,有12个未知量，需要6对已知对应点,代入上式,同样有,解方程组，得到ai, bi. 12个系数。,记做向量矩阵形式:,通常实际应用中，会取多余的对应点对，这时A不是方阵而是高矩阵，这时矩阵的逆用广义逆矩阵来求解。高矩阵的广义逆矩阵为,在广告制作和计算机动画中常常要使物体变形。,内插方法除了零级内插和双线性内插外，还有B样条插值和sinc插值函数内插,几何畸变复原的一套方法也可以用于使图像畸变的工作中。,clear f=checkerboard(24); subplot(2,2,1),imshow(f) title(a) 原始图像,FontSiz

14、e,10) s=0.7; theta=pi/6; T=s*cos(theta) s*sin(theta) 0; -s*sin(theta) s*cos(theta) 0; 0 0 1; tform=maketform(affine,T); g1=imtransform(f,tform,nearest); subplot(2,2,2),imshow(g1) title(b) 最邻近插值变换,FontSize,10) g2=imtransform(f,tform); subplot(2,2,3),imshow(g2) title(c) 双线性插值变换,FontSize,10) g3=imtrans

15、form(f,tform,Fillvalue,0.5); subplot(2,2,4),imshow(g3) title(e) 修改c的背景为灰色,FontSize,10),clear f=imread(Fig0515(a)(base-with-control-points).tif); subplot(2,2,1),imshow(f) title(a) 原始图像,FontSize,10) inputpoints=83 81;450 56;43 293;249 392;436 442; outputpoints=68 66;375 47;42 286;275 434;523 532; tform=cp2tform(inputpoints,outputpoints,projective); g1=imtransform(f,tform); subplot(2,2,2),imshow(g1) title(b) 几何失真图像,FontSize,10) tform=cp2tform(outputpoints,inputpoints,projective); g2=imtransform(g1,tf