图像退化与复原

内容简介
1.图像退化与复原的原理
- 1.1 图像退化的数学模型
- 1.2 图像退化的原理
- 1.3 图像复原的原理
2. 图像去噪
- 2.1 噪声模型
- 2.2 噪声参数的估计
- 2.3 针对噪声的图像复原——空间、频域滤波
3. 退化函数的估计
- 3.1 图像观察估计
- 3.2 试验估计
- 3.3 数学建模估计
4. 图像复原方法
- 4.1 逆滤波
- 4.2 最小均方误差滤波(维纳滤波)
- 4.3 约束最小二乘方滤波
- 4.4 几何均值滤波
- 4.5 度量图像复原效果的指标——信噪比(SNR)

内容简介

图像复原和图像增强两者有较大的重合部分。但图像增强更多的是主观过程，是人们希望通过某种方法加强图像特定信息或细节的操作。而图像复原可以说是客观过程，是因为在使用设备拍照过程中，由于某些客观原因所导致拍摄出来的图像与设备拍摄的理想图像有差异，所以希望通过某些操作来消除这些客观因素给图像带来的干扰，从而呈现出更真实的图像过程。
本篇内容主要知识来自冈萨雷斯第三版的《数字图像处理》，如对细节内容感兴趣的读者可以查阅相关书籍进行参考。本篇博文的内容与原书内容的顺序有所出入，如有不足欢迎留言指正，谢谢。

1.图像退化与复原的原理

1.1 图像退化的数学模型

简单来说，图像的退化是由于某种原因，图像从理想图像转变为实际我们看到的有瑕疵图像的过程。而图像复原，就是通过某种方法，将退化后的图像进行改善，尽量使复原后的图像接近理想图像的过程。那么，整个退化和复原的过程可以用如下图表示：

f(x,y)

退化函数h(x,y)

噪声η(x,y)

g(x,y)

复原滤波器

复原后的图像

这里，我们将理想图像用f(x,y)f(x,y)f(x,y)表示，将退化后的实际图像用g(x,y)g(x,y)g(x,y)表示，退化过程可以分为退化函数h(x,y)h(x,y)h(x,y)以及加载在图像上的噪声η(x,y)\eta(x,y)η(x,y)表示，复原后的图像可以用f^(x,y)\hat f(x,y)f^(x,y)表示。因此，整个图像退化的数学模型即为：
g(x,y)=h(x,y)⋆f(x,y)+η(x,y)g(x,y)=h(x,y)\star f(x,y)+\eta(x,y)g(x,y)=h(x,y)⋆f(x,y)+η(x,y)
频率域表示为：G(u,v)=H(u,v)F(u,v)+N(u,v)频率域表示为：G(u,v)=H(u,v)F(u,v)+N(u,v)频率域表示为：G(u,v)=H(u,v)F(u,v)+N(u,v)
其中，⋆\star⋆是卷积的符号，大写字母G(u,v)表示g(x,y)傅里叶变换后的函数，其他大写函数含义相同。
因此，图像的退化过程可以理解为，经过了一次未知的退化函数卷积，并夹杂着噪声η(x,y)\eta(x,y)η(x,y)的过程。

1.2 图像退化的原理

我们假设，图像的退化过程为：
g(x,y)=H[f(x,y)]+η(x,y)g(x,y)=H[f(x,y)]+\eta(x,y)g(x,y)=H[f(x,y)]+η(x,y)
其中，H是一个算子，具有线性(线性算子具有加性和均匀性)和位置不变性的特点。注意，上述假设是本篇博文内所有内容的前提假设，非常重要。
一幅图像可以当做是对一个连续二维函数进行离散采样得到的，每个像素点上都通过一个冲激函数(也叫脉冲函数，即公式中的δ(x,y)\delta(x,y)δ(x,y))对连续函数进行采样，因此，一幅二维图像可以表示为：
f(x,y)=∫−∞∞∫−∞∞f(α,β)δ(x−α,y−β)dαdβf(x,y)=\int_{-∞}^∞\int_{-∞}^∞f(\alpha,\beta)\delta(x-\alpha,y-\beta)d\alpha d\betaf(x,y)=∫−∞∞∫−∞∞f(α,β)δ(x−α,y−β)dαdβ
先不考虑加性噪声η(x,y)\eta(x,y)η(x,y)的影响，即η(x,y)=0\eta(x,y)=0η(x,y)=0。所以退化过程变为：
g(x,y)=H[f(x,y)]=H[∫−∞∞∫−∞∞f(α,β)δ(x−α,y−β)dαdβ]g(x,y)=H[f(x,y)]=H\bigg[\int_{-∞}^∞\int_{-∞}^∞f(\alpha,\beta)\delta(x-\alpha,y-\beta)d\alpha d\beta\bigg]g(x,y)=H[f(x,y)]=H[∫−∞∞∫−∞∞f(α,β)δ(x−α,y−β)dαdβ]
由于算子H具有线性特征，因此可以拓展到积分上(积分相当于求和)。另外，由于f(α,β)f(\alpha,\beta)f(α,β)和x,y均无关，所以有：
g(x,y)=∫−∞∞∫−∞∞f(α,β)H[δ(x−α,y−β)]dαdβg(x,y)=\int_{-∞}^∞\int_{-∞}^∞f(\alpha,\beta)H[\delta(x-\alpha,y-\beta)]d\alpha d\betag(x,y)=∫−∞∞∫−∞∞f(α,β)H[δ(x−α,y−β)]dαdβ
令h(x,α,y,β)=H[δ(x−α,y−β)]h(x,\alpha,y,\beta)=H[\delta(x-\alpha,y-\beta)]h(x,α,y,β)=H[δ(x−α,y−β)]称为系统H的冲击响应。又因为H算子具有位置不变性，再把忽略的噪声η(x,y)\eta(x,y)η(x,y)加上后有：
g(x,y)=∫−∞∞∫−∞∞f(α,β)h(x−α,y−β)dαdβ+η(x,y)g(x,y)=\int_{-∞}^∞\int_{-∞}^∞f(\alpha,\beta)h(x-\alpha,y-\beta)d\alpha d\beta+\eta(x,y)g(x,y)=∫−∞∞∫−∞∞f(α,β)h(x−α,y−β)dαdβ+η(x,y)
其中，∫−∞∞∫−∞∞f(α,β)h(x−α,y−β)dαdβ\int_{-∞}^∞\int_{-∞}^∞f(\alpha,\beta)h(x-\alpha,y-\beta)d\alpha d\beta∫−∞∞∫−∞∞f(α,β)h(x−α,y−β)dαdβ就是卷积的定义式，因此，退化过程可以表达为
g(x,y)=h(x,y)⋆f(x,y)+η(x,y)g(x,y)=h(x,y)\star f(x,y)+\eta(x,y)g(x,y)=h(x,y)⋆f(x,y)+η(x,y)

