冈萨雷斯Matlab版第五章（图像复原）知识点要点记录

前言

复原的目的是在预定义的意义上改善给定的图像。尽管图像增强与图像复原之间有重叠的部分，但前者主要是主观的处理。图像复原大部分是客观的处理。复原通过使用退化现象的先验知识试图重建或回复一幅退化的图像。因此，复原技术趋向于将退化模型化并用相反的处理来回复原图像，考虑用去模糊函数来消除图像的模糊。

5.1 图像退化/复原处理的模型

本章中用退化函数把退化过程模型化，它和加性噪声项一起，作用于输入图像f(x, y)，产生一幅退化的图像g(x, y):
g(x, y) = H[f(x, y)] + η(x, y)
给定g(x, y)，一些关于退化函数H的知识以及一些关于加性噪声η(x, y)的知识，复原的目标就是得到原图像的一个估计f‘，并使得这个估计尽可能接近原始的输入图像。

若H是线性、空间不变的过程，则退化图像在空间域通过下式给出。
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)
其中，大写字母表示的项是卷积方程式中对应项的傅里叶变换。
退化函数H(u, v)有时称为光学传递函数（OTF），在空间域，h(x, y)称为点扩散函数（PSF）。对于任何种类的输入，让h(x, y)作用于光源上的一个点来得到退化的特征，点扩散函数就是来源于此。MATLAB工具箱提供两个函数：otf2psf和psf2otf用于OTF和PSF之间的转换。有时退化过程称为“用PSF或OTF对图像进行卷积”，复原过程称为反卷积。

5.2 噪声模型

模拟噪声的行为和影响的能力是图像复原的核心。

5.2.1 使用函数`imnoise`添加噪声

工具箱采用imnoise来使用噪声污染一幅图像，基本语法为：
g = imnoise(f, type, parameters)
其中，f是输入图像，type和parameters使用如下。
函数imnoise在给图像添加噪声之前，将它转换为范围[0, 1]内的double类图像。

g = imnoise(f, 'gaussian', m, var)将均值为m，方差为var的高斯噪声加到图像f上，默认值m = 0, var = 0.01。
g = imnoise(f, 'localvar', V)将均值为0、局部方差为V的高斯噪声添加到图像f上，其中V是与f大小相同的一个数组，它包含了每个点的理想方差值。
g = imnoise(f, 'localvar', image_intensity, var)将均值为0的高斯噪声添加到图像f上，其中噪声的局部方差var是图像f的亮度值的函数，参量image_tensity和var是大小相同的向量，plot(image_intensity, var)绘制出噪声方差与图像亮度的函数关系，向量iamge_intensity必须包含范围在[0, 1]内的归一化亮度值。
g = imnoise(f, 'salt & pepper', d)用椒盐噪声污染图像f，其中d是噪声密度（即包含噪声值的图像区域的百分比）。因此，大约有d * numel(f)个像素受到了影响，默认噪声密度为0.65。
g = imnoise(f, 'speckle', var)用方程g = f + n * f将乘性噪声添加到图像f上，其中n是均值为0、方差为var的均匀分布的随机噪声。var的默认值是0.04。
g = imnoise(f, 'poisson')从数据中生成泊松噪声，而不是将人工的噪声添加到数据中。为了遵守泊松统计，uint8类和uint16类图像的亮度必须和光子的数量相符合。当每像素的光子数量大于65535但小于10¹²时，就要使用双精度图像。亮度值在0~1之间变化，并且对应于光子的数量除以10¹²。

5.2.2 使用指定的分布产生空间随机噪声

先验知识：若w是一个在区间(0, 1)内均匀分布的随机变量，则可以通过求解下面的方程来得到具有指定CDF（累积分布函数）和F_z的一个随机变量z：
z = F_z^-1(w)
即，寻找方程F_z(z) = w的解。
例5.1 使用均匀随机数生成具有批定分布的随机数
假设我们有一个其均匀随机数w在区间(0, 1)内的随机数生成器，并假设我们要用它来生成具有瑞利CDF的随机数z，它具有如下形式：

为得到z，可解方程

或

由于平方根项是非负的，所以可以确定不会产生值小于a的z，正如瑞利CDF所定义要求的那样，因此，我们的生成器产生的均匀随机数w可以用在前面的方程中，以便用来生成一个参数为a和b的瑞利分布的随机变量z。
在MATLAB中很容易推广到一个M * N的随机数数组R：
R = a + sqrt(b * log(1 - rand(M, N)))
其中，log是自然对数，rand生成一个大小为M * N的数组，数组的元素为（0， 1）区间内均匀分布的随机数。
TIPS：

