数字图像处理 第5章-图像复原与重建

数字图像处理 第5章-图像复原与重建

正像图像增强那样，图像复原技术的主要目的是以预先确定的目标来改善图像。尽管两者有相覆盖的领域，但图像增强主要是一个主观过程，而图像复原则大部分是一个客观过程。图像复原试图利用退化现象的某种先验知识来复原被退化的图像。因而，复原技术是面向退化模型的，并且采用相反的过程进行处理，以便恢复出原图像。

相比之下，图像增强技术基本上是一个探索性过程,即根据人类视觉系统的生理特点来设计一种改善图像的方法。

例如。对比度拉伸被认为是一种增强技术，因为它主要基于给观看者提供其喜欢接受的图像，而通过去模糊函数去除图像模糊则被认为是一种图像复原技术。

5.1 图像退化/复原过程的模型

退化过程被建模为一个退化函数和一个加性噪声项，对一幅输入图像$f(x,y)$进行处理，产生一幅退化后的图像 $g(x,y)$。给定$g(x,y)$和关于退化函数$H$的一些知识以及关于加性噪声项$\eta (x,y)$的一些知识后，图像复原的目的就是获得原始图像的一个估计$\hat{f}(x,y)$。通常，我们希望这一估计尽可能地接近原始输入图像，并且$H$和$\eta$的信息知道得越多，所得到的$\hat{f}(x,y)$就会越接近$f(x,y)$。

如果$H$是一个线性的、位置不变的过程，那么空间域中的退化图像可由下式(5.1-1)给出：
$$g(x,y)=h(x,y)★f(x,y)+\eta (x,y)$$
其中，$h(x,y)$是退化函数的空间表示，符号“★”表示空间卷积。

空间域中的卷积等同于频率域中的乘积，因此可以把式(5.1-1)中的模型写成等价的频率域表示：
$$G(u,v)=H(u,v)F(u,v)+N(u,v)$$
其中的大写字母项是式(5.1-1)中相应项的傅里叶变换。
下图为图像退化/复原过程的模型

5.2 噪声模型

数字图像中，噪声主要来源于图像的获取和/或传输过程。成像传感器的性能受各种因素的影响如图像获取过程中的环境条件和传感元器件自身的质量。例如，在使用CCD摄像机获取图像时，光照水平和传感器温度是影响结果图像中噪声数量的主要因素。图像在传输中被污染主要是由于传输信道中的干扰。例如，使用无线网络传输的图像可能会因为光照或其他大气因素而污染。

5.2.1 噪声的空间和频率特性

与我们的讨论相关的是定义噪声空间特性的参数，以及噪声是否与图像相关。频率特性是指傅里叶域中噪声的频率内容(即相对于电磁波谱的频率)。例如，当噪声的傅里叶谱是常量时，噪声通常称为白噪声。这个术语是从白光的物理特性派生出来的，它以相等的比例包含可见光谱中的几乎所有频率。

5.2.2 一些重要的噪声概率密度函数

基于前一节的假设，我们关心的空间噪声描述子就是图5.1中模型的噪声分量中灰度值的统计特性。可以认为它们是由概率密度函数(PDF)表征的随机变量。下面是在图像处理应用中最常见的 PDF。

高斯噪声

在空间域和频率域中，由于高斯噪声在数学上的易处理性，故实践中常用这种噪声(也称为正态噪声)模型。事实上，这种易处理性非常方便，以至于高斯模型常常应用于在一定程度上导致最好结果的场合。
高斯随机变量$z$的 PDF 由下式给出：
$$p(z)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-(z-\bar{z}{)}^2/2{\sigma }^2}$$
其中，z表示灰度值，之表示$\bar{z}$的均值，$\sigma$表示$z$的标准差。标准差的平方${\sigma }^2$称为$z$的方差。高斯函数的曲线如图5.2(a)所示。当z服从式(5.2-1)的分布时，其值有大约70%落在范围[($\bar{z}$-σ),($\bar{z}$+σ)]内，有大约 95%落在范围[($\bar{z}$-2σ),($\bar{z}$+2σ)]内。

瑞利噪声

