一、图像退化/复原过程的模型
图像退化过程可以理解为将原始图片 f ( x , y ) f(x,y) f(x,y)经过退化函数 H H H的处理，在加上一个噪声项从而获得退化后的图像 g ( x , y ) g(x,y) g(x,y)。而复原过程即为结合给定的退化函数 H H H与噪声 η ( x , y ) \eta(x,y) η(x,y)重构原始图像的估计结果 f ^ ( x , y ) \hat f(x,y) f^(x,y)，并在复原过程希望 f ^ ( x , y ) \hat f(x,y) f^(x,y)尽可能与 f ( x , y ) f(x,y) f(x,y)相似。
退化函数通常遵循如下范式：(“ ⋆ \star ⋆”表示空间卷积)
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)
上述为课本中的理论介绍，而在实际应用场景中模糊、散焦、运动模糊等等都属于退化。而复原即为将获取到的图片复原成原场景下的图片。
二、噪声模型
数字图像的噪声主要来源于图像的获取或传输过程。
2.1噪声的空间和频率特性
频率特性是指傅里叶域中噪声的频率内容（即相对于电磁披谱的频率）。 例如，当噪声的傅里叶谱是常量时，通常称为白噪声。之所以称之为白噪声是因为物理中白光等比例包含可见光谱中几乎所有频率，而以相同比例包含所有频率的函数的傅里叶谱是常量。
2.2一些重要的噪声概率密度函数
2.2.1高斯噪声
高斯噪声的概率密度函数(PDF)由下式给出：
p ( z ) = 1 2 π σ e − ( z − z ‾ ) 2 / 2 σ 2 p(z)=\frac{1}{\sqrt{2\pi }\sigma }e^{-(z-\overline{z})^{2}/2\sigma ^{2}} p(z)=2π σ1e−(z−z)2/2σ2
上述式子中， z z z表示灰度值， z ‾ \overline{z} z表示 z z z的均值， σ \sigma σ表示 z z z的标准差，标准差的平方 σ 2 \sigma^{2} σ2称为 z z z的方差。 z z z的值有70%落在 [ z ‾ − σ , z ‾ + σ ] [\overline{z} -\sigma ,\overline{z} +\sigma ] [z−σ,z+σ]范围内，有95%落在 [ z ‾ − 2 σ , z ‾ + 2 σ ] [\overline{z} -2\sigma ,\overline{z} +2\sigma ] [z−2σ,z+2σ]范围内。高斯函数的曲线及加噪效果如图所示：
2.2.2瑞利噪声
瑞利噪声的PDF由下式给出：
p ( z ) = { 2 b ( z − a ) e − ( z − a ) 2 / b z ≥ a 0 z < a p(z)=\left\{\begin{matrix} \frac{2}{b}(z-a)e^{-(z-a)^{2}/b} & z\geq a\\ 0 & z< a \end{matrix}\right. p(z)={b2(z−a)e−(z−a)2/b0z≥az<a
概率密度的均值和方差如下所示：
μ = a + π b / 4 \mu =a+\sqrt{\pi b/4} μ=a+πb/4
σ 2 = b ( 4 − π ) 4 \sigma ^{2}=\frac{b(4-\pi )}{4} σ2=4b(4−π)
距原点的位移和密度的基本形状向右变形了这一事实。瑞利密度对于近似歪斜的直方图十分适用。瑞利密度的曲线及加噪效果如下所示：
2.2.3爱尔兰（伽马）噪声
伽马噪声的PDF由下式给出：
p ( z ) = { a b z b − 1 ( b − 1 ) ! e − a z z ≥ a 0 z < a p(z)=\left\{\begin{matrix} \frac{a^{b}z^{b-1}}{(b-1)!}e^{-az}&z\geq a \\ 0 & z< a \end{matrix}\right. p(z)={(b−1)!abzb−1e−az0z≥az<a
其概率密度的均值和方差如下所示：
μ = b a , σ 2 = b a 2 \mu =\frac{b}{a},\sigma ^{2}=\frac{b}{a^{2}} μ=ab,σ2=a2b
尽管通常将该概率密度函数称为伽马密度，但只有分母为伽马函数 Γ ( b ) \Gamma(b) Γ(b)才是正确的，其名称应为爱尔兰密度。该密度函数的曲线及加噪效果如下所示：
2.2.4指数噪声
指数噪声的PDF由下式给出（a>0）：
p ( z ) = { a e − a z z ≥ 0 0 z < 0 p(z)=\left\{\begin{matrix}ae^{-az} & z\geq 0\\ 0 & z< 0\end{matrix}\right. p(z)={ae−az0z≥0z<0
该概率密度函数的均值和方差如下所示：
μ = 1 a , σ 2 = 1 a 2 \mu =\frac{1}{a},\sigma ^{2}=\frac{1}{a^{2}} μ=a1,σ2=a21
指数分布的PDF是当b=1时爱尔兰分布的特殊情况。该密度函数的曲线及加噪效果如下所示：
2.2.5均匀噪声
均匀分布噪声的PDF由下式给出：
p ( z ) = { 1 b − a a ≤ z ≤ b 0 其他 p(z)=\left\{\begin{matrix}\frac{1}{b-a} & a\leq z\leq b\\ 0 & 其他\end{matrix}\right. p(z)={b−a10a≤z≤b其他
该概率密度函数的均值和方差如下所示：
μ = a + b 2 , σ 2 = ( b − a ) 2 12 \mu =\frac{a+b}{2},\sigma ^{2}=\frac{(b-a)^{2}}{12} μ=2a+b,σ2=12(b−a)2
该密度函数的曲线及加噪效果如下所示：
2.2.6脉冲（椒盐）噪声
脉冲噪声的PDF由下式给出：
p ( z ) = { P a z = a P b z = b 1 − P a − P b 其他 p(z)=\left\{\begin{matrix} P_{a} &z=a \\P_{b} &z=b \\ 1-P_{a}-P_{b} & 其他 \end{matrix}\right. p(z)=⎩ ⎨ ⎧PaPb1−Pa−Pbz=az=b其他
如果 P a P_{a} Pa或 P b P_{b} Pb为零，则脉冲噪声称为单极脉冲；如果 P a P_{a} Pa或 P b P_{b} Pb均不为零，则脉冲噪声称为双极脉冲噪声或椒盐噪声。
脉冲噪声可以为正，也可为负。标定以后，脉冲噪声总是数字化为最大值（纯黑或纯白）。通常，负脉冲以黑点（胡椒点）出现，正脉冲以白点（盐点）出现。该密度函数的曲线及加噪效果如下所示：
2.3周期噪声
周期噪声是在图像获取中从电力或机电干扰中产生，周期噪声可以通过频率域滤波显著减少。
三、只存在噪声的复原----空间滤波
当唯一退化是噪声时吗，退化函数变为如下形式：
g ( x , y ) = f ( x , y ) + η ( x , y ) g(x,y)=f(x,y)+\eta (x,y) g(x,y)=f(x,y)+η(x,y)
噪声项通常是未知的，可以选择空间滤波方法进行图像复原。
3.1均值滤波器
3.1.1算术均值滤波器
这是最简单的均值滤波器，令 S x y S_{xy} Sxy表示中心在 ( x , y ) (x,y) (x,y)，尺寸为 m × n m×n m×n的矩形窗口，平滑了一幅图像的局部变化，在模糊了结果的同时减少了噪声。
f ^ ( x , y ) = 1 m n ∑ ( s , t ) ∈ S x y g ( s , t ) \widehat{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)
3.1.2几何均值滤波器
几何均值滤波器所达到的平滑度可以与算术 均值滤波器相比，但几何均值滤波器在滤波过程中，与算术均值滤波器相比，会丢失更少的图像细节。
f ^ ( x , y ) = [ ∏ ( s , t ) ∈ S x y g ( s , t ) ] 1 m n \widehat{f}(x,y)=[\prod _{(s,t)\in S_{xy}}g(s,t)]^{\frac{1}{mn}} f (x,y)=[(s,t)∈Sxy∏g(s,t)]mn1
3.1.3谐波均值滤波器
谐波均值滤波器对于“盐”噪声效果好，但不适用于“胡椒”噪声，善于处理高斯噪声。
f ^ ( x , y ) = m n ∑ ( s , t ) ∈ S x y 1 g ( s , t ) \widehat{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
3.1.4逆谐波均值滤波器
f ^ ( x , y ) = ∑ ( s , t ) ∈ S x y g ( s , t ) Q + 1 ∑ ( s , t ) ∈ S x y g ( s , t ) Q \widehat{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，逆谐波均值滤波器等价于谐波均值滤波器。
3.2统计排序滤波器
3.2.1中值滤波器
f ^ ( x , y ) = m e d i a n ( s , t ) ∈ S x y g ( s , t ) \widehat{f}(x,y)=median_{(s,t)\in S_{xy}}g(s,t) f (x,y)=median(s,t)∈Sxyg(s,t)
在相同尺寸下，比起均值滤波器引起的模糊少，对单极或双极脉冲噪声非常有效。
3.2.2最大值和最小值滤波器
最大值滤波器如下所示：
f ^ ( x , y ) = m a x ( s , t ) ∈ S x y g ( s , t ) \widehat{f}(x,y)=max_{(s,t)\in S_{xy}}g(s,t) f (x,y)=max(s,t)∈Sxyg(s,t)
;用于发现图像中的最亮点，可以有效过滤胡椒噪声）。
最小值滤波器如下所示：
f ^ ( x , y ) = m i n ( s , t ) ∈ S x y g ( s , t ) \widehat{f}(x,y)=min_{(s,t)\in S_{xy}}g(s,t) f (x,y)=min(s,t)∈Sxyg(s,t)
用于发现图像中的最暗点，可以有效过滤盐粒噪声。
3.2.3中点滤波器
f ^ ( x , y ) = 1 2 [ m a x ( s , t ) ∈ S x y g ( s , t ) + m i n ( s , t ) ∈ S x y g ( s , t ) ] \widehat{f}(x,y)=\frac{1}{2}[max_{(s,t)\in S_{xy}}g(s,t)+min_{(s,t)\in S_{xy}}g(s,t)] f (x,y)=21[max(s,t)∈Sxyg(s,t)+min(s,t)∈Sxyg(s,t)]
结合了顺序统计和求平均，对于高斯和均匀随机分布这类噪声有最好的效果。
3.2.4修正的阿尔法均值滤波器
假设在 S x y S_{xy} Sxy邻域内去掉 g ( s , t ) g(s,t) g(s,t)最高灰度值的 d / 2 d/2 d/2和最低灰度值的 d / 2 d/2 d/2。令 g r ( s , t ) g_{r}(s,t) gr(s,t)代表剩余的 m n − d mn-d mn−d个像素，当 d = 0 d=0 d=0，等价于算术均值滤波器；当 d = ( m n − 1 ) / 2 d=(mn-1)/2 d=(mn−1)/2，等价于中值滤波器；当 d d d取其它值时，适用于包括多种噪声的情况下，例如高斯噪声和椒盐噪声混合的情况。
f ^ ( x , y ) = 1 m n − d ∑ ( s , t ) ∈ S ( x y ) g r ( s , t ) \widehat{f}(x,y)=\frac{1}{mn-d}\sum _{(s,t)\in S_(xy)}g_{r}(s,t) f (x,y)=mn−d1(s,t)∈S(xy)∑gr(s,t)
3.3自适应滤波器
3.3.1自适应局部降低噪声滤波器
滤波器作用于局部区域 S x y S_{xy} Sxy，滤波器在该区域中心任何一点 ( x , y ) (x,y) (x,y)的响应基于以下4个量：a. g ( x , y ) g(x,y) g(x,y)，带噪图像在点 ( x , y ) (x,y) (x,y)上的值；b. σ n 2 \sigma _{n}^{2} σn2，污染 f ( x , y ) f(x,y) f(x,y)以形成 g ( x , y ) g(x,y) g(x,y)的噪声方差；c. m L m_{L} mL， S x y S_{xy} Sxy中像素的局部均值; d. σ L 2 \sigma _{L}^{2} σL2， S x y S_{xy} Sxy中像素的局部方差。
滤波器的预期性能如下：
1、如果 σ n = 0 \sigma _{n}=0 σn=0，滤波器返回g(x,y)的值。因为在g(x,y)下零噪声的情况等同于f(x,y)。
2、如果局部方差 σ L 2 \sigma _{L}^{2} σL2与 σ n 2 \sigma _{n}^{2} σn2高相关，滤波器返回一个g(x,y)的近似值。
3、如果 σ L 2 = σ n 2 \sigma _{L}^{2}=\sigma _{n}^{2} σL2=σn2，滤波器返回区域 S x y S_{xy} Sxy上像素的算术均值。这样局部噪声用求平均 m L m_{L} mL来降低。
基于这些假设得到的 f ^ ( x , y ) \widehat{f}(x,y) f (x,y)的自适应表达式可以写成：
f ^ ( x , y ) = g ( x , y ) − σ η 2 σ L 2 [ g ( x , y ) − m L ] \widehat{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]
唯一需要知道或估计的未知量是噪声方差 σ n 2 \sigma _{n}^{2} σn2，其它参数可以从 S x y S_{xy} Sxy中的像素计算出来。
3.3.2自适应中值滤波器
传统中值滤波器只能处理空间密度不大的冲激噪声（ p a p_{a} pa， p b p_{b} pb<0.2），而自适应中值滤波器可以处理具有更大概率的冲激噪声，可以在平滑非冲激噪声时保存细节，而传统中值滤波器无法做到。
考虑如下符号：
1、 z m i n = S x y z_{min}=S_{xy} zmin=Sxy中灰度级的最小值
2、 z m a x = S x y z_{max}=S_{xy} zmax=Sxy中灰度级的最大值
3、 z m e d = S x y z_{med}=S_{xy} zmed=Sxy中灰度级的中值
4、 z x y z_{xy} zxy＝在坐标（x,y）上的灰度级
5、 S m a x = S x y S_{max}=S_{xy} Smax=Sxy允许的最大尺寸
进程A：
A 1 = z m e d − z m i n A_{1}=z_{med}-z_{min} A1=zmed−zmin
A 2 = z m e d − z m a x A_{2}=z_{med}-z_{max} A2=zmed−zmax
如果 A 1 ＞ 0 A_{1}＞0 A1＞0且 A 2 ＜ 0 A_{2}＜0 A2＜0，则转到进程B，否则增大窗口尺寸
如果窗口尺寸≤ S m a x S_{max} Smax，则重复进程A
否则输出 z m e d z_{med} zmed
进程B：
B 1 = z x y − z m i n B_{1}=z_{xy}-z_{min} B1=zxy−zmin
B 1 = z x y − z m a x B_{1}=z_{xy}-z_{max} B1=zxy−zmax
B 1 ＞ 0 B_{1}＞0 B1＞0且 B 2 ＜ 0 B_{2}＜0 B2＜0，则输出 z x y z_{xy} zxy
否则输出 z m e d z_{med} zmed
算法主要目的：除去“椒盐”噪声（冲激噪声）；平滑其它非冲激噪声；减少物体边界细化或粗化等失真。
四、用频率域滤波消除周期噪声
4.1带阻滤波器
阻止一定频率范围内的信号通过而允许其它频率范围内的信号通过，消除或衰减傅里叶变换原点处的频段。
4.2带通滤波器
允许一定频率范围内的信 号通过而阻止其它频率范围内的信号通过。
H B P ( u , v ) = 1 − H B R ( u , v ) H_{BP}(u,v)=1-H_{BR}(u,v) HBP(u,v)=1−HBR(u,v)
H B P ( u , v ) H_{BP}(u,v) HBP(u,v)表示带通滤波器， H B R ( u , v ) H_{BR}(u,v) HBR(u,v)表示相应的带阻滤波器。
4.3陷波滤波器
阻止或通过事先定义的中心频率邻域内的频率。由于傅里叶变换是对称的，陷波滤波器必须以关于原点对称的形式出现。如果陷波滤波器位于原点处，则以它本身形式出现。
传递函数由下式给出：
H N P ( u , v ) = 1 − H N R ( u , v ) H_{NP}(u,v)=1-H_{NR}(u,v) HNP(u,v)=1−HNR(u,v)
式中， H N P ( U , V ) H_{NP}(U,V) HNP(U,V)是陷波滤波器的传递函数，这个陷波滤波器与传递函数为 H N R ( U , V ) H_{NR}(U,V) HNR(U,V)的陷波滤波器相对应。
五、退化函数
当噪声被处理掉之后，我们将图像复原这一过程近似的看为如下方式：
H s ( u , v ) = G s ( u , v ) F ^ s ( u , v ) H_{s}(u,v) = \frac{G_{s}(u,v)}{\hat F_{s}(u,v)} Hs(u,v)=F^s(u,v)Gs(u,v)
六、维纳滤波
通过尽可能使重建图像和原图像之间的均方根最小来实现图像的复原。均方误差的计算方式如下：
e 2 = E { ( f − f ^ ) 2 } e^2=E\{(f-\hat f)^2 \} e2=E{(f−f^)2}
其中， E { ⋅ } E\{\cdot\} E{⋅}是参数的期望值。
因此误差函数的最小值在频率域中为如下表示：
F ^ ( u , v ) = [ H ∗ ( u , v ) ∣ H ( u , v ) ∣ 2 + S η ( u , v ) / S f ( u , v ) ] G ( u , v ) \hat F(u,v)=[\frac{H^*(u,v)}{|H(u,v)|^2+S_\eta(u,v)/S_f(u,v)}]G(u,v) F^(u,v)=[∣H(u,v)∣2+Sη(u,v)/Sf(u,v)H∗(u,v)]G(u,v)
H ∗ ( u , v ) H^*(u,v) H∗(u,v)是 H ( u , v ) H(u,v) H(u,v)的复共轭， ∣ H ( u , v ) ∣ 2 = H ∗ ( u , v ) H ( u , v ) |H(u,v)|^2=H^*(u,v)H(u,v) ∣H(u,v)∣2=H∗(u,v)H(u,v)。
比例式 S η ( u , v ) / S f ( u , v ) S_\eta(u,v)/S_f(u,v) Sη(u,v)/Sf(u,v)被称为信噪功率比。在频率域中用下式来表示：
S N R = ∑ u = 0 M − 1 ∑ v = 0 N − 1 ∣ F ( u , v ) ∣ 2 ∑ u = 0 M − 1 ∑ v = 0 N − 1 ∣ N ( u , v ) ∣ 2 SNR = \frac{\sum_{u=0}^{M-1}\sum_{v=0}^{N-1}|F(u,v)|^2}{\sum_{u=0}^{M-1}\sum_{v=0}^{N-1}|N(u,v)|^2} SNR=∑u=0M−1∑v=0N−1∣N(u,v)∣2∑u=0M−1∑v=0N−1∣F(u,v)∣2
七、小结
图像复原过程可以分为两步，第一步根据可能存在的噪声类型采用相应的滤波器减弱噪声的影响。第二步，采用与退化函数相对应的复原函数实现图像的复原。