[数字图像处理笔记] 第六章 图像复原

1. 图像复原及退化模型

1.1 图像复原

图像复原：根据退化原因，建立相应的数学模型，从被污染或畸变的图像信号中提取所需要的信息，沿着使图像降质的逆过程恢复图像本来面貌。

根据不同指标进行分类：

• 在给定退化模型条件下：无约束 和 有约束

• 根据是否需要外界干预：自动 和 交互

• 根据处理所在域：频率域 和 空间域

图像增强不需要考虑图像是如何退化的，而是试图采用各种技术来增强图像的视觉效果。

而图像复原需知道图像退化的机制和过程等先验知识，据此找出一种相应的逆处理方法，从而得到复原的图像。

如果图像已退化，应先作复原处理，再作增强处理。二者的目的都是为了改善图像的质量。

1.2 图像退化

图像退化：图像在形成传输过程中，由于成像系统、传输介质和设备的不完善，使图像的质量变坏。

原因一般有以下几点：

• 成像系统的像差、畸变、带宽有限导致的 图像失真。

• 成像器件拍摄姿态和扫描非线性引起的图像 几何失真。

• 成像传感器与被拍摄景物之间的相对运动，引起所成图像的 运动模糊。

• 灰度失真。光学系统或成像传感器本身特性不均匀，造成同样亮度景物成像灰度不同。

• 辐射失真。由于场景能量传输通道中的介质特性如大气湍流效应、大气成分变化引起图像失真。

• 图像在成像、数字化、采集和处理过程中引入的 噪声 等。

退化模型如下：

\[g(x, y) = H[f(x, y)] + n(x, y) \]

场景辐射能量在物平面上分布用 \(f(x, y)\) 描述，在通过成像系统 \(H\) 时，在像平面所得图像为 \(H[f(x, y)]\) ，一般称 \(H\) 为 退化函数。一般还有加性噪声 \(n(x, y)\)。

性质：

• 线性性

  \[H(k_1f_1(x, y) + k_2f_2(x, y)) = k_1H(f_1(x, y)) + k_2H(f_2(x, y)) \]

• 空间不变性

  \[H(f(x - a, y - b)) = g(x - a, y - b) \]

如果系统 \(H\) 是一个线性、位移不变性的过程，其点扩散函数用 \(h(x, y)\) 表示，\(f(x, y)\) 表示理想的、没有退化的图像。

那么在空间域中给出的退化图像可由下式给出：

\[g(x, y) = h(x, y) * f(x, y) + n(x, y) \]

其中， \(h(x, y)\) 是退化函数的空间描述，* 表示空间卷积。等价的频域描述：

\[G(u, v) = H(u, v)F(u, v) + N(u, v) \]

许多种退化都可以用线性位移不变模型来近似，这样线性系统中的许多数学工具如 线性代数，能用于求解图像复原问题，从而使运算方法简捷和快速。

当 退化不太严重 时，一般用线性位移不变系统模型 来复原图像，在很多应用中有较好的复原结果，且计算大为简化。

尽管实际 非线性和位移可变 的情况能更加准确而普遍地反映图像复原问题的本质，但在 数学上求解困难。只有在要求很精确的情况下才用位移可变的模型去求解，其求解也常以位移不变的解法为基础加以修改而成。

2. 噪声模型

2.1 概率密度函数

• 高斯噪声

  \[p(z) = \frac{1}{\sqrt{2 \pi}\sigma}\text{exp}(-\frac{(\mu - z)^2}{2 \sigma ^2}) \]

  其中 \(z\) 表示灰度值，\(\mu\) 表示灰度值的均值，\(\sigma\) 表示灰度值的标准差。

• 均匀分布

  \[p(z) = \begin{cases} \frac{1}{b - a}, & a \le z \le b \\ 0, & else \end{cases} \]

• 脉冲噪声（椒盐噪声）

  \[p(z) = \begin{cases} P_a, & z = a \\ P_b, & z = b \\ 0, & else \end{cases} \]

如图所示，分别为原始图像以及添加了三种噪声后的图像，以及对应的直方图。

相关代码：

I = imread('3.png');
I = rgb2gray(I);
I = im2double(I);

subplot(2, 4, 1);
imshow(I), title('initial');
subplot(2, 4, 5);
hist(I, 10);

I1 = imnoise(I, 'gaussian', 0.1);
subplot(2, 4, 2);
imshow(I1), title('gaussian');
subplot(2, 4, 6);
hist(I1, 10);

I2 = imnoise(I, 'speckle', 0.1);
subplot(2, 4, 3);
imshow(I2), title('speckle');
subplot(2, 4, 7);
hist(I2, 10);

I3 = imnoise(I, 'salt & pepper', 0.1);
subplot(2, 4, 4);
imshow(I3), title('salt & pepper');
subplot(2, 4, 8);
hist(I3, 10);

2.2 噪声参数估计

(1) 典型的周期噪声参数通过检测图像的 傅里叶谱 来估计。

(2) 噪声概率密度函数参数从 传感器的技术说明 中获取。

(3) 噪声概率密度函数参数由 *截取一组“平坦”环境的图像 估计。