随机数生成器方差
表达式

有时称为随机数生成器方程，因为它确定了如何生成需要的随机数。
randn(M, N)
生成一个大小为M*N的数组，它的元素是零均值、单位方差的正态（高斯）数。
find(A)
I = find(A)返回数组A中所有非零元素的索引，若一个也没有，返回空矩阵，以格式A(:)处理数组A，所以I是一个列向量。
[r, c] = find(A)返回矩阵A的非零元素的行和列索引。
[r, c, v] = find(A)除了返回行索引和列索引外，还以列向量v返回A的非零值。
随机变量的生成表
imnoise2函数
R = imnoise2(type, M, N, a, b)
产生一个大小为M * N的噪声数组R，用户可直接指定需要的噪声参数，注意，由椒盐噪声产生的噪声数组有三个值：相应于胡椒噪声的0，相应于盐粒噪声的1，以及相应于无噪声的0.5。
源码见【MATLAB】图像处理相关函数imnoise2函数

5.2.3 周期噪声

图像的周期噪声一般产生于图像采集过程中的电气和/或电机干扰。这是本章唯一考虑的一种空间从属噪声。我们的周期噪声的模型是二维正弦波，它有如下所示的方程：

其中，A是振幅，u₀和v₀分别关于x轴和y轴确定正弦频率。B_x和B_y是关于原点的相移。该方程的M*N DFT为

我们所看到的是一对分别位于（u + u₀, v + v₀）和（u - u₀, v - v₀）的复共轭冲击。

imnoise3函数
接受冲击位置（频率坐标）的一个任意数，每个冲击位置都有自己的振幅、频率和相移参数，并且以前一段中所描述的正弦和的形式计算r(x, y)。这个函数也输出正弦波之和的傅里叶变换R(u, v)以及R(u, v)的谱。正弦波通过逆DFT从给定的冲击位置中产生。这使得它更加直观，而且也简化了空间噪声模式中的频率显示。确定一个冲击的位置仅需要一对坐标。这个程序将产生共轭对称冲击。
源码见【MATLAB】图像处理相关函数imnoise3函数。

5.2.4 估计噪声参数

假设z_i是用来表示一幅图像的灰度级的一个离散随机变量，令p（z_i）， i = 0, 1,2 ···, L - 1 是相应的归一化直方图，其中L是可能的亮度值的数目。直方图分量p（z_j）是亮度值z_j出现的概率的一个估计，该直方图也可以看成是亮度PDF的近似。
描述直方图分布形状的一种主要方法是通过它的中心距（即均值的矩），它定义为

其中，n是矩的阶，m是均值。

均值m

因为假设直方图已归一化，其所有分量之和为1，所以由前面的方程可知，μ₀ = 1和μ₁ = 0.
方差μ₂
二阶矩

是方差。
函数statmoments
函数statmoments计算均值和n阶中心矩，并返回行向量v。因为μ₀ = 1和μ₁ = 0. 所以statmoments忽略这两个矩，改为令v(1) = m, v(k) = μk， k = 2, 3, ···， n，语法如下：
[v, unv] = statmoments(p, n)
其中， p为直方图向量，n是计算的矩的数量。要求对于uint8类图像，p的分量数等于2⁸;对于uint16类图像，p的分量数等于2¹⁶；对于double类图像，p的分量数等于2⁸或2¹⁶。输出向量v包含了以随机变量值为基础的归一化矩，而随机变量已被标度区间[0, 1]内，所以，所有的矩也在这个区间内。向量unv包含了与v相同的矩，但用位于原始值区间内的数值计算。例如，若length(p) = 256, v(1) = 0.5, so unv(1) = v(1) * （length(p) - 1） = 127.5
源码见
函数roipoly
噪声参数通常直接由带噪声的图像或一组带噪声的图像来估计，在这种情况下，所选择的方法是，在一幅图像中选择一个尽可能与背景一样无特色的区域，以便使该区域的亮度值的可变性主要由噪声产生。在MATLAB中，我们使用函数roipoly来选择一个感兴趣的区域（ROI），该函数将产生一个多边形的ROI。语法为：
B = roipoly(f, c, r)
其中，f是我们感兴趣的图像，c和r是多边形顶点的相应列坐标和行坐标（列坐标先被指定），输出B是一个二值图像，它是一个与f大小相同且在感兴趣的区域外为0、在感兴趣区域内为1的二值图像。图像B是一个用来将操作限制在感兴趣区域内的掩模。
为了交互式地指定一个多边形，可以使用以下语法：
B = roipoly(f)
它将把图像f显示到屏幕上，让用户使用鼠标来指定多边形。若省略f，则函数将在最后显示的图像上操作
按Backspace键或者Delete键可以删除先前选中的顶点；
按Shift键加上单击、右击、或双击鼠标可为选区添加最后一个顶点，然后开始用1 来填充多边形区域；
按Return键则会结束选区而不添加任何顶点。
[B, c, r] = roipoly(...)可以得到二值图像和多边形的顶点列表，其中roipoly(...)表示该函数任何有效语法。这种格式在交互式地指定ROI时特别有用。
histroi函数
用于计算一幅图像在多边形区域内的直方图。源码见

