图像处理-Ch5-图像复原与重建

Ch5 图像复原

图像退化与复原(Image Degradation and Restoration)

图像复原目的：以某种预定义的方式改善给定图像。

Q: 图像增强 v.s. 图像复原?

A: 图像增强是一个主观的过程：自行选定不同的工具，比如低通、高通滤波，使得图像主观上看上去比较美观(因个人审美不同而不同)。

而图像复原的大部分过程是客观的。原先是什么样子、恢复原来的样子。

$$\begin{gathered} g(x,y)=H[f(x,y)]+\eta (x,y) \\ g(x,y)=h(x,y)\ast f(x,y)+\eta (x,y) \\ G(x,y)=H(u,v)F(u,v)+N(u,v) \end{gathered}$$

噪声模型(Noise Models)

假设$H(u,v)=1$: $g(x,y)=f(x,y)+\eta (x,y)$

i.i.d.空间随机噪声(Generating Spatial Random Noise with a Specified Distribution)

如果$w$是在区间$(0,1)$上均匀分布的随机变量，那么可以通过求解方程$z=F_z^{-1}(w)$得到一个具有特定累积分布函数（CDF）$F_z$的随机变量$z$。
$$z=F_z^{-1}(w)$$

例：目标是生成具有瑞利（Rayleigh）累积分布函数（CDF）的随机数$z$。

