第五章 图像的恢复与重构

1. 第五章 图像的恢复与重构 • 什么是图像退化：图像的质量变坏叫做退化。退化的形式有图像模糊、图像有干扰等。 • 图像退化的处理方法：无论是由光学、光电或电子方法获得的图像都会有不同程度的退化；退化的形式多种多样。如传感器噪声、摄像机未聚焦、物体与摄像设备之间的相对移动、随机大气湍流、光学系统的相差、成像光源或射线的散射等；如果我们对退化的类型、机制和过程都十分清楚，那么就可以利用其反过程来复原图像。 典型的图像复原方法是根据图像退化的先验知识建立一个退化模型，以此模型为基础，采用滤波等手段进行处理，使得复原后的图像符合一定的准则，达到改善图像质量的目的。

2. 一、图像退化模型 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没有关系

3. 称 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)为噪声图像。

4. 二、离散图像退化模型 对于图像降质过程进行数学建模，设： f(i, j)为原始图像；y(i, j)为降质图像；h(i, j; k, l)为点扩散函数；图像为M×N维。有 假设为空间不移变h(i, j; k, l)，则：

5. 线性位移不变的图像退化模型则表示为： g(x,y)=f(x,y)*h(x,y)+n(x,y) f(x,y) g(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)]

6. f(x) g(x) h(x) + n(x) 三、循环矩阵及傅立叶化 一个一维离散序列通过一个系统发生失真的过程可用下图表示 如果考虑加性噪声，根据离散序列的卷积定理，有 用矩阵表示，可以写成 扩展为周期为M的序列 f= H-1[g- n]

7. A=5 B=5 M=9 M=9 由于离散卷积的周期性，有he(x)=he(x+M)，H可以写成 H是一个循环阵。 结论：离散卷积都可以写成：输入矩阵×循环矩阵！ 对数字图象的二维离散函数也是如此。 对图像退化模型而言，有

8. = + MN×1 MN×MN MN×1 MN×1 用矩阵形式表示上式：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：

9. 1、一维信号序列循环矩阵的对角化和傅立叶化1、一维信号序列循环矩阵的对角化和傅立叶化 解矩阵方程： f= H-1[g- n] 最简单的计算方法就是对角化，H → H-1也是对角阵。 对角化H的方法——求取其特征值和特征矢量。 对循环矩阵而言，设：其有M个特征值和特征矢量。 由于w(k)是由傅立叶系数构成的，因此w(k)彼此是正交的。所以，由w(k)构成的变换矩阵是可逆的。

10. 用特征矢量组成的矩阵：W=[w(0),w(1)…w(M-1)] 生成对角矩阵D：D=W -1HW；且D(k,k)=λ(k)。而 凑一个常数 根据周期性， M – i = - i 所以 H=WDW -1

11. 对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)。上面的过程称之为循环矩阵的傅立叶化。

12. = + 2、二维信号序列块循环矩阵的对角化和傅立叶化 如果图像的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)]

13. 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)有固定的或近似的数学模型。

14. 四、常见的线性移不变降质算子 • 运动模糊：通常在拍摄过程中，相机或物体移动造成的运动模糊可以用一维均匀邻域像素灰度的平均值来表示： • 大气扰动模糊：这种模糊经常出现在遥感和航空摄影中，由于曝光时间过长引起的模糊可用高斯点扩散函数来表示： 式中K是一个归一化常数，保证模糊的大小为单位值，σ2可以决定模糊的程度。

15. 均匀不聚焦模糊 这是由于相机聚焦不准确引起的，虽然不聚焦由许多参数决定，如相机的焦距、相机光圈的大小、形状、物体和相机之间的距离等，但在研究中为了简单起见，我们用下列函数表示聚焦不准引起的模糊： • 均匀二维模糊 这是最常见的一种模糊，可以用来近似聚焦不准引起的模糊： 其中L是奇数。