5.3 仅有噪声的复原：空间滤波

若出现的退化仅仅是噪声，在这种情况下，所选择的降低噪声的方法是空间滤波。本节中，我们将总结和实现几种降低噪声的空间滤波器。

5.3.1 空间噪声滤波器

表5.2列出了本节所感兴趣的空间滤波器，其中S_xy表示输入的噪声图像g的m*n子图像（区域）。S的小标表示在子图像中心的坐标（x， y），f’(x, y)（f的一个估值）表示滤波器在这些坐标处的响应。

spfilt函数
spfilt函数与表5.2中列出的任何滤波器在空间域执行滤波，注意表中提到的函数imlincomb的用途是计算输入的线性组合，该函数的语法是：
B = imlincomb(c1, A1, c2, A2, ..., ck, Ak)
即实现等式B = c1 * A1 + c2 * A2 + ...+ ck * Ak
其中，c是double类实标量，A是有着相同类和大小的数组。
最后测试得到

=源码见

5.3.2 自适应空间滤波器

一个自适应中值滤波器的详细算法如下所示

（1）令
z_min表示S_xy中的最小亮度值
z_max表示S_xy中的最大亮度值
z_med表示S_xy中的亮度中值
z_xy表示坐标（x，y）处的亮度值
这个自适应中值滤波算法工作在两个层面，表示为level A和level B：
（2）Level A：若z_min< z_med<z_max则转向Level B
否则增加窗口尺寸
若窗口尺寸≤S_max，重复Level A
否则输出z_med
（3）Level B：若z_min< z_med<z_max，则输出z_xy
否则输出z_med

其中，S_max表示允许的最大自适应滤波器窗口的大小，Level A最后一步的另一种选择是输出z_xy来代替中值。

adpmedian函数实现算法
f = adpmedian(g, Smax)
其中，g是将被滤波的图像，Smax是允许的最大自适应滤波器窗口的大小。
最后测试得到

源码见【MATLAB】图像处理相关函数adpmedian函数

5.4 通过频域滤波来降低周期噪声

n阶巴特沃斯陷波滤波器
传递函数为：

其中，

``(u0, v0)和(-u-, -v0)是“陷波”的位置，D₀是它们的半径。注意，滤波器指定了频率矩形的中心，所以它必须在使用前用函数fftshift进行预处理。

5.5 退化函数建模

图像模糊
在图像获取时传感器和场景之间的均匀线性运动产生的图像模糊。

对图像模糊建模

fspecial函数
PSF = fspecial('motion', len, theta); 近似于由有着len个像素的摄像机的线性移动的效果。参数theta以度为单位，以顺时针方向对正水平轴度量。len默认值是9， theta默认值是0，在水平方向上它对应于9个像素的移动。
imfilter函数
g = imfilter(f, PSF, 'circular'); 创建一个已知PSF的退化图像，其中，‘circular’用来减少边界效应，然后通过添加适当的噪声来构造退化的图像模型：
g = g + noise;
其中noise是一幅与g大小相同的随机噪声图像。
checkerboard函数
C = checkerboard(NP, M, N); 产生测试图案，且图案的大小可以被缩放而不影响它的主要特征。其中，NP是每一个正方形一边的像素值，M是行数，N是列数，若省略N，则默认N = M；若省略M和N，则将产生一个一面为8个正方形的正方形测试板；若省略NP，默认NP= 10.测试板左半部分的亮正方形是白色，右半部分的亮正方形是灰色，命令：
K = im2double(checkerboard(NP, M, N)) >0.5;
产生一个亮正方形都是白色的测试板。
pixeldup函数
B = pixeldup(A, m, n)将A的每个像素在垂直方向上总共复制了m次，在水平方向上总共复制了n次，若省略n，默认n = m。
例5.7 模糊噪声图像的建模
源码见github
测试得到