瑞利分布的累积分布函数（CDF）为:
$$F_z(z)=\{\begin{array}{ll}1-e^{-(z-a{)}^2/b}& \text{if }z\ge a\\ 0& \text{if }z<a\end{array}$$
通过求解方程$1-e^{-(z-a{)}^2/b}=w$，可以得到$z$的值：
$$z=a+\sqrt{-b\ln (1-w)}$$

周期噪声(Periodic Noise)

在图像中，周期噪声是在图像获取时从电力或机电干扰中产生的。这是唯一一种空间依赖型噪声。

周期噪声信号(Periodic Noise Signal) ：
$$r(x,y)=A\sin (2\pi u_0(x+B_x)/M+2\pi v_0(y+B_y)/N)$$
对于傅里叶变换：（实际上无法推出这个，但是大概是这个形式）
$$R(u,v)=\frac{A}{j2}\left[\exp \left(\frac{j2\pi u_0{B}_x}{M}\right)\delta (u+u_0,v+v_0)-\exp \left(\frac{j2\pi v_0{B}_y}{N}\right)\delta (u-u_0,v-v_0)\right]$$

估计噪声参数(Estimating Noise Parameters)

假如没有被噪声污染，那么图中有三个灰度级，直方图中应该有三根线。但是有噪声，现在每根线向两边扩展。

图l对应椒盐噪声的直方图：可以看见一共有4条线（两边分别有两条），椒盐噪声产生黑白、黑色估计和原有的背景色重合，于是直方图中最左的线特别的长（大概的解释）。

估计噪声参数：

1. 周期噪声：对图像的傅里叶谱审视，可以直接发现：比如在除中心的亮点之外，还有其他的亮点。

   

2. 具有先验知识：事先知道这种图像中会存在噪声、会出现在哪个频段。

3. 在平坦的图像区域，噪声会暴露出来：比如说一块纯色的地板，应该灰度级是差不多的。

Q: 怎么做？

1. 最简单的方法是利用图像中的采样数据来估计噪声的均值和方差。
2. 通过直方图的形状来辨识最接近的概率密度函数(PDF)的匹配。
3. 假如形状类似高斯，那么高斯只需要均值与方差。

统计矩与中心矩：

统计矩:
$${\mu }_n=\sum _{i=0}^{L-1}(z_i-m{)}^{n}p(z_i)\text{where}m=\sum _{i=0}^{L-1}{z}_{i}p(z_i)$$
中心距：

n	对应中心距
0	${\mu }_0={\sum }_{i=0}^{L-1}(z_i-m{)}^{0}p(z_i)=1$
1	${\mu }_1={\sum }_{i=0}^{L-1}(z_i-m)p(z_i)={\sum }_{i=0}^{L-1}{z}_{i}p(z_i)-m{\sum }_{i=0}^{L-1}p(z_i)=0$
2	${\mu }_2={\sum }_{i=0}^{L-1}(z_i-m{)}^{2}p(z_i)=\text{Var}(z)$

在仅有噪声情况下图像复原-空域滤波

如果只有噪声, 图像退化模型可以表示为：
$$\begin{gathered} g(x,y)=f(x,y)+\eta (x,y) \\ G(u,v)=F(u,v)+N(u,v) \end{gathered}$$

均值滤波器(Mean Filters)

均值相当于：给出一组数，经过一系列运算，得到一个新的值。这个值，比这组数中最小值大，比最大值小，那么调和均值滤波器、反调和均值滤波器就算均值滤波器。

调和均值滤波器能够很好地去除盐噪声，但不能去除椒噪声。对于高斯噪声很好。

反调和均值滤波器非常适合消除椒盐噪声。当$Q$为正值时，该滤波器消除胡椒噪声。当$Q$为负值时，它消除盐噪声。它不能同时消除两种噪声。

• 当Q=0时：算术均值是反调和均值的一个特例。
  $$\hat{y}(x,y)=\frac{{\sum }_{(s,t)\in S_{xy}}g(s,t)}{{\sum }_{(s,t)\in S_{xy}}1}=\frac{1}{mn}\sum _{(s,t)\in S_{xy}}g(s,t)$$

• 当Q=-1时：是调和均值滤波器。
  $$\hat{f}(x,y)=\frac{{\sum }_{(s,t)\in S_{xy}1}}{{\sum }_{(s,t)\in S_{xy}}g(s,t)}=\frac{mn}{{\sum }_{(s,t)\in S_{xy}}g(s,t)}$$

统计排序滤波器(Order-statistics Filters)

椒盐噪声反复使用中值滤波器，总能完全去除；同时边缘信息也会丢失。

用频域滤波器消除周期噪声(Periodic Noise Reduction by Frequency Domain Filtering)

关于以下三种滤波器的奇妙比喻：来自王伟强老师。

入冬时存储的土豆，在开春时会发芽。

• 带阻滤波器：连皮肉带芽一起削掉。
• 带通滤波器：与带阻合为1体。
• 陷波滤波器：只把有芽的部分去除。

带阻滤波器(Bandreject Filters)

带通滤波器(Bandpass Filters)

$$H_{bp}=1-H_{br}(u,v)$$

带通滤波器执行与带阻滤波器相反的操作。得到噪声信号。

带通滤波器与带阻滤波器是一对儿；低通滤波器与高通滤波器是一对儿。

陷波滤波器(Notch filters)

成对出现。(a)是原图像，可以看见存在一些横纹(噪声)。然会得到(b)是(a)的傅里叶谱；通过©中的陷波滤波器的传递函数，可以去除掉噪声。(d)是图像复原后的图；(e)是提取出来的噪声。

公式中还定义了$D_1(u,v)$和$D_2(u,v)$的具体形式：
$$\begin{gathered} D_1(u,v)={\left[{\left(u-\frac{M}{2}-u_0\right)}^2+{\left(v-\frac{N}{2}-v_0\right)}^2\right]}^{\frac{1}{2}} \\ {D}_2(u,v)={\left[{\left(u-\frac{M}{2}+u_0\right)}^2+{\left(v-\frac{N}{2}+v_0\right)}^2\right]}^{\frac{1}{2}} \end{gathered}$$
这里$M$和$N$是图像的尺寸，$u_0$和$v_0$是滤波器的中心坐标。 这些高通滤波器在图像处理中常用于增强图像的高频细节，例如边缘和纹理。

不同的滤波器适用于不同的应用场景，第一种滤波器是理想高通滤波器，具有最尖锐的截止，但会导致振铃效应；Butterworth和Gaussian高通滤波器则提供了更平滑的过渡，减少了振铃效应。

最佳陷波滤波算法(Optimum Notch Filtering)

当存在多个干扰分量时，之前提到的方法并不总是可行的，因为在滤波过程中会去除过多的图像信息。这里讨论的方法是最优的，因为它使恢复估计值$\hat{f}(x,y)$的局部方差最小化。

首先，通过以下方式获得噪声的初始估计：
$$\begin{gathered} N(u,v)=F_N(u,v)G(u,v) \\ \eta (x,y)={\mathfrak{J}}^{-1}(F_N(u,v)G(u,v)) \end{gathered}$$
其中$F_N(u,v)$被构建为仅通过与干扰模式相关的分量。

令：
$$\hat{f}(x,y)=g(x,y)-w(x,y)\eta (x,y)$$
$\eta (x,y)$估计已知，我们将确定调制函数$w(x,y)$，以使$\hat{f}(x,y)$的局部方差最小化。

目标：局部很好的平滑性，相邻的像素之间比较相似。
$$\begin{gathered} \min {\sigma }^2(x,y)=\frac{1}{(2a+1)(2b+1)}\sum _{s=-a}^a\sum _{t=-b}^b[\hat{f}(x+s,y+t)-\bar{f}{]}^2 \\ \bar{f}=\frac{1}{(2a+1)(2b+1)}\sum _{s=-a}^a\sum _{t=-b}^b\hat{f}(x+s,y+t) \end{gathered}$$
对上式展开：
$$\begin{gathered} {\sigma }^2(x,y)=\frac{1}{(2a+1)(2b+1)}\sum _{s=-a}^a\sum _{t=-b}^b \\ {\left[[g(x+s,y+t)-w(x+s,y+t)\cdot \eta (x+s,y+t)]-[\overline{g(x,y)}-\overline{w(x,y)\eta (x,y)}]\right]}^2 \end{gathered}$$
这个式子中只有$w(x+s,y+t)$是变量，变量一共有$(2a+1)(2b+1)w(x+s,y+t)$个变量。可以看到每一个$(x,y)$位置都能建立这样的目标函数、有这么多的变量去求解，实在是难以计算。因此进行简化：将原来的M*N个位置的多目标优化问题，转化为了M*N个i.i.d,优化问题。

令$w(x+s,y+t)=w(x,y)$, (将$(2a+1)(2b+1)$个变元简化为$1$ )我们有单变元函数：
$$\begin{gathered} {\sigma }^2(x,y)=\frac{1}{(2a+1)(2b+1)}\sum _{s=-a}^a\sum _{t=-b}^b \\ {\left([g(x+s,y+t)-w(x,y)\cdot \eta (x+s,y+t)]-[\overline{g(x,y)}-\overline{w(x,y)\eta (x,y)}]\right)}^2 \end{gathered}$$
为了最小化${\sigma }^2$，我们求解：
$$\begin{array}{rl}\frac{\partial {\sigma }^2}{\partial w(x,y)}& =0\\ 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)}\end{array}$$
$\overline{{\eta }^2}(x,y)$是平方的均值，${\bar{\eta }}^2(x,y)$是均值的平方。