1.3 图像复原的原理

由于退化过程被建模为卷积的结果，因此，图像复原就是需要找到具有相反过程的卷积核，所以图像复原通常也叫做“图像去卷积”，所使用的复原滤波器也叫做“去卷积滤波器”。
根据上述分析可知，想要将实际图像恢复成理想图像，主要要完成两个工作，一个是消除图像的噪声干扰，称为图像去噪；另一个则是找到复原滤波器，称为图像去卷积。

2. 图像去噪

2.1 噪声模型

俗话说的好，知己知彼，百战不殆。要去除图像的噪声，首先要了解噪声的种类和性质。通常，图像的噪声模型就是噪声在图像中的像素值的统计特性，可以用概率统计中的概率密度函数(PDF)来表征，表示的是噪声在某个灰度级z产生的概率，一旦发生，则在图像中某个位置会出现一个灰度级为z的像素代替原像素。以下给出常见的噪声模型，其中，z表示噪声的灰度值，zˉ\bar zzˉ表示噪声灰度值的均值，σ\sigmaσ表示z的标准差。：
1).高斯噪声
p(z)=12πσe−(z−zˉ)2/2σ2p(z)=\frac{1}{\sqrt{2\pi}\sigma}e^{-(z-\bar z)^2/2\sigma^2}p(z)=2πσ1e−(z−zˉ)2/2σ2
高斯噪声通常源于诸如电子电路噪声及低照明度或高温带来的传感器噪声。高斯噪声的概率密度图像如下：