5.6 直接逆滤波

构成如下形式的估计的这种方法称为逆滤波。

5.7 维纳滤波

是最早也是人们最熟知的线性图像复原方法。

算法
维纳滤波器寻找一个使统计误差函数

最小的估计f‘，其中，E是期望值操作符，f是未退化的图像。该表达式在频域可表示为：

其中

比率S_η（u， v）/ S_f(u, v)称为噪信功率比。若关于u和v的所有相关值，噪声功率谱为零，则这个比率就变为零，且维纳滤波器就成为逆滤波器。
我们感兴趣的两个量是平均噪声功率和平均图像功率（标准常量），分别定义为：

其中，M和N分别表示图像和噪声数组的垂直和水平大小，比率为：

也是一个标量。
实现方法
deconvwnr函数，语法如下：
1、fr = deconvwnr(g, PSF); g表示退化图像，fr表示复原图像，信噪比为零。
2、fr = deconvwnr(g, PSF, NSPR); g表示退化图像，fr表示复原图像，信噪比已知，或是个常量或是个数组，在这个情况下，NSPR可以是一个交互的标量输入。
3、fr = deconvwnr(g, PSF， NACORR， FACORR); g表示退化图像，fr表示复原图像，假设噪声和未退化图像的自相关函数NACORR和FACORR是已知的。
若复原图像呈现出由算法中使用的离散傅里叶变换所引入的振铃，则往往在调用deconvwnr函数之后要调用edgetaper函数，语法是：
J = edgetaper(I, PSF)
该函数利用点扩散函数PSF模糊了输入图像I的边缘。这样，输出图像J就是图像I与I的模糊版本的加权和。
举例
例5.8 使用deconvwnr函数复原模糊噪声图像
得到的结果如下：

源码见github

5.8 约束的最小二乘方（正则）滤波

实现线性复原。
这个二维离散卷积的定义如下：

利用这个公式，我们可以用向量矩阵的形式将线性退化模型g(x, y) = h(x, y) * f(x, y) + η(x, y)表示为：
g = Hf + η
假设g(x, y)的大小是MN，则我们可以用g(x, y)的第一行的图像元素形成向量g的前N个元素，用第二行的图像元素形成向量g的下N个元素，以此类推。最终的向量维数为MN1.
准则函数C定义为：

该函数的约束条件为：

其中，||w||^2 = w^T w是欧几里得向量范数，f‘是非退化图像的估计。这个最优化问题的频域解决办法由下式得到：

其中，γ是一个必须加以调整的参量，以便约束条件得到满足，P(u, v)是函数（称为拉普拉斯算子）

的傅里叶变换。

deconvreg函数
语法为：
fr = deconvreg(g, PSF, NOISEPOWER, RANGE)
其中，g是被污染的图像，fr是复原的图像，NOISEPOWER与||η||^2成比例，RANGE为值的范围，在求解γ时，该算法受中国范围的限制。默认范围是[10^-9, 10⁹]，若将上述两个参数排除在参量之外，则函数deconvreg就会产生一个逆滤波方案。NOISEPOWER的较好初始估计为MN[σ_η² + m_η²]，其中M和N代表图像的维数，括号中的参量代表噪声方差和噪声均值的平方。该估计只是一个简单地开始。
例5.9 使用函数deconvreg复原模糊噪声图像
源码见github
得到结果

5.9 使用Lucy-Richardson算法的迭代非线性复原

L-R算法是从最大似然公式中引出来的，在这个方程中，图像是用泊松统计加以模型化。当下面这个迭代收敛时，模型的最大似然函数可以得到一个令人满意的方程：

其中，*代表卷积， f’是未退化图像的估计。

deconvlucy函数
fr = deconvlucy(g, PSF, NUMIT, DAMPAR, WEIGHT);
其中，fr代表复原的图像，g代表退化的图像，PSF是点扩散函数， NUMIT是迭代的次数（默认值10），DAMPAR（默认值0，表示无衰减）是一个标量，指定了结果图像与原图像g之间的偏离阈值。当像素值偏离原值的范围在DAMPAR之内时，就不用再迭代。这既抑制了这些像素上的噪声，又保存了必要的图像细节。
WEIGHT是一个与g同样大小的数组，它为每一个像素分配一个权重来反映其重量。例如从某个有缺陷的成像数组得出的一个不良像素最终会被赋以权重值零，从而排除该像素来求解；也可以根据平坦区域修正量调整像素的权重。当用一个指定的PSF来模拟模糊时，WEIGHT可以从计算像素中剔除那些来自图像边界的像素点，因此，PSF造成的模糊是不同的。若PSF的大小为n*n，则在WEIGHT中用到的零边界的宽度是ceil(n/2)。默认值是同输入图像g同等大小的一个单位数组。
例5.10 使用deconvlucy复原模糊噪声图像

