图像的恢复与重构课件

第五章图像的恢复与重构

什么是图像退化：图像的质量变坏叫做退化。退化的形式有图像模糊、图像有干扰等。图像退化的处理方法：无论是由光学、光电或电子方法获得的图像都会有不同程度的退化；退化的形式多种多样。如传感器噪声、摄像机未聚焦、物体与摄像设备之间的相对移动、随机大气湍流、光学系统的相差、成像光源或射线的散射等；如果我们对退化的类型、机制和过程都十分清楚，那么就可以利用其反过程来复原图像。典型的图像复原方法是根据图像退化的先验知识建立一个退化模型，以此模型为基础，采用滤波等手段进行处理，使得复原后的图像符合一定的准则，达到改善图像质量的目的。第五章图像的恢复与重构什么是图像退化：图像的质量变1一、图像退化模型f(i,j)：原始图像g(i,j)：降质图像T(·)：成像系统的作用，则：g(x,y)=T[f(x,y)]设T是线性移不变的。一幅连续的图像f(x,y)可以用抽样函数的二维卷积表示：因此，令h(x,α;y,β)=T[δ(x-α,y-β)]，则有：定义于不在原点的二维δ函数由于f(α,β)与x，y没有关系一、图像退化模型f(i,j)：原始图像因此，定义于不在原点2称h(x,α;y,β)为点扩散函数(PSF)或系统冲击响应。多数情况下它表现为时不变的，反映在图像中为位移不变的，则h(x,α;y,β)可以表示为h(x-α,y-β)其中*表示卷积运算。如果T(·)是一个h可分离系统，即：h(x,α;y,β)=h1(x,α)h2(y,β)则二维运算可以分解为列和行两次一维运算来代替。在加性噪声情况下，图像退化模型可以表示为：g(x,y)=f(x,y)*h(x,y)+n(n,y)其中n(x,y)为噪声图像。

称h(x,α;y,β)为点扩散函数(PSF)或系统冲击响3二、离散图像退化模型对于图像降质过程进行数学建模，设：f(i,j)为原始图像；y(i,j)为降质图像；h(i,j;k,l)为点扩散函数；图像为M×N维。有假设为空间不移变h(i,j;k,l)，则：二、离散图像退化模型对于图像降质过程进行数学建模，设：假设为4线性位移不变的图像退化模型则表示为：

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)+线性位移不变的图像退化模型则表示为：结论：如果已知g(x,y5三、循环矩阵及傅立叶化一个一维离散序列通过一个系统发生失真的过程可用下图表示f=H-1[g-n]用矩阵表示，可以写成

g(x)f(x)n(x)h(x)+如果考虑加性噪声，根据离散序列的卷积定理，有扩展为周期为M的序列三、循环矩阵及傅立叶化一个一维离散序列通过一个系统发生失真的6由于离散卷积的周期性，有he(x)=he(x+M)，H可以写成H是一个循环阵。结论：离散卷积都可以写成：输入矩阵×循环矩阵！对数字图象的二维离散函数也是如此。对图像退化模型而言，有A=5B=5M=9M=9由于离散卷积的周期性，有he(x)=he(x+M)，H可以写7用矩阵形式表示上式：g=Hf+n

g、f和n分别表示M×N的函数矩阵ge(i,j)、fe(i,j)和ne(i,j)的各行前后相连而成的列矢量(堆叠矢量)。如果假设原始图像是M×N维矩阵，则H是MN×MN循环矩阵，且H是一个分块(M×M个)循环矩阵：每一个子矩阵Hi自身也是循环矩阵N×N：=+MN×1MN×MNMN×1MN×1用矩阵形式表示上式：g=Hf+n每一个子矩阵Hi自身也是循81、一维信号序列循环矩阵的对角化和傅立叶化解矩阵方程：f=H-1[g-n]最简单的计算方法就是对角化，H→H-1也是对角阵。对角化H的方法——求取其特征值和特征矢量。对循环矩阵而言，设：其有M个特征值和特征矢量。由于w(k)是由傅立叶系数构成的，因此w(k)彼此是正交的。所以，由w(k)构成的变换矩阵是可逆的。1、一维信号序列循环矩阵的对角化和傅立叶化解矩阵方程：f=9用特征矢量组成的矩阵：W=[w(0),w(1)…w(M-1)]生成对角矩阵D：D=W-1HW；且D(k,k)=λ(k)。而根据周期性，M–i=-i凑一个常数所以H=WDW-1用特征矢量组成的矩阵：W=[w(0),w(1)…w(M-1)10对g=Hf+n而言，可以写成g=WDW-1f+n，有W－1g=DW-1f+W－1n；其中对

W-1g列矢量的每一行G(k)而言，有对W-1f列矢量的每一行F(k)而言对W-1n而言有同样的结果。所以对W－1g=DW-1f+W－1n而言，G(u)=MH(u)F(u)+N(u)。上面的过程称之为循环矩阵的傅立叶化。对g=Hf+n而言，可以写成g=WDW-1f+n，有W－111如果图像的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、二维信号序列块循环矩阵的对角化和傅立叶化=+如果图像的g、f、n采用堆叠矢量的方法构成，g=Hf+n。123、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)有固定的或近似的数学模型。3、H(u,v)的获取要知道一个图象降质系统13四、常见的线性移不变降质算子运动模糊：通常在拍摄过程中，相机或物体移动造成的运动模糊可以用一维均匀邻域像素灰度的平均值来表示：

大气扰动模糊：这种模糊经常出现在遥感和航空摄影中，由于曝光时间过长引起的模糊可用高斯点扩散函数来表示：式中K是一个归一化常数，保证模糊的大小为单位值，σ2可以决定模糊的程度。四、常见的线性移不变降质算子运动模糊：通常在拍摄过程中，相机14均匀不聚焦模糊这是由于相机聚焦不准确引起的，虽然不聚焦由许多参数决定，如相机的焦距、相机光圈的大小、形状、物体和相机之间的距离等，但在研究中为了简单起见，我们用下列函数表示聚焦不准引起的模糊：均匀二维模糊这是最常见的一种模糊，可以用来近似聚焦不准引起的模糊：其中L是奇数。均匀不聚焦模糊15五、无约束恢复逆滤波对于图像退化模型：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)。五、无约束恢复逆滤波16什么是无约束恢复当对噪声一无所知时，使n无约束的小。由于g=Hf+n，假设通过恢复可以得到一个不错的f的估计f’。显然，f’应满足关系g－Hf’＝n，我们希望n尽可能的小。于是问题转化为f’在什么情况下n最小——对矩阵而言就是它的迹(对角线之和)的平方最小。什么是无约束恢复当对噪声一无所知时，使n无约束的小。17实际情况中，噪声是不可避免的，因而只能求F(u,v)的估计值：