2).瑞利噪声
p(z)={2b(z−a)e−(z−a)2/b,z≥a0,z<a,其中,zˉ=a+πb/4,σ2=b(4−π)4p(z)=\begin{cases} \frac{2}{b}(z-a)e^{-(z-a)^2/b},z\ge a\\ 0,z<a \end{cases},其中,\bar z=a+\sqrt{\pi b/4},\sigma^2=\frac{b(4-\pi)}{4}p(z)={b2(z−a)e−(z−a)2/b,z≥a0,z<a,其中,zˉ=a+πb/4,σ2=4b(4−π)
瑞利噪声常见于深度成像，其概率密度图像如下：

3).伽马噪声
p(z)={abzb−1(b−1)!e−az,z≥a0,z<a,其中,zˉ=ba,σ2=ba2p(z)=\begin{cases} \frac{a^bz^{b-1}}{(b-1)!}e^{-az},z\ge a\\ 0,z<a \end{cases},其中,\bar z=\frac{b}{a},\sigma^2=\frac{b}{a^2}p(z)={(b−1)!abzb−1e−az,z≥a0,z<a,其中,zˉ=ab,σ2=a2b
伽马噪声在激光成像中比较常见，其概率密度图像如下：

4).指数噪声
p(z)={ae−az,z≥a0,z<a,其中,zˉ=1a,σ2=1a2p(z)=\begin{cases} ae^{-az},z\ge a\\ 0,z<a \end{cases},其中,\bar z=\frac{1}{a},\sigma^2=\frac{1}{a^2}p(z)={ae−az,z≥a0,z<a,其中,zˉ=a1,σ2=a21
指数噪声是伽马噪声当b=1的特殊情况，其概率密度图像如下：

5).均匀噪声
p(z)={1b−a,a≤z≤b0,其他,其中,zˉ=a+b2,σ2=(b−a)212p(z)=\begin{cases} \frac{1}{b-a},a\le z\le b\\ 0,其他 \end{cases},其中,\bar z=\frac{a+b}{2},\sigma^2=\frac{(b-a)^2}{12}p(z)={b−a1,a≤z≤b0,其他,其中,zˉ=2a+b,σ2=12(b−a)2
均匀噪声的概率密度图像如下：

6).椒盐(脉冲)噪声
p(z)={Pa,z=aPb,z=b1−Pa−Pb,其他p(z)=\begin{cases} P_a,z=a\\ P_b,z=b\\ 1-P_a-P_b,其他 \end{cases}p(z)=⎩⎪⎨⎪⎧Pa,z=aPb,z=b1−Pa−Pb,其他
如果b>a，则灰度级b在图像中将显示为一个亮点；同时，灰度级a在图像中呈现为一个暗点。如果Pa=0P_a=0Pa=0，则图像中没有暗点，只有亮点噪声，此时称为盐粒噪声；如果Pb=0P_b=0Pb=0，则图像中只有暗点，此时称为胡椒噪声。胡椒噪声和盐粒噪声都是椒盐噪声的特例，也称为单极脉冲噪声。椒盐噪声通常发生在设备快速过渡时期，其概率密度图像如下：

7).周期噪声
以上几种是常见的与图像空间坐标无关的噪声。周期噪声是一种空间型噪声，其在图像中通常呈现的样子是，图像中有等间隔的条纹等干扰(如下图)。这说明图中存在周期噪声，周期噪声通常是在图像获取期间由电力或机电干扰产生的。

周期噪声可以通过频率域的滤波来消除，因为周期性噪声通常在频率域会呈现出成对的亮点或线段。

2.2 噪声参数的估计

如前文所述，周期噪声通常通过傅里叶变换到频率域来观测，因此不另做介绍，详情可以参见傅里叶变换的相关内容。
这里主要讲述其他噪声的未知参数如何确定。如果拍照设备和拍照环境已知，那么最好的方式是：

构建同样的条件对纯色的平板状物体进行拍照，得到理论上纯色的图像；
对图像进行直方图分析，与各个噪声模型的概率密度函数图像进行匹配，确定噪声类型；
计算图像的像素均值zˉ\bar zzˉ和标准差σ\sigmaσ，根据这两个参数计算不同噪声模型中的未知参数a,b。

如果无法还原拍照条件或拍照设备未知，那么可以取原图像中恒定灰度值的一小部分替代第一步中的纯色图像，然后再进行后续操作。但是这样做有点以偏概全的意思，需要慎重使用该方法。

2.3 针对噪声的图像复原——空间、频域滤波

理论上来说，只需要实际图像减去噪声函数η(x,y)\eta(x,y)η(x,y)即可消去噪声干扰。但是实际中，噪声函数通常是完全无法测量的，且噪声本身具有随机性，即使确切知道噪声的种类和表达式，也无法确定图像每个位置的噪声情况。因此，直接减去噪声函数来去噪是不可行的。
但是没关系，只针对噪声的情况下，可以用空间滤波的方式进行去噪。当然，对于周期噪声，还是需要使用频域滤波来解决。以下介绍常见的滤波器。
符号说明：

SxyS_{xy}Sxy：表示滤波器的卷积核移动到某位置后与图像重合的区域；
m,n：表示滤波器卷积核的尺寸，m×nm×nm×n

III.空间滤波
1).均值滤波器
a.算术均值滤波器
f^(x,y)=1mn∑(s,t)∈Sxyg(s,t)\hat f(x,y)=\frac{1}{mn}\sum_{(s,t)\in S_{xy}}g(s,t)f^(x,y)=mn1(s,t)∈Sxy∑g(s,t)
算术均值滤波器可以降低噪声，但是对椒盐噪声效果不好，原因很好理解，因为均值滤波将区域内所有像素进行求和取均值，这些像素中就包含了椒盐过亮或过暗的噪声，对均值的影响较大，所以滤波效果不好。另外，均值滤波器由于对某个区域内像素取均值来代替当前位置像素，所以图像会更加平滑，也就是会模糊化。

b.几何均值滤波器
f^(x,y)=[∏s,t∈Sxyg(s,t)]1mn\hat f(x,y)=\bigg[ \prod_{s,t\in S_{xy}}g(s,t) \bigg]^{\frac{1}{mn}}f^(x,y)=[s,t∈Sxy∏g(s,t)]mn1
几何均值滤波器相比于算术均值滤波器来说，丢失的图像信息更少。

c.谐波均值滤波器
f^(x,y)=mn∑(s,t)∈Sxy1g(s,t)\hat f(x,y)=\frac{mn}{\sum_{(s,t)\in S_{xy}}\frac{1}{g(s,t)}}f^(x,y)=∑(s,t)∈Sxyg(s,t)1mn
谐波均值滤波器对于盐粒噪声效果较好，但不适用于胡椒噪声。同时，它也适合处理高斯噪声那类的其他噪声。

d.逆谐波均值滤波器
f^(x,y)=∑(s,t)∈Sxyg(s,t)Q+1∑(s,t)∈Sxyg(s,t)Q\hat f(x,y)=\frac{\sum_{(s,t)\in S_{xy}}g(s,t)^{Q+1}}{\sum_{(s,t)\in S_{xy}}g(s,t)^Q}f^(x,y)=∑(s,t)∈Sxyg(s,t)Q∑(s,t)∈Sxyg(s,t)Q+1
其中，Q称为滤波器的阶数，当Q为正时，适合消除胡椒噪声；当Q为负时，适合消除盐粒噪声。但无法同时消除两种噪声。当Q=0时，退化为算术均值滤波器，当Q=-1时，退化为谐波均值滤波器。

2).统计排序滤波器

中值滤波器
f^(x,y)=median(s,t)∈Sxy{g(s,t)}\hat f(x,y)=median_{(s,t)\in S_{xy}}\{g(s,t)\}f^(x,y)=median(s,t)∈Sxy{g(s,t)}
即取滤波器覆盖范围内所有像素的中位数。中值滤波器对于椒盐噪声或单极椒盐噪声非常有效！
最大值、最小值滤波器
f^(x,y)=max⁡(s,t)∈Sxy{g(s,t)}orf^(x,y)=min⁡(s,t)∈Sxy{g(s,t)}\hat f(x,y)=\max\limits_{(s,t)\in S_{xy}}\{g(s,t)\} \quad or\quad \hat f(x,y)=\min\limits_{(s,t)\in S_{xy}}\{g(s,t)\}f^(x,y)=(s,t)∈Sxymax{g(s,t)}orf^(x,y)=(s,t)∈Sxymin{g(s,t)}
最大最小值滤波器分别适合处理胡椒噪声和盐粒噪声。
中点滤波器
f^(x,y)=12[max⁡(s,t)∈Sxy{g(s,t)}+min⁡(s,t)∈Sxy{g(s,t)}]\hat f(x,y)=\frac{1}{2}[\max\limits_{(s,t)\in S_{xy}}\{g(s,t)\}+\min\limits_{(s,t)\in S_{xy}}\{g(s,t)\}]f^(x,y)=21[(s,t)∈Sxymax{g(s,t)}+(s,t)∈Sxymin{g(s,t)}]
中点滤波器适用于处理随机分布的噪声，如高斯噪声或均匀噪声。不太适合处理椒盐噪声。
修正的阿尔法均值滤波器
f^(x,y)=1mn−d∑(s,t)∈Sxygr(s,t)\hat f(x,y)=\frac{1}{mn-d}\sum_{(s,t)\in S_{xy}}g_r(s,t)f^(x,y)=mn−d1(s,t)∈Sxy∑gr(s,t)
假设在领域SxyS_{xy}Sxy内去掉g(s,t)最小的d/2个灰度值和最大的d/2个灰度值后，剩下的灰度值用gr(s,t)g_r(s,t)gr(s,t)来代替，即可得到上式。其中d的取值范围为[0,mn-1]。当d=mn-1时，退化成中值滤波器。当d为其他值时，适合处理混合多种噪声的情况，如高斯噪声和椒盐噪声混合的情况。

3).自适应滤波器
自适应滤波器的性能通常要优于上述几种滤波器，但是滤波器的复杂度明显提高。

自适应局部降低噪声滤波器
f^(x,y)=g(x,y)−ση2σL2[g(x,y)−mL]\hat f(x,y)=g(x,y)-\frac{\sigma_{\eta}^2}{\sigma_L^2}[g(x,y)-m_L]f^(x,y)=g(x,y)−σL2ση2[g(x,y)−mL]
其中，σL\sigma_LσL和mLm_LmL是SxyS_{xy}Sxy中像素的局部标准差和均值。g(x,y)是实际图像中的像素值。ση\sigma_{\eta}ση是噪声在图像中的标准差，该值通常是未知的，因此通常需要人为设置，但是注意，该值通常都满足ση<σL\sigma_{\eta}<\sigma_Lση<σL的条件，如果不满足，需要在实际操作中强行将比值设置为1，否则会出现负灰度的情况。当然也可以出现负灰度后最后对图像进行重新标准化，但是这样会损失图像的动态范围，所以不建议这么做。
自适应中值滤波器
符号说明：
zminz_{min}zmin：SxyS_{xy}Sxy中的最小灰度值；
zmaxz_{max}zmax：SxyS_{xy}Sxy中的最大灰度值；
zmedz_{med}zmed：SxyS_{xy}Sxy中灰度值的中位数；
zxyz_{xy}zxy：坐标(x,y)处的灰度值；
SmaxS_{max}Smax：SxyS_{xy}Sxy允许的最大尺寸；
对于滤波器窗口移动到图像中的某个位置后，执行以下步骤：
步骤1：
A1=zmed−zmin,A2=zmed−zmaxA_1=z_{med}-z_{min},\quad A_2=z_{med}-z_{max}A1=zmed−zmin,A2=zmed−zmax
如果A1A_1A1>0且A2<0A_2<0A2<0，那么转到步骤二；
否则，增大窗口尺寸；如果窗口尺寸小于SmaxS_{max}Smax，则重复步骤1；
否则，输出zmedz_{med}zmed。
步骤2：
B1=zxy−zmin,B2=zxy−zmaxB_1=z_{xy}-z_{min},\quad B_2=z_{xy}-z_{max}B1=zxy−zmin,B2=zxy−zmax
如果B1>0B_1>0B1>0且B2<0B_2<0B2<0，则输出zxyz_{xy}zxy；
否则输出zmedz_{med}zmed

从该算法的第一步可以看出，首先要确认当前正在处理的图像区域的中值zmedz_{med}zmed是不是极端值，如果是极端值，很可能是椒盐噪声，因此不可取当前中值，需要扩大滤波器窗口，选择其他中值。如果扩大到一定程度都是极端值，那么就只能接受，同时也说明在较大范围内，极端值是比较均衡的，或许不是噪声。
第二步，如果zmedz_{med}zmed不是极端值，那么进而检查当前考察的像素zxyz_{xy}zxy是不是极端值，如果不是极端值，说明当前像素很可能不是椒盐噪声，可以取原值；反之，如果是极端值，那么很可能当前像素是噪声，因此用中值替代。
通过上述操作，自适应中值滤波器可以消除椒盐噪声的同时，减少图像本身的细节损失，尽可能地保留原图像的信息，只处理噪声点。另外，在滤波器移动到下一个位置时，可以仅仅用新加入的像素来更新中值，不需要所有像素重新计算一遍，可以提高算法效率。

IIIIII.频域滤波
频域滤波需要先用傅里叶变换将图像转换到频域，得到图像的频谱，然后对频谱图像进行处理。处理后的频谱再进行反傅里叶变换得到处理后的图像。关键就在于对频谱的处理方式。
1).带阻滤波器