源码见Github

5.10 盲去卷积

在盲去卷积中，最优化问题用规定的约束条件并假定收敛时通过迭代来求解，得到的最大f(x, y)好h(x, y)就是还原的图像和PSF。

deconvblind函数
[f, PSFe] = deconvnlind(g, INITPSF);
其中，g表示退化图像，INITPSF是点扩散函数的初始估计，PSFe是这个函数最终计算得到的估计值，fr是利用估计的PSF复原的图像。PSF估计受其初始推测尺寸的巨大影响，而很少受其他值的影响。同样，这个函数也包含附加参数来控制迭代的次数和复原的其他特性，如：
[fr, PSFe] = deconvblind(g, INTPSF, NUMIT, DAMPAR, WEIGHT);
例5.11 使用函数deconvblind估计PSF

源码见github

5.11 几何变换与图像配准

几何变换常常用来实现图像配准。图像配准就是取两张相同场景的图像并加以对准，从而可以将它们合并，以便目测或者定量比较。

5.11.1 空间几何变换

假设有一幅定义在(w, z)坐标系上的图像f经过几何变形后产生了定义在(x, y)坐标系上的图像g。(注意：y表示行，x表示列)。这个变换可以表示为
(x, y) = T{(w, z)}
仿射变换是最常用的一种空间变换形式，仿射变换可以用如下矩阵形式表示：

该变换可以按比例缩放、旋转、平移或者剪切，具体取决于T的元素的取值，具体如下：

maketform函数
tform = maketform(transform_type, transform_parameters);
第一个输入变量transform_type是'affine', 'projective', 'box', 'composite', 'custom'字符串之一，transform_parameters是依赖于变换类型的附加参量，tform结构表示空间变换。
tformfwd函数
XY = tformfwd(WZ, tform);
计算正变换T(w, z)，其中wz是一个大小为P2的点矩阵。WZ的每一行都包含一个点的w坐标和z坐标，类似地，XY是一个大小为P2的点矩阵，每一行都包含变换点的x和y坐标。
tforminv函数
WZ = tforminv(XY, tform);
计算逆变换T^-1{(x, y)}，含义同上。
自定义函数vistformfwd函数
vistformfwd(tform, wdata, zdata, N);
建立点的栅格结构，使用函数tformfwd对栅格进行变换，然后将原来的栅格和变换后的栅格紧挨着画在一起，以便对比。效果如下：
缩放

水平剪切

旋转

函数源码见：【MATLAB】图像处理相关函数tformfwd函数
函数使用例程见github

5.11.2 对图像应用空间变换

主要有两类：

正向映射
按顺序扫描每个输入像素，由T{(w, z)}决定的位置上直接将值复制给输出图像上的同一个像素点。
反向映射
按顺序扫描每个输出像素，用T^-1{(x, y)}计算出输入图像上的对应位置，并将其插值到最近的输入像素之中，以决定输出像素的值。
imtransform函数，使用反向传播
g = imtransform(f, tform, interp);
其中， interp是一个字符串，用来指明输入像素怎样插值来获取输出像素，取值如下：'nearest', ''bilinear', 'bicubic'，默认值'bilinear'。
例5.13 图像空间变换
主要讨论图像变换的几种情况，仿射变换矩阵有如下形式：

效果如下：
原图

图像变换

其中最后一个得到和原图一模一样的结果，这是由于在默认情况下在边框只进行反向映射变换，这种方法可以有效撤销平移，通过详细指定参量XData和YData的值，XData是一个包含两个元素的向量，指定输出图像的最左列和最右列的位置；YData也是一个包含两个元素的向量，指定输出图像的最顶部行和最底部行的位置。
所以使用以下语句得到平移效果
g5 = imtransform(f, tform2, 'XData', [1 400], 'YData', [1 400], 'FillValue', 0.5);
平移

源码见github

5.11.3 图像配准

图像配准法是寻求将两幅相同场景的图像加以对准。工具箱支持以控制点为基础的图像配准，控制点也称联结点，即在两幅图像中的位置已知或可以交互地选择的像素的一个子集。

cp2tform函数和maketform函数的变换类型