如果H(u,v)有许多零点，必然使得复原的结果受到极大影响。或者如果H(u,v)不为零但是有非常小的值，也即病态条件，也会使复原效果受到影响。解决这个问题的方法是避开H(u,v)的零点。幸好一般的H(u,v)在低频附近的有限区域内不为零。因此逆滤波可以在原点附近进行，相当于在频域乘上一低通窗口函数W(u,v)。

实际情况中，噪声是不可避免的，因而只能求F(u,v)的估计值18为了防止随着u、v的增大H(u,v)的迅速减小而增设一些条件由于截断地原因，被恢复的图象振铃较大。一种改进的方法是取(?)其中，d、k均为小于1的常数。逆滤波的应用条件：退化图像g(x,y)是信噪比较高的图像。设：为了防止随着u、v的增大H(u,v)的迅速减小而增设一些条件19对降质图像g=Hf+n而言，所谓恢复可以看成是对g通过某种线性变换L而得到原图像的估计值f’。其中，L是H的某种逆过程。如果f’与f之间的差满足某些条件——约束条件，例如：满足给定的均方误差，则称该图像被恢复了。六、有约束恢复(某种准则)的一般公式噪声是独立的，与图像不相关。求其期望应为0。期望值越小，其迹也越小对降质图像g=Hf+n而言，所谓恢复可以看成是对g通过某种线20因此对求迹而言，式中第2项和第3项的相等，有问题转化为：L等于什么时，均方误差最小。由此可得最佳L：0000将Rf，Rn看成是一个循环块矩阵，对其进行对角化，有：图像的相关度约为20－30个象素合理的假设因此对求迹而言，式中第2项和第3项的相等，有问题转化为：L等21由最佳L：可得两边乘以W和W-1：又因为当D为对角阵时有所以当L为分块循环矩阵时，WLW-1也是一对角阵。其对角线上元素的值为循环序列(“主值”序列)的傅立叶变换值。同理，D，D*，A，B都为对角阵，其对角线上的元素的值为对应循环序列(“主值”序列)的傅立叶变换值。于是上式可以写成标量形式合理的假设f’(x,y)=g(x,y)*l(x,y)由最佳L：可得两边乘以W和W-1：又因为当D为对角阵时有所以22七、维纳滤波恢复基本原理研究发现，逆滤波复原方法对噪声极为敏感，要求信噪比较高(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’。七、维纳滤波恢复基本原理23用拉格朗日法求解设Rf和Rn为f和n的相关矩阵，根据前述分析结果：若QTQ用Rf－1Rn来代替，有：拉格朗日系数用拉格朗日法求解设Rf和Rn为f和n的相关矩阵，根据前述分析24当D为对角阵，分块循环矩阵为因此：写成标量形式为：统计量的得到是非常困难的一件事。其中Sf(u,v),Sn(u,v)分别是f(x,y)和n(x,y)的功率谱。当D为对角阵，分块循环矩阵为因此：统计量的得到是非常困难的一25在实际应用中，当噪声情况未知时八、最小平方(最大平滑)恢复与维纳滤波恢复的不同在于QTQ的选择。维纳滤波中QTQ用Rf－1Rn，根据前面的推导，其依据的准则为均方误差最小。而最小平方恢复依据的准则为换言之：在约束条件||g-Hf’||2=||n||2条件下，恢复出来的图像的梯度最小(或称最大平滑)。根据此条件选择Q。在实际应用中，当噪声情况未知时八、最小平方(最大平滑)恢复换26由图像增强一章可知，对图像求梯度可用如下操作如果f(x,y)尺寸是A×B，则图像与p(x,y)相卷的结果后的图像为M×N，M≥A+3-1，N≥B+3-1。构造p(x,y)、f(x,y)的扩展周期图像pe(x,y)与fe(x,y)做卷积。*************************0-10-14-10-10由图像增强一章可知，对图像求梯度可用如下操作如果f(x,y)27如果图像的g、f采用堆叠矢量的方法构成，有g=CfMN×1MN×MNMN×1=对C进行对角化，有E=W-1CW在约束条件||g-Hf’||2=||n||2条件下，恢复出来的图像的梯度最小(或称最大平滑)。根据此条件选择Q。如果图像的g、f采用堆叠矢量的方法构成，有g=CfMN×128令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，继续计算；如果||r||2>||n||2+a，减少s，继续计算。点对点计算令Q＝C根据29九、几何失真校正在不同条件下拍摄的图像，一个物体的图像常会发生几何失真，出现歪斜变形的现象。从太空中宇航器拍摄的地球上的等距平行线，其图像会变为歪斜或不等距；用光学和电子扫描仪摄取的图像常会有桶形失真和枕形失真；用普通的光学摄影与测试雷达拍摄的同一地区的景物二者在几何形状上有较大的差异。以一幅图像为基准，去校正另一种方式摄入的图像，以期校正其几何失真，就叫做图像的几何失真复原或者几何失真校正。九、几何失真校正在不同条件下拍摄的图像，一个物体的图像常会发301、空间变换几何基准图像的坐标系统用f(x,y)来表示需要校正的图像的坐标系统用g(x’,y’)表示设两个图像坐标系统之间的关系用解析式表示x’=s(x,y)y’=t(x,y)通常s(x,y)和t(x,y)用多项式来表示：1、空间变换几何基准图像的坐标系统用f(x,y)来表示31

通常用线性失真来近似较小的几何失真x’=a0+a1x+a2yy’=b0+b1x+b2y更精确一些可以用二次型来近似x’=a0+a1x+a2y+a3x2+a4xy+a5y2y’=b0+b1x+b2y+b3x2+b4xy+b5y2若基准图像为f(x,y)，失真图像为g(x’,y’)，对于景物上的同一个点，假定其灰度不变，则f(x,y)=g(x’,y’)，可利用已知点的对应点的坐标构造方程组求取ai、bj。几何失真复原的一套方法也可以用于使图像失真的工作中：在广告制作和计算机动画中常常要使物体变形。假设的变形关系通常用线性失真来近似较小的几何失真假设的变形关系322、灰度插值有时会遇到知道两点灰度值，要计算两点之间点的灰度值问题，如加大图像尺寸等。两点之间点的灰度值问题——灰度差值问题。双线性差值法g(E)=(xE-xA)[g(B)-g(A)]+g(A)g(F)=(xF-xC)[g(D)-g(C)]+g(C)g(H)=(yH-yA)[g(C)-g(A)]+g(A)g(I)=(yI-yB)[g(D)-g(B)]+g(B)x3×37×7ABCDEFyIH02、灰度插值有时会遇到知道两点灰度值，要计算两点之间点的灰度33十、图像质量的优劣的客观指标图像质量的优劣既可以通过人眼主观视觉效果来判断，也可以通过客观指标来衡量。1、均方误差(MSE)：2、峰值信噪比(PSNR)：其中M、N分别是x方向、y方向图像像素点的个数，f(i,j)和f’(i,j)分别是原始图像和测试图像在(i,j)点上的取值，L是图像中灰度取值的范围，对8比特的灰度图像而言，L=255。十、图像质量的优劣的客观指标图像质量的优劣既可以通过人眼主观343、模糊信噪比(BSNR,theBlurredSignal-to-NoiseRatio)

表示由模糊和叠加噪声引起的降质程度。3、模糊信噪比(BSNR,theBlurredSig35

ISNR只是评价图像恢复算法好坏的一个客观指标，ISNR高并不一定主观视觉效果好。4、信噪比的改善度量ISNR(theImprovementinSNR)ISNR只是评价图像恢复算法好坏的一个客观指标，ISNR36十一、图象的重构——从投影数据重构图象在我们观察的图象中，有一部分是通过数据人为构造的图象，如：CT图象、地址断面图象、安全检查设备输出的图象等。其中一类是由投影信息重构的图象——本节的重点。1917年，奥地利数学家J.Radon证明了二维或三维物体可以从许多投影来重构其内部的数据。1963年，美国科学家A.M.Cormack首先将这一理论用于医学图象重构；1972年，英国科学家G.N.Homsfield设计出第一台X射线扫描仪ComputerTomgraph——CT。二人于1979年双双获诺贝尔医学奖。十一、图象的重构——从投影数据重构图象在我们观察的图象中，有371、投影数据平行光按传输方向投射到一个物体后，在垂直于平行光传输方向的平面上生成影象称为投影；用某种传感器取得的投影影象的数据称为投影数据。不同物质组合体的视图自然光投影视图X射线投影视图片状X射线投影视图最后一幅图象的数据就是重构断面图象的投影数据。显然，仅凭这组数据不能得到断面图象(里面物体的方园无法判断)，我们需要更多的数据。1、投影数据平行光按传输方向投射到一个物体后，在垂直于平行光382、物质对X射线的吸收研究表明：当强度为I0的X射线通过吸收率为μ(x)的物体时，有下面的关系：μ(x)xabI0III00计算机采样示意图基准检测器检测器物体+I0lnI0lnI+-A/D2、物质对X射线的吸收研究表明：当强度为I0的X射线通过吸收393、由投影数据重构图像——简单方法对应一个吸收率不均匀物质，μ(x,y)不为常数。可否根据得到的投影值求得μ(x,y)？——结论：如果μ(x,y)已知，既可以重构图象。方法一：解联立方程法设每个网格中的吸收系数为x1～

xN，第i条射线与第xj个象素的相交长度为aij(可以计算得出)，代表第j个像素沿第i条射线的贡献的权值。如果用pi(θk)表示沿射线方向在射线角度为θk时的总吸收测量值，则可通过解方程组的方法得到x1～

xN。长度aij可以计算出来，只要找出63个独立方程即可。1第i条射线第j个像素接收器23N553、由投影数据重构图像——简单方法对应一个吸收率不均匀物质，40方法二：叠代法被照射物质内的某一确切位置的吸收函数值是固定的。构造叠代公式进行叠代计算，其收敛值应是吸收函数值。上次求得的μk(i,j)值某角度为θ时的投影数据值某角度为θ时射线经过的像元数之和同一投影线路上的上次μk(i,j)值之和2.5+(3-5)/2=1.52.5+(7-5)/2=3.5给定叠代结束条件：μ1μ1μ1μ12.52.52.52.53.53.51.51.53.752.752.251.253.5321.542.752.25143213725364154给定初值2.5方法二：叠代法被照射物质内的某一确切位置的吸收函数值是固定的414、投影定理和傅立叶重构xysstp(s,θ)θf(x,y)X射线方向1)、构造两个坐标系，相交成θ角。根据解析几何的知识有根据投影数据重构图像的理论进行图像重构。4、投影定理和傅立叶重构xysstp(s,θ)θf(x,y)422)、寻找p(s,θ)与f(s,t)之间的关系根据X射线与物体吸收率的关系有：3)、求投影p(s,θ)在某一θ时对s的傅立叶变换：其中，T=0。可以写成F(S,T)=F(S,0)角度θ固定后，p(s,θ)的傅立叶变换＝s轴上各点的傅立叶变换。xysstp(s,θ)θf(x,y)X射线方向L沿L方向的微分2)、寻找p(s,θ)与f(s,t)之间的关系根据X射线与物434)、寻找f(x,y)与p(s,θ)、f(s,t)之间的关系由于吸收值与座标系统无关(仅差坐标变换系数)，有：对f(x,y)和f(s,t)的傅立叶变换为：将(1)代入(2)有(雅可比行列式＝1，dxdy=dsdt)：雅可比行列式4)、寻找f(x,y)与p(s,θ)、f(s,t)之间的关44观察傅立叶变换公式，与(3)比较，恰好满足旋转公式这个变换是频域的一个旋转变换，于是有投影定理(以前图所示)：某角度投影数据对s的傅立叶变换＝s轴上对应点的傅立叶变换。对比观察傅立叶变换公式，与(3)比较，恰好满足旋转公式这个变换是455、傅立叶重建原理如果已知F(X,Y)，则该式即为重构图象的生成公式。雅可比行列式5、傅立叶重建原理如果已知F(X,Y)，则该式即为重构图象的46实际应用中往往采用下面的方法公式F(S,T)=P(S,θ)=F(S,0)T=0告诉我们：通过采样数据得到的傅立叶变换值，实际上是x-y座标系统旋转θ角度后在s-t座标系统中s轴上数据的傅立叶变换，如下图所示。如果要得到真正图象的傅立叶变换，需对下图做适当的修改。XYYX实际应用中往往采用下面的方法公式F(S,T)=P(S,θ)47一般用线性内插公式计算(X,Y)座标点上的值。式中：C为直角坐标定位点；P1，P2，P3，P4为与C点距离最近之四点之值；d1，d2，d3，d4为C点与P1，P2，P3，P4四点的距离。对于直角坐标定位点在旋转座标轴上之点，可采用零级内插公式得到：CXY得到频域图象后，对其进行反变换，即可得到f(x,y)。最近点的傅立叶变换值一般用线性内插公式计算(X,Y)座标点上的值。式中：C为直角486、卷积反投影重建原理对空域图象而言，有空域、频域座标对应关系为xystθf(x,y)φXYSTθF(X,Y)r6、卷积反投影重建原理对空域图象而言，有空域、频域座标对应关49根据傅立叶重建公式及座标变换式，有P(S,θ)乘以|S|等效为两个时域函数的卷积。P(S,θ)的傅立叶反变换是投影数据p(s,θ)；而其余部分可设为h(s)，即于是根据傅立叶重建公式及座标变换式，有P(S,θ)乘以|S|等效507、核函数是什么样子根据|S|=(X2+Y2)1/2，可知XY|S|sh(s)h(s)sh(s)sX|S||S|X|S|X|S|Xh(s)s7、核函数是什么样子根据|S|=(X2+Y2)1/2，可知51第五章图像的恢复与重构