如上图所示，带阻滤波器就是将黑色范围的频谱归零，只留下白色部分的频谱的处理方式。其中，黑色圆环的半径是截止频率D0D_0D0，宽度是带宽W。
2).带通滤波器
功能与带阻滤波器恰恰相反。由于带通滤波器会消除图像大量信息，通常不会直接使用，而是将带通滤波器作为提取工具，即提取白色圆环范围内的图像信息，单独进行查看。通常，提取的信息就是噪声模板，即可以将噪声图像单独提取出来。

3).陷波滤波器
陷波滤波器通常是设计一对关于频域中心对称的区域进行过滤。对称是因为傅里叶变换具有对称性，要获得有效的滤波效果，滤波器需要以对称的形式出现。陷波滤波器的设计比较自由，可以根据情况设计任何形状、任何对数的区域。下图是陷波滤波器的示意图，上方是理想的滤波器，下方两种分别是布特沃斯和高通陷波滤波器，它们更接近于实际的硬件滤波器，因为硬件的滤波是无法达到理想滤波器这种截断效果的。
4).最佳陷波滤波
步骤：

用带通滤波器或类似思想，提取图像中的噪声部分得到噪声的傅里叶频谱N(u,v)；
对N(u,v)作傅里叶反变换得到噪声函数η(x,y)\eta(x,y)η(x,y)；
通过设置加权函数w(x,y)w(x,y)w(x,y)对噪声函数进行加权，然后从真实图像中减去加权后的噪声函数，得到理想图像的估计：
f^(x,y)=g(x,y)−w(x,y)η(x,y)\hat f(x,y)=g(x,y)-w(x,y)\eta(x,y)f^(x,y)=g(x,y)−w(x,y)η(x,y)
加权函数的设计遵循以下原则：使估计值f^(x,y)\hat f(x,y)f^(x,y)在每一点(x,y)的指定领域上的方差最小。最终得到加权函数的计算公式如下：
w(x,y)=g(x,y)η(x,y)‾−gˉ(x,y)ηˉ(x,y)η2‾(x,y)−ηˉ2(x,y)w(x,y)=\frac{\overline{g(x,y)\eta(x,y)}-\bar g(x,y)\bar \eta(x,y)}{\overline{\eta^2}(x,y)-\bar\eta^2(x,y)}w(x,y)=η2(x,y)−ηˉ2(x,y)g(x,y)η(x,y)−gˉ(x,y)ηˉ(x,y)

