Matlab数字图像处理——图像复原

图像复原技术

图像在采集、传送、转换过程中会混入噪声，造成图像模糊、失真、有噪声，图像复原技术的目的是使图像尽可能恢复到原本的样子。

一、常见的噪声类型

• 高斯噪声
• 椒盐噪声
• 均匀分布噪声
• 指数分布噪声
• 伽马分布噪声

在Matlab中，函数 imnoise 可以给图像加入高斯噪声、椒盐噪声、泊松分布噪声和乘性噪声

二、空域内的滤波复原方法

1、均值滤波

(1)、算术均值滤波

令 S_xy 表示中心在点(x, y)处、大小为 m * n的矩形子窗口区域(邻域)的一组坐标，然后计算这一组坐标所在像素值的平均值，将其赋值给点 (x, y)。

通俗点讲，就是选图像上的一块区域，把这块区域的所有像素点相加求平均值再赋值给这块区域中心的那个像素点，后边的几种方法只是算法不同而已。

公式为：$$\hat{f}(x,y)=\frac{1}{mn}\sum _{(s,t)\in S_{xy}}g(s,t)$$

(2)、几何均值滤波

算术均值滤波是把 S_xy 所代表点的像素值求和再取平均值，而几何均值滤波是把像素值求乘再求像素的数量次幂。几何均值滤波实现的平滑可与算术均值滤波相比，且丢失的图像细节更少
公式为：$$\hat{f}(x,y)={\left[\prod _{(s,t)\in S_{xy}}g(s,t)\right]}^{\frac{1}{mn}}$$

(3)、谐波均值滤波器

谐波均值滤波器对于盐粒噪声效果很好，但不适用于胡椒噪声，擅于处理高斯噪声
公式为：$$\hat{f}(x,y)=\frac{mn}{{\sum }_{(s,t)\in S_{xy}}\frac{1}{g(s,t)}}$$

(4)、逆谐波均值滤波器

Q 称为滤波器的阶数，适用于减少或在实际中消除椒盐噪声干扰的影响。当 Q 为正时，该滤波器可消除胡椒噪声，当 Q 为负值时该滤波器可消除盐粒噪声，当 Q 为 0 时，公式可化简为算术均值滤波，当 Q = -1 时则为谐波均值滤波器。
公式为：$$\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}$$

下列代码为算术均值滤波和几何均值滤波的实例代码及结果：

I = imread('火影9.jpg');
I = im2double(I);
I = imnoise(I, 'gaussian', 0.05);   %加入高斯噪声
%算术均值滤波
PSF = fspecial('average', 3);   %创建预定义的二维滤波器
J = imfilter(I, PSF);
%几何均值滤波
K = exp(imfilter(log(I), PSF));

figure;
subplot(131);imshow(I);title('加入高斯噪声');
subplot(132);imshow(J);title('算术均值滤波');
subplot(133);imshow(K);title('几何均值滤波');

冈萨雷斯大佬的《数字图像处理》那本书上确实写了几何均值滤波比算术均值滤波要保留更多细节，我tm咋感觉是假的呢 o(一︿一+)o，结果如下：

2、统计排序滤波

其实我觉得和均值滤波的原理差不多，都是选取一块区域，把这块区域里的像素们(一长串数字)进行各种运算操作。

(1)、中值滤波器

中值滤波函数：medfilt2

正如其名字一样，使用一个像素邻域中的灰度级的中值来替代该像素的值。该滤波器针对单极或双极脉冲噪声尤其有效
公式为： f ^ ( x , y ) = m e d i a n ( s , t ) ∈ S x y { g ( s , t ) } \hat{f}(x, y) = median_{(s, t) \in S_{xy}} \left\{ g(s, t) \right\} f^​(x,y)=median(s,t)∈Sxy​​{g(s,t)}

I = imread('火影9.jpg');
I = im2double(I);
I = imnoise(I, 'salt & pepper', 0.05);   %加入椒盐噪声

J = medfilt2(I, [3, 3]);    %中值滤波

figure;
subplot(121);imshow(I);title('加入椒盐噪声');
subplot(122);imshow(J);title('中值滤波');

从结果可以看出，图像很好的去除了椒盐噪声，但是图像有些模糊了：

(2)、最大值滤波器

就是对邻域中的所有像素点取最大值赋值给该像素点，最大值滤波器对于发现图像中的最亮点非常有用！另外它会从黑色物体边缘去除一些黑色像素
公式为： f ^ ( x , y ) = m a x ( s , t ) ∈ S x y { g ( s , t ) } \hat{f}(x, y) = max_{(s, t) \in S_{xy}} \left\{ g(s, t) \right\} f^​(x,y)=max(s,t)∈Sxy​​{g(s,t)}