什么是图像退化：图像的质量变坏叫做退化。退化的形式有图像模糊、图像有干扰等。图像退化的处理方法：无论是由光学、光电或电子方法获得的图像都会有不同程度的退化；退化的形式多种多样。如传感器噪声、摄像机未聚焦、物体与摄像设备之间的相对移动、随机大气湍流、光学系统的相差、成像光源或射线的散射等；如果我们对退化的类型、机制和过程都十分清楚，那么就可以利用其反过程来复原图像。典型的图像复原方法是根据图像退化的先验知识建立一个退化模型，以此模型为基础，采用滤波等手段进行处理，使得复原后的图像符合一定的准则，达到改善图像质量的目的。第五章图像的恢复与重构什么是图像退化：图像的质量变52一、图像退化模型f(i,j)：原始图像g(i,j)：降质图像T(·)：成像系统的作用，则：g(x,y)=T[f(x,y)]设T是线性移不变的。一幅连续的图像f(x,y)可以用抽样函数的二维卷积表示：因此，令h(x,α;y,β)=T[δ(x-α,y-β)]，则有：定义于不在原点的二维δ函数由于f(α,β)与x，y没有关系一、图像退化模型f(i,j)：原始图像因此，定义于不在原点53称h(x,α;y,β)为点扩散函数(PSF)或系统冲击响应。多数情况下它表现为时不变的，反映在图像中为位移不变的，则h(x,α;y,β)可以表示为h(x-α,y-β)其中*表示卷积运算。如果T(·)是一个h可分离系统，即：h(x,α;y,β)=h1(x,α)h2(y,β)则二维运算可以分解为列和行两次一维运算来代替。在加性噪声情况下，图像退化模型可以表示为：g(x,y)=f(x,y)*h(x,y)+n(n,y)其中n(x,y)为噪声图像。

