[TOC]
数字图像在获取的过程中，由于光学系统的像差、 光学成像衍射、 成像系统的非线性畸变、 摄影胶片的感光的非线性、 成像过程的相对运动、 大气的湍流效应、环境随机噪声等原因， 图像会产生一定程度的退化。因此，必须采取一定的方法尽可能地减少或消除图像质量的下降，恢复图像的本来面目，这就是图像复原， 也称为图像恢复。
图像复原与图像增强有类似的地方， 都是为了改善图像。但是它们又有着明显的不同。图像复原是试图利用退化过程的先验知识使已退化的图像恢复本来面目，即根据退化的原因， 分析引起退化的环境因素，建立相应的数学模型， 并沿着使图像降质的逆过程恢复图像。从图像质量评价的角度来看， 图像复原就是提高图像的可理解性。而图像增强的目的是提高视感质量，图像增强的过程基本上是一个探索的过程， 它利用人的心理状态和视觉系统去控制图像质量， 直到人们的视觉系统满意为止。
图像复原是利用退化现象的某种先验知识，建立退化现象的数学模型，再根据模型进行反向的推演运算，以恢复原来的景图像。因而，图像复原可以理解为图像降质过程的反向过程。建立图像复原的反向过程的数学模型，就是图像复原的主要务。经过反向过程的数学模型的运算，要想恢复全真的景物图像比较困难。所以， 图像复原本身往往需要有一个质量准， 即衡量接近全真景物图像的程度，或者说，对原图像的估计是否到达最佳的程度。
由于引起退化的因素众多而且性质不同，为了描述图像退化过程所建立的数学模型往往多种多样，而恢复的质量标准也往存在差异性，因此图像复原是一个复杂的数学过程，图像复原的方法、技术也各不相同。
图像复原处理的关键问题在于建立退化模型。输入图像f(x, y)经过某个退化系统后输出的是一幅退化的图像。为了讨论方便， 把噪声引起的退化即噪声对图像的影响一般作为加性噪声考虑， 这也与许多实际应用情况一致，如图像数字化时的量化噪声、 随机噪声等就可以作为加性噪声，即使不是加性噪声而是乘性噪声， 也可以用对数方式将其转化为相加形式。
原始图像 $f(x, y)$ 经过一个退化算子或退化系统 $H(x, y)$ 的作用, 再和噪声 $n(x, y)$ 进行叠加, 形成退化后的图像 $g(x, y)$ 。图11-1表示退化过程的输入和输出的关系, 其中 $H(x, y)$ 概括了退化系统的物理过程, 就是所要寻找的退化数学模型。
此处是进行了简化，因为乘法噪声可以转化为加法噪声
数字图像的图像恢复问题可看作是: 根据退化图像 $g(x, y)$ 和退化算子 $H(x, y)$ 的形式, 沿着反向过程去求解原始图像 $f(x, y)$, 或者说是逆向地寻找原始图像的最佳近似估计。图像退化的过程可以用数学表达式写成如下的形式:
$$
g(x, y)=H[f(x, y)]+n(x, y)\tag{5-1}
$$
在这里, $n(x, y)$ 是一种统计性质的信息。在实际应用中, 往往假设噪声是白噪声, 即它的频谱密度为常数, 并且与图像不相关。
在图像复原处理中， 尽管非线性、 时变和空间变化的系统模型更具有普遍性和准确性，更与复杂的退化环境相接近，但它给实际处理工作带来了巨大的困难， 常常找不到解或者很难用计算机来处理。因此，在图像复原处理中， 往往用线性系统和空间不变系统模型来加以近似。这种近似的优点使得线性系统中的许多理论可直接用于解决图像复原问题，同时又不失可用性。
一幅连续图像 $f(x, y)$ 可以看作是由一系列点源组成的。因此, $f(x, y)$ 可以通过点源函数的卷积来表示。即
$$
f(x, y)=\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} f(\alpha, \beta) \delta(x-\alpha, y-\beta) d \alpha d \beta\tag{5-2}
$$
式中, $\delta$ 函数为点源函数, 表示空间上的点脉冲。
在不考虑噪声的一般情况下, 连续图像经过退化系统 $H$ 后的输出为
$$
\quad g(x, y)=H[f(x, y)]\tag{5-3}
$$
把式(5-2)代入式(5-3)得
$$
\begin{gathered}
g(x, y)=H[f(x, y)]=H\left[\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} f(\alpha, \beta) \delta(x-\alpha, y-\beta) d \alpha d \beta\right]
\end{gathered}\tag{5-4}
$$
在线性和空间不变系统的情况下, 退化算子 $H$ 具有如下性质:
(1) 线性: 设 $f_1(x, y)$ 和 $f_2(x, y)$ 为两幅输入图像, $k_1$ 和 $k_2$ 为常数, 则
$$
H\left[k_1 f_1(x, y)+k_2 f_2(x, y)\right]=k_1 H\left[f_1(x, y)\right]+k_2 H\left[f_2(x, y)\right]\tag{5-5}
$$
由该性质还可推出下面两个结论:
(1) 当 $k_1=k_2=1$ 时, 式(11-5)变为
$$
H\left[f_1(x, y)+f_2(x, y)\right]=H\left[f_1(x, y)\right]+H\left[f_2(x, y)\right]\tag{5-6}
$$
(2) 如果 $f_2(x, y)=0$, 则式 (11-5) 变为
$$
H\left[k_1 f_1(x, y)\right]=k_1 H\left[f_1(x, y)\right]\tag{5-7}
$$
(2) 空间不变性: 如果对任意 $f(x, y)$ 以及 $a$ 和 $b$, 有
$$
H[f(x-a, y-b)]=g(x-a, y-b)\tag{5-8}
$$
则对于线性空间不变系统, 输入图像经退化后的输出为
$$
\begin{aligned}
g(x, y) & =H[f(x, y)]=H\left[\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} f(\alpha, \beta) \delta(x-\alpha, y-\beta) d \alpha d \beta\right] \\
& =\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} f(\alpha, \beta) H[\delta(x-\alpha, y-\beta)] d \alpha d \beta \\
& =\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} f(\alpha, \beta) h(x-\alpha, y-\beta) d \alpha d \beta
\end{aligned}\tag{5-9}
$$
式中, $h(x-\alpha, y-\beta)$ 为该退化系统的点扩展函数, 或叫系统的冲激响应函数。它表示系统对坐标为 $(a, \beta)$ 处的冲激函数 $\delta(x-\alpha, y-\beta)$ 的响应。也就是说, 只要系统对冲激函数的响应为已知, 那么就可以清楚图像退化是如何形成的。因为对于任一输入 $f(a, \beta)$ 的响应, 都可以通过上式计算出来。
此时, 退化系统的输出就是输入图像信号 $f(x, y)$ 与点扩展函数 $h(x, y)$ 的卷积, 即
$$
g(x, y)=\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} f(\alpha, \beta) h(x-\alpha, y-\beta) d \alpha d \beta=f(x, y) * h(x, y)\tag{5-10}
$$
图像退化除了受到成像系统本身的影响外, 有时还要受到噪声的影响。假设噪声 $n(x, y)$ 是加性白噪声, 这时上式可写成
$$
\begin{aligned}
g(x, y) & =\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} f(\alpha, \beta) h(x-\alpha, y-\beta) d \alpha d \beta+n(x, y) \\
& =f(x, y) * h(x, y)+n(x, y)
\end{aligned}\tag{5-11}
$$
在频域上, 式(11-11)可以写成
$$
G(u, v)=F(u, v) H(u, v)+N(u, v)\tag{5-12}
$$
其中, $G(u, v) 、 F(u, v) 、 N(u, v)$ 分别是退化图像 $g(x, y)$ 、原图像 $f(x, y)$、噪声信号 $n(x, y)$ 的傅立叶变换; $H(u, v)$ 是系统的点冲激响应函数 $h(x, y)$ 的傅立叶变换, 称为系统在频率域上的传递函数。
式(5-11)和式(5-12)就是连续函数的退化模型。可见, 图像复原实际上就是已知 $g(x, y)$ 求 $f(x, y)$ 的问题或已知 $G(u, v)$ 求 $F(u, v)$ 的问题, 它们的不同之处在于一个是在空域,一个是在频域。
显然, 进行图像复原的关键问题是寻找降质系统在空间域上的冲激响应函数 $h(x, y)$, 或者降质系统在频率域上的传递函数 $H(u$, $v)$ 。一般来说, 传递函数比较容易求得。因此, 在进行图像复原之前, 一般应设法求得完全的或近似的降质系统传递函数, 要想得到 $h(x, y)$, 只需对 $H(u, v)$ 求傅立叶逆变换即可。
- 一维离散退化模型
设 $f(x)$ 为具有 $A$ 个采样值的离散输入函数, $h(x)$ 为具有 $B$ 个采样值的退化系统的冲激响应函数, 则经退化系统后的离散输出函数 $g(x)$ 为输入 $f(x)$ 和冲激响应 $h(x)$ 的卷积, 即
$$
g(x)=f(x)^* h(x)
$$
为了避免上述卷积所产生的各个周期重叠（设每个采样函数的周期为 $M)$ ，分别对 $f(x)$ 和 $h(x)$ 用添零延伸的方法扩展成周期 $M=A+B-1$ 的周期函数, 即
$$
f_e(x)=\left{\begin{array}{cc}
f(x) & 0 \leq x \leq A-1 \
0 & A \leq x \leq M-1
\end{array}\right. \tag{5-13a}
$$
$$
h_e(x)=\left{\begin{array}{cc}
h(x) & 0 \leq x \leq B-1 \
0 & B \leq x \leq M-1
\end{array}\right.\tag{5-13b}
$$
输出为
$$
g_e(x)=f_e(x) * h_e(x)=\sum_{m=0}^{M-1} f_e(m) h_e(x-m)\tag{5-14}
$$
式中, $x=0,1,2, \ldots, M-1$ 。
因为 $f_{e}(x)$ 和 $h_{e}(x)$ 已扩展成周期函数, 故 $g_{e}(x)$ 也是周期性函数, 用矩阵表示为
$$
\left[\begin{array}{c}
g(0) \
g(1) \
g(2) \
M \
g(M-1)
\end{array}\right]=\left[\begin{array}{cccc}
h_e(0) & h_e(-1) & \Lambda & h_e(-M+1) \
h_e(1) & h_e(0) & \Lambda & h_e(-M+2) \
h_e(2) & h_e(1) & \Lambda & h_e(-M+3) \
M & M & M & M \
h_e(M-1) & h_e(M-2) & \Lambda & h_e(0)
\end{array}\right]\left[\begin{array}{c}
f_e(0) \
f_e(1) \
f_e(2) \
M \
f_e(M-1)
\end{array}\right]\tag{5-15}
$$
因为 $h_{e}(x)$ 的周期为 $M$, 所以 $h_{e}(x)=h_{e}(x+M)$, 即
$$
\begin{gathered}
h_e(-1)=h_e(M-1) \
h_e(-2)=h_e(M-2) \
h_e(-3)=h_e(M-3) \
M \
h_e(-M+1)=h_e(1)
\end{gathered}
$$
$M \times M$ 阶矩阵 $\boldsymbol{H}$ 可写为
$$
H=\left[\begin{array}{cccc}
h_e(0) & h_e(-1) & \Lambda & h_e(1) \
h_e(1) & h_e(0) & \Lambda & h_e(2) \
h_e(2) & h_e(1) & \Lambda & h_e(3) \
M & M & \Lambda & \Lambda \
h_e(M-1) & h_e(M-2) & h_e(0)
\end{array}\right]\tag{5-16}
$$
可将式（5-15）写成更简洁的形式，即
$$
g=H f\tag{5-17}
$$
式中, $\boldsymbol{g} 、 \boldsymbol{f}$ 都是 $M$ 维列向量, $\boldsymbol{H}$ 是 $M \times M$ 阶矩阵, 矩阵中的每一行元素均相同, 只是每行以循环方式右移一位, 因此矩阵 $\boldsymbol{H}$ 是循环矩阵。循环矩阵相加或相乘得到的还是循环矩阵。
- 二维离散模型
设输入的数字图像 $f(x, y)$ 大小为 $A \times B$, 点扩展函数 $h(x, y)$ 被均匀采样为 $C \times D$ 大小。为避免交叠误差, 仍用添零扩展的方法, 将它们扩展成 $M=A+C-1$ 和 $N=B+D-1$ 个元素的周期函数。
$$
f_e(x, y)=\left{\begin{array}{cc}
f(x, y) & 0 \leq x \leq A-1 \text { 且 } 0 \leq y \leq B-1 \\
0 & \text { 其他 }
\end{array}\right. \tag{5-18a}
$$
$$
h_e(x, y)=\left{\begin{array}{cc}
h(x, y) & 0 \leq x \leq C-1 \text { 且 } 0 \leq y \leq D-1 \\
0 & \text { 其他 }
\end{array}\right.\tag{5-18b}
$$
则输出的降质数字图像为
$$
g_e(x, y)=\sum_{m=0}^{M-1} \sum_{n=0}^{N-1} f_e(m, n) h_e(x-m, y-n)=f(x, y)^* h(x, y)\tag{5-19}
$$
式中: $x=0,1,2, \ldots, M-1 ; y=0,1,2, \ldots, N-1$ 。
式(5-19)的二维离散退化模型同样可以用式 (5-17) 所示的矩阵形式表示, 即
$$
g=H f
$$
式中, $g 、 f$ 是 $M N \times 1$ 维列向量, $H$ 是 $M N \times M N$ 维矩阵。其方法是将 $g(x, y)$ 和 $f(x, y)$ 中的元素排成列向量。
$$
f_1=\left[\underbrace{f_1(0,0),f_1(0,1),\cdots,f_1(0,N-1)}{第0行元素},\underbrace{f_1(1,0),f_1(1,1),\cdots,f_1(1,N-1)}{第1行元素},\cdots,\underbrace{f_1(M-1,0),f_1(M-1,1),\cdots,f_1(M-1,N-1)}{第M-1行元素}\right]^T\
g_1=\left[\underbrace{g_1(0,0),g_1(0,1),\cdots,g_1(0,N-1)}{第0行元素},\underbrace{g_1(1,0),g_1(1,1),\cdots,g_1(1,N-1)}{第1行元素},.\cdots,\underbrace{g_1(M-1,0),g_1(M-1,1),\cdots,g_1(M-1,N-1)}{第M-1行元素}\right]^T
\tag{5-20}
$$
$$
H=\left[\begin{array}{cccc}
H_0 & H_{M-1} & H_{M-2} & \cdots& H_1 \
H_1 & H_0 & H_{M-1} & \cdots& H_2 \
\vdots & \vdots & \vdots & \ddots&\vdots \
H_{M-1} & H_{M-2} & H_{M-3} & \cdots &H_0
\end{array}\right]
\tag{5-21}
$$
$H_i(i=0,1,2, \ldots, M-1)$ 为子矩阵, 大小为 $N \times N$, 即 $H$ 矩阵由 $M \times M$ 个大小为 $N \times N$ 的子矩阵组成, 称为分块循环矩阵。分块矩阵是由延拓函数 $h_{e}(x, y)$ 的第 $j$ 行构成的, 构成方法如下:
$$
H=\left[\begin{array}{cccc}
h_e(j, 0) & h_e(j, N-1) & h_e(j, N-2) & \cdots& h_e(j, 1) \
h_e(j, 1) & h_e(j, 0) & h_e(j, N-1) & \cdots &h_e(j, 0) \
\vdots & \vdots & \vdots &\ddots &\vdots \
h_e(j, N-1) & h_e(j, N-2) & h_e(j, N-3) & \cdots &h_e(j, 0)
\end{array}\right]
\tag{5-22}
$$
若把噪声考虑进去, 则离散图像退化模型为
$$
g_e(x, y)=\sum_{m=0}^{M-1} \sum_{n=0}^{N-1} f_e(m, n) h_e(x-m, y-n)+n_e(x, y)
\tag{5-23}
$$
写成矩阵形式为
$$
g=H f+n\tag{5-24}
$$
上述线性空间不变退化模型表明, 在给定了 $g(x, y)$, 并且知道退化系统的点扩展函数 $h(x, y)$ 和噪声分布 $n(x, y)$ 的情况下, 可估计出原始图像 $f(x, y)$ 。
假设图像大小 $M=N=512$, 相应矩阵 $M H$ 的大小为 $M N \times M N=262144 \times 262144$, 这意味着要解出 $f(x, y)$ 需要解262 144 个联立方程组, 其计算量十分惊人。
图像恢复目的是在给定的 $g(x, y), h(x, y), n(x, y)$ 情况下, 求 $f(x, y)$. 由于 $G , F$ 维数很大, 通过解方程组不可能进行图像恢复.
- 代数复原方法:寻找一个估计 $f(x, y)$.
由式(5-24)可得
$$
n=g-H f\tag{5-25}
$$
在对 $n$ 没有先验知识的情况下, 可以依据这样的最优准则, 即寻找一个 $\hat{f}$, 使得 $H \hat{f}$ 在最小二乘方误差的意义下最按近 $g$, 即要使 $n$ 的模或范数 (norm) 最小:
$$
|n|^2=n^T n=|g-H \hat{f}|^2=(g-H \hat{f})^T(g-\hat{H f})\tag{5-26}
$$
即优化目标为:
$$
L(\hat{f})=|g-H \hat{f}|^2\tag{5-27}
$$
如果我们在求最小值的过程中, 不做任何约束, 称这种复原为非约束复原。由极值条件
$$
\frac{\partial L(\hat{f})}{\partial \hat{f}}=0 \Rightarrow H^T(g-H \hat{f})=0\tag{5-28}
$$
解出 $\hat{f}$ 为
$$
\hat{f}=\left(H^T H\right)^{-1} H^T g\tag{5-29}
$$
非约束复原是指除了使准则函数 $L(\hat{f})=|g-\hat{H f}|^2$ 最小外, 再没有其他的约束条件。因此只需了解降质系统的传递函数或点扩展函数, 就能利用如前所述的方法进行复原。但是由于传递函数存在病态问题, 复原只能局限在靠近原点的有限区域内进行, 这使得非约束图像复原具有相当大的局限性。
最小二乘类约束复原是指除了要求了解关于退化系统的传递函数之外, 还需要知道某些噪声的统计特性或噪声与图像的某些相关情况。根据所了解的噪声的先验知识的不同, 采用不同的约束条件, 可得到不同的图像复原技术。在最小二乘类约束复原中, 要设法寻找一个最优估计 $\hat{f}$, 使得形式为 $|Q \hat{f}|^2=|n|^2$ 的函数最小化。求这类问题的最小化, 常采用拉格朗日乘子算法。也就说, 要寻找一个 $\hat{f}$, 使得优化目标函数:
$$
J(\hat{f})=|\hat{O \hat{f}}|^2+\alpha\left(|g-H \hat{f}|^2-|n|^2\right)\tag{5-30}
$$
为最小。式中, $Q$ 为 $\hat{f}$ 的线性算子, $\alpha$ 为一常数, 称为拉格朗日乘子。对式(5-30)求导得
$$
\begin{aligned}
& \frac{\partial J(\hat{f})}{\partial f}=0 \
& Q^T Q \hat{f}-\alpha H^T(g-H \hat{f})=0
\end{aligned}
$$
求解 $\hat{f}$ 得到
$$
\hat{f}=\left(H^T H+\gamma Q^T Q\right)^{-1} H^T g\tag{5-31}
$$
只能接近，无法完全恢复
1.逆滤波基本原理
由退化图像的付里叶变换求得图像的付里叶变换的估计方法称为逆滤波.
逆滤波
$$
\begin{aligned}
\overline{g(x, y)}= & \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x, A) h(x-\imath v, y-A) d x d, \partial+n(x, y) \
& G(u, v)=F(u, v) H(u, v)+N(u, v)
\end{aligned}
$$
不考虑噪声:
$$
\begin{aligned} F(u, v) & =\frac{G(u, v)}{H(u, v)} \ f(x, y) & =\mathcal{G}^{-1}[F(u, v)]\end{aligned}
$$
有噪声:
$$
\begin{aligned}
& F(u, v)=\frac{G(u, v)}{H(u, v)}-\frac{N(u, v)}{H(u, v)} \\
& f(x, y)=\xi^{-1}[F(u, v)]
\end{aligned}
$$
图5.2 逆滤波
逆滤波
2.去除由均匀直线运动引起的模糊
- 模糊模型
$$
g(x, y)=\int_0^T f\left[x-x_0(t), y-y_0(t)\right] d t
$$
$f(x, y)$ : 原图像; $g(x, y)$ : 模糊图像; $T$ : 暴光时间
- 图像恢复 $G(u, v)=H(u, v) F(u, v)$
其中: $H(u, v)=\int_0^T \exp \left{-j 2 \pi\left[u x_0(t)+v y_0(t)\right]\right} d t$
在实验室中，可以反复做实验的情景
重点
知道图像，知道有关分布特征
1.最小二乘方滤波
设原始图像, 退化图像和噪声分别为 $f(x, y), g(x, y)$ 和 $n(x, y)$, 则 :
$$
g(x, y)=H[f(x, y)]=H\left[\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} f(\alpha, \beta) \delta(x-\alpha, y-\beta) d \alpha d \beta\right]+n(x, y)$ 在给定的 $g(x, y)$ 时, 找出 $f(x, y)$ 的估计 $f(x, y)
$$
使：
$$
e^2=E\left{[f(x, y)-\hat{f}(x, y)]^2\right}
$$
最小.
设
$$
\quad \hat{f}(x, y)=\iint m(x, y, \alpha, \beta) g(\alpha, \beta) d \alpha d \beta
$$
2. 求解
如果 $f(x, y), g(x, y)$ 和 $n(x, y)$ 是均匀随机场, 则 :
$$
\hat{f}(x, y)=\iint m(x-\alpha, y-\beta) g(\alpha, \beta) d \alpha d \beta
$$
因此图像恢复就是确定m $m, y)$. 通过求解得:
$$
\begin{aligned}
M(u, v) & =\left[\frac{H^*(u, v)}{|H(u, v)|^2+\left[S_n(u, v) / S_f(u, v)\right]}\right] \
& =\left[\frac{1}{H(u, v)} \cdot \frac{|H(u, v)|^2}{\left|H(u, v)^2\right|+\left[S_n(u, v) / S_f(u, v)\right]}\right]
\end{aligned}
$$
$H(u, v)$ 是的h $(x, y)$ 付里叶变换.
$S_f, S_n$ 分别是 $f(x, y), n(x, y)$ 的谱密度
于是:
$$
\begin{aligned}
\hat{F}(u, v) & =\left[\frac{H^(u, v)}{|H(u, v)|^2+\left[S_n(u, v) / S_f(u, v)\right]}\right] G(u, v) \
& =\left[\frac{1}{H(u, v)} \cdot \frac{|H(u, v)|^2}{\left|H(u, v)^2\right|+\left[S_n(u, v) / S_f(u, v)\right]}\right] G(u, v)
\end{aligned}\tag{5-32}
$$
式中，
$$
u, v=0,1,2, \ldots, N-1,|H(u, v)|^2=H^(u, v) H(u, v)
$$
- 讨论
-
$S_n=0$ 时, 是理想的逆滤波器
-若随机过程的统计特性不知道, 可采用:
$$
M(u, v)=\frac{1}{H(u, v)} \frac{|H(u, v)|^2}{|H(u, v)|^2+\Gamma}
$$
式中 $\Gamma$ 是噪声对信号的功率密度比.
-
维纳滤波器是假设线性系统, 实际上图像的记录和评价图像的视觉系统都是非线性系统.
-
维纳滤波器是根据最小均方差设计的，这个准则不见得与人的视觉判决准则相符.
-
维纳滤波器是基于平稳随机过程的模型,实际图像不一定符合这一条件.
（仅供了解）
1.图像功率谱滤波器
$$
H_{R}(u, v)=\left[\frac{W_{FI}(u, v)}{\left|H_{D}(u, v)\right|^2 W_{F 1}(u, v)+W_{N}(u, v)}\right]^{\frac{1}{2}}
$$
式中 $H_{D}(u, v)$ 是图像退化的传递函数。 $W_{F 1}(u, v)$ 代表滤波器输出功率频谱,
$$
W_{FI}(u, v)=\left|H_{D}(u, v)\right|^2 W_{FO}(u, v)
$$
式中 $W_{FO}(u, v)$ 代表观测的功率频谱。它与理想图像的功率频谱的关系是
$$
W_{FO}(u, v)=\left|H_{D}(u, v)\right|^2 \cdot W_{FI}(u, v)+W_{N}(u, v)
$$
式中 $W_{N}(u, v)$ 是噪声功率频谱。
2.几何平均滤波器.
几何平均滤波器的传递函数由下式表示:
$$
H_{R}(u, v)=\left[H_{D}(u, v)\right]^{-s}\left[\frac{H_{D}^*(u, v) W_{FI}(u, v)}{\left|H_{D}(u, v)\right|^2 W_{FI}(u, v)+W_{N}(u, v)}\right]^{1-s}
$$
S是一个参数
3.约束最小方滤波器
约束最小方滤波器的传递函数如下:
$$
H_{\mathrm{R}}(u, v)=\frac{H_{\mathrm{D}}^*(u, v)}{\left|H_{\mathrm{D}}(u, v)\right|^2+r|L(u, v)|^2}
$$
$r$是一个参数，$L(u,v)$是一个频率参数。$r=1$为维纳滤波器
最大后验复原是一种统计方法, 它把原图像 $f(x, y)$ 和退化图像 $g(x, y)$ 都作为随机场, 在已知 $g(x, y)$ 的前提下, 求出后验条件概率密度函数 $P(f(x, y) / g(x, y))$ 。若 $\hat{f}(x, y)$ 使式
$$
\max _{f} p(f \mid g)=\max _f p(g \mid f) p(f) / p(g)=\max _f p(g \mid f) p(f)\tag{5-33}
$$
最大, 则 $\hat{f}(x, y)$ 就代表已知退化图像 $g(x, y)$ 时, 最可能的原始图像 $f(x, y)$ 。这种图像方法称之为最大后验图像复原方法。
最大后验图像复原方法把图像看作是非平稳随机场, 把图像模型表示成一个平稳随机过程对于一个不平稳的均值作零均值Gauss起伏, 可得出求解迭代序列:
$$
\hat{f}_{k+1}=\hat{f}_k-h^* \sigma_n^{-2}\left(g-h^* \hat{f}_k\right)-\sigma_n^{-2}\left(\hat{f}_k-\bar{f}\right)\tag{5-34}
$$
其中, $k$ 为迭代次数, $\sigma_n^{-2}$ 和 $\sigma_f^{-2}$ 分别为 $f$ 和 $n$ 的方差的倒数, $f$ 是随空间而变的均值, 它是一个常数, 但要经过多次迭代才能收敛到最后的解。
最大熵复原方法是通过最大化某种反映图像平滑性的准则函数来作约束条件, 以解决图像复原中反向滤波法存在的病态问题。
首先简单介绍一下熵的概念。
熵的定义为
$$
H=-\int_{-\infty}^{+\infty} P(x) \ln P(x) d x\tag{5-35}
$$
式中, $P(x)$ 为随机变量 $x$ 的概率密度。
对于离散信号, 熵的定义为
$$
H=-\sum_{k=1}^M P(k) \ln P(k)\tag{5-36}
$$
熵是表征随机变量集合的随机程度的量度。当所有随机变量等可能性时, 也就是说 $P_1=P_2=\ldots=P_{m}$ 时熵最大, 为 $H=\ln M$ 。由于概率 $P(k)$ 介于 $0 \sim 1$ 之间, 因此最大熵的范围为 $0 \sim \ln M, H$ 不可能出现负值。
在二维数字图像中, 熵的定义为
$$
H_f=-\sum_{m=1}^M \sum_{n=1}^N f(m, n) \ln f(m, n)\tag{5-37}
$$
最大熵复原的原理是将 $f(x, y)$ 写成随机变量的统计模型，然后在一定的约束条件下，找出用随机变量形式表示的熵的表达式, 运用求极大值的方法, 求得最优估计解 $\hat{f}(x, y)$ 。最大熵复原的含义是对 $\hat{f}(x, y)$ 的最大平滑估计。 最大熵复原常用Friend和Burg两种方法, 这两种方法基本原理相同, 这里仅介绍Friend法。
首先定义一幅大小为 $M \times N$ 的图像 $f(x, y)$, 显然 $f(x, y)$ 非负。 图像的总能量 $E$ 和熵分别为
$$
E=\sum_{i=j}^M \sum_{j=1}^N f\left(x_i, y_j\right)\tag{5-38}
$$
和
$$
H_f=\sum_{i=j}^M \sum_{j=1}^N f\left(x_i, y_j\right) \operatorname{lnf}\left(x_i, y_j\right)\tag{5-39}
$$
类似地可定义噪声的熵 $H_n$ :
$$
H_n=-\sum_{i=1}^M \sum_{j=1}^N n^{\prime}\left(x_i, y_j\right) 1 n n^{\prime}\left(x_i, y_j\right)\tag{5-40}
$$
式中: $n^{\prime}(x, y)=n(x, y)+B, B$ 为最大噪声负值。
恢复就是在满足式(5-38) 和图像退化模型的约束条件下, 使恢复后的图像熵和噪声熵达到最大。熵通常取决于 $f$ 的形状, 当图像具有均匀的灰度时熵最大。因此用最大熵恢复图像具有某种平滑性。
引入拉格朗日 (Lagrange) 函数
$$
\begin{aligned}
R & =H_f+\rho H_n+\sum_{m=1}^N \sum_{n=1}^N \lambda_{m m}\left{\sum_{x=1}^N \sum_{n=1}^N h(m-x, n-y) f(x, y)+n_T(m, n)-B-g(m, n)\right. \
& +\beta\left{\sum_{x=1}^N \sum_{y=1}^N f(x, y)-E\right}
\end{aligned}\tag{5-41}
$$
式中: $\lambda_{m n}(m, n=1,2, \ldots, N)$ 和 $\beta$ 是拉格朗日乘子， $\rho$ 是加权因子，表示 $H_f$ 和 $H_n$ 相互之间的权重。
$\hat{f}(x, y)$ 和 $\hat{n}(m, n)$ 分别表示 $f(x, y)$ 和 $n^{\prime}(x, y)$ 的估计值, 则有
$$
\frac{\partial R}{\partial f(x, y)}=0 \tag{5-42}
$$
$$
\frac{\partial R}{\partial n^{\prime}(x, y)}=0\tag{5-43}
$$
把式(5-41)分别代入式(5-42)和式(5-43), 可得
$$
\hat{f}(x, y)=\exp \left[-1+\beta+\sum_{m=1}^N \sum_{n=1}^N \lambda_{m m} h(m-x, n-y)\right] \quad x, y=1,2, \ldots, N \tag{5-44}
$$
$$
\hat{n}(m, n)=\exp \left(-1+\frac{\lambda_{m n}}{\rho}\right) \quad \quad m, n=1,2, \ldots, N\tag{5-45}
$$
并且$\hat{f}(x,y)$和$\hat{n}(m,n)$满足下列约束条件
$$
\sum_{x=1}^N \sum_{y=1}^N f(x, y)=E\tag{5-46}
$$
$$
\sum_{x=1}^N \sum_{y=1}^N h(m-x, n-y) f(x, y)+\hat{n}(m, n)-\beta=g(m, n)\tag{5-47}
$$
式(5-44)为图像恢复函数。把式(5-44)和式(5-45)代入式(5$46)$ 和式(5-47)可得 $\left(N^2+1\right)$ 个方程。由此联立方程组可解得 $\left(N^2+1\right)$ 个末知数解 $\lambda_{m n}(m, n=1,2, \ldots, N)$, 解上述方程组可求得 $\hat{f}(x, y)$ 的值。
投影复原法是用代数方程组来描述线性和非线性退化系统的。该系统可用下式描述：
$$
g(x, y)=H[f(x, y)]+n(x, y)
$$
其中: $f(x, y)$ 是原始图像, $g(x, y)$ 是退化图像, $n(x, y)$ 是系统噪声, $H$ 是退化算子，表示对图像进行某种运算。
在使用投影复原法进行图像复原时, 引进一些先验信息附加的约束条件, 可改善图像复原效果。
$$
g(x, y)=H[f(x, y)]+n(x, y)
$$
投影复原方法就是迭代方法之一。设初始估值 $f^0(x, y)$, 然后进行迭代运算.
设 $f^0(x, y)=g(x, y)$
设 $H$ 是非线性的, 并忽略噪声, 则上式写成:
$$
\begin{array}{ll}
a_{11} f_1+a_{12} f_2+\ldots & a_{1 N} f_N=g_1 \
a_{21} f_1+a_{22} f_2+\ldots & a_{2 N} f_N=g_2 \
\ldots & \
a_{M 1} f_1+a_{M 2} f_2+\ldots & a_{M N} f_N=g_M
\end{array}\tag{5-48}
$$
$$
f=\left[f_1 \ldots f_N\right]
$$
取
$$
f^{(0)}=\left[f_1^{(0)}, f_2^{(0)}, \ldots f_N^{(0)}\right]
$$
$$
f^{(0)}=\left[g_1, g_2, \ldots g_N\right]
$$
$f^{(1)}$ 取 $f^{(0)}$ 在第一个超平面$a_{11} f_1+a_{12} f_2+\ldots \quad a_{1 N} f_{N}=g_1$上的投影。
即:
$$
f^{(1)}=f^{(0)}-\frac{f^{(0)} \cdot a_1-g_1}{a_1 \bullet a_2} a_1\tag{5-49}
$$
$f^{(2)}$ 取 $f^{(1)}$ 在第二个超平面$a_{21} f_1+a_{22} f_2+\ldots \quad a_{2 N} f_{N}=g_2$ 上的投影。
依次继续下去, 直到满足 $( 5 . 4 8 )$ 式中最后一个方程式。这就实现了迭代的第一个循环然后再从 (5.48)式中第一个方程式开始第二次迭代。即取 $f^{(M)}$ 在第一个超平面$a_{11} f_1+a_{12} f_2+\ldots a_{1 N} f_{N}=g_1$上的投影。并称之为 $f^{(M+1)}$ 。依次继续, 直到
$$
f^{(k+1)}=f^{(k)}\tag{5-50}
$$
则:
$$
f^{(k)}=f\tag{5-51}
$$
数字图像在获取过程中，由于成像系统的非线性，成像后的图像与原景物图像相比，会产生比例失调，甚至扭曲，我们把这类图像退化现象称之为几何畸变。典型的几何失真如图5-4所示。
![image-20230426164500568]()
图5-4 几种典型的几何失真
(a) 原图像； (b) 梯形失真； (c) 枕形失真； (d) 桶形失真
一般，几何畸变校正要对失真的图像进行精确的几何校正， 通常是先确定一幅图像为基准，然后去校正另一幅图像的几何形状。因此，几何畸变校正一般分两步来做：第一步是图像空间坐标的变换；第二步是重新确定在校正空间各像素点的取值。
按照一幅标准图像 $g(u, v)$ 或一组基准点去校正另一幅几何失真图像 $f(x, y)$, 称之为空间几何坐标变换。根据两幅图像的一些 已知对应点对(也称为控制点对)建立起函数关系式, 将失真图像的 $x-y$ 坐标系变换到标准图像 $u-v$ 坐标系, 从而实现失真图像按标准图像的几何位置校正, 使 $f(x, y)$ 中的每一像点都可在 $g(u, v)$ 中找到对应像点。
图像的几何失真虽然是非线性的，但在一个局部小区域内可近似认为是线性的，基于这一假设，将标准图像和被校正图像之间的对应点对划分成一系列小三角形区域，三角形顶点为三个控制点， 在三角形区内满足以下线性关系：
$$
\left{\begin{array}{l}
x=a u+b v+c \
y=d u+e v+f
\end{array}\right.\tag{5-52}
$$
若三对控制点在两个坐标系中的位置分别为 $\left(x_1, y_1\right) 、\left(x_2, y_2\right)$ 、 $\left(x_3, y_3\right)$ 和 $\left(u_1, v_1\right) 、\left(u_2, v_2\right) 、\left(u_3, v_3\right)$, 则可建立两级方程组:
$$
\left{\begin{array}{l}
x_1=a u_1+b v_1+c \
x_2=a u_2+b v_2+c \
x_3=a u_3+b v_3+c
\end{array}\right.\tag{5-53a}
$$
$$
\left{\begin{array}{l}
y_1=d u_1+e v_1+f \\
y_2=d u_2+e v_2+f \\
y_3=d u_3+e v_3+f
\end{array}\right.\tag{5-53b}
$$
解方程组(5-53)，可求出a, b, c, d, e, f六个系数。用式(5-53)可实现该三角形区内其他像点的坐标变换。对于不同的三角形控制区域，这六个系数的值是不同的。
三角形线性法简单，能满足一定的精度要求，这是因它是以局部范围内的线性失真去处理大范围内的非线性失真， 所以选择的控制点对越多，分布越均匀，三角形区域的面积越小， 则变换的精度越高。但是控制点过多又会导致计算量的增加， 因此需要综合考虑。
图像经几何位置校正后，在校正空间中各像点的灰度值等于被校正图像对应点的灰度值。一般校正后的图像某些像素点可能挤压在一起，或者分散开，不会恰好落在坐标点上，因此常采用内插法来求得这些像素点的灰度值。经常使用的方法有如下两种。
- 最近邻点法
最近邻点法是取与像素点相邻的4个点中距离最近的邻点灰度值作为该点的灰度值。如图5-5所示。显然，最近邻点法计算简单，但精度不高，同时校正后的图像亮度有明显的不连续性。
图5-5 最近邻点法
- 内插法
用像素点周围４个邻点的灰度值加权内插作为g(u0, v0)，如图5-6所示，设像素点(x0, y0)周围４个点为$(x_1, y_1), (x_1+1, y_1), (x_1, y_1+1), (x_1+1, y_1+1)$，则校正值为
$$
\begin{aligned}
g\left(u_0, v_0\right)= & (1-\alpha)(1-\beta) f\left(x_1, y_1\right)+\alpha(1-\beta) f\left(x_1+1, y_1\right) \\
& +(1-\alpha) \beta f\left(x_1, y_1+1\right)+\alpha \beta f\left(x_1+1, y_1+1\right)
\end{aligned}\tag{5-54}
$$
式中：
$$
\alpha=\left|x_0-x_1\right|, \beta=\left|y_0-y_1\right|
$$
图5-6 内插法几何校正
盲目图像复原法是在没有图像退化先验知识的情况下，对观察目标的多幅图像以某种方式抽出退化信息， 从而进行图像复原的方法。
加性噪声的模糊图像复原方法一般有两种： 直接测量法和间接估计法。
用直接测量法复原图像时，需要测量图像的模糊脉冲响应和噪声功率谱或协方差函数。在所观察的景物中，点光源能量往往直接指示出冲激响应；另外，图像边缘是否陡峭也能用来推测模糊冲激响应。在背景亮度相对恒定的区域内测量图像的协方差，可以估计出观测图像的噪声协方差函数。
间接估计法复原图像类似于多图像平均法处理。如对一个景物连续拍摄 $N$ 次, 每一次获取的图像用下式表示:
$$
g_i(x, y)=f(x, y)+n_i(x, y) \quad i=1,2, \ldots, N\tag{5-55}
$$
式中: $f(x, y)$ 是原始图像, $g_i(x, y)$ 是第次获取的图像, $n_i(x, y)$ 是第 $i$ 次的噪声函数。
$$
f(x, y)=\frac{1}{N} \sum_{i=1}^N g_i(x, y)-\frac{1}{N} \sum_{i=1}^N n_i(x, y)\tag{5-56}
$$
当 $N$ 很大时, 上式右边的噪声项的值趋于它的数学期望 $E{n(x$, $y)}$ 。一般情况下白色高斯噪声在所有 $(x, y)$ 上的数学期望等于零。 因此，原始图像为
$$
f(x, y)=\frac{1}{N} \sum_{i=1}^N g_i(x, y)\tag{5-57}
$$
盲目图像复原
若 $g_{i}(x, y)=f_{i}(x, y)^* h_{i}(x, y)$
则 $G_i(n, v)=F_i(u, v) H_{i}(u, v)$
$\ln G_{i}(n, v)=\ln \left[F_{i}(u, v)\right]+\ln \left[H_{i}(u, v)\right]$
如果帧间退化冲激响应不相关, 则
$$
\sum_{i=1}^M \ln \left[G_i(u, v)\right]=M \ln \left[F_i(u, v)\right]+\sum_{i=1}^M \ln \left[H_i(u, v)\right]
$$
当M很大时, 传递函数的对数和接近一恒定值
$$
\begin{gathered}
K_H=\lim {M \rightarrow \infty} \sum{i=1}^M \ln \left[H_i(u, v)\right] \quad \hat{F}i(u, v)=\exp \left{\frac{k_H(u, v)}{M} \prod{i=1}^M\left[G_i(u, v)\right]^{\frac{1}{M}}\right} \
\hat{f}(x, y)=J^{-1}\left[\hat{F_G}(u, v)\right]
\end{gathered}
$$
根据前面章节的讨论， 图像退化的模型为
$$
g(x, y)=f(x, y) * h(x, y)+n(x, y)
$$
那么，图像复原的过程可认为是已知$g(x,y),h(x,y),n(x,y)$的一些先验知识，求出$f(x,y)$。对于不同的退化函数和噪声性质， 可推导出前面介绍的一些图像复原方法。本节介绍实际中经常用到的逆滤波复原图像法。
图像的点扩展函数$h$为
$$
h=\frac{1}{49}
\left[
\begin{array}{lllllll}
1 &1& 1& 1& 1& 1& 1\
1 &1& 1& 1& 1& 1& 1\
1 &1& 1& 1& 1& 1& 1\
1 &1& 1& 1& 1& 1& 1\
1 &1& 1& 1& 1& 1& 1\
1 &1& 1& 1& 1& 1& 1\
1 &1& 1& 1& 1& 1& 1\
\end{array}
\right]
$$
首先使用退化函数$h$对图像进行模糊操作， 生成一幅退化的图像，退化系统为
$$
g=F^{-1}(F(h) \cdot F(f))
$$
上式的计算过程是，先求出图像和退化函数的傅立叶变换， 在频域相乘后，再按下式求逆傅立叶变换：
$$
\hat{f}=F^{-1}\left(\frac{G}{H}\right)=F^{-1}\left(\frac{F(g)}{F(h)}\right)
$$
图像复原的具体实现过程如下。
图11-7 图像复原菜单
在文档类中添加两个消息映射函数， 分别是OnRestoreBlur()和OnRestoreInvfilt()。OnRestorBlur()调用图像复原类CRestore的ConvBlur()函数生成一幅退化图像。OnRestoreInvfilt()调用图像复原类CRestore的InvFilter()函数对退化图像进行复原。代码如下：
void CDipDoc: : OnRestoreBlur()
{
// TODO: Add your command handler code here
//判断当前是否载入图像
if( m_pImageObject == NULL ) return; D
//定义一个图像复原对象
CRestore Restore( m_pImageObject );
//调用卷积模糊函数
Restore.ConvBlur();
//更新视图
UpdateAllViews(NULL);
}
OnRestoreInvfilt()函数主要代码如下， 其他代码与OnRestoreBlur()相同。
void CDipDoc: : OnRestoreInvfilt()
{
//调用逆滤波复原函数
Restore.InvFilter();
}
功能： 实现图像复原
函数列表：
ConvBlur() 对图像进行卷积模糊处理， 生成一幅待复原的退化图像；
InvFilter() 逆滤波复原法， 对由卷积模糊所造成的退化图像进行复原；
ConvNoiseBlur() 对图像进行卷积噪声模糊处理， 生成一幅待复原的退化图像；
Wiener() 维纳滤波复原法， 对由卷积噪声所造成的退化图像进行复原。
ConvBlur()函数关键代码
{
//定义卷积模板， 大小为7×7
Mm m_matFilter;
m_matFilter=zeros(nHeight, nWidth);
int i, j;
for(i=1; i<7; i++)
{
for(j=1; j<7; j++)
{
m_matFilter.r(i, j)=0.02;
}
}
// m_matBits是一个二维矩阵， 存储的是图像数据
//对原图像进行快速傅立叶变换， 将图像从空域变换到频域
m_matBits=fft2(m_matBits);
//对模板进行快速傅立叶变换
m_matFilter=fft2(m_matFilter);
//频率相乘
for(i=1; i<=nHeight; i++)
{
for(j=1; j<=nWidth; j++)
{
m_matBits.r(i, j)=m_matBits.r(i, j)*m_matFilter.r(i, j)-m_matBits.i(i, j)*m_matFilter.i(i, j);
m_matBits.i(i, j)=m_matBits.r(i, j)*m_matFilter.i(i, j)+m_matBits.i(i, j)*m_matFilter.r(i, j);
}
}
//对图像做逆傅立叶变换， 将图像从频域变换到空域
m_matBits=ifft2(m_matBits);
}
InvFilter()函数关键代码
{
//对原图像进行快速傅立叶变换， 将图像从空域变换到频域
m_matBits=fft2(m_matBits);
//对模板进行快速傅立叶变换
m_matFilter=fft2(m_matFilter);
//频率相乘
double a, b, c, d;
for(i=1; i<=nHeight; i++)
{
for(j=1; j<=nWidth; j++)
{
a=m_matBits.r(i, j);
b=m_matBits.i(i, j);
c=m_matFilter.r(i, j);
d=m_matFilter.i(i, j);
if(c*c+d*d>0.003)
{
m_matBits.r(i, j)=(a*c+b*d)/(c*c+d*d);
m_matBits.i(i, j)=(b*c-a*d)/(c*c+d*d);
}
}
}
//对图像进行快速逆傅立叶变换， 将图像从频域变换到空域
m_matBits=ifft2(m_matBits);
}
图11-8 逆滤波法复原图像
(a) 原始图像； (b) 退化图像； (c)复原图像
