一、实验目的
利用逆滤波和维纳滤波，对Lena加噪运动模糊降质图像进行复原，比较不同参数选择对复原结果的影响。
二、实验内容
1) 输入Lena图像，对图像进行运动降质；降质模型：
2) 对图像叠加高斯白噪声；
3) 寻找最佳逆滤波半径r；
4) 逆滤波；
5) IFFT，展示结果；
6) 再寻找最佳维纳滤波K值；
7) 维纳滤波；
8) IFFT，展示结果。
三、实验过程和数学原理
1、逆滤波
在无噪情况下，逆滤波是完美的：
然而实际情况都有噪声：
如果H(u,v)存在零点，那么在H(u,v)零点附近进行复原，会导致第二项变得很大很大，复原效果很差。
理论上说，我们应该找出H(u,v)的所有零点，然后规避这些零点进行逆滤波。然而，降质模型零点十分分散（实际上是一条条斜线，实验中会看到），并且要作无数个圆区域，编程非常麻烦。
实际编程中，我们采用如下思路：既然fft后的频谱中，信号频谱主要集中在低频分量，那么，我们用fftshift，将频谱移到中心；以频谱中心为圆心，规定一个圆区域，在圆内正常逆滤波；而圆外是大量的较小的噪声分量，不给任何机会，直接赋值为0。这样，我们就不需要考虑H(u,v)的零点影响了。
此外，再通过PSNR最大准则，寻找最佳滤波半径R。实验结果如下：
实验代码：
function Imagerestoration %~~~~~~~~~~~~~~~~~~~~~~~~~~实验说明~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~% % 1、由于T=5,a=b=1效果太差太差，几乎无法看到复原现象，因此本实验采用T=1,a=b=0.02降质模型。 % 2、作为实验讲义的补充，本实验加入均值为0、方差为1e-3的AGWN模型。 % 若方差过大，逆滤波效果也不理想。 % 3、一定要对fftshift后的频谱进行运动模糊处理！！！否则现象都是错误的！！！ % 2018-5-19 by XING %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~% clear;close all;clc; fprintf('-------------------------------逆滤波实验-------------------------------\n'); I = im2double(imread('lena512.bmp'));% [0,1] [M,~] = size(I);% square % Display the original image. figure; subplot(1,3,1), imshow(I); title('\fontsize{20}原始图像'); %% Simulate a Motion Blur：H(u,v) T=1;a=0.02;b=0.02; v=[-M/2:M/2-1];u=v'; A=repmat(a.*u,1,M)+repmat(b.*v,M,1); H=T/pi./A.*sin(pi.*A).*exp(-1i*pi.*A); H(A==0)=T;% replace NAN %% Get the blurred Image % Warning: fftshift should be written F=fftshift(fft2(I)); FBlurred=F.*H; % Display the blurred image IBlurred =real(ifft2(ifftshift(FBlurred))); subplot(1,3,2), imshow(uint8(255.*mat2gray(IBlurred))); title('\fontsize{20}运动模糊图像'); %% Deblur perfectly without Noise FDeblurred=FBlurred./H; IDeblurred=real(ifft2(ifftshift(FDeblurred))); subplot(1,3,3), imshow(uint8(255.*mat2gray(IDeblurred))); title('\fontsize{20}无噪情况下直接逆滤波'); %% Simulate Noise Model noise_mean = 0; noise_var = 1e-3; noise=imnoise(zeros(M),'gaussian', noise_mean,noise_var); FNoise=fftshift(fft2(noise)); %% Get the Blurred_Noised Image FBlurred_Noised=FNoise+FBlurred; % Display the blurred_noised image IBlurred_Noised=real(ifft2(ifftshift(FBlurred_Noised))); figure; subplot(1,3,1), imshow(uint8(255.*mat2gray(IBlurred_Noised))); title('\fontsize{20}加噪运动模糊图像'); %% Deblur when Ignoring Noise FDeblurred2=FBlurred_Noised./H; FH1=abs(FDeblurred2); IDeblurred2=real(ifft2(ifftshift(FDeblurred2))); subplot(1,3,2), imshow(uint8(255.*mat2gray(IDeblurred2))); title ('\fontsize{20}有噪情况下直接逆滤波'); %% Find out the best Radius maxPSNR=0; bestRadius=0; tic; for Radius=33:1e-2:34 % 预实验bestr约为33.8左右 FDeblurred2=zeros(M); for a=1:M for b=1:M if sqrt((a-M/2).^2+(b-M/2).^2)<Radius FDeblurred2(a,b)=FBlurred_Noised(a,b)./H(a,b); end end end % Calculate PSNR and compare with the best IDeblurred2=real(ifft2(ifftshift(FDeblurred2))); PSNR=PSNRcal(IDeblurred2,I); if PSNR>maxPSNR maxPSNR=PSNR; bestRadius=Radius; end end fprintf(' 最佳滤波半径: %.1f\n', bestRadius); fprintf(' 最大PSNR： %d dB\n', round(maxPSNR)); fprintf(' 寻找最佳半径耗时: %.1f s\n', toc); % Displace the best Restoration FDeblurred2=zeros(M); for a=1:M for b=1:M if sqrt((a-M/2).^2+(b-M/2).^2)<bestRadius FDeblurred2(a,b)= FBlurred_Noised(a,b)./H(a,b); end end end FH2=abs(FDeblurred2); IDeblurred2=real(ifft2(ifftshift(FDeblurred2))); subplot(1,3,3), imshow(uint8(255.*mat2gray(IDeblurred2))); title(['\fontsize{20}最佳半径为 ', num2str(bestRadius),'的圆内逆滤波']); fprintf('\n 半径逆滤波原理: \n'); fprintf(' 1：通过fft2，使有效信号集中在低频区域，噪声集中在高频区域；\n'); fprintf(' 2：通过fftshift，把频谱移到矩阵中心；\n'); fprintf(' 3：以矩阵中心为圆心，规定一个半径为r的圆域；圆外直接赋0，抑制噪声；圆内正常逆滤波。\n'); fprintf(' 这比寻找、避开H(u,v)的零点滤波，编程要简单很多，效果也不错！\n'); fprintf('\n 实验结果说明: \n'); fprintf(' 1：对于无噪运动模糊图像，逆滤波复原几乎完美 ^ ^\n'); fprintf(' 2：对于有噪运动模糊图像，直接逆滤波是灾难，取半径滤波效果尚可。\n'); fprintf('\n 同时说明: PSNR是多么不靠谱！\n'); figure; subplot(1,2,1),imshow(im2double(uint8(FH1))); title ('\fontsize{20}有噪情况下直接逆滤波得到的图像频谱'); subplot(1,2,2),imshow(im2double(uint8(FH2))); title ('\fontsize{20}圆内逆滤波得到的图像频谱'); fprintf('\n~~~~~~~程序已暂停；按任意键进行维纳滤波实验~~~~~~~~\n'); pause;
2、维纳滤波
理想维纳滤波器为：
实际应用中，NSR难以被精确计算。因此，我们常常设为K，并寻找最佳K值。
实验结果：
实验代码：
%% Deblur Image Using Wiener Filter fprintf('\n-------------------------------维纳滤波实验-------------------------------\n'); fprintf(' 根据大量实验，我发现PSNR无法作为寻找最佳K值的标准：\n'); fprintf(' 1：参数K将会停留在K=0处，即逆滤波；\n'); fprintf(' 2：最大PSNR达到70dB以上，实际图像质量极差！\n'); fprintf('\n 因此，我通过观察，选择了最佳K=0.05~\n'); fprintf(' ps. 程序中保留了根据PSNR寻找K值的代码，有兴趣可以尝试 ^ ^\n'); % Display the blurred_noised image again figure(); subplot(1,3,1); imshow(uint8(255.*mat2gray(IBlurred_Noised))); title('\fontsize{20}加噪运动模糊图像'); % Deblur with theoretic NSR buf=(abs(H)).^2; % Notice '.' !!!!!!!! NSR=FNoise./F; FDeblurred3=FBlurred_Noised./H.*buf./(buf+NSR); IDeblurred3=real(ifft2(ifftshift(FDeblurred3))); subplot(1,3,2), imshow(uint8(255.*mat2gray(IDeblurred3))); title('\fontsize{20}K=NSR的理论维纳滤波'); % Find out the best K % tic; % maxPSNR=0; % beskK=0; % for K=0:1e-2:1 % FDeblurred2=zeros(M); % FDeblurred3=FBlurred_Noised./H.*buf./(buf+bestK); % IDeblurred3=real(ifft2(ifftshift(FDeblurred3))); % % % Calculate PSNR and compare with the best % PSNR=PSNRcal(IDeblurred3,I); % if PSNR>maxPSNR % maxPSNR=PSNR; % bestK=K; % end % end % % fprintf(' 最佳K值: %.2f\n', bestK); % fprintf(' 最大PSNR: %d dB\n', round(maxPSNR)); % fprintf(' 寻找最佳K值耗时: %.1f s\n', toc); % Deblur with best K bestK=0.05; FDeblurred3=FBlurred_Noised./H.*buf./(buf+bestK); IDeblurred3=real(ifft2(ifftshift(FDeblurred3))); % Display the best restored Image subplot(1,3,3), imshow(uint8(255.*mat2gray(IDeblurred3))); title(['\fontsize{20}实际最佳K= ', num2str(bestK),'的维纳滤波']); fprintf(' Written by XING\n'); fprintf(' 2018-5-19 Beijing\n'); end
实验中用到的PSNR计算子程序：
function PSNR=PSNRcal(I,I2) h=512;w=512; B=8;% 编码一个像素用8个二进制位 MAX=2^B-1;% 图像有多少灰度级 MES=sum(sum((I-I2).^2))/(h*w);% 均方差 PSNR=20*log10(MAX/sqrt(MES));% 峰值信噪比 end