16. 五、无约束恢复 逆滤波 对于图像退化模型： 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)。

17. 什么是无约束恢复 当对噪声一无所知时，使n无约束的小。 由于g=Hf+n，假设通过恢复可以得到一个不错的f的估计f’。显然，f’应满足关系g－Hf’＝n，我们希望n尽可能的小。于是问题转化为f’在什么情况下n最小——对矩阵而言就是它的迹(对角线之和)的平方最小。

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

19. 设： 为了防止随着u、v的增大H(u,v)的迅速减小而增设一些条件 由于截断地原因，被恢复的图象振铃较大。一种改进的方法是取(?) 其中，d、k均为小于1的常数。逆滤波的应用条件：退化图像g(x,y)是信噪比较高的图像。

20. 六、有约束恢复(某种准则)的一般公式 对降质图像g=Hf+n而言，所谓恢复可以看成是对g通过某种线性变换L而得到原图像的估计值f’。其中，L是H的某种逆过程。如果f’与f之间的差满足某些条件——约束条件，例如：满足给定的均方误差，则称该图像被恢复了。 期望值越小，其迹也越小 噪声是独立的，与图像不相关。求其期望应为0。

21. 0 0 因此 对求迹而言，式中第2项和第3项的相等，有 问题转化为：L等于什么时，均方误差最小。 合理的假设 由此可得最佳L： 将Rf，Rn看成是一个循环块矩阵，对其进行对角化，有： 图像的相关度约为20－30个象素 0 0

22. 可得 由最佳L： 两边乘以W和W-1： 又因为当D为对角阵时有 所以 合理的假设f’(x,y)=g(x,y)*l(x,y) 当L为分块循环矩阵时，WL W-1也是一对角阵。其对角线上元素的值为循环序列(“主值”序列)的傅立叶变换值。同理，D，D*，A，B都为对角阵，其对角线上的元素的值为对应循环序列(“主值”序列)的傅立叶变换值。于是上式 可以写成标量形式

23. 七、维纳滤波恢复 基本原理 研究发现，逆滤波复原方法对噪声极为敏感，要求信噪比较高(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’。

24. 用拉格朗日法求解 拉格朗日系数 设Rf和Rn为f和n的相关矩阵，根据前述分析结果： 若QTQ用Rf－1Rn来代替，有：

25. 当D为对角阵，分块循环矩阵为 因此： 写成标量形式为： 统计量的得到是非常困难的一件事。 其中Sf(u,v), Sn(u,v)分别是f(x,y)和n(x,y)的功率谱。

26. 在实际应用中，当噪声情况未知时 八、最小平方(最大平滑)恢复 与维纳滤波恢复的不同在于QTQ的选择。维纳滤波中QTQ用Rf－1Rn，根据前面的推导，其依据的准则为均方误差最小。而最小平方恢复依据的准则为 换言之：在约束条件||g-Hf’||2=||n||2条件下，恢复出来的图像的梯度最小(或称最大平滑)。根据此条件选择Q。

27. 由图像增强一章可知，对图像求梯度可用如下操作由图像增强一章可知，对图像求梯度可用如下操作 如果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 -1 0 -1 4 -1 0 -1 0

28. = MN×1 MN×MN MN×1 如果图像的g、f采用堆叠矢量的方法构成， 有g=Cf 对C进行对角化，有 E=W-1CW 在约束条件||g-Hf’||2=||n||2条件下，恢复出来的图像的梯度最小(或称最大平滑)。根据此条件选择Q。

29. 令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，继续计算。

30. 九、几何失真校正 在不同条件下拍摄的图像，一个物体的图像常会发生几何失真，出现歪斜变形的现象。 从太空中宇航器拍摄的地球上的等距平行线，其图像会变为歪斜或不等距； 用光学和电子扫描仪摄取的图像常会有桶形失真和枕形失真； 用普通的光学摄影与测试雷达拍摄的同一地区的景物二者在几何形状上有较大的差异。 以一幅图像为基准，去校正另一种方式摄入的图像，以期校正其几何失真，就叫做图像的几何失真复原或者几何失真校正。

31. 1、空间变换 几何基准图像的坐标系统用f(x, y)来表示 需要校正的图像的坐标系统用g(x’, y’)表示 设两个图像坐标系统之间的关系用解析式表示 x’=s(x, y) y’=t(x,y) 通常s(x,y)和t(x,y)用多项式来表示：

32. 通常用线性失真来近似较小的几何失真 x ’ =a0+a1x+a2y y ’ =b0+b1x+b2y 更精确一些可以用二次型来近似 x ’ =a0+a1x+a2y+a3x2+a4xy+a5y2 y ’ =b0+b1x+b2y+b3x2+b4xy+b5y2 若基准图像为f(x,y)，失真图像为g(x ’,y ’)，对于景物上的同一个点，假定其灰度不变，则f(x, y)=g(x’, y’)，可利用已知点的对应点的坐标构造方程组求取ai、bj。 几何失真复原的一套方法也可以用于使图像失真的工作中：在广告制作和计算机动画中常常要使物体变形。 假设的变形关系

33. 3×3 7×7 y C F D H I A E B 2、灰度插值 有时会遇到知道两点灰度值，要计算两点之间点的灰度值问题，如加大图像尺寸等。 两点之间点的灰度值问题 ——灰度差值问题。 双线性差值法 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) 0 x

