图像的恢复与重构

单击此处编辑母版标题样式,*,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,第五章 图像的恢复与重构,什么是图像退化：图像的质量变坏叫做退化。退化的形式有图像模糊、图像有干扰等。,图像退化的处理方法：无论是由光学、光电或电子方法获得的图像都会有不同程度的退化；退化的形式多种多样。如传感器噪声、摄像机未聚焦、物体与摄像设备之间的相对移动、随机大气湍流、光学系统的相差、成像光源或射线的散射等；如果我们对退化的类型、机制和过程都十分清楚，那么就可以利用其反过程来复原图像。,典型的图像复原方法是根据图像退化的先验知识建立一个退化模型，以此模型为基础，采用滤波等手段进行处理，使得复原后的图像符合一定的准则，达到改善图像质量的目的。,一、图像退化模型,f(i, j)：,原始图像,g(i, j)：,降质图像,T(,)：,成像系统的作用，则：,g(x,y)=Tf(x,y),设,T,是线性移不变的。,一幅连续的图像,f(x,y),可以用抽样函数的二维卷积,表示：,因此，,令,h(x,;y,),=T,(x-,y-) ，,则有：,定义于不在原点的二维,函数,由于,f(, ),与,x，y,没有关系,称,h(x,;y,),为点扩散函数(,PSF),或系统冲击响应。多数情况下它表现为时不变的，反映在图像中为位移不变的，则,h(x,;y,),可以表示为,h(x-,y-),其中*表示卷积运算。如果,T(,),是一个,h,可分离系统，即：,h(x,;y,)=h,1,(x,) h,2,(y, ),则二维运算可以分解为列和行两次一维运算来代替。,在加性噪声情况下，图像退化模型可以表示为：,g(x,y)=f(x,y)*h(x,y)+n(n,y),其中,n(x, y),为噪声图像。,二、离散图像退化模型,对于图像降质过程进行数学建模，设：,f(i, j),为原始图像；,y(i, j),为降质图像；,h(i, j; k, l),为点扩散函数；图像为,M,N,维。有,假设为空间不移变,h(i, j; k, l),，,则：,线性位移不变的图像退化模型则表示为：,g(x,y)=f(x,y)*h(x,y)+n(x,y),结论：如果已知,g(x,y)、 n(x,y)、 h(x,y)，,则,f(x,y),可以计算出来。对等式两端取傅立叶变换有：,G(u,v)=F(u,v)H(u,v)+N(u,v),F(u,v)= (G(u,v)- N(u,v)/ H(u,v),f(x,y)= F,-1,F(u,v),g(x,y),f(x,y),n(x,y),h(x,y),+,三、循环矩阵及傅立叶化,一个一维离散序列通过一个系统发生失真的过程可用下图表示,f= H,-1,g- n,用矩阵表示，可以写成,g(x),f(x),n(x),h(x),+,如果考虑加性噪声，根据离散序列的卷积定理，有,扩展为周期为,M,的,序列,由于离散卷积的周期性，有,h,e,(x)=h,e,(x+M)，H,可以写成,H,是一个循环阵。,结论：,离散卷积都可以写成：输入矩阵循环矩阵！,对数字图象的二维离散函数也是如此。,对图像退化模型而言，有,A=5,B=5,M=9,M=9,用矩阵形式表示上式,：,g=,Hf,+n,g、f,和,n,分别表示,MN,的函数矩阵,g,e,(i, j)、,f,e,(i, j),和,n,e,(i, j),的各行前后相连而成的列矢量(堆叠矢量)。,如果假设原始图像是,M,N,维矩阵，则,H,是,MN,MN,循环矩阵，且,H,是一个分块(,M,M,个)循环矩阵：,每一个子矩阵,H,i,自身也是循环矩阵,N,N：,=,+,MN1,MNMN MN1,MN1,1、一维信号序列循环矩阵的对角化和傅立叶化,解,矩阵方程：,f= H,-1,g- n,最简单的计算方法就是对角化，,H,H,-1,也是对角阵。,对角化,H,的方法,求取其特征值和特征矢量。,对循环矩阵而言，设：其有,M,个特征值和特征矢量。,由于,w,(k),是由傅立叶系数构成的，因此,w,(k),彼此是正交的。所以，由,w,(k),构成的变换矩阵是可逆的。,用特征矢量组成的矩阵：,W,=,w,(0),w,(1),w,(M-1),生成对角矩阵,D,：,D=W,-1,HW,；,且,D(k,k)=(k)。,而,根据周期性，,M i = - i,凑,一个常数,所以,H=WDW,-1,对,g=,Hf,+n,而言，可以写成,g=WDW,-1,f+n,，,有,W,1,g=DW,-1,f+ W,1,n,；,其中对,W,-1,g,列矢量的每一行,G(k),而言，有,对,W,-1,f,列矢量的每一行,F(k),而言,对,W,-1,n,而言有同样的结果。所以对,W,1,g=DW,-1,f+ W,1,n,而言，,G(u)=MH(u)F(u)+N(u)。,上面的过程称之为循环矩阵的傅立叶化。,如果图像的,g,、,f,、,n,采用堆叠矢量的方法构成，,g=,Hf,+n,。,同一维的情况类似，不同的地方是,H,为块循环矩阵，以及其中的傅立叶变换是二维的，但最后结论是一样的。,G(u)=MH(u)F(u)+N(u) ,一维情况下的结论,G(u,v)=MNH(u,v)F(u,v)+N(u,v),二维情况下的结论,在实际应用中，认为,f(x,y)、h(x,y)、g(x,y)、n(x,y),的维数是相等的。,F(u,v) =G(u,v)-N(u,v)/ MNH(u,v) f(x,y),=F,-1,F(u,v),2、二维信号序列块循环矩阵的对角化和傅立叶化,=,+,3、,H(u,v),的获取,要知道一个图象降质系统的,H(u,v),是一件非常困难的事情。但因为,f(x,y)*h(x,y)=g(x,y)，,有,F(u,v)H(u,v)=G(u,v),如果,f(x,y)=(x,y)，F(x,y)=1,则,H(u,v)=G(u,v),所以，可以用实验的方法得到,h(x,y),和,H(u,v)；H(u,v),可用点源的输出图像的傅立叶变换来近似。,另外，有一些图象降质系统的,H(u,v),有固定的或近似的数学模型。,四、常见的线性移不变降质算子,运动模糊：通常在拍摄过程中，相机或物体移动造成的运动模糊可以用一维均匀邻域像素灰度的平均值来表示：,大气扰动模糊：这种模糊经常出现在遥感和航空摄影中，由于曝光时间过长引起的模糊可用高斯点扩散函数来表示：,式中,K,是一个归一化常数，保证模糊的大小为单位值，,2,可以决定模糊的程度。,均匀不聚焦模糊,这是由于相机聚焦不准确引起的，虽然不聚焦由许多参数决定，如相机的焦距、相机光圈的大小、形状、物体和相机之间的距离等，但在研究中为了简单起见，我们用下列函数表示聚焦不准引起的模糊：,均匀二维模糊,这是最常见的一种模糊，可以用来近似聚焦不准引起的模糊：,其中,L,是奇数。,五、无约束恢复,逆,滤波,对于图像退化模型：,g(x,y)=f(x,y)*h(x,y)+n(x,y),两边取傅立叶变换：,G(u,v)= F(u,v) H(u,v) +N(u,v),H(u,v),又称为系统的转移函数(或滤波函数)，它使图像退化。,在无噪声的情况下，上式可以简化为：,G(u,v)= F(u,v) H(u,v),F(u,v) =G(u,v)/H(u,v),这种1/,H(u,v),的形式称为逆滤波,。再进行傅立叶逆变换就可以得到,f(x,y)。,什么是无约束恢复,当对,噪声一无所知时，使,n,无约束的小。,由于,g=,Hf,+n,，,假设通过恢复可以得到一个不错的,f,的估计,f,。,显然，,f,应,满足关系,g,Hf,n,，,我们希望,n,尽可能的小。于是问题转化为,f,在什么情况下,n,最小,对矩阵而言就是它的迹(对角线之和)的平方最小。,实际情况中，噪声是不可避免的，因而只能求,F(u,v),的估计值：,如果,H(u,v),有许多零点，必然使得复原的结果受到极大影响。 或者如果,H(u,v),不为零但是有非常小的值，也即病态条件，也会使复原效果受到影响。,解决这个问题的方法是避开,H(u,v),的零点。 幸好一般的,H(u,v),在低频附近的有限区域内不为零。,因此逆滤波可以在原点附近进行，相当于在频域乘上一低通窗口函数,W(u,v)。,为了防止随着,u、v,的增大,H(u,v),的迅速减小而增设一些条件,由于截断地原因，被恢复的图象振铃较大。,一种改进的方法是取(?),其中，,d、k,均为小于1的常数。,逆滤波的应用条件：退化图像,g(x,y),是信噪比较高的图像。,设：,对,降质图像,g=,Hf,+n,而言，所谓恢复可以看成是,对,g,通过某种线性变换,L,而得到原图像的估计值,f,。,其中，,L,是,H,的某种逆过程。如果,f,与,f,之间的差满足某些条件,约束条件，例如：满足给定的均方误差，则称该图像被恢复了。,六、有约束恢复(某种准则,),的一般公式,噪声是独立的，与图像不相关。求其期望应为0。,期望值越小，其迹也越小,因此,对求迹而言，式中第2项和第3项的相等，有,问题转化为：,L,等于什么时，均方误差最小。,由此可得最佳,L,：,0,0,0,0,将,R,f,，,R,n,看成是一个循环块矩阵，,对其进行对角化，有：,图像的相关度约为2030个象素,合理的假设,由最佳,L,：,可得,两边乘以,W,和,W,-1,：,又,因为当,D,为对角阵时有,所以,当,L,为分块循环矩阵时,，,WL W,-1,也是一对角阵。其对角线上元素的值为循环序列(“主值”序列)的傅立叶变换值。同理，,D,，,D,*,，,A,，,B,都为对角阵，其对角线上的元素的值为对应循环序列(“主值”序列)的傅立叶变换值。于是上式,可以写成标量形式,合理的假设,f(x,y)=g(x,y)*l(x,y),七、维纳滤波恢复,基本原理,研究发现，逆滤波复原方法对噪声极为敏感，要求信噪比较高(100以上,)，,通常不满足该条件。,因此希望找到一种方法，在有噪声条件下，从退化图像,g(x,y),复原出,f(x,y),的估计值,f,(x,y) ，,该估计值应符合一定的准则。,用向量,f,g,n,来表示,f(x,y), g(x, y), n(x,y)，,Q,为对,f,的线性算子，在约束条件,|,g,-,Hf,|,2,=|,n,|,2,下求,Qf,的,最小化|,Qf,|,2,而得到,f,的最佳估计,f,。,用拉格朗日法求解,设,R,f,和,R,n,为,f,和,n,的相关矩阵，根据前述分析结果：,若,Q,T,Q,用,R,f,1,R,n,来代替，有：,拉格朗日系数,当,D,为对角阵，分块循环矩阵为,因此：,写成标量形式为：,统计量的得到是非常困难的一件事。,其中,S,f,(u,v),S,n,(u,v),分别是,f(x,y),和,n(x,y),的功率谱。,在实际应用中，当噪声情况未知时,八、最小平方,(,最大平滑),恢复,与维纳滤波恢复的不同在于,Q,T,Q,的选择。维纳滤波中,Q,T,Q,用,R,f,1,R,n,，,根据前面的推导，其依据的准则为,均方误差最小。而,最小平方恢复依据的准则为,换言之：,在约束条件|,g,-,Hf,|,2,=|,n,|,2,条件下，恢复出来的图像的梯度最小(或称最大平滑)。根据此条件选择,Q,。,由图像增强一章可知，对图像求梯度可用如下操作,如果,f(x,y),尺寸是,AB，,则图像与,p(x,y),相卷的结果后的图像为,MN，MA+3-1， NB+3-1。,构造,p(x,y)、f(x,y),的扩展周期图像,p,e,(x,y),与,f,e,(x,y),做卷积。,* * * * *,* * * * *,* * * * *,* * * * *,* * * * *,0 -1 0,-1 4 -1,0 -1 0,如果图像的,g,、,f,采用堆叠矢量的方法构成， 有,g=,Cf,MN1,MNMN,MN1,=,对,C,进行对角化，有,E=W,-1,CW,在约束条件|,g,-,Hf,|,2,=|,n,|,2,条件下，恢复出来的图像的梯度最小(或称最大平滑)。根据此条件选择,Q,。,令,Q,C,根据,有：,左乘,W,-1,，,有,根据对角化的讨论，有,与维纳滤波不同此值为已知，不是统计量。,剩下的问题是对,s,的估测。构造残差矢量,r,=,g,-,Hf,，,调节,s,使其满足|,r,|,2,=|,n,|,2,a(a,是一个准确度系数)。,赋给,s,某个初始值；计算,f,和|,r,|,2,。,如果满足|,r,|,2,=|,n,|,2,a，,停止计算；如果|,r,|,2,|,n,|,2,+a，,减少,s，,继续计算。,点对点,计算,九、几何失真校正,在不同条件下拍摄的图像，一个物体的图像常会发生几何失真，出现歪斜变形的现象。,从太空中宇航器拍摄的地球上的等距平行线，其图像会变为歪斜或不等距；,用光学和电子扫描仪摄取的图像常会有桶形失真和枕形失真；,用普通的光学摄影与测试雷达拍摄的同一地区的景物二者在几何形状上有较大的差异。,以一幅图像为基准，去校正另一种方式摄入的图像，以期校正其几何失真，就叫做图像的几何失真复原或者几何失真校正。,1、空间变换,几何基准图像的坐标系统用,f(x, y),来表示,需要校正的图像的坐标系统用,g(x, y,),表示,设,两个图像坐标系统之间的关系,用解析式表示,x,=s,(x, y) y,=t,(x,y),通常,s(x,y),和,t(x,y),用多项式来表示：,通常用线性失真来近似较小的几何失真,x,=a,0,+a,1,x+a,2,y y,=b,0,+b,1,x+b,2,y,更精确一些可以用二次型来近似,x,=a,0,+a,1,x+a,2,y+a,3,x,2,+a,4,xy+a,5,y,2,y,=b,0,+b,1,x+b,2,y+b,3,x,2,+b,4,xy+b,5,y,2,若基准图像为,f(x,y)，,失真图像为,g(x,y,)，,对于景物上的同一个点，,假定其灰度不变,，则,f(x, y)=g(x, y,)，,可利用已知点的对应点的坐标构造方程组求取,a,i,、,b,j,。,几何失真复原的一套方法也可以用于使图像失真的工作中：在广告制作和计算机动画中常常要使物体变形。,假设的变形关系,2、灰度插值,有时会遇到知道两点灰度值，要计算两点之间点的灰度值问题，如加大图像尺寸等。,两点之间点的灰度值问题,灰度差值问题。,双线性差值法,g(E)=(,x,E,-,x,A,)g(B)-g(A)+g(A),g(F)= (,x,F,-,x,C,)g(D)-g(C)+g(C),g(H)=(,y,H,-,y,A,)g(C)-g(A)+g(A),g(I) = (,y,I,-,y,B,)g(D)-g(B)+g(B),x,33,77,A,B,C,D,E,F,y,I,H,0,十、图像质量的优劣,的,客观指标,图像质量的优劣既可以通过人眼主观视觉效果来判断，也可以通过客观指标来衡量。,1、均方误差(,MSE)：,2、,峰值信噪比(,PSNR)：,其中,M 、N,分别是,x,方向、,y,方向图像像素点的个数，,f(i,j),和,f,(i,j),分别是原始图像和测试图像在(,i, j),点上的取值，,L,是图像中灰度取值的范围，对8比特的灰度图像而言，,L=255。,3、模糊信噪比(,BSNR, the Blurred Signal-to-Noise Ratio,),表示由模糊和叠加噪声引起的降质程度。,ISNR,只是评价图像恢复算法好坏的一个客观指标，,ISNR,高并不一定主观视觉效果好。,4、,信噪比的改善度量,ISNR(the Improvement in SNR),十一、图象的重构从投影数据重构图象,在,我们观察的图象中，有一部分是通过数据人为构造的图象，如：,CT,图象、地址断面图象、安全检查设备输出的图象等。其中一类是由投影信息重构的图象本节的重点。,1917,年，奥地利数学家,J.Radon,证明了二维或三维物体可以从许多投影来重构其内部的数据。1963年，美国科学家,A.M.,Cormack,首先将这一理论用于医学图象重构；1972年，英国科学家,G.N.,Homsfield,设计出第一台,X,射线扫描仪,Computer,Tomgraph,CT。,二人于1979年双双获诺贝尔医学奖。,1、投影数据,平行光按传输方向投射到一个物体后，在垂直于平行光传输方向的平面上生成影象称为投影；用某种传感器取得的投影影象的数据称为投影数据。,不同物质组合体的视图,自然光投影视图,X,射线投影视图,片状,X,射线投影视图,最后一幅图象的数据就是重构断面图象的投影数据。显然，仅凭这组数据不能得到断面图象(,里面物体的方园无法判断,)，,我们需要更多的数据。,2、物质对,X,射线的吸收,研究表明：当强度为,I,0,的,X,射线通过吸收率为,(x),的物体时，有下面的关系：,(x),x,a,b,I,0,I,I,I,0,0,计算机采样示意图,基准检测器,检测器,物体,+,I,0,lnI,0,lnI,+,-,A/D,3、由投影数据重构图像简单方法,对应一个吸收率不均匀物质，,(x,y),不为常数。可否根据得到的投影值求得,(x,y)？,结论：如果,(x,y),已知，既可以重构图象。,方法一：解联立方程法,设每个网格中的吸收系数为,x,1,x,N,，,第,i,条射线与第,x,j,个象素的相交长度为,a,ij,(,可以计算得出)，代表第,j,个像素沿第,i,条射线的贡献的权值。,如果用,p,i,(,k,),表示沿射线方向在射线角度为,k,时,的总吸收测量值，则可通过解方程组的方法得到,x,1,x,N,。,长度,a,ij,可以计算出来，只要找出63个独立方程即可。,1,第,i,条射线,第,j,个,像素,接收器,2,3,N,55,方法二：叠代法,被照射物质内的某一确切位置的吸收函数值是固定的。构造叠代公式进行叠代计算，其收敛值应是吸收函数值。,上次求得的,k,(i,j),值,某,角度为,时的投影数据值,某,角度为,时射线经过的像元数之和,同一投影线路上的上次,k,(i,j),值之和,2.5+(3-5)/2=1.5,2.5+(7-5)/2=3.5,给定叠代结束条件：,1,1,1,1,2.5,2.5,2.5,2.5,3.5,3.5,1.5,1.5,3.75,2.75,2.25,1.25,3.5,3,2,1.5,4,2.75,2.25,1,4,3,2,1,3,7,2,5,3,6,4,1,5,4,给定初值,2.5,4、投影定理和傅立叶重构,x,y,s,s,t,p(s,),f(x,y),X,射线方向,1)、构造两个坐标系，相交成,角。根据解析几何的知识有,根据投影数据重构图像的理论进行图像重构。,2)、寻找,p(s,),与,f(s,t),之间的关系,根据,X,射线与物体吸收率的关系有：,3)、求投影,p(s,),在某一,时对,s,的傅立叶变换：,其中，,T=0。,可以写成,F(S,T)=F(S,0),角度,固定后，,p(s,),的傅立叶变换,s,轴上各点的傅立叶变换。,x,y,s,s,t,p(s,),f(x,y),X,射线方向,L,沿,L,方向的微分,4)、寻找,f(x,y),与,p(s,)、f(s,t),之间的关系,由于吸收值与座标系统无关(仅差坐标变换系数,)，,有：,对,f(x,y),和,f(s,t),的傅立叶变换为：,将(1)代入(2)有,(雅可比行列式1，,dxdy,=,dsdt,)：,雅可比行列式,观察傅立叶变换公式，与(3)比较，恰好满足旋转公式,这个变换是频域的一个旋转变换，于是有,投影定理(以前图所示)：,某角度投影数据对,s,的傅立叶变换,s,轴上对应点的傅立叶变换。,对比,5、傅立叶重建原理,如果已知,F(X,Y)，,则,该式即为重构图象的生成公式。,雅可比行列式,实际应用中往往采用下面的方法,公式,F(S,T)= P(S,)=F(S,0),T=0,告诉我们：通过采样数据得到的傅立叶变换值，实际上是,x-y,座标系统旋转,角度后在,s-t,座标系统中,s,轴上数据的傅立叶变换，如下图所示。如果要得到真正图象的傅立叶变换，需对下图做适当的修改。,X,Y,Y,X,一般用线性内插公式计算(,X,Y),座标点上的值,。,式中：,C,为直角坐标定位点；,P,1,，P,2,，P,3，,P,4,为与,C,点距离最近之四点之值；,d,1,，d,2,，d,3，,d,4,为,C,点与,P,1,，P,2,，P,3,，P,4,四点的距离。对于直角坐标定位点在旋转座标轴上之点，可采用零级内插公式得到：,C,X,Y,得到频域图象后，对其进行反变换，即可得到,f(x,y)。,最近点的傅立叶变换值,6、卷积反投影重建原理,对,空域图象而言，有,空域、频域座标对应关系为,x,y,s,t,f(x,y),X,Y,S,T,F(X,Y),r,根据傅立叶重建公式及座标变换式，有,P(S,),乘以|,S|,等效为两个时域函数的卷积。,P(S,),的傅立叶反变换是投影数据,p(s,) ；,而其余部分可设为,h(s)，,即,于是,7、核函数是什么样子,根据|,S|= (X,2,+Y,2,),1/2,，,可知,X,Y,|S|,s,h(s),h(s),s,h(s),s,X,|S|,|S|,X,|S|,X,|S|,X,h(s),s,