一、实验目的
利用逆滤波和维纳滤波,对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(-1ipi.*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);
<span class="k" style="color:rgb(0,112,32);font-weight:bold;">for</span> <span class="n">a</span><span class="p">=</span><span class="mi" style="color:rgb(64,160,112);">1</span><span class="p">:</span><span class="n">M</span>
<span class="k" style="color:rgb(0,112,32);font-weight:bold;">for</span> <span class="n">b</span><span class="p">=</span><span class="mi" style="color:rgb(64,160,112);">1</span><span class="p">:</span><span class="n">M</span>
<span class="k" style="color:rgb(0,112,32);font-weight:bold;">if</span> <span class="nb" style="color:rgb(0,112,32);">sqrt</span><span class="p">((</span><span class="n">a</span><span class="o" style="color:rgb(102,102,102);">-</span><span class="n">M</span><span class="o" style="color:rgb(102,102,102);">/</span><span class="mi" style="color:rgb(64,160,112);">2</span><span class="p">)</span><span class="o" style="color:rgb(102,102,102);">.^</span><span class="mi" style="color:rgb(64,160,112);">2</span><span class="o" style="color:rgb(102,102,102);">+</span><span class="p">(</span><span class="n">b</span><span class="o" style="color:rgb(102,102,102);">-</span><span class="n">M</span><span class="o" style="color:rgb(102,102,102);">/</span><span class="mi" style="color:rgb(64,160,112);">2</span><span class="p">)</span><span class="o" style="color:rgb(102,102,102);">.^</span><span class="mi" style="color:rgb(64,160,112);">2</span><span class="p">)</span><span class="o" style="color:rgb(102,102,102);"><</span><span class="n">Radius</span>
<span class="n">FDeblurred2</span><span class="p">(</span><span class="n">a</span><span class="p">,</span><span class="n">b</span><span class="p">)=</span><span class="n">FBlurred_Noised</span><span class="p">(</span><span class="n">a</span><span class="p">,</span><span class="n">b</span><span class="p">)</span><span class="o" style="color:rgb(102,102,102);">./</span><span class="n">H</span><span class="p">(</span><span class="n">a</span><span class="p">,</span><span class="n">b</span><span class="p">);</span>
<span class="k" style="color:rgb(0,112,32);font-weight:bold;">end</span>
<span class="k" style="color:rgb(0,112,32);font-weight:bold;">end</span>
<span class="k" style="color:rgb(0,112,32);font-weight:bold;">end</span>
<span class="c" style="color:rgb(96,160,176);font-style:italic;">% Calculate PSNR and compare with the best</span>
<span class="n">IDeblurred2</span><span class="p">=</span><span class="nb" style="color:rgb(0,112,32);">real</span><span class="p">(</span><span class="n">ifft2</span><span class="p">(</span><span class="n">ifftshift</span><span class="p">(</span><span class="n">FDeblurred2</span><span class="p">)));</span>
<span class="n">PSNR</span><span class="p">=</span><span class="n">PSNRcal</span><span class="p">(</span><span class="n">IDeblurred2</span><span class="p">,</span><span class="n">I</span><span class="p">);</span>
<span class="k" style="color:rgb(0,112,32);font-weight:bold;">if</span> <span class="n">PSNR</span><span class="o" style="color:rgb(102,102,102);">></span><span class="n">maxPSNR</span>
<span class="n">maxPSNR</span><span class="p">=</span><span class="n">PSNR</span><span class="p">;</span>
<span class="n">bestRadius</span><span class="p">=</span><span class="n">Radius</span><span class="p">;</span>
<span class="k" style="color:rgb(0,112,32);font-weight:bold;">end</span>
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值。
实验结果:
![](https://i-blog.csdnimg.cn/blog_migrate/71e6b6847b75deccd91123aefb2aec37.png)
%% 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))/(hw);% 均方差
PSNR=20log10(MAX/sqrt(MES));% 峰值信噪比
end