其PDF如下：
$$p(z)==\{\begin{array}{ll}\frac{2}{b}(z-a)e^{-(z-a{)}^2/b}& z\ge a\\ 0& z<a\end{array}$$

爱尔兰（伽马）噪声

爱尔兰噪声的 PDF由下式给出
$$p(z)==\{\begin{array}{ll}\frac{a^b{z}^{b-1}}{(b-1)!e^{-az}}& z\ge a\\ 0& z<a\end{array}$$

指数噪声

指数噪声的 PDF由下式给出
$$p(z)==\{\begin{array}{ll}a{e}^{-az}& t\ge 0\\ 0& z<0\end{array}$$

均匀噪声

均匀噪声的 PDF由下式给出
$$p(z)==\{\begin{array}{ll}\frac{1}{b-a}& a\le z\le b\\ 0& 其他\end{array}$$

脉冲（椒盐）噪声

$$p(z)==\{\begin{array}{ll}{P}_a& z=a\\ P_b& z=b\\ 1-P_a-P_b& 其他\end{array}$$
如果$b>a$，则灰度级$b$在图像中将显示为一个亮点；反之，灰度级$a$在图像中将显示为一个暗点。若$P_a$或$P_b$为零，则脉冲噪声称为单极脉冲。如果$P_a$和$P_b$两者均不可能为零，尤其是它们近似相等时，则脉冲噪声值将类似于在图像上随机分布的胡椒和盐粉微粒。由于这个原因，双极脉冲噪声也称为椒盐噪声。这种类型的噪声也可以使用散粒噪声和尖峰噪声来称呼。

5.2.3 周期噪声

一幅图像中的周期噪声是在图像获取期间由电力或机电干扰产生的。
周期噪声可通过频率域滤波来显著地减少。

5.2.4 噪声参数的估计

典型地，周期噪声的参数是通过检测图像的傅里叶谱来估计的。正如前几节提及的那样，周期声趋向于产生频率尖峰，甚至通过视觉分析也经常可以检测到这些尖峰。另一种方法是试图直接由图像推断出噪声分量的周期性，但这仅在非常简单的情况下才是可能的。在噪声尖峰格外显著或可以使用关于干扰的频率分量一般位置的某些知识时，自动分析是可能的。

5.3 只存在噪声的复原——空间滤波

当一幅图像中唯一存在的退化是噪声时，式(5.1-1)和式(5.1-2)变成：
$$g(x,y)=f(x,y)+\eta (x,y)$$
和
$$G(u,v)=F(u,v)+N(u,v)$$
噪声项是未知的,故从$g(x,y)$或$G(u,v)$中减去它们不是一个现实的选择。在周期噪声的情况下，由$G(u,v)$的谱来估计$N(u,v)$通常是可能的。在这种情况下，从 $G(u,v)$中减去 $N(u,v)$可得到原图像的一个估计。然而，这种类型的知识通常只是例外而不是规律。

5.3.1 均值滤波器

算术均值滤波器

这是最简单的均值滤波器。令$S_{xy}$表示中心在点(x,y)处、大小为 $m\ast n$ 的矩形子图像窗口(邻域)的一组坐标。算术均值滤波器在$S_{xy}$定义的区域中计算被污染图像$g(x,y)$的平均值。在点(x,y)处复原图像$\hat{f}$的值，就是简单地使用$S_{xy}$定义的区域中的像素计算出的算术均值，即
$$\hat{f}(x,y)=\frac{1}{mn}\sum _{(s,t)\in S_{xy}}g(s,t)$$
这个操作可以使用大小为$m\ast n$的一个空间滤波器来实现，其所有的系数均为其值的 $1/mn$。均值滤波平滑一幅图像中的局部变化，虽然模糊了结果，但降低了噪声。

几何均值滤波器

使用几何均值滤波器复原的一幅图像由如下表达式给出：
$$\hat{f}(x,y)=[\prod _{(s,t)\in S_{xy})}g(s,t){]}^{\frac{1}{mn}}$$
其中，每个复原的像素由子图像窗口中像素的乘积的$1/mn$ 次幂给出。

谐波均值滤波器

谐波均值滤波操作由如下表达式给出：
$$\hat{f}(x,y)=\frac{mn}{{\sum }_{(s,t)\in S_{xy}}\frac{1}{g(s,t)}}$$
谐波均值滤波器对于盐粒噪声效果较好，但不适用于胡椒噪声。它善于处理像高斯噪声那样的其他噪声。

逆谐波均值滤波器

逆谐波均值滤波器基于如下表达式产生一幅复原的图像：
$$\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}$$
其中 $Q$称为滤波器的阶数。这种滤波器适合减少或在实际中消除椒盐噪声的影响。当$Q$值为正时，该滤波器消除胡椒噪声;当$Q$值为负时，该滤波器消除盐粒噪声。但它不能同时消除这两种噪声。注意，当$Q$=0时，逆谐波均值滤波器简化为算术均值滤波器；而当$Q$=-1时，则为谐波均值滤波器。

5.3.2 统计排序滤波器

中值滤波器

最著名的统计排序滤波器是中值滤波器，如其名称所暗示的那样，它使用一个像素邻域中的灰度级的中值来替代该像素的值，即
f ^ ( x , y ) = m e d i a n ( s , t ) ∈ S x y { g ( s , t ) } {\hat f (x,y)}=median_{(s,t)∈S_{xy}}\{g(s,t)\} f^​(x,y)=median(s,t)∈Sxy​​{g(s,t)}
在(x,y)处的像素值是计算的中值。中值滤波器的应用非常普遍，因为对于某些类型的随机噪声，它们可提供良好的去噪能力，且比相同尺寸的线性平滑滤波器引起的模糊更少。在存在单极或双极脉冲噪声的情况下，中值滤波器尤其有效。

最大值和最小值滤波器

尽管中值滤波器是目前为止图像处理中最常用的一种统计排序滤波器，但它绝不是唯一的一种中值相当于顺序排列的数组中间的那个数，但是您会从基本的统计学回想到排列本身还有很多其他的可能性。例如，可以使用序列中的最后一个数值，称为最大值滤波器，由下式给出
f ^ ( x , y ) = m a x ( s , t ) ∈ S x y { g ( s , t ) } {\hat f (x,y)}=max_{(s,t)∈S_{xy}}\{g(s,t)\} f^​(x,y)=max(s,t)∈Sxy​​{g(s,t)}

这种滤波器对于发现图像中的最亮点非常有用。同样，因为胡椒噪声的值非常低，作为子图像区域S中这种最大值选择过程的结果，可以用这种滤波器降低它。选择起始值的滤波器称为最小值滤波器，它由下式给出：
f ^ ( x , y ) = m i n ( s , t ) ∈ S x y { g ( s , t ) } {\hat f (x,y)}=min_{(s,t)∈S_{xy}}\{g(s,t)\} f^​(x,y)=min(s,t)∈Sxy​​{g(s,t)}
这种滤波器对于发现图像中的最暗点非常有用。同样,作为最小值操作的结果，它可以降低盐粒噪声。

中点滤波器

中点滤波器简单地计算滤波器包围区域中最大值和最小值之间的中点，即
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 ) } ] {\hat f (x,y)}={1\over 2}[max_{(s,t)∈S_{xy}}\{g(s,t)\}+min_{(s,t)∈S_{xy}}\{g(s,t)\}] f^​(x,y)=21​[max(s,t)∈Sxy​​{g(s,t)}+min(s,t)∈Sxy​​{g(s,t)}]
注意，这种滤波器结合了统计排序和求平均。它对于随机分布噪声工作得最好，如高斯噪声或均匀噪声。