其中，是对应函数在指定滤波器领域范围内的像素均值。比如g(x,y)η(x,y)‾\overline{g(x,y)\eta(x,y)}g(x,y)η(x,y)是函数g(x,y)η(x,y)g(x,y)\eta(x,y)g(x,y)η(x,y)在滤波器指定范围内的像素均值。通常，w(x,y)w(x,y)w(x,y)不用计算每个x,y的数值，通过移动滤波器，使得滤波器范围每次移动都不重叠，可以假设在滤波器范围内的每个w(x,y)w(x,y)w(x,y)的数值都相等。这样可以简化计算。

3. 退化函数的估计

3.1 图像观察估计

根据图像退化函数H算子的位置不变性，可以选择图像中一个具有很强信号内容的区域(高对比度区域)，然后忽略该区域中噪声的影响。最后对该区域的图像进行单独、精致地图像处理，使得该区域图像尽可能地恢复到理想图像的样子(甚至可以用手动方式处理)。假设gs(x,y)g_s(x,y)gs(x,y)是要观察的子图像，f^s(x,y)\hat f_s(x,y)f^s(x,y)是处理过的子图像，即我们处理子图像后认为它是理想图像的一个估计。由于忽略了噪声影响，所以可以得到：
Hs(u,v)=Gs(u,v)F^(u,v)H_s(u,v)=\frac{G_s(u,v)}{\hat F(u,v)}Hs(u,v)=F^(u,v)Gs(u,v)
根据图像退化函数H算子的位置不变性，就可以假定H(u,v)=Hs(u,v)H(u,v)=H_s(u,v)H(u,v)=Hs(u,v)。但是，其实这种方式非常繁琐，效果也不一定好，需要反复试错来调整效果，因此不是非常好的办法。

3.2 试验估计

前提是可以获得与拍摄退化图像相同的设备，然后利用该设备进行实验，通过对实验图像进行处理得到退化函数。最简单的实验设计是，利用一个冲激亮点照在纯色背景上，然后用设备进行拍照，得到的图像就是gs(x,y)g_s(x,y)gs(x,y)，由于冲激函数傅里叶变换后是常数，所以实验估计的退化函数是
Hs(u,v)=Gs(u,v)AH_s(u,v)=\frac{G_s(u,v)}{A}Hs(u,v)=AGs(u,v)

3.3 数学建模估计