称h(x,α;y,β)为点扩散函数(PSF)或系统冲击响54二、离散图像退化模型对于图像降质过程进行数学建模，设：f(i,j)为原始图像；y(i,j)为降质图像；h(i,j;k,l)为点扩散函数；图像为M×N维。有假设为空间不移变h(i,j;k,l)，则：二、离散图像退化模型对于图像降质过程进行数学建模，设：假设为55线性位移不变的图像退化模型则表示为：

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)+线性位移不变的图像退化模型则表示为：结论：如果已知g(x,y56三、循环矩阵及傅立叶化一个一维离散序列通过一个系统发生失真的过程可用下图表示f=H-1[g-n]用矩阵表示，可以写成

g(x)f(x)n(x)h(x)+如果考虑加性噪声，根据离散序列的卷积定理，有扩展为周期为M的序列三、循环矩阵及傅立叶化一个一维离散序列通过一个系统发生失真的57由于离散卷积的周期性，有he(x)=he(x+M)，H可以写成H是一个循环阵。结论：离散卷积都可以写成：输入矩阵×循环矩阵！对数字图象的二维离散函数也是如此。对图像退化模型而言，有A=5B=5M=9M=9由于离散卷积的周期性，有he(x)=he(x+M)，H可以写58用矩阵形式表示上式：g=Hf+n

g、f和n分别表示M×N的函数矩阵ge(i,j)、fe(i,j)和ne(i,j)的各行前后相连而成的列矢量(堆叠矢量)。如果假设原始图像是M×N维矩阵，则H是MN×MN循环矩阵，且H是一个分块(M×M个)循环矩阵：每一个子矩阵Hi自身也是循环矩阵N×N：=+MN×1MN×MNMN×1MN×1用矩阵形式表示上式：g=Hf+n每一个子矩阵Hi自身也是循591、一维信号序列循环矩阵的对角化和傅立叶化解矩阵方程：f=H-1[g-n]最简单的计算方法就是对角化，H→H-1也是对角阵。对角化H的方法——求取其特征值和特征矢量。对循环矩阵而言，设：其有M个特征值和特征矢量。由于w(k)是由傅立叶系数构成的，因此w(k)彼此是正交的。所以，由w(k)构成的变换矩阵是可逆的。1、一维信号序列循环矩阵的对角化和傅立叶化解矩阵方程：f=60用特征矢量组成的矩阵：W=[w(0),w(1)…w(M-1)]生成对角矩阵D：D=W-1HW；且D(k,k)=λ(k)。而根据周期性，M–i=-i凑一个常数所以H=WDW-1用特征矢量组成的矩阵：W=[w(0),w(1)…w(M-1)61对g=Hf+n而言，可以写成g=WDW-1f+n，有W－1g=DW-1f+W－1n；其中对

W-1g列矢量的每一行G(k)而言，有对W-1f列矢量的每一行F(k)而言对W-1n而言有同样的结果。所以对W－1g=DW-1f+W－1n而言，G(u)=MH(u)F(u)+N(u)。上面的过程称之为循环矩阵的傅立叶化。对g=Hf+n而言，可以写成g=WDW-1f+n，有W－162如果图像的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、二维信号序列块循环矩阵的对角化和傅立叶化=+如果图像的g、f、n采用堆叠矢量的方法构成，g=Hf+n。633、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)有固定的或近似的数学模型。3、H(u,v)的获取要知道一个图象降质系统64四、常见的线性移不变降质算子运动模糊：通常在拍摄过程中，相机或物体移动造成的运动模糊可以用一维均匀邻域像素灰度的平均值来表示：

大气扰动模糊：这种模糊经常出现在遥感和航空摄影中，由于曝光时间过长引起的模糊可用高斯点扩散函数来表示：式中K是一个归一化常数，保证模糊的大小为单位值，σ2可以决定模糊的程度。四、常见的线性移不变降质算子运动模糊：通常在拍摄过程中，相机65均匀不聚焦模糊这是由于相机聚焦不准确引起的，虽然不聚焦由许多参数决定，如相机的焦距、相机光圈的大小、形状、物体和相机之间的距离等，但在研究中为了简单起见，我们用下列函数表示聚焦不准引起的模糊：均匀二维模糊这是最常见的一种模糊，可以用来近似聚焦不准引起的模糊：其中L是奇数。均匀不聚焦模糊66五、无约束恢复逆滤波对于图像退化模型：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)。五、无约束恢复逆滤波67什么是无约束恢复当对噪声一无所知时，使n无约束的小。由于g=Hf+n，假设通过恢复可以得到一个不错的f的估计f’。显然，f’应满足关系g－Hf’＝n，我们希望n尽可能的小。于是问题转化为f’在什么情况下n最小——对矩阵而言就是它的迹(对角线之和)的平方最小。什么是无约束恢复当对噪声一无所知时，使n无约束的小。68实际情况中，噪声是不可避免的，因而只能求F(u,v)的估计值：