修正的阿尔法均值滤波器

假设在邻域$S$内去掉 $g(s,t)$最低灰度值的 $d/2$和最高灰度值的 $d/2$。令$g_r(s,t)$代表剩下的 $mn-d$ 个像素。由这些剩余像素的平均值形成的滤波器称为修正的阿尔法均值滤波器
$$\hat{f}(x,y)=\frac{1}{mn-d}\sum _{(s,t)\in S_{xy}}{g}_r(s,t)$$
其中，$d$的取值范围可为0到$mn-1$。当$d=0$时，修正的阿尔法均值滤波器退化为前一节讨论的算术均值滤波器。如果选择$d=mn-1$，则修正的阿尔法均值滤波器将退化为中值滤波器。当$d$取其他值时，修正的阿尔法均值滤波器在包括多种噪声的情况下很有用，例如高斯噪声和椒盐噪声混合的情况下。

5.3.3 自适应滤波器

滤波器作用于局部区域$S_{xy}$。滤波器在该区域中心任意一点(x,y)上的响应基于以下4个量：
(a)$g(x,y)$，带噪图像在点(x,y)上的值
(b)${\sigma }_{\eta }^2$，污染 f(x,y)以形成 g(x,y)的噪声的方差
©$m_L$，S中像素的局部均值
(d)${\sigma }_L^2$，S中像素的局部方差。
我们希望滤波器的性能如下：
1、如果${\sigma }_{\eta }^2$为零，则滤波器应该简单地返回 g(x,y)的值。这无关紧要，在零噪声情况下 g(x,y)等于 f(x,y)。
2、如果局部方差与${\sigma }_{\eta }^2$是高度相关的，则滤波器返回g(x,y)的一个近似值。典型地，高局部方差与边缘相关，并且应该保护这些边缘。
3、如果两个方差相等，我们则希望滤波器返回$S_{xy}$中像素的算术均值。这种情况发生在局部区域与整个图像有相同特性的条件下，并且局部噪声将通过简单地求平均来降低。
基于这些假设得到的$\hat{f}(x,y)$的自适应表达式可以写成
$$\hat{f}(x,y)=g(x,y)-\frac{{\sigma }_{\eta }^2}{{\sigma }_L^2}[g(x,y)-m_L]$$