这样我们就能求解出每个位置的weight: $w(x,y)$，需要求M*N次。

估计退化函数(Estimating the Degradation Function)

该节分为以下两个部分进行讨论：

1. 假设没有噪声的情况：$\eta (x,y)=0,g(x,y)=h(x,y)\ast f(x,y)$
2. 假设存在噪声的情况，如何处理？

进行图像复原，第一步要做的，就是找到$H(u,v)$. 有以下三种方法：

运动引起的图像模糊(Image Blur due to Motion)

现在，假设一幅图像由于均匀线性运动而变得模糊。
$$g(x,y)={\int }_0^{T}f\left(x-x_0(t),y-y_0(t)\right)dt$$
$T$是快门时间。然后，它的傅里叶变换是:
$$\begin{array}{rl}G(u,v)& ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }g(x,y)\exp \left[-j2\pi (ux+vy)\right]dxdy\\ & ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\left[{\int }_0^{T}f\left(x-x_0(t),y-y_0(t)\right)dt\right]\exp \left[-j2\pi (ux+vy)\right]dxdy\\ & ={\int }_0^T\left[{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f\left(x-x_0(t),y-y_0(t)\right)\exp \left[-j2\pi (ux+vy)\right]dxdy\right]dt\\ & ={\int }_0^{T}F(u,v)\exp \left[-j2\pi (u{x}_0(t)+v{y}_0(t))\right]dt\\ & =F(u,v){\int }_0^T\exp \left[-j2\pi (u{x}_0(t)+v{y}_0(t))\right]\\ & \\ H(u,v)& ={\int }_0^T\exp \left[-j2\pi (u{x}_0(t)+v{y}_0(t))\right]dtdt\end{array}$$

如果$x_0(t)=at/T$，$y_0(t)=0$（垂直方向上没有发生运动，只是在水平方向上抖动了以下），那么：
$$H(u,v)={\int }_0^T\exp \left(-j2\pi uat/T\right)dt=\frac{T}{\pi ua}\sin (\pi ua)\exp \left(-j\pi ua\right)$$
如果$y_0(t)=bt/T$而不是0(水平、垂直方向上都发生了运动)，那么：
$$H(u,v)=\frac{T}{\pi (ua+vb)}\sin (\pi (ua+vb))\exp \left(-j\pi (ua+vb)\right)$$

直接逆滤波(Direct Inverse Filtering)

对于被退化函数$H$退化的图像，最简单的恢复方法是直接逆滤波:
$$\hat{F}(u,v)=\frac{G(u,v)}{H(u,v)}$$
可以看到，这个式子中假设$N(u,v)=0$，且我们需要已知$G(u,v),H(u,v)$.

那要是$N(u,v)\ne 0$，我们可以得到:
$$\hat{F}(u,v)=F(u,v)+\frac{N(u,v)}{H(u,v)}$$
如果退化函数$H(u,v)$在某一个频率点有零值或非常小的值(一般都是离原点比较远的地方存在)，那么$N(u,v)/H(u,v)$的比值会很大、可能会压制$F(u,v)$的估计，造成$F(u,v)$的估计存在较大误差。

一种解决方法是利用原点附近退化函数的值来进行逆滤波的计算。比如设置一个半径：$u^2+v^2\le R^2$，足够小的时候来计算$\hat{F}(u,v)$; $u^2+v^2>R^2$的地方，做低通滤波。

可以推断出来噪声的频率大于70.（85把噪声包含进去，此时就是$N(u,v)/H(u,v)$的比值会很大、可能会压制$F(u,v)$的估计，造成$F(u,v)$的估计存在较大误差。得到了依托）

维纳滤波器(Wiener Filtering)

理论上最优的一种图像复原方法。

维纳滤波器寻求使统计误差函数最小化的估计值$\hat{f}$
$$e^2=E\{(f-\hat{f}{)}^2\}$$
$E$是期望算子，$f$是未退化图像。

此表达式在频域中的解为:
$$\hat{F}(u,v)=\left[\frac{1}{H(u,v)}\frac{|H(u,v){|}^2}{|H(u,v){|}^2+S_n(u,v)/S_f(u,v)}\right]G(u,v)$$
其中 :

| $H(u,v)$ | $|H(u,v){|}^2=H^{\ast }(u,v)H(u,v)$ | $H^{\ast }(u,v)$ | $S_n(u,v)=|N(u,v){|}^2$ | $S_f(u,v)=|F(u,v){|}^2$ |
| --------- | -------------------------------- | ------------------- | ------------------------- | ------------------------- |
| 退化函数 | — | 是$H(u,v)$的复共轭 | 是噪声的功率谱 | 未退化图像的功率谱 |

实际上很难在现实中求解：主要是因为$S_n(u,v)/S_f(u,v)$、很难求噪声的功率谱以及未退化图像的功率谱。通常使用近似来应用：$S_n(u,v)/S_f(u,v)=K$
$$\hat{F}(u,v)=\left[\frac{1}{H(u,v)}\frac{|H(u,v){|}^2}{|H(u,v){|}^2+K}\right]G(u,v)$$

约束最小二乘图像复原算法(Constrained Least Squares Filtering)

了解退化函数$H$的问题 - 在本章讨论的所有方法中，都存在需要了解退化函数$H$的问题。

维纳算法不好的是：得事先知道噪声功率谱和未退化图像谱，但这个很难知道。

接下来介绍约束最小二乘滤波方法，这个方法只需要知道噪声的均值和方差。

如果我们用矩阵来模拟退化过程，有：
$$\begin{gathered} g(x,y)=h(x,y)\ast f(x,y)+\eta (x,y) \\ g=Hf+\eta \end{gathered}$$
将其转化成矩阵形式，相当于是做$vec()$操作拉伸。比如原先的$f(x,y):M\times N$，现在的$f$是$MN\times 1$.

该方法的核心是$H$对噪声的敏感性问题。一种缓解该问题的方法是基于平滑度的度量进行恢复。

拉普拉斯（图像的二阶导数:${\nabla }^{2}f=\frac{{\partial }^{2}f}{\partial x^2}+\frac{{\partial }^{2}f}{\partial y^2}$）似乎是一个很好的选择。

一阶导反映变化快慢、二阶导反映了平滑程度。

需要最小化的代价函数$C$为 :
$$\begin{gathered} \min C=\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}[{\nabla }^{2}f(x,y){]}^2 \\ s.t.\Vert \mathbf{g}-\mathbf{H}\hat{\mathbf{f}}\Vert =\Vert \eta \Vert \end{gathered}$$