3. 空间域滤波复原

均值滤波器：算术均值滤波器、几何均值滤波器、谐波均值滤波器 、逆谐波均值滤波器。

• 算术均值滤波器

  \[\hat f(x, y) = \frac{1}{mn}\sum_{(s, t) \in S_{xy}}g(s, t) \]

• 几何均值滤波器

  \[\hat f(x, y) = [\prod_{(s, t)\in S_{xy}}g(s, t)]^{\frac{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 > 0\)，滤波器用于消除 胡椒 噪声；

  当 \(Q < 0\)，滤波器用于消除 盐 噪声；

  当 \(Q = 0\)，逆谐波均值滤波器退化为算数均值滤波器

  当 \(Q = -1\)，逆谐波均值滤波器退化为谐波均值滤波器

算术均值滤波器 和 几何均值滤波器（尤其是后者）更适合处理 高斯噪声 或 均匀噪声 等随机噪声；

逆谐波均值滤波器 更适合处理 脉冲噪声，但必须知道噪声是暗噪声还是亮噪声，以便选择合适的 \(Q\) 符号。

如图所示：

相关代码：

I = imread('3.png');
I = rgb2gray(I);
I = im2double(I);

subplot(2, 3, 1);
imshow(I), title('initial img');

I_noise = imnoise(I, 'gaussian', 0.2);
subplot(2, 3, 2);
imshow(I_noise), title('noised img');

I_mean = imfilter(I_noise, fspecial('average', 3));
subplot(2, 3, 3);
imshow(I_mean), title('ave filtered');

I_mean = exp(imfilter(log(I_noise), fspecial('average', 3)));
subplot(2, 3, 4);
imshow(I_mean), title('jihe filtered');

Q = -1.5;
I_mean = imfilter(I_noise .^ (Q + 1), fspecial('average', 3)) ./ imfilter(I_noise .^ Q, fspecial('average', 3));
subplot(2, 3, 5);
imshow(I_mean), title('nixiebo filtered 1');

Q = 1.5;
I_mean = imfilter(I_noise .^ (Q + 1), fspecial('average', 3)) ./ imfilter(I_noise .^ Q, fspecial('average', 3));
subplot(2, 3, 6);
imshow(I_mean), title('nixiebo filtered 2');

顺序统计滤波器：中值滤波器、二维中值滤波器、修正后的阿尔法均值滤波器、最大值/最小值滤波器、中点滤波器。

• 中值滤波器

  • 一维

    设包围某点的一维数据集 \(A\)，将所有数据按照大小升序排序。对该点的中值滤波为：

    \[y = \text{Med}(x_1, x_2, ... , x_n) = \begin{cases} x^{'}_{\frac{n}{2}}, & n 为奇数 \\ \frac{1}{2}[x^{'}_{\frac{n}{2}} + x^{'}_{\frac{n}{2} + 1}], & n 为偶数 \end{cases} \]

  • 二维

    中值滤波可以去掉椒盐噪声，平滑结果优于均值滤波，在抑制随机噪声的同时保持图像边缘少受模糊。

    \[\hat f(x, y) = \text{Med}_{(s,t) \in S_{xy}}\{g(s, t)\} \]

  • 修正后的 \(\alpha\)

    在 \(S_{xy}\) 邻域内 去掉 \(g(s, t)\) 最高的 \(\frac{d}{2}\) 个灰度值 和 最低的 \(\frac{d}{2}\) 个灰度值，用 \(g(s, t)\) 来代替剩余的 \(mn - d\) 个像素。

    由这些像素点的均值形成的滤波器称为 修正后 \(\alpha\) 均值滤波器 。

    \[\hat f(x, y) = \frac{1}{mn - d}\sum_{(s,t) \in S_{xy}}g_r(s, t) \]

    当 \(d = 0\)，退化为算数均值滤波器；

    当 \(d = mn - 1\)，退化为中值滤波器；

    适用于高斯噪声和椒盐噪声的混合情况。

如图所示：

可以看出对于 椒盐噪声，中值滤波好于均值滤波 。

相关代码：

I = imread('3.png');
I = rgb2gray(I);
I = im2double(I);

subplot(2, 2, 1);
imshow(I), title('initial img');

I_noise = imnoise(I, 'salt & pepper', 0.1);
subplot(2, 2, 2);
imshow(I_noise), title('noised img');

I1 = imfilter(I_noise, fspecial('average', 5));
subplot(2, 2, 3);
imshow(I1), title('ave filtered img');

I1 = medfilt2(I_noise);
subplot(2, 2, 4);
imshow(I1), title('med filtered img');

• 最大/最小滤波器

  • 最大

    \[\hat f(x, y) = \max_{(s, t)\in S_{xy}}\{g(s, t)\} \]

    在发现最亮点时比较有用，胡椒噪声是比较低的值，此滤波器可以很好地消除胡椒噪声。

  • 最小

    \[\hat f(x, y) = \min_{(s, t)\in S_{xy}}\{g(s, t)\} \]

    在发现最暗点时比较有用，盐噪声是比较高的值，此滤波器可以很好地消除盐噪声。

• 中点滤波器

  \[\hat f(x, y) = \frac{1}{2}[\max_{(s, t) \in S_{xy}}\{g(s, t)\} + \min_{(s,t) \in S_{xy}}\{g(s, t)\}] \]

  结合了顺序统计和求平均的操作，对于高斯噪声和均匀随机分布这类噪声具有最好的滤波效果。

4. 频率域滤波复原

带阻、带通和陷波滤波器，它们能削减或消除 周期性噪声。

原理：空间域卷积相当于频率域乘积。可以在频率域中直接设计滤波器，对图像进行恢复处理。

4.1 带阻滤波器

• 理想带阻滤波器

  \[H(u, v) = \begin{cases} 1, & D(u, v) < D_0 - \frac{W}{2} \\ 0, & D_0 - \frac{W}{2} \le D(u, v) \le D_0 + \frac{W}{2} \\ 1, & D(u, v) > D_0 + \frac{W}{2} \end{cases} \]

  其中 \(W\) 是频带宽度，\(D_0\) 是频带中心半径，\(D(u, v)\) 是到中心化频率矩形原点的距离。

• \(n\) 阶巴特沃斯带阻滤波器

  \[H(u, v) = \frac{1}{1 + [\frac{D(u, v)W}{D^2(u, v) - D_0^2}]^{2n}} \]

• 高斯带阻滤波器

  \[H(u, v) = 1 - \text{exp}\{- \frac{1}{2}[\frac{D^2(u, v) - D_0^2}{D(u, v)W}]\} \]

4.2 带通滤波器

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

\[H_{bp}(u, v) = 1 - H_{br}(u, v) \]

4.3 陷波滤波器

5. 估计退化函数

5.1 图像观察估计法

使用目标和背景的样本灰度级，构建一个不模糊的图像，该图像和看到的子图像有相同的大小和特性。

用 \(g_s(x, y)\) 定义观察的子图像，用 \(\hat f(x, y)\) 表示构建的子图像。

\[H_s(u, v) = \frac{G_s(u, v)}{\hat F_s(u, v)} \]

5.2 试验估计法

\[H(u, v) = \frac{G(u, v)}{A} \]

其中 \(G(u, v)\) 是观察图像的傅里叶变换，\(A\) 是一个常量，表示 冲激强度。

5.3 模型估计法

• 大气湍流模型

  \[H(u, v) = \text{exp}[-k(u^2 + v^2)^{\frac{5}{6}}] \]

  其中 \(k\) 是常数。

• 运动模糊模型

  \[g(x, y) = \int_0^T f[x - x_0(t), y - y_0(t)] \text{d}t \]

  其中 \(g(x, y)\) 为模糊的图像。

  如图所示：

  

  相关代码：

  I = imread('3.png');
  I = rgb2gray(I);
  I = im2double(I);
  
  subplot(1, 2, 1);
  imshow(I), title('initial img');
  
  len = 25;
  theta = 11;
  
  I_psf = imfilter(I, fspecial('motion', len, theta), 'circular', 'conv');
  subplot(1, 2, 2);
  imshow(I_psf), title('blurred img');
  

6. 逆滤波

逆滤波（Inverse Filtering）是图片处理中常用的一种去除模糊的技术。在许多情况下，图片会因为多种原因（如相机抖动，雾气等）模糊。若可得知导致模糊的原因（数学上称为点扩散函数，或者说冲激响应），那么我们可以构造一种“逆滤波”去除原始图像的模糊。

I = imread('blurred_image.jpg');
I = rgb2gray(I);
subplot(1,2,1); imshow(I);

PSF = fspecial('motion', 21, 11);
filtered = deconvwnr(I, PSF, 0); 
subplot(1,2,2); imshow(filtered);

7. 最小均方误差滤波-维纳滤波

为了解决高噪声下的图像问题，最小均方滤波器来解决。

I = imread('blurred_image.jpg');
I = rgb2gray(I);

% 定义点扩散函数 PSF（Point Spread Function）
PSF = fspecial('motion', 21, 11);

% 定义噪声水平
noise_var = 0.001;

% Wiener滤波
filtered = deconvwnr(I, PSF, noise_var); 

subplot(1,2,1); imshow(I);;
subplot(1,2,2); imshow(filtered);

8. 几何失真校正

在图像的获取或显示过程中往往会产生 几何失真。

例如成像系统有一定的几何非线性。这主要是由于视像管摄像机及阴极射线管显示器的扫描偏转系统有一定的非线性，因此会造成如图所示的枕形失真或桶形失真。

8.1 空间变换

假设一幅图像为 \(f(x, y)\) ，经过几何失真变成了 \(g(u, v)\) ，这里的 \((u, v)\) 表示失真图像的坐标，此时已不是原图像的坐标。

\[\begin{split} u = r(x, y) \\\\ v = s(x, y) \end{split} \]

其中，\(r(x, y)\) 和 \(s(x, y)\) 是空间变换，产生了几何失真图像 \(g(u, v)\) 。在未知情况下，通常 \(r(x, y)\) 和 \(s(x, y)\) 可用多项式来近似。