逆滤波与维纳滤波（一）

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

 

 

 　　退化过程可以被模型化为一个退化函数和一个加性噪声项，处理一幅输入图像 f(x, y) 产生一幅退化图像 g(x, y)。给定 g(x, y) 和关于退化函数 H 的一些知识以及外加噪声项 η(x, y)，图像复原的目的是获得关于原始图像的近似估计。

噪声模型

 

 

 图像复原

 

 

　　在退化复原模型中，输入输出关系可以表示为

　　        

　　如果系统 H 是一个线性系统，那么 H 应该满足可加性和均匀性。 H[f1(x, y) + f2(x, y)] = H[f1(x, y)] + H[f2(x, y)] H[a f1(x, y)] = a H[f1(x, y)] 这里，a 是比例常数，f1(x, y) 和 f2(x, y) 是任意两幅输入图像。

　　　要想要复原图像 需要知道退化函数，一般我们是不知道退化函数的，所以需要提前估计退化函数：

　大致有三种方法：

　　1.观察法

　　2.试验法

　　3.数学建模法

　　我们着重介绍获得 退化函数之后的操作：

　　1.逆滤波

　　2.维纳滤波

逆滤波：

　　在该方法中，用退化函数除退化图像的傅里叶变换来计算原始图像的傅里叶变换估计：

      　　　

　　考虑到噪声的影响

　　　

 弊端：

　　如果退化是零或非常小的值，N(u, v)/H(u, v) 之比很容易决定估计值。 一种解决退化是零或者很小值问题的途径，是限制滤波的频率使其接近原点值。

一：对指定的一幅灰度图像，先用3*3均值滤波器进行模糊处理，形成退化图像1；再叠加椒盐噪声，形成退化图像2；再对上述退化图像1和2采用逆滤波进行复原，给出复原结果图像。分析对比在对H零点问题采用不同处理方法下的复原结果。

1.1图像退化（均值滤波+椒盐噪声）

%---------------------9x9均值滤波模糊处理-------------------
clc;         %清空控制台
clear;       %清空工作区
close all;   %关闭已打开的figure图像窗口
color_pic=imread('lena.jpg');  %读取彩色图像
gray_pic=rgb2gray(color_pic);    %将彩色图转换成灰度图
double_gray_pic=im2double(gray_pic);   %将uint8转成im2double型便于后期计算
[width,height]=size(double_gray_pic);
H=fspecial('average',9);  %生成9x9均值滤波器,图像更模糊，3x3几乎看不出差别
degrade_img1=imfilter(double_gray_pic, H, 'conv', 'circular');  %使用卷积滤波，默认是相关滤波

%--------------------在均值滤波模糊图像基础上添加椒盐噪声------------------
degrade_img2=imnoise(degrade_img1,'salt & pepper',0.05);  %给退化图像1添加噪声密度为0.05的椒盐噪声(退化图像2)
figure('name','退化图像');
subplot(2,2,1);imshow(color_pic,[]);title('原彩色图');
subplot(2,2,2);imshow(double_gray_pic,[]);title('原灰度图');
subplot(2,2,3);imshow(degrade_img1,[]);title('退化图像1');
subplot(2,2,4);imshow(degrade_img2,[]);title('退化图像2');

 1.2直接逆滤波还原图像

%-------------------------对退化图像1逆滤波复原-------------------------
fourier_H=fft2(H,width,height);  %注意此处必须得让H从9x9变成与原图像一样的大小此处为512x512，否则ifft2 ./部分会报错矩阵不匹配
fourier_degrade_img1=fft2(degrade_img1);    %相当于 G(u,v)=H(u,v)F(u,v)，已知G(u,v),H(u,v)，求F(u,v)
restore_one=ifft2(fourier_degrade_img1./fourier_H);  %因为是矩阵相除要用./  
figure('name','退化图像1逆滤波复原');
subplot(2,2,1);imshow(im2uint8(degrade_img1),[]);title('退化图像1');
subplot(2,2,2);imshow(im2uint8(restore_one),[]);title('复原图像1');

%-------------------------对退化图像2直接逆滤波复原-------------------------
fourier_degrade_img2=fft2(degrade_img2); %相当于 G(u,v)=H(u,v)F(u,v)+N(u,v)
restore_two=ifft2(fourier_degrade_img2./fourier_H);
figure('name','退化图像1逆滤波复原');
subplot(2,2,3);imshow(im2uint8(degrade_img2),[]);title('退化图像2');
subplot(2,2,4);imshow(im2uint8(restore_two),[]);title('复原图像2');

 

　　在叠加椒盐噪声后，使得退化图像的方差变大了，所以噪声对图像复原的影响较大。

 

 