这种方法是最科学，也是最难的。需要对拍照退化模型进行建模，以及对拍照环境、条件进行建模，最终直接得到真实图像的函数、退化函数和理想图像函数之间的关系。
比如，拍照时拍照设备进行了x,y方向的匀速运动，那么模型就可以构建为：g(x,y)=∫0Tf[x−x0(t),y−y0(t)]dtg(x,y)=\int_0^T f[x-x_0(t),y-y_0(t)]dtg(x,y)=∫0Tf[x−x0(t),y−y0(t)]dt，对上式傅里叶变换后即可得到：
G(u,v)=F(u,v)∫0Te−j2π[ux0(t)+vy0(t)]dt,∴H(u,v)=∫0Te−j2π[ux0(t)+vy0(t)]dtG(u,v)=F(u,v)\int_0^T e^{-j2\pi[ux_0(t)+vy_0(t)]}dt,∴H(u,v)=\int_0^T e^{-j2\pi[ux_0(t)+vy_0(t)]}dtG(u,v)=F(u,v)∫0Te−j2π[ux0(t)+vy0(t)]dt,∴H(u,v)=∫0Te−j2π[ux0(t)+vy0(t)]dt

4. 图像复原方法

注意，下述方法的前提均是已知退化函数h(x,y)或H(u,v)，且噪声变量与图像变量不相关。在此前提下，才可以对退化后的图像进行复原。
为什么不能直接根据F^(u,v)=G(u,v)/H(u,v)\hat F(u,v)=G(u,v)/H(u,v)F^(u,v)=G(u,v)/H(u,v)，然后傅里叶逆变换对图像复原？因为实际情况下，图像都是同时存在噪声函数和退化函数的，上式可以变换为：F^(u,v)=G(u,v)/H(u,v)=F(u,v)+N(u,v)/H(u,v)\hat F(u,v)=G(u,v)/H(u,v)=F(u,v)+N(u,v)/H(u,v)F^(u,v)=G(u,v)/H(u,v)=F(u,v)+N(u,v)/H(u,v)，其中N(u,v)是噪声函数的频域表达。而噪声函数通常是无法确切知道的，我们之前也只是用空间滤波的方式来减弱噪声，并没有真正得到噪声函数。除此之外，如果退化函数H(u,v)的数值很小，很明显，上式加号的后一项会很容易支配真实图像F(u,v)的数值，从而导致图像在转换回空间域后失真的情况。所以，直接采取傅里叶逆变换复原图像是很难得到比较好的效果的。
本节中的滤波均指的是复原滤波。

4.1 逆滤波

假设不考虑噪声函数的影响，为了解决上文中提到的H(u,v)数值很小产生的图像失真问题，我们可以采取频率限制的办法进行逆滤波。
通常，高频范围内的H(u,v)数值较小，低频范围通常不会出现H(u,v)过小的问题，因此，只需要对H(u,v)进行低通滤波，截断高频信息，再进行图像逆滤波，则可以一定程度解决图像失真的问题。而这个截断频率则需要人为选择和调试。但就算参数调的比较好，也无法解决噪声函数干扰的问题。

4.2 最小均方误差滤波(维纳滤波)

该方法建立在图像和噪声都是随机变量的基础上，且噪声和图像不相关，其中一个或另一个均值为0。原则是找到一幅理想图像的估计f^\hat ff^，使得理想图像和估计图像之间的均方误差最小：min⁡e2=E{(f−f^)2}\min e^2=E\{(f-\hat f)^2\}mine2=E{(f−f^)2}。
具体推导在此省略，详情可以参阅N.Wiener[1942]年提出的论文或其他文献。基于上述原则，可以得到如下估计图像的频域表达式：
F^(u,v)=[1H(u,v)∣H(u,v)∣2∣H(u,v)∣2+Sη(u,v)/Sf(u,v)]⇒F^(u,v)=[1H(u,v)∣H(u,v)∣2∣H(u,v)∣2+K]\hat F(u,v)=\bigg[ \frac{1}{H(u,v)}\frac{|H(u,v)|^2}{|H(u,v)|^2+S_{\eta}(u,v)/S_f(u,v)} \bigg]\Rightarrow \hat F(u,v)=\bigg[ \frac{1}{H(u,v)}\frac{|H(u,v)|^2}{|H(u,v)|^2+K} \bigg]F^(u,v)=[H(u,v)1∣H(u,v)∣2+Sη(u,v)/Sf(u,v)∣H(u,v)∣2]⇒F^(u,v)=[H(u,v)1∣H(u,v)∣2+K∣H(u,v)∣2]
其中，Sη(u,v)=∣N(u,v)∣2S_{\eta}(u,v)=|N(u,v)|^2Sη(u,v)=∣N(u,v)∣2是噪声的功率谱；Sf(u,v)=∣F(u,v)∣2S_f(u,v)=|F(u,v)|^2Sf(u,v)=∣F(u,v)∣2是未退化图像的功率谱。通常，除了白噪声的功率谱是常数以外，其他噪声的功率谱是未知的，所以通常使用将对应的项用常数K来代替，简化上述计算。K是可以人为选择和调整的超参数。

4.3 约束最小二乘方滤波