已知：$g,H,\eta$。该问题在频域的解为 :
$$\hat{F}(u,v)=\left[\frac{H^{\ast }(u,v)}{|H(u,v){|}^2+\gamma |P(u,v){|}^2}\right]G(u,v)$$
其中$P(x,y)$是拉普拉斯算子在空域中的一个形式(傅里叶变换)。
$$P(x,y)=\left(\begin{array}{ccc}0& -1& 0\\ -1& 4& -1\\ 0& -1& 0\end{array}\right)$$

令$r=g-H\hat{f}$，且$\phi (\gamma )=r^{T}r=\Vert r{\Vert }^2$。

可以证明$\phi (\gamma )$是$\gamma$的单调递增函数。因此，我们可以调整$\gamma$，使得 :
$$\Vert r{\Vert }^2=\Vert \eta {\Vert }^2\pm \alpha$$

当$\alpha$足够小，那么$\Vert r{\Vert }^2\sim \Vert \eta {\Vert }^2$. 回看频域的解：$\hat{F}(u,v)=\left[\frac{H^{\ast }(u,v)}{|H(u,v){|}^2+\gamma |P(u,v){|}^2}\right]G(u,v)$.

已知的条件有：$H(u,v),P(u,v),G(u,v)$,那么当$\gamma$确定时：$\gamma \Rightarrow \hat{F}\Rightarrow \hat{f}\Rightarrow g-H\hat{f}=r$