5.4 用频率域滤波消除周期噪声

用频率域技术可以有效地分析并滤除周期噪声。其基本概念是在傅里叶变换中，周期噪声在对应于周期干扰的频率处,以集中的能量脉冲形式出现。

5.4.1 带阻滤波器

在第四章中介绍的理想、布特沃斯和高斯带阻滤波器的传递函数总结中。图5.15 显示了这些滤波器的透视图，并且，下面的例子说明了使用一个带阻滤波器降低周期噪声的效果。

5.4.2 带通滤波器

带通滤波器执行与带阻滤波器相反的操作

5.4.3 陷波滤波器

陷波滤波器阻止(或通过)事先定义的中心频率的邻域内的频率。

5.5 线性、位置不变的退化

图像复原前，图 5.1中的输入/输出关系可以表示为
$$g(x,y)=H[f(x,y)]+\eta (x,y)$$
现在，让我们假设$n(x,y)=0$，则 $g(x,y)=H[f(x,y)]$。如果
$$H[a{f}_1(x,y)+b{f}_2(x,y)]=aH[f_1(x,y)]+bH[f_2(x,y)]$$
则系统$H$是一个线性系统，其中a和b是标量，$f_1(x,y)$和$f_2(x,y)$是两幅输入图像。
若a=b=1，则变为
$$H[f_1(x,y)+f_2(x,y)]=H[f_1(x,y)]+H[f_2(x,y)]$$
这就是所谓的加性。这一特性简单地表明，如果 H为线性算子，那么两个输入之和的响应等于两个响应之和。
如果£(x,y)=0，则式(5.5-2)变为
$$H[a{f}_1(x,y)]=aH[f_1(x,y)]$$
我们称之为均匀性。它表明任何与常数相乘的输人的响应等于该输人响应乘以相同的常数，即一个线性算子具有加性和均匀性。
一个算子对于任意$f(x,y)$，$\alpha$和$\alpha$，如果
$$H[f(x-\alpha ,y-\beta )]=g(x-\alpha ,y-\beta )$$
则输人/输出关系具有 $g(x,y)=H[f(x,y)]$的算子称为位置(或空间)不变系统。该定义说明图像中任意点处的响应只取决于该点处的输入值，而与该点的位置无关。

$h(x,\alpha ,y,\beta )=H[\delta (x-,y-\beta )]$为系统$H$的冲激响应。

如果H是位置不变的，线性系统H完全可由其冲激响应来表征：
$$g(x,y)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(\alpha ,\beta )h(x-\alpha ,y-\beta )d\alpha d\beta$$

5.6 估计退化函数

5.6.1 图像观察估计

假设给我们一幅退化图像，而没有关于退化函数 H的任何知识。说白了就是基于图像本身来收集信息。

