第四章 图像恢复和重建
图像恢复技术的主要目的是提前确定目标来改善图像，其中大部分属于客观处理，面向退化模型，并采用相反的过程来恢复原始图像。
图像增强技术基本上是一个探索性的过程，即根据人类视觉系统的生理和心理特征来设计改善图像的方法。
4.1 图像退化/恢复处理的模型
本章将退化过程建模为退化函数和加性噪声项f(x, y)退化后的图像处理g(x, y)：g(x, y) = H[f(x, y)] η(x, y)。给定g(x, y)一些关于退化函数H和加性噪声的知识η(x, y)经过一些知识，图像恢复的目的是获得源图像的估计 f ^ ( x , y ) \hat f (x,y) f^(x,y)。通常，我们希望这一估计尽可能地接近原始输入图像，并且H和η你知道的信息越多，你得到的信息就越多。 f ^ ( x , y ) \hat f (x,y) f^(x,y)越接近f(x, y)。
若H是一个线性、空间不变的过程，则空间域中的退化图像可由式给出： g ( x , y ) = h ( x , y ) ★ f ( x , y ) η ( x , y ) g(x, y)=h(x, y)★f(x, y) η(x, y) g(x,y)=h(x,y)★f(x,y) η(x,y)。式中，h(x, y)是该退化函数的空间表示。
空间域的卷积和频率域的乘法组成了一个傅里叶变换对，所以可以用等价的频率域来写出前面的模型：G(u, v)=H(u, v)F(u, v)+N(u, v)。式中，用大写字母表示的是空间域中相应项的傅里叶变换。退化函数F(u, v)有时称为光传递函数(OTF)（源于光学系统的傅里叶分析）。在空间域中，h(x, y)称为点扩散函数(PSF)（来源：对于任何类型的输入，让h(x, y)作用域一个点光源来得到退化的特征）。OTF和PSF是一个傅里叶变换对，工具箱函数otf2psf和psf2otf用于在OTF和PSF之间相互转换。
由于线性、空间不变的退化函数H可以由卷积来建模，所以退化处理有时也被称为PSF与图像卷积。类似的，复原处理有时也被称为去卷积。
4.2 噪声模型
模拟噪声特性和影响的能力是图像复原的核心。
4.2.1 使用函数imnoise对图像添加噪声
函数imnoise可以用噪声污染一幅图像，语法：g=imnoise(f, type, parameters)，其中f是输入图像。该函数在给图像添加噪声之前将图像转换为区间[0, 1]上的double类。语法形式如下表：
| 语法 | 注释 |
|---|---|
| g=imnoise(f, ‘gaussian’, m, var) | 将均值为吗、方差为var的高斯噪声添加到图像f上。默认值是均值为0、方差0.01的噪声。 |
| g=imnoise(f, ‘localvar’, V) | 将均值为0、局部方差为V的高斯噪声添加到图像f上，其中V是大小与f相同的数组，它在每个点包含了期望的方差值。 |
| g=imnoise(f, ‘localvar’, image_intensity, var) | 将均值为0的高斯噪声添加到图像f上，其中噪声的局部方差var是图像f的灰度值的函数。自变量image_intensity和var是大小相同的向量，plot(image_intensity, var)绘制出噪声方差和图像灰度间的函数关系。向量image_intensity必须是[0, 1]上的归一化灰度值。 |
| g=imnoise(f, ‘salt%pepper’， d) | 用椒盐噪声污染图像f，其中d是噪声密度（包含噪声值的图像区域的百分比），大约有d*numel(f)个像素受到影响，默认噪声密度0.05。 |
| g=imnoise(f, ‘speckle’, var) | 使用方程g=f+n.*f将乘性噪声添加到图像f上，其中n是均值为0、方差为var的均匀分布的随机噪声，var默认值0.04。 |
| g=imnoise(f, ‘poisson’) | 由数据生成泊松噪声，而不将人为噪声添加到数据中。为服从泊松统计，uint8、uint16类图像的灰度必须对应于光子数（或其他量子信息）。当每个像素中的光子数大于65535（但小于 1 0 12 10^{12} 1012)时需要使用双精度图像，灰度值在0~1之间变化，并且对应于光子的数量除以 1 0 12 10^{12} 1012。 |
4.2.2 使用规定分布生成空间随机噪声
空间噪声是随机数，由概率密度函数(PDF)或等效的累积分布函数(CDF)来表征。
许多随机数生成器以区间(0, 1)上的一个均匀CDF来表达随机数的生成问题为基础。有些情形下所选的基本随机数生成器是均值为0、方差为1的高斯随机数生成器。
本节的基本方法是概率论的一个著名结果：若w是在区间(0,1)上均匀分布的一个随机变量，则通过求解方程 z = F − 1 ( w ) z=F^{-1}(w) z=F−1(w)可得到一个具有规定CDF和F的随机变量z，该结果也可以等效表述为求方程 F ( z ) = w F(z)=w F(z)=w的解。
使用均匀随机数生成具有规定分布的随机数
假设有一个区间(0, 1)上的均匀随机数w的生成器，并假设要用该生成器生成具有瑞利CDF的随机数z，有如下形式： F ( z ) = { 1 − e − ( z − a ) 2 b , z ≥ a 0 , z < a ( b > 0 ) F(z)=\begin{cases} 1-e^{-\frac{(z-a)^2}{b}},z\geq a \\0,z<a \end{cases}(b>0) F(z)={ 1−e−b(z−a)2,z≥a0,z<a(b>0)。求出z需要解方程 1 − e − ( z − a ) 2 b = w 1-e^{-\frac{(z-a)^2}{b}}=w 1−e−b(z−a)2=w或 z = a + − b l n ( 1 − w ) z=a+\sqrt{-bln(1-w)} z=a+−bln(1−w) 。由于平方根项非负，因此可以确定不会生成小于a的z值，与瑞利CDF的定义所要求的一致。该生成器生成的均匀随机数w可以用在前面的方程中，以便生成一个具有瑞利分布且带有参数a、b的随机变量z。
语句*R=a+sqrt(b×log(1-rand(M, N)));*可以生成一个随机数数组R，rand可以生成区间(0, 1)上的均匀分布随机数。若M=N=1，则该语句由一个具有瑞利分布且以参数a和b来表征的随机变量生成一个单一值。
表达式 z = a + − b l n ( 1 − w ) z=a+\sqrt{-bln(1-w)} z=a+−bln(1−w) 有时称为随机数生成方程，因为它确定了生成期望随机数的方式。此情形下可以求出一个解析解，但并不总是可能的，其输出近似于具有给定CDF的随机数。
表4.1 随机变量的生成
| 名称 | 均值和方差 | CDF | 生成器* | |
|---|---|---|---|---|
| 均匀 | p ( z ) = { 1 b − a , 0 ≤ z ≤ b 0 , 其 他 p(z)=\begin{cases} \frac{1}{b-a},0\leq z \leq b \\ 0,其他 \end{cases} p(z)={ b−a1,0≤z≤b0,其他 | m = a + b 2 , σ 2 = ( b − a ) 2 12 m=\frac{a+b}{2},\sigma^2=\frac{(b-a)^2}{12} m=2a+b,σ2=12(b−a)2 | F ( z ) = { 0 , z < a z − a b − a , a ≤ z ≤ b 1 , z > b F(z)=\begin{cases} 0,z<a \\ \frac{z-a}{b-a},a\leq z \leq b \\ 1,z>b\end{cases} F(z)=⎩⎪⎨⎪⎧0,z<ab−az−a,a≤z≤b1,z>b | MATLAB函数rand |
| 高斯 | p ( z ) = 1 2 π b e − ( z − a ) 2 2 b 2 , − ∞ < z < ∞ p(z)=\frac{1}{\sqrt{2\pi}b}e^{-\frac{(z-a)^2}{2b^2}},-\infty<z<\infty p(z)=2π b1e−2b2(z−a)2,−∞<z<∞ | m = a , σ 2 = b 2 m=a,\sigma^2=b^2 m=a,σ2=b2 | F ( z ) = ∫ − ∞ z p ( v ) d v F(z)=\int_{-\infty}^{z} p(v)dv F(z)=∫−∞zp(v)dv | MATLAB函数randn |
| 对数正态 | p ( z ) = 1 2 π b z e [ l n ( z ) − a ] 2 2 b 2 , z > 0 p(z)=\frac{1}{\sqrt{2\pi}bz}e^{\frac{[ln(z)-a]^2}{2b^2}},z>0 p(z)=2π bz1e2b2[ln(z)−a]2,z>0 | m = e a + ( b 2 2 ) , σ 2 = [ e b 2 − 1 ] e 2 a + b 2 m=e^{a+(\frac{b^2}{2})},\sigma^2=[e^{b^2}-1]e^{2a+b^2} m=ea+(2b2),σ2=[eb2−1]e2a+b2 | F ( z ) = ∫ 0 z p ( v ) d v F(z)=\int_0^z p(v)dv F(z)=∫0zp(v)dv |