逆濾波與維納濾波（一）

圖像退化/復原過程的模型

 

 

 　　退化過程可以被模型化為一個退化函數和一個加性噪聲項，處理一幅輸入圖像 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

　　維納濾波請見下篇博客。