5.6.2 试验估计

是使用与获取退化图像的设备相似的装置，得到一个准确的退化估计。

5.6.3 建模估计

是使用退化建模来解决图像复原的问题

5.7 逆滤波

最简单的复原方法是直接做逆滤波，在这里，我们用退化函数除退化图像的傅里叶变换 $G(u,v)$来计算原始图像傅里叶变换的估计$\hat{F}(x,y)$，即
$$\hat{F}(u,v)=\frac{G(u,v)}{H(u,v)}$$
该除法是阵列操作。

根据之前的讨论，进一步我们还能得到：
$$\hat{F}(u,v)=F(u,v)+\frac{N(u,v)}{H(u,v)}$$

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

这一方法综合了退化函数和噪声统计特征进行复原处理的方法。该方法建立在图像和噪声都是随机变量的基础上，目标是找到未污染图像$f$的一个估计$\hat{f}$使它们之间的均方误差最小。

5.9 约束最小二乘方滤波

维纳滤波却存在一些其他的困难：未退化图像和噪声的功率谱必须是已知的。在前面的章节给出的近似方法能得到很好的结果。然而，功率谱比的常数估计并不总是一个合适的解。。正如之前讨论的那样，这些参数通常可从一幅给定的退化图像计算出来，因此这是一个很重要的优点。另一个不同是，维纳滤波建立在最小化统计准则的基础之上，因此在平均意义上它是最优的。

5.10 几何均值滤波

可对5.8节中讨论的维纳滤波器稍加推广。这种推广就是所谓的几何均值滤波器的形式：

$$\hat{F}(u,v)=[\frac{H^{\ast }(u,v)}{|H(u,v){|}^2}{]}^{\alpha }[\frac{H^{\ast }(u,v)}{|H(u,v){|}^2+\beta [\frac{S_{\eta }(u,v)}{S_f(u,v)}]}{]}^{1-\alpha }G(u,v)$$
其中，$\alpha$和$\beta$是正的实常数。几何均值滤波器由两个括号内的幂次分别为$\alpha$和$1-\alpha$的表达式组成

5.11 由投影重建图像

在本章的前几节中，我们讨论了各种复原退化图像的技术。在这一节，我们研究从一系列投影重建一幅图像的问题，集中讨论X射线计算机断层(CT)。这是最早且应用最广泛的CT类型，也是当前数字图像处理在医学中的主要应用之一。

5.11.1 计算机断层(CT)原理

简单来说，就是使用射线从许多不同的方向穿过物体而得到该物体内部结构。
即使用X射线来穿透物体，由于物体会吸收X射线的能量，而在另一端我们可以使用X射线吸收检测器来收集相关的信息，然后根据相关信息来成像。

5.11.2 傅里叶切片定理

式子如下：
$$G(\omega ,\theta )=[F(u,v){]}_{u=\omega cos\theta ;v=\omega sin\theta }==F(\omega cos\theta ,\omega sin\theta )$$
它说明一个投影的傅里叶变换是得到投影区域的二维傅里叶变换的一个切片。

5.11.3 使用平行射线束滤波反投影的重建

使用平行射线束滤波反投影的重建是为了解决直接使用反投影时会生成不可接受的模糊结果而出现的
其步骤如下：

1、计算每一个投影的一维傅里叶变换。
2、如上面说明的那样，用滤波函数$|\omega |$乘以每一个傅里叶变换，就是乘以一个合适的窗(如汉明窗)。
3、得到每一个滤波后的变换的一维反傅里叶变换。
4、对步骤3得到的所有一维反变换积分(求和)。

5.11.4 使用扇形射线束滤波反投影的重建

其公式为：

$$f(r,\phi )={\int }_0^{2\pi }\frac{1}{R^2}[{\int }_{-{\alpha }_m}^{{\alpha }_m}q(\alpha ,\beta )h({\alpha }^{\prime }-\alpha )]d\phi$$
其中
$$h(\alpha )=\frac{1}{2}(\frac{\alpha }{sin\alpha }{)}^{2}s(\alpha )$$
和
$$q(\alpha ,\beta )=p(\alpha ,\beta )Dcos\alpha$$