如果H(u,v)有许多零点，必然使得复原的结果受到极大影响。或者如果H(u,v)不为零但是有非常小的值，也即病态条件，也会使复原效果受到影响。解决这个问题的方法是避开H(u,v)的零点。幸好一般的H(u,v)在低频附近的有限区域内不为零。因此逆滤波可以在原点附近进行，相当于在频域乘上一低通窗口函数W(u,v)。

实际情况中，噪声是不可避免的，因而只能求F(u,v)的估计值69为了防止随着u、v的增大H(u,v)的迅速减小而增设一些条件由于截断地原因，被恢复的图象振铃较大。一种改进的方法是取(?)其中，d、k均为小于1的常数。逆滤波的应用条件：退化图像g(x,y)是信噪比较高的图像。设：为了防止随着u、v的增大H(u,v)的迅速减小而增设一些条件70对降质图像g=Hf+n而言，所谓恢复可以看成是对g通过某种线性变换L而得到原图像的估计值f’。其中，L是H的某种逆过程。如果f’与f之间的差满足某些条件——约束条件，例如：满足给定的均方误差，则称该图像被恢复了。六、有约束恢复(某种准则)的一般公式噪声是独立的，与图像不相关。求其期望应为0。期望值越小，其迹也越小对降质图像g=Hf+n而言，所谓恢复可以看成是对g通过某种线71因此对求迹而言，式中第2项和第3项的相等，有问题转化为：L等于什么时，均方误差最小。由此可得最佳L：0000将Rf，Rn看成是一个循环块矩阵，对其进行对角化，有：图像的相关度约为20－30个象素合理的假设因此对求迹而言，式中第2项和第3项的相等，有问题转化为：L等72由最佳L：可得两边乘以W和W-1：又因为当D为对角阵时有所以当L为分块循环矩阵时，WLW-1也是一对角阵。其对角线上元素的值为循环序列(“主值”序列)的傅立叶变换值。同理，D，D*，A，B都为对角阵，其对角线上的元素的值为对应循环序列(“主值”序列)的傅立叶变换值。于是上式可以写成标量形式合理的假设f’(x,y)=g(x,y)*l(x,y)由最佳L：可得两边乘以W和W-1：又因为当D为对角阵时有所以73七、维纳滤波恢复基本原理研究发现，逆滤波复原方法对噪声极为敏感，要求信噪比较高(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’。七、维纳滤波恢复基本原理74用拉格朗日法求解设Rf和Rn为f和n的相关矩阵，根据前述分析结果：若QTQ用Rf－1Rn来代替，有：拉格朗日系数用拉格朗日法求解设Rf和Rn为f和n的相关矩阵，根据前述分析75当D为对角阵，分块循环矩阵为因此：写成标量形式为：统计量的得到是非常困难的一件事。其中Sf(u,v),Sn(u,v)分别是f(x,y)和n(x,y)的功率谱。当D为对角阵，分块循环矩阵为因此：统计量的得到是非常困难的一76在实际应用中，当噪声情况未知时八、最小平方(最大平滑)恢复与维纳滤波恢复的不同在于QTQ的选择。维纳滤波中QTQ用Rf－1Rn，根据前面的推导，其依据的准则为均方误差最小。而最小平方恢复依据的准则为换言之：在约束条件||g-Hf’||2=||n||2条件下，恢复出来的图像的梯度最小(或称最大平滑)。根据此条件选择Q。在实际应用中，当噪声情况未知时八、最小平方(最大平滑)恢复换77由图像增强一章可知，对图像求梯度可用如下操作如果f(x,y)尺寸是A×B，则图像与p(x,y)相卷的结果后的图像为M×N，M≥A+3-1，N≥B+3-1。构造p(x,y)、f(x,y)的扩展周期图像pe(x,y)与fe(x,y)做卷积。*************************0-10-14-10-10由图像增强一章可知，对图像求梯度可用如下操作如果f(x,y)78如果图像的g、f采用堆叠矢量的方法构成，有g=CfMN×1MN×MNMN×1=对C进行对角化，有E=W-1CW在约束条件||g-Hf’||2=||n||2条件下，恢复出来的图像的梯度最小(或称最大平滑)。根据此条件选择Q。如果图像的g、f采用堆叠矢量的方法构成，有g=CfMN×179令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，继续计算；如果||r||2>||n||2+a，减少s，继续计算。点对点计算令Q＝C根据80九、几何失真校正在不同条件下拍摄的图像，一个物体的图像常会发生几何失真，出现歪斜变形的现象。从太空中宇航器拍摄的地球上的等距平行线，其图像会变为歪斜或不等距；用光学和电子扫描仪摄取的图像常会有桶形失真和枕形失真；用普通的光学摄影与测试雷达拍摄的同一地区的景物二者在几何形状上有较大的差异。以一幅图像为基准，去校正另一种方式摄入的图像，以期校正其几何失真，就叫做图像的几何失真复原或者几何失真校正。九、几何失真校正在不同条件下拍摄的图像，一个物体的图像常会发811、空间变换几何基准图像的坐标系统用f(x,y)来表示需要校正的图像的坐标系统用g(x’,y’)表示设两个图像坐标系统之间的关系用解析式表示x’=s(x,y)y’=t(x,y)通常s(x,y)和t(x,y)用多项式来表示：1、空间变换几何基准图像的坐标系统用f(x,y)来表示82

通常用线性失真来近似较小的几何失真x’=a0+a1x+a2yy’=b0+b1x+b2y更精确一些可以用二次型来近似x’=a0+a1x+a2y+a3x2+a4xy+a5y2y’=b0+b1x+b2y+b3x2+b4xy+b5y2若基准图像为f(x,y)，失真图像为g(x’,y’)，对于景物上的同一个点，假定其灰度不变，则f(x,y)=g(x’,y’)，可利用已知点的对应点的坐标构造方程组求取ai、bj。几何失真复原的一套方法也可以用于使图像失真的工作中：在广告制作和计算机动画中常常要使物体变形。假设的变形关系通常用线性失真来近似较小的几何失真假设的变形关系832、灰度插值有时会遇到知道两点灰度值，要计算两点之间点的灰度值问题，如加大图像尺寸等。两点之间点的灰度值问题——灰度差值问题。双线性差值法g(E)=(xE-xA)[g(B)-g(A)]+g(A)g(F)=(xF-xC)[g(D)-g(C)]+g(C)g(H)=(yH-yA)[g(C)-g(A)]+g(A)g(I)=(yI-yB)[g(D)-g(B)]+g(B)x3×37×7ABCDEFyIH02、灰度插值有时会遇到知道两点灰度值，要计算两点之间点的灰度84十、图像质量的优劣的客观指标图像质量的优劣既可以通过人眼主观视觉效果来判断，也可以通过客观指标来衡量。1、均方误差(MSE)：2、峰值信噪比(PSNR)：其中M、N分别是x方向、y方向图像像素点的个数，f(i,j)和f’(i,j)分别是原始图像和测试图像在(i,j)点上的取值，L是图像中灰度取值的范围，对8比特的灰度图像而言，L=255。十、图像质量的优劣的客观指标图像质量的优劣既可以通过人眼主观853、模糊信噪比(BSNR,theBlurredSignal-to-NoiseRatio)

表示由模糊和叠加噪声引起的降质程度。3、模糊信噪比(BSNR,theBlurredSig86

ISNR只是评价图像恢复算法好坏的一个客观指标，ISNR高并不一定主观视觉效果好。4、信噪比的改善度量ISNR(theImprovementinSNR)ISNR只是评价图像恢复算法好坏的一个客观指标，ISNR87十一、图象的重构——从投影数据重构图象在我们观察的图象中，有一部分是通过数据人为构造的图象，如：CT图象、地址断面图象、安全检查设备输出的图象等。其中一类是由投影信息重构的图象——本节的重点。1917年，奥地利数学家J.Radon证明了二维或三维物体可以从许多投影来重构其内部的数据。1963年，美国科学家A.M.Cormack首先将这一理论用于医学图象重构；1972年，英国科学家G.N.Homsfield设计出第一台X射线扫描仪ComputerTomgraph——CT。二人于1979年双双获诺贝尔医学奖。十一、图象的重构——从投影数据重构图象在我们观察的图象中，有881、投影数据平行光按传输方向投射到一个物体后，在垂直于平行光传输方向的平面上生成影象称为投影；用某种传感器取得的投影影象的数据称为投影数据。不同物质组合体的视图自然光投影视图X射线投影视图片状X射线投影视图最后一幅图象的数据就是重构断面图象的投影数据。显然，仅凭这组数据不能得到断面图象(里面物体的方园无法判断)，我们需要更多的数据。1、投影数据平行光按传输方向投射到一个物体后，在垂直于平行光892、物质对X射线的吸收研究表明：当强度为I0的X射线通过吸收率为μ(x)的物体时，有下面的关系：μ(x)xabI0III00计算机采样示意图基准检测器检测器物体+I0lnI0lnI+-A/D2、物质