然后根据$\Vert r{\Vert }^2=\Vert \eta {\Vert }^2\pm \alpha$确定$\Vert r{\Vert }^2$是否落入区间，不在区间的话，可以不断地调整$\gamma$,使$\Vert r{\Vert }^2<\Vert \eta {\Vert }^2\pm \alpha$落入邻域。此时，所估计的$\hat{f}$是问题的解。

Q: 前面假设噪声已知，那实际中如何计算$\Vert \eta {\Vert }^2$？
$$\begin{gathered} \Vert \eta {\Vert }^2=MN[{\sigma }_{\eta }^2+m_{\eta }^2] \\ \{\begin{array}{l}{\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)\end{array} \end{gathered}$$
${\eta }^2,m_{\eta }^2$是噪声的方差和期望。

$$\begin{array}{rl}{\sigma }_{\eta }^2& =\frac{1}{MN}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}[\eta (x,y)-m_{\eta }{]}^2\\ & =\frac{1}{MN}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}[{\eta }^2(x,y)-2m_{\eta }\eta (x,y)+m_{\eta }^2]\\ & =\frac{1}{MN}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}{\eta }^2(x,y)-2m_{\eta }\frac{1}{MN}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}\eta (x,y)+m_{\eta }^2\\ & =\frac{1}{MN}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}{\eta }^2(x,y)-2m_{\eta }^2+m_{\eta }^2\\ & =\frac{1}{MN}\sum _{x=0}^{M-1}\sum _{y=0}^{N-1}{\eta }^2(x,y)-m_{\eta }^2\\ \sum _{x=0}^{M-1}\sum _{y=0}^{N-1}{\eta }^2(x,y)& =MN[{\sigma }_{\eta }^2+m_{\eta }^2]\\ \Vert \eta {\Vert }^2& =MN[{\sigma }_{\eta }^2+m_{\eta }^2]\end{array}$$

下图都展示了经过约束最小二乘滤波处理后的图像。可以看到效果比维纳滤波好一点。

)+m_{\eta}^2\ &=\frac{1}{MN}\sum_{x = 0}^{M-1}\sum_{y = 0}^{N-1}\eta2(x,y)-2m_{\eta}^2+m_{\eta}2\ &=\frac{1}{MN}\sum_{x = 0}^{M-1}\sum_{y = 0}^{N-1}\eta2(x,y)-m_{\eta}^2 \

\sum_{x = 0}^{M-1}\sum_{y = 0}^{N-1}\eta2(x,y)&=MN[\sigma_{\eta}^2 + m_{\eta}^2]\
|\mathbf{\eta}|^2 &= MN[\sigma_{\eta}^2+m_{\eta}2]
\end{align*}
$$

下图都展示了经过约束最小二乘滤波处理后的图像。可以看到效果比维纳滤波好一点。

[外链图片转存中…(img-UnOHPG8p-1735022584536)]