该方法与维纳方法相比，最好的地方在于，除了退化函数外，它只需要得到噪声方差和均值，就可以进行图像复原。而这两个物理量相对于维纳滤波所需的两个功率谱来说，是更加容易得到的。
根据公式F^(u,v)=F(u,v)+N(u,v)/H(u,v)\hat F(u,v)=F(u,v)+N(u,v)/H(u,v)F^(u,v)=F(u,v)+N(u,v)/H(u,v)可知，退化函数H相对于未知的噪声是非常敏感的。为了减小H对噪声的敏感度，最好的办法就是改变图像复原最优化原则，使用一种以平滑度度量为主的优化原则，如使用图像的拉普拉斯算子作为优化基础：
最小准则函数C=∑x=0M−1∑y=0N−1[▽2f(x,y)]2，其约束条件为：∥g−Hf^∥=∥η∥2最小准则函数C=\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}[\bigtriangledown^2 f(x,y)]^2，其约束条件为：\|g-H\hat f\|=\|\eta\|^2最小准则函数C=x=0∑M−1y=0∑N−1[▽2f(x,y)]2，其约束条件为：∥g−Hf^∥=∥η∥2
该问题在频率域中求解得到:
F^(u,v)=[1H(u,v)∣H(u,v)∣2∣H(u,v)∣2+γ∣P(u,v)∣2]\hat F(u,v)=\bigg[ \frac{1}{H(u,v)}\frac{|H(u,v)|^2}{|H(u,v)|^2+\gamma|P(u,v)|^2} \bigg]F^(u,v)=[H(u,v)1∣H(u,v)∣2+γ∣P(u,v)∣2∣H(u,v)∣2]
其中，P(u,v)P(u,v)P(u,v)是拉普拉斯滤波器p(x,y)=[0−10−14−10−10]p(x,y)=\begin{bmatrix}0&-1&0\\-1&4&-1\\0&-1&0\end{bmatrix}p(x,y)=⎣⎡0−10−14−10−10⎦⎤的傅里叶变换。
γ\gammaγ是后续需要求解的系数，注意到，当γ=0\gamma=0γ=0时，上述解退化成逆滤波过程。
如何求解γ\gammaγ，使得上述解能尽量满足最优问题的约束条件呢？
定义残差向量r=g−Hf^r=g-H\hat fr=g−Hf^，注意，此处的残差向量是图像中所有像素值的残差。由于f^\hat ff^是γ\gammaγ的函数，所以残差r也是关于γ\gammaγ的函数，且可以证明：
ϕ(γ)=rTr=∥r∥2\phi(\gamma)=r^Tr=\|r\|^2ϕ(γ)=rTr=∥r∥2
是γ\gammaγ的单调递增函数。理论上，要满足约束条件，就有∥r∥2=∥η∥2\|r\|^2=\|\eta\|^2∥r∥2=∥η∥2，由于∥r∥2\|r\|^2∥r∥2的单调性，所以可以设置我们能接受的精度范围aaa，然后进行如下循环：

设定初始的γ\gammaγ值；
计算R(u,v)=G(u,v)−H(u,v)F^(u,v)R(u,v)=G(u,v)-H(u,v)\hat F(u,v)R(u,v)=G(u,v)−H(u,v)F^(u,v)，傅里叶逆变换后得到r(x,y)，再计算∥r∥2\|r\|^2∥r∥2；
如果满足精度条件∥r∥2∈[∥η∥2−a,∥η∥2+a]\|r\|^2\in\big[\|\eta\|^2-a,\|\eta\|^2+a\big]∥r∥2∈[∥η∥2−a,∥η∥2+a]，那么停止循环，返回当前γ\gammaγ；否则，如果∥r∥2<∥η∥2−a\|r\|^2<\|\eta\|^2-a∥r∥2<∥η∥2−a，那么增大γ\gammaγ；如果∥r∥2>∥η∥2+a\|r\|^2>\|\eta\|^2+a∥r∥2>∥η∥2+a，那么减小γ\gammaγ，重新返回步骤2。
其中，∥η∥2=MN[ση2+mη2]\|\eta\|^2=MN[\sigma_{\eta}^2+m_{\eta}^2]∥η∥2=MN[ση2+mη2]，ση,mη\sigma_{\eta},m_{\eta}ση,mη可以通过 2.2节噪声参数的估计 所述方法进行计算。

∥η∥2\|\eta\|^2∥η∥2求解的推导过程：
已知：
{ση2=1MN∑x=0M−1∑y=0N−1[η(x,y)−mη]2mη=1MN∑x=0M−1∑y=0N−1η(x,y)∥η∥2=sumx=0M−1∑y=0N−1η(x,y)\begin{cases} \sigma_{\eta}^2=\frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}[\eta(x,y)-m_{\eta}]^2\\ m_{\eta}=\frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}\eta(x,y)\\ \|\eta\|^2=sum_{x=0}^{M-1}\sum_{y=0}^{N-1}\eta(x,y) \end{cases}⎩⎪⎨⎪⎧ση2=MN1∑x=0M−1∑y=0N−1[η(x,y)−mη]2mη=MN1∑x=0M−1∑y=0N−1η(x,y)∥η∥2=sumx=0M−1∑y=0N−1η(x,y)
∴ση2=1MN∑x=0M−1∑y=0N−1[η2(x,y)−2η(x,y)mη+mη2]=1MN∥η∥2−2mη2+mη2∴\sigma_{\eta}^2=\frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}[\eta^2(x,y)-2\eta(x,y)m_{\eta}+m_{\eta}^2]=\frac{1}{MN}\|\eta\|^2-2m_{\eta}^2+m_{\eta}^2∴ση2=MN1x=0∑M−1y=0∑N−1[η2(x,y)−2η(x,y)mη+mη2]=MN1∥η∥2−2mη2+mη2
∴∥η∥2=MN[ση2+mη2]∴\|\eta\|^2=MN[\sigma_{\eta}^2+m_{\eta}^2]∴∥η∥2=MN[ση2+mη2]