34. 十、图像质量的优劣的客观指标 图像质量的优劣既可以通过人眼主观视觉效果来判断，也可以通过客观指标来衡量。 1、均方误差(MSE)： 2、峰值信噪比(PSNR)： 其中M 、N分别是x方向、y方向图像像素点的个数，f(i,j)和f’(i,j)分别是原始图像和测试图像在(i, j)点上的取值，L是图像中灰度取值的范围，对8比特的灰度图像而言，L=255。

35. 3、模糊信噪比(BSNR, the Blurred Signal-to-Noise Ratio )表示由模糊和叠加噪声引起的降质程度。

36. 4、信噪比的改善度量ISNR(the Improvement in SNR) ISNR只是评价图像恢复算法好坏的一个客观指标，ISNR高并不一定主观视觉效果好。

37. 十一、图象的重构——从投影数据重构图象 在我们观察的图象中，有一部分是通过数据人为构造的图象，如：CT图象、地址断面图象、安全检查设备输出的图象等。其中一类是由投影信息重构的图象——本节的重点。 1917年，奥地利数学家J.Radon证明了二维或三维物体可以从许多投影来重构其内部的数据。1963年，美国科学家A.M.Cormack首先将这一理论用于医学图象重构；1972年，英国科学家G.N.Homsfield设计出第一台X射线扫描仪Computer Tomgraph——CT。 二人于1979年双双获诺贝尔医学奖。

38. 不同物质组合体的视图 自然光投影视图 X射线投影视图 片状X射线投影视图 1、投影数据 平行光按传输方向投射到一个物体后，在垂直于平行光传输方向的平面上生成影象称为投影；用某种传感器取得的投影影象的数据称为投影数据。 最后一幅图象的数据就是重构断面图象的投影数据。显然，仅凭这组数据不能得到断面图象(里面物体的方园无法判断)，我们需要更多的数据。

39. a b μ(x) I0 I I0 I 基准检测器 物体 检测器 lnI lnI0 I0 x + - + 0 A/D 2、物质对X射线的吸收 研究表明：当强度为I0的X射线通过吸收率为μ(x) 的物体时，有下面的关系： 计算机采样示意图