(3)、最小值滤波器

就是对邻域中的所有像素点取最小值赋值给该像素点，最小值滤波器对于发现图像中的最暗点非常有用！另外它会从白色物体边缘去除一些白色像素
公式为： f ^ ( x , y ) = m i n ( s , t ) ∈ S x y { g ( s , t ) } \hat{f}(x, y) = min_{(s, t) \in S_{xy}} \left\{ g(s, t) \right\} f^​(x,y)=min(s,t)∈Sxy​​{g(s,t)}

(4)、中点滤波器

选取最大值与最小值的中点(就是最大最小值的平均)替换该点的像素值，结合了统计排序和求平均，对于处理随机噪声表现很好，比如高斯噪声或均匀噪声。
公式为： 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) = \frac{1}{2} [ max_{(s, t) \in S_{xy}} \left\{ g(s, t) \right\} + min_{(s, t) \in S_{xy}} \left\{ g(s, t) \right\} ] f^​(x,y)=21​[max(s,t)∈Sxy​​{g(s,t)}+min(s,t)∈Sxy​​{g(s,t)}]

• 二维顺序统计滤波函数：ordfilt2，可以实现中值、最大最小值滤波

I = imread('火影9.jpg');
I = rgb2gray(I);
I = im2double(I);
I = imnoise(I, 'salt & pepper', 0.1);   %加入椒盐噪声

domain = [1 1 1 1; 1 1 1 1; 1 1 1 1; 1 1 1 1];
J1 = ordfilt2(I, 1, domain);
J2 = ordfilt2(I, 8, domain);
J3 = ordfilt2(I, 16, domain);

figure;
subplot(221);imshow(I);title('加入椒盐噪声');
subplot(222);imshow(J1);title('最小值滤波');
subplot(223);imshow(J2);title('中值滤波');
subplot(224);imshow(J3);title('最大值滤波');

结果展示：

3、自适应滤波

二维自适用去噪滤波：wiener2，自适应滤波器的性能比目前所有滤波器的性能都好

I = imread('火影9.jpg');
I = rgb2gray(I);

%I = imcrop(I, [100, 100, 200, 200]);    %图像裁剪

I = imnoise(I, 'gaussian', 0, 0.003);   %加入高斯噪声
[K, noise] = wiener2(I, [5, 5]);

figure;
subplot(121);imshow(I);title('加入高斯噪声');
subplot(122);imshow(K);title('自适应滤波');

可以看出，效果还不错：

三、图像复原方法

1、逆滤波

大意就是我们知道了原图受到了什么算法(叫做退化函数 H )带来的噪声，那么我们从噪声图反着用这个算法就可以推算出原图。
不过即使知道了退化函数，也还是不能完全恢复到原图。

I = imread('cameraman.tif');
I = im2double(I);
[m, n] = size(I);
M = 2 * m;
N = 2 * n;

u = -m / 2:(m / 2 - 1);
v = -n / 2:(n / 2 - 1);
[U, V] = meshgrid(u, v);           %基于向量 u 和 v 中包含的坐标返回二维网格坐标
D = sqrt(U.^2 + V.^2);
D0 = 130;                           %滤波器截止频率为20
H = exp(-(D.^2)./(2 * (D0^2)));
N = 0.01 * ones(size(I, 1), size(I, 2));
N = imnoise(N, 'gaussian', 0, 0.001);
J = fftfilter(I, H) + N;

HC = zeros(m, n);
M1 = H > 0.1;
HC(M1) = 1 ./ H(M1);
K = fftfilter(J, HC);
HC = zeros(m, n);
M2 = H > 0.01;
HC(M2) = 1 ./ H(M2);
L = fftfilter(J, HC);

figure;
subplot(221);imshow(J);
subplot(222);imshow(J, []);
subplot(223);imshow(K, []);
subplot(224);imshow(L, []);

其中，fftfilter函数如下：

function Z = fftfilter(X, H)
F = fft2(X, size(H, 1), size(H, 2));
Z = H .* F;
Z = ifftshift(Z);
Z = abs(ifft2(Z));
Z = Z(1:size(X, 1), 1:size(X, 2));
end

完整目录

Matlab数字图像处理——图像处理工具箱Image Processing Toolbox
Matlab数字图像处理——图像类型的转换
Matlab数字图像处理——图像文件的读取
Matlab数字图像处理——图像文件的显示
Matlab数字图像处理——视频文件的读写
Matlab数字图像处理——图像的像素运算(灰度变换)
Matlab数字图像处理——图像的空间变换
Matlab数字图像处理——图像的平移、邻域操作、区域选取
Matlab数字图像处理——图像增强
Matlab数字图像处理——图像复原