4.4 几何均值滤波

其实就是在维纳滤波的基础上，增加几何均值的框架，公式如下：
F^(u,v)=[1H(u,v)]α[1H(u,v)∣H(u,v)∣2∣H(u,v)∣2+βSη(u,v)/Sf(u,v)]1−α\hat F(u,v)=\bigg[ \frac{1}{H(u,v)} \bigg]^{\alpha}\bigg[ \frac{1}{H(u,v)}\frac{|H(u,v)|^2}{|H(u,v)|^2+\beta S_{\eta}(u,v)/S_f(u,v)} \bigg]^{1-\alpha}F^(u,v)=[H(u,v)1]α[H(u,v)1∣H(u,v)∣2+βSη(u,v)/Sf(u,v)∣H(u,v)∣2]1−α
其中α,β\alpha,\betaα,β都是人为设定的参数。当α=0.5\alpha=0.5α=0.5时，就是几何均值的形式。当α=0.5，β=1\alpha=0.5，\beta=1α=0.5，β=1时，称为谱均衡滤波器。上式其实根据系数的不同表示了不同的滤波器，可以看作是一个滤波器族，是复原滤波器的整合数学表达。

4.5 度量图像复原效果的指标——信噪比(SNR)

注意，通常度量图像效果的数学指标和人眼观测的实际效果有细微的出入，所以不代表指标效果越好，所看到的图像效果越好。
信噪比原理上是SNR=未退化图像的功率谱噪声的功率谱SNR=\frac{未退化图像的功率谱}{噪声的功率谱}SNR=噪声的功率谱未退化图像的功率谱，但是由于公式中两者都无法准确得到，因此通常用如下形式来近似：
SNR=∑x=0M−1∑y=0N−1f^(x,y)2∑x=0M−1∑y=0N−1[f(x,y)−f^(x,y)]2SNR=\frac{\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}\hat f(x,y)^2}{\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}[f(x,y)-\hat f(x,y)]^2}SNR=∑x=0M−1∑y=0N−1[f(x,y)−f^(x,y)]2∑x=0M−1∑y=0N−1f^(x,y)2

数字图像处理——图像退化与复原相关推荐

数字图像处理——图像退化（大气湍流模型与运动模糊模型）与图像复原（逆滤波与维纳滤波）
一.图像退化一般来说,图像的退化模型可以表示为其中g(x,y) 表示退化后的图像,h(x,y)表示退化模型,f(x,y)表示原图像,n(x,y)表示噪声. 在频域上面可以表示为下面介绍常见的两种 ...
数字图像学笔记 —— 16. 图像退化与复原（自适应滤波之「最小均方差滤波」）
文章目录图像恢复的一般运算过程什么是「最小均方差滤波」实现步骤实现代码最后的结果图像恢复的一般运算过程我们从前几章的基本理论出发,退化信号恢复成原始信号的步骤,可以概括成两步基本公式.对 ...
数字图像处理图像反转的实现_使用8086微处理器反转16位数字
数字图像处理图像反转的实现 Problem statement: 问题陈述: Write an assembly language program in 8086 microprocessor to ...
数字图像处理图像反转的实现_反转8位数字| 8085微处理器
数字图像处理图像反转的实现 Problem statement: 问题陈述: To reverse 8 bits number using 8085 microprocessors. 使用8085微处 ...
数字图像处理课设图像的锐化_数字图像处理图像锐化处理.ppt
数字图像处理图像锐化处理 4.7.2 灰度级到彩色转换灰度级到彩色转换(例) 在HSI彩色空间的直方图均衡强度均衡处理没有改变图像的色调和饱和度值,但它的确影响了整体图像的彩色感观. 向量分量可以用 ...
matlab 求其骨架,数字图像处理图像的骨架生成和提取（Matlab）三种方法
[实例简介] 数字图像处理图像的骨架生成和提取(Matlab),有三种方法,推荐给大家! [实例截图] [核心代码] Programe ├── Programe1 │ ├── 00.JPG │ ...
数字图像学笔记——13. 图像退化与复原（退化函数的评估方法：观察法、实验法、数学建模法与湍流导致的退化）
在对受到多种原因影响的图像进行复原时,我们经常需要先行评估对图像质量产生影响的退化函数,有时甚至需要尝试建模.通过这些手段,能够最大程度上恢复图像上的噪音,并重建高清的图像细节. 文章目录线性位置不 ...
数字图像学笔记——14. 图像退化与复原（线性退化）
文章目录运动导致的退化(线性退化) 水平运动导致的退化垂直运动导致的退化运动导致的退化(线性退化) 在上一章 <数字图像学笔记--13. 图像退化与复原(退化函数的评估方法:观察法.实验法 ...
Matlab数字图像处理——图像的空间变换
Matlab空间变换函数 imtransform Matlab空间变换函数 imtransform 可以实现图像仿射变换(如平移.旋转.剪切.缩放).投影变换, 该函数可与 maketform 配合 ...

数字图像处理——图像退化与复原