40. 1 2 3 接收器 第 i 条射线 第j个像素 55 N 3、由投影数据重构图像——简单方法 方法一：解联立方程法 设每个网格中的吸收系数为x1～xN ，第i条射线与第xj个象素的相交长度为aij(可以计算得出)，代表第j个像素沿第i条射线的贡献的权值。 对应一个吸收率不均匀物质，μ(x,y)不为常数。可否根据得到的投影值求得μ(x,y)？——结论：如果μ(x,y)已知，既可以重构图象。 如果用pi(θk)表示沿射线方向在射线角度为θk时的总吸收测量值，则可通过解方程组的方法得到x1～xN 。 长度aij可以计算出来，只要找出63个独立方程即可。

41. μ1 μ1 2.5 2.5 1.5 1.5 μ1 μ1 2.5 2.5 3.5 3.5 3 2 7 1.5 2 5 3 3.5 3 1 2 4 1.25 6 2.25 3 4 2.75 3.75 1 1 2.25 5 2.75 4 4 被照射物质内的某一确切位置的吸收函数值是固定的。构造叠代公式进行叠代计算，其收敛值应是吸收函数值。 方法二：叠代法 上次求得的μk(i,j)值 某角度为θ时射线经过的像元数之和 某角度为θ时的投影数据值 同一投影线路上的上次μk(i,j)值之和 给定初值 2.5 2.5+(3-5)/2=1.5 2.5+(7-5)/2=3.5 给定叠代结束条件：

42. p(s,θ) s y s t θ x f(x,y) X射线方向 4、投影定理和傅立叶重构 1)、构造两个坐标系，相交成θ角。根据解析几何的知识有 根据投影数据重构图像的理论进行图像重构。

43. p(s,θ) 沿L方向的微分 s y s t θ x L f(x,y) X射线方向 2)、寻找p(s,θ)与f(s,t)之间的关系 3)、求投影p(s,θ)在某一θ时对s的傅立叶变换： 根据X射线与物体吸收率的关系有： 其中，T=0。可以写成F(S,T)=F(S,0) 角度θ固定后，p(s,θ)的傅立叶变换＝s轴上各点的傅立叶变换。

44. 4)、寻找f(x,y)与p(s,θ)、f(s,t) 之间的关系 由于吸收值与座标系统无关(仅差坐标变换系数)，有： 对f(x,y)和f(s,t)的傅立叶变换为： 雅可比行列式 将(1)代入(2)有(雅可比行列式＝1，dxdy=dsdt)：

45. 观察傅立叶变换公式，与(3)比较，恰好满足旋转公式观察傅立叶变换公式，与(3)比较，恰好满足旋转公式 对比 这个变换是频域的一个旋转变换，于是有 投影定理(以前图所示)： 某角度投影数据对s的傅立叶变换＝s轴上对应点的傅立叶变换。

46. 5、傅立叶重建原理 如果已知F(X,Y)，则 雅可比行列式 该式即为重构图象的生成公式。

47. 实际应用中往往采用下面的方法 公式F(S,T)= P(S,θ)=F(S,0) T=0告诉我们：通过采样数据得到的傅立叶变换值，实际上是x-y座标系统旋转θ角度后在s-t座标系统中s轴上数据的傅立叶变换，如下图所示。如果要得到真正图象的傅立叶变换，需对下图做适当的修改。 Y Y X X

48. Y C X 一般用线性内插公式计算(X,Y)座标点上的值。 式中：C为直角坐标定位点；P1，P2，P3，P4为与C点距离最近之四点之值； d1，d2，d3，d4为C点与P1，P2，P3，P4四点的距离。对于直角坐标定位点在旋转座标轴上之点，可采用零级内插公式得到： 最近点的傅立叶变换值 得到频域图象后，对其进行反变换，即可得到f(x,y)。

49. f(x,y) y s t φ θ Y T x S F(X,Y) θ X 6、卷积反投影重建原理 对空域图象而言，有 空域、频域座标对应关系为 r

50. 根据傅立叶重建公式及座标变换式，有 P(S,θ)乘以|S|等效为两个时域函数的卷积。 P(S,θ)的傅立叶反变换是投影数据p(s,θ) ；而其余部分可设为h(s)，即 于是