數字影像處理-運動模糊&逆濾波&維納濾波(Matlab)

1、對指定的一幅灰度影像，先用3*3均值濾波器進行模糊處理，形成退化影像1；再疊加椒鹽噪聲，形成退化影像2；再對上述退化影像1和2采用逆濾波進行復原，給出復原結果影像，分析對比在對H零點問題采用不同處理方法下的復原結果，

1-1 影像退化(均值濾波+椒鹽噪聲)

• Matlab代碼塊：

%---------------------9x9均值濾波模糊處理-------------------
clc;         %清空控制臺
clear;       %清空作業區
close all;   %關閉已打開的figure影像視窗
color_pic=imread('lena512color.bmp');  %讀取彩色影像
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 直接逆濾波還原影像

• Matlab代碼塊：

%-------------------------對退化影像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(1,2,1);imshow(im2uint8(degrade_img1),[]);title('退化影像1');
subplot(1,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);

• 運行效果：

1-3 去掉噪聲分量逆濾波還原影像

• Matlab代碼塊：

%------------------------去掉噪聲分量逆濾波復原-----------------------
noise=degrade_img2-degrade_img1;   %提取噪聲分量
fourier_noise=fft2(noise);   %對噪聲進行傅里葉變換
restore_three=ifft2((fourier_degrade_img2-fourier_noise)./fourier_H);  %G(u,v)=H(u,v)F(u,v)+N(u,v),解得F(u,v)=[G(u,v)-N(u,v)]/H(u,v)
figure('name','退化影像2逆濾波復原');
subplot(2,2,1);imshow(double_gray_pic,[]);title('原灰度圖');
subplot(2,2,2);imshow(im2uint8(degrade_img2),[]);title('退化影像2');
subplot(2,2,3);imshow(im2uint8(restore_two),[]);title('直接逆濾波復原');
subplot(2,2,4);imshow(im2uint8(restore_three),[]);title('去掉噪聲分量逆濾波復原');

• 運行效果：

• 總結:逆濾波依賴噪聲屬性
  當有噪聲時，直接進行逆濾波,此時輸出為： G ( u , v ) H ( u , v ) \dfrac{G(u,v)}{H(u,v)} H(u,v)G(u,v)?

G ( u , v ) = H ( u , v ) F ( u , v ) + N ( u , v ) G(u,v)=H(u,v)F(u,v)+N(u,v) G(u,v)=H(u,v)F(u,v)+N(u,v)

G ( u , v ) H ( u , v ) = F ( u , v ) + N ( u , v ) H ( u , v ) \dfrac{G(u,v)}{H(u,v)}=F(u,v)+\dfrac{N(u,v)}{H(u,v)} H(u,v)G(u,v)?=F(u,v)+H(u,v)N(u,v)?
??從上式可見，如果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 影像退化(運動模糊+高斯噪聲)

• Matlab代碼塊：

% ----------------2、對一幅灰度影像進行運動模糊并疊加高斯噪聲，并采用維納濾波進行復原-------------
clc;         %清空控制臺
clear;       %清空作業區
close all;   %關閉已打開的figure影像視窗
color_pic=imread('lena512color.bmp');  %讀取彩色影像
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('運動模糊添加噪聲');

• 運行效果：

2-2 維納濾波(Matlab自帶函式：deconvwnr)

• 函式引數介紹：

J = deconvwnr(I,PSF,NSR) deconvolves image I using the Wiener filter algorithm, returning deblurred image J. Image I can be an N-dimensional array. PSF is the point-spread function with which I was convolved. NSR is the noise-to-signal power ratio of the additive noise. NSR can be a scalar or a spectral-domain array of the same size as I. Specifying 0 for the NSR is equivalent to creating an ideal inverse filter.

• Matlab代碼塊：

%----------------------------維納濾波matlab自帶函式deconvwnr-------------------------
restore_ignore_noise = deconvwnr(motion_blur_noise, H_motion, 0);  %nsr=0,忽視噪聲
signal_var=var(double_gray_pic(:));
estimate_nsr=noise_var/signal_var;   %噪信比估值
restore_with_noise=deconvwnr(motion_blur_noise,H_motion,estimate_nsr);  %信號的功率譜使用影像的方差近似估計
figure('name','函式法維納濾波');
subplot(1,2,1);imshow(im2uint8(restore_ignore_noise),[]);title('忽視噪聲直接維納濾波(nsr=0)，相當于逆濾波');
subplot(1,2,2);imshow(im2uint8(restore_with_noise),[]);title('考慮噪聲維納濾波');

• 運行效果：

2-3 根據公式撰寫維納濾波

• 公式原理：

• Matlab代碼塊:

%---------------------------公式法----------------------------
fourier_H_motion=fft2(H_motion,width,height);  %H(u,v)
pow_H_motion=abs(fourier_H_motion).^2;   %|H(u,v)|^2
noise=motion_blur_noise-motion_blur;    %提取噪聲分量
fourier_noise=fft2(noise);    % N(u,v)  噪聲傅里葉變換
fourier_double_gray_pic=fft2(double_gray_pic);  %F(u,v)為未經過退化的圖片
nsr=abs(fourier_noise).^2./abs(fourier_double_gray_pic).^2;   %噪信比=|N(u,v)|^2/|F(u,v)|^2
H_w=1./fourier_H_motion.*pow_H_motion./(pow_H_motion+nsr); %H_w(u,v)=1/H(u,v)*|H(u,v)|^2/[|H(u,v)|^2+NSR] 
fourier_motion_blur_noise=fft2(motion_blur_noise);  %G(u,v)
restore_with_noise=ifft2(fourier_motion_blur_noise.*H_w);  %輸出頻域=G(u,v)H_w(u,v)，時域為頻域傅里葉逆變換
figure('name','公式法維也納濾波');
subplot(1,2,1);imshow(motion_blur_noise,[]);title('運動模糊添加噪聲')
subplot(1,2,2);imshow(restore_with_noise,[]);title('維也納濾波')

• 運行效果：