逆滤波和维纳滤波

最新推荐文章于 2024-03-04 15:28:15 发布

梅津太郎

最新推荐文章于 2024-03-04 15:28:15 发布

阅读量7.3k

点赞数 8

分类专栏：图像处理

本文链接：https://blog.csdn.net/gaocui883/article/details/106643316

版权

图像处理专栏收录该内容

88 篇文章 1 订阅

订阅专栏

一、实验目的

利用逆滤波和维纳滤波，对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);">&lt;</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);">&gt;</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值。

实验结果：

实验代码：

%% 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

梅津太郎

关注

8
点赞
踩
28

收藏

觉得还不错? 一键收藏
0
评论
逆滤波和维纳滤波

一、实验目的        利用逆滤波和维纳滤波，对Lena加噪运动模糊降质图像进行复原，比较不同参数选择对复原结果的影响。二、实验内容1)     输入Lena图像，对图像进行运动降质；降质模型：2) ...
复制链接

扫一扫

专栏目录