　　当有噪声时，直接进行逆滤波,此时输出为：
　　　　　　

                

  从上式可见，如果H(u,v)足够小，则输出结果被噪声所支配，复原出来的图像便是一幅充满椒盐噪声的图。
  所以当我们已知道所加的噪声信号时，还原出的信号为F ( u , v ) = G ( u , v ) − N ( u , v ) H ( u , v ) F(u,v)=\dfrac{G(u,v)-N(u,v)}{H(u,v)}F(u,v)=H(u,v)G(u,v)−N(u,v)​，即如果不知道噪声信号，则清晰还原困难。

2.对一幅灰度图像进行运动模糊并叠加高斯噪声，并采用维纳滤波进行复原

2.1图像退化(运动模糊+高斯噪声)

% ----------------2、对一幅灰度图像进行运动模糊并叠加高斯噪声，并采用维纳滤波进行复原-------------
clc;         %清空控制台
clear;       %清空工作区
close all;   %关闭已打开的figure图像窗口
color_pic=imread('lena.jpg');  %读取彩色图像
gray_pic=rgb2gray(color_pic);    %将彩色图转换成灰度图
double_gray_pic=im2double(gray_pic);   %将uint8转成im2double型便于后期计算
[width,height]=size(double_gray_pic);

%-------------------------添加运动模糊----------------------
H_motion = fspecial('motion', 18, 90);%运动长度为18，逆时针运动角度为90°
motion_blur = imfilter(double_gray_pic, H_motion, 'conv', 'circular');%卷积滤波
noise_mean=0;  %添加均值为0
noise_var=0.001; %方差为0.001的高斯噪声
motion_blur_noise=imnoise(motion_blur,'gaussian',noise_mean,noise_var);%添加均值为0，方差为0.001的高斯噪声
figure('name','运动模糊加噪');
subplot(1,2,1);imshow(motion_blur,[]);title('运动模糊');
subplot(1,2,2);imshow(motion_blur_noise,[]);title('运动模糊添加噪声');

　　对上述分别进行逆滤波操作：

fourier_H=fft2(H_motion,width,height);  %注意此处必须得让H从9x9变成与原图像一样的大小此处为512x512，否则ifft2 ./部分会报错矩阵不匹配
fourier_degrade_img1=fft2(motion_blur);    %相当于 G(u,v)=H(u,v)F(u,v)，已知G(u,v),H(u,v)，求F(u,v)
restore_one=ifft2(fourier_degrade_img1./fourier_H);  %因为是矩阵相除要用./  
figure('name','运动模糊逆滤波复原');
subplot(2,2,1);imshow(im2uint8(motion_blur),[]);title('退化图像1');
subplot(2,2,2);imshow(im2uint8(restore_one),[]);title('复原图像1');

%-------------------------对退化图像2直接逆滤波复原-------------------------
fourier_degrade_img2=fft2(motion_blur_noise); %相当于 G(u,v)=H(u,v)F(u,v)+N(u,v)
restore_two=ifft2(fourier_degrade_img2./fourier_H);
figure('name','运动模糊+高斯噪声逆滤波复原');
subplot(2,2,3);imshow(im2uint8(motion_blur_noise),[]);title('退化图像2');
subplot(2,2,4);imshow(im2uint8(restore_two),[]);title('复原图像2');

　　运动模糊

 

　　运动模糊+高斯噪声

 

 　上述复原图像中出现较大误差，已经完全看不出来原有图像，可能是因为在偏离中心处H(u,v) 比较小，容易受噪声的影响，所以可以采用低通滤波器。

fourier_H=fft2(H_motion,width,height);  %注意此处必须得让H从9x9变成与原图像一样的大小此处为512x512，否则ifft2 ./部分会报错矩阵不匹配
fourier_degrade_img1=fft2(motion_blur);    %相当于 G(u,v)=H(u,v)F(u,v)，已知G(u,v),H(u,v)，求F(u,v)
restore_one=ifft2(fourier_degrade_img1./fourier_H);  %因为是矩阵相除要用./  
h=fspecial('gaussian',width,5);
restore_three=imfilter(restore_one,h,'conv','circular');
figure('name','运动模糊逆滤波复原');
subplot(2,2,1);imshow(im2uint8(motion_blur),[]);title('退化图像1');
subplot(2,2,2);imshow(im2uint8(restore_three),[]);title('复原图像1');

%-------------------------对退化图像2直接逆滤波复原-------------------------
fourier_degrade_img2=fft2(motion_blur_noise); %相当于 G(u,v)=H(u,v)F(u,v)+N(u,v)
restore_two=ifft2(fourier_degrade_img2./fourier_H);
restore_four=imfilter(restore_two,h,'conv','circular');
figure('name','运动模糊+高斯噪声逆滤波复原');
subplot(2,2,3);imshow(im2uint8(motion_blur_noise),[]);title('退化图像2');
subplot(2,2,4);imshow(im2uint8(restore_four),[]);title('复原图像2');

 

 

　　参考：https://blog.csdn.net/qq_23023937/article/details/109518677

　　维纳滤波请见下篇博客。