[数字图像处理]图像复原--逆滤波
1.逆滤波的问题点
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
2.两个退化函数的模型
2.1 大气湍流模型
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
2.2 运动模糊模型
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
3.图像的逆滤波
3.1 实验步骤与实验用图像
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
3.2 直接逆滤波
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
3.3 维纳滤波器
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
![](/assets/blank.gif)
3.4 约束最小二乘方滤波
![](/assets/blank.gif)
![](/assets/blank.gif)
4 实验代码
close all;
clear all;
clc;
%% ----------init-----------------------------
f = imread('./original_DIP.tif');
f = mat2gray(f,[0 255]);
f_original = f;
[M,N] = size(f);
P = 2*M;
Q = 2*N;
fc = zeros(M,N);
for x = 1:1:M
for y = 1:1:N
fc(x,y) = f(x,y) * (-1)^(x+y);
end
end
F_I = fft2(fc,P,Q);
figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Original Image');
subplot(1,2,2);
imshow(log(1 + abs(F_I)),[ ]);
xlabel('b).Fourier spectrum of a).');
%% ------motion blur------------------
H = zeros(P,Q);
a = 0.02;
b = 0.02;
T = 1;
for x = (-P/2):1:(P/2)-1
for y = (-Q/2):1:(Q/2)-1
R = (x*a + y*b)*pi;
if(R == 0)
H(x+(P/2)+1,y+(Q/2)+1) = T;
else H(x+(P/2)+1,y+(Q/2)+1) = (T/R)*(sin(R))*exp(-1i*R);
end
end
end
%% ------the atmospheric turbulence modle------------------
H_1 = zeros(P,Q);
k = 0.0025;
for x = (-P/2):1:(P/2)-1
for y = (-Q/2):1:(Q/2)-1
D = (x^2 + y^2)^(5/6);
D_0 = 60;
H_1(x+(P/2)+1,y+(Q/2)+1) = exp(-k*D);
end
end
%% -----------noise------------------
a = 0;
b = 0.2;
n_gaussian = a + b .* randn(M,N);
Noise = fft2(n_gaussian,P,Q);
figure();
subplot(1,2,1);
imshow(n_gaussian,[-1 1]);
xlabel('a).Gaussian noise');
subplot(1,2,2);
imshow(log(1 + abs(Noise)),[ ]);
xlabel('b).Fourier spectrum of a).');
%%
G = H .* F_I + Noise;
% G = H_1 .* F_I + Noise;
gc = ifft2(G);
gc = gc(1:1:M+27,1:1:N+27);
for x = 1:1:(M+27)
for y = 1:1:(N+27)
g(x,y) = gc(x,y) .* (-1)^(x+y);
end
end
gc = gc(1:1:M,1:1:N);
for x = 1:1:(M)
for y = 1:1:(N)
g(x,y) = gc(x,y) .* (-1)^(x+y);
end
end
figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Original Image');
subplot(1,2,2);
imshow(log(1 + abs(F_I)),[ ]);
xlabel('b).Fourier spectrum of a).');
figure();
subplot(1,2,1);
imshow(abs(H),[ ]);
xlabel('c).The motion modle H(u,v)(a=0.02,b=0.02,T=1)');
subplot(1,2,2);
n = 1:1:P;
plot(n,abs(H(400,:)));
axis([0 P 0 1]);grid;
xlabel('H(n,400)');
ylabel('|H(u,v)|');
figure();
subplot(1,2,1);
imshow(real(g),[0 1]);
xlabel('d).Result image');
subplot(1,2,2);
imshow(log(1 + abs(G)),[ ]);
xlabel('e).Fourier spectrum of d). ');
%% --------------inverse_filtering---------------------
%F = G ./ H;
%F = G ./ H_1;
for x = (-P/2):1:(P/2)-1
for y = (-Q/2):1:(Q/2)-1
D = (x^2 + y^2)^(0.5);
if(D < 258)
F(x+(P/2)+1,y+(Q/2)+1) = G(x+(P/2)+1,y+(Q/2)+1) ./ H_1(x+(P/2)+1,y+(Q/2)+1);
% no noise D < 188
% noise D < 56
else F(x+(P/2)+1,y+(Q/2)+1) = G(x+(P/2)+1,y+(Q/2)+1);
end
end
end
% Butterworth_Lowpass_Filters
H_B = zeros(P,Q);
D_0 = 70;
for x = (-P/2):1:(P/2)-1
for y = (-Q/2):1:(Q/2)-1
D = (x^2 + y^2)^(0.5);
%if(D < 200) H_B(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D/D_0)^100);end
H_B(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D/D_0)^20);
end
end
F = F .* H_B;
f = real(ifft2(F));
f = f(1:1:M,1:1:N);
for x = 1:1:(M)
for y = 1:1:(N)
f(x,y) = f(x,y) * (-1)^(x+y);
end
end
%% ------show Result------------------
figure();
subplot(1,2,1);
imshow(f,[0 1]);
xlabel('a).Result image');
subplot(1,2,2);
imshow(log(1 + abs(F)),[ ]);
xlabel('b).Fourier spectrum of a).');
figure();
n = 1:1:P;
plot(n,abs(F(400,:)),'r-',n,abs(F(400,:)),'b-');
axis([0 P 0 1000]);grid;
xlabel('Number of rows(400th column)');
ylabel('Fourier amplitude spectrum');
legend('F_{limit}(u,v)','F(u,v)');
figure();
n = 1:1:P;
plot(n,abs(H(400,:)),'g-');
axis([0 P 0 1]);grid;
xlabel('H''_{s}(n,400)');
ylabel('|H''_{s}(u,v)|');
%% ----------Wiener filters-----------
% K = 0.000014;
K = 0.02;
%H_Wiener = ((abs(H_1).^2)./((abs(H_1).^2)+K)).*(1./H_1);
H_Wiener = ((abs(H).^2)./((abs(H).^2)+K)).*(1./H);
F_Wiener = H_Wiener .* G;
f_Wiener = real(ifft2(F_Wiener));
f_Wiener = f_Wiener(1:1:M,1:1:N);
for x = 1:1:(M)
for y = 1:1:(N)
f_Wiener(x,y) = f_Wiener(x,y) * (-1)^(x+y);
end
end
[SSIM_Wiener mssim] = ssim_index(f_Wiener,f_original,[0.01 0.03],ones(8),1);
SSIM_Wiener
%% ------show Result------------------
figure();
subplot(1,2,1);
%imshow(f_Wiener(1:128,1:128),[0 1]);
imshow(f_Wiener,[0 1]);
xlabel('d).Result image by Wiener filter');
subplot(1,2,2);
imshow(log(1+abs(F_Wiener)),[ ]);
xlabel('c).Fourier spectrum of c).');
% subplot(1,2,2);
% %imshow(f(1:128,1:128),[0 1]);
% imshow(f,[0 1]);
% xlabel('e).Result image by inverse filter');
figure();
n = 1:1:P;
plot(n,abs(F(400,:)),'r-',n,abs(F_Wiener(400,:)),'b-');
axis([0 P 0 500]);grid;
xlabel('Number of rows(400th column)');
ylabel('Fourier amplitude spectrum');
legend('F(u,v)','F_{Wiener}(u,v)');
figure();
subplot(1,2,1);
imshow(log(1 + abs(H_Wiener)),[ ]);
xlabel('a).F_{Wiener}(u,v).');
subplot(1,2,2);
n = 1:1:P;
plot(n,abs(H_Wiener(400,:)));
axis([0 P 0 80]);grid;
xlabel('Number of rows(400th column)');
ylabel('Amplitude spectrum');
%% ------------Constrained_least_squares_filtering---------
p_laplacian = zeros(M,N);
Laplacian = [ 0 -1 0;
-1 4 -1;
0 -1 0];
p_laplacian(1:3,1:3) = Laplacian;
P = 2*M;
Q = 2*N;
for x = 1:1:M
for y = 1:1:N
p_laplacian(x,y) = p_laplacian(x,y) * (-1)^(x+y);
end
end
P_laplacian = fft2(p_laplacian,P,Q);
F_C = zeros(P,Q);
r = 0.2;
H_clsf = ((H')./((abs(H).^2)+r.*P_laplacian));
F_C = H_clsf .* G;
f_c = real(ifft2(F_C));
f_c = f_c(1:1:M,1:1:N);
for x = 1:1:(M)
for y = 1:1:(N)
f_c(x,y) = f_c(x,y) * (-1)^(x+y);
end
end
%%
figure();
subplot(1,2,1);
imshow(f_c,[0 1]);
xlabel('e).Result image by constrained least squares filter (r = 0.2)');
subplot(1,2,2);
imshow(log(1 + abs(F_C)),[ ]);
xlabel('f).Fourier spectrum of c).');
[SSIM_CLSF mssim] = ssim_index(f_c,f_original,[0.01 0.03],ones(8),1);
figure();
subplot(1,2,1);
imshow(log(1 + abs(H_clsf)),[ ]);
xlabel('a).F_{clsf}(u,v).');
subplot(1,2,2);
n = 1:1:P;
plot(n,abs(H_clsf(400,:)));
axis([0 P 0 80]);grid;
xlabel('Number of rows(400th column)');
ylabel('Amplitude spectrum');
原文发于博客:http://blog.csdn.net/thnh169/
=============更新日志===================
NULL
[数字图像处理]图像复原--逆滤波相关推荐
- 数字图像处理--傅里叶(逆)变换
数字图像处理–傅里叶(逆)变换 主要内容 (1)对一副图像进行缩放,显示原始图像和缩放后的图像,分别对其进行傅里叶变换,显示变换后结果: (2)对一副图像进行旋转,显示原始图像和旋转后的图像,分别对其 ...
- 领域平均法matlab代码实验,数字图像处理邻域平均法滤波实验报告matlab实现.doc...
数字图像处理邻域平均法滤波实验报告matlab实现 数字图像处理 实验报告 实验三 邻域平均法滤波 学号 姓名 实验三 邻域平均法滤波 一.实验内容 选取噪声较明显的图像,分别采用3*3.5*5.7* ...
- 数字图像处理之平滑滤波
数字图像处理之平滑滤波 by方阳 版权声明:本文为博主原创文章,转载请标出转载地址 http://www.cnblog ...
- 数字图像处理实验(五)|图像复原{逆滤波和伪逆滤波、维纳滤波deconvwnr、大气湍流扰动模型、运动模糊处理fspecial}(附matlab实验代码和截图)
文章目录 一.实验目的 二.实验仪器 三.实验原理 四.实验内容 1.逆滤波:选择MATLAB文件夹中的foggy图像作为实验图像. (1)生成退化函数: (2)复原 (a)直接逆滤波 (b)修正函数 ...
- matlab 数字图像滤波,数字图像处理 (基于Matlab) 滤波
<数字图像处理> 实验报告 一.实验目的(不少于200字) 一.第一个实验用的是各种空间域的方式来滤波,也就是直接把图像和空间滤波器的模板做卷积,当 然图像处理很重要的一个部分还有频域的处 ...
- 数字图像处理-频率域滤波原理
from:https://blog.csdn.net/forrest02/article/details/55510711?locationNum=15&fps=1 写在前面的话 作者是一名在 ...
- 数字图像处理之频域滤波
前段时间看了很多的概念和知识,发现因为是走马观花的过了一遍,所以看得稀里糊涂的,然后许多地方混淆了概念,特别是关于图像频率域的部分的理解(包括图像频率域滤波之类的),所以下面总结一下这段时间重新看&l ...
- 系统学习数字图像处理之频域滤波
最近在看模板匹配,虽然很简单,但还是想认真过下基础,因此把信号处理频域相关的内容,接着图像处理再过一遍. 理论上,对连续变量t的连续函数f(t)的傅里叶变换为F(u),利用f(t)取样后的函数重建f( ...
- 数字图像处理——中值滤波及其改进算法
一.算法介绍 中值滤波器是非线性滤波器的一个例子,它在保留图像特征方面非常有效. 但是,滤波器的窗口大小直接影响中值滤波器的性能. 较小的窗口保留了特征,但会导致噪声抑制的减少. 在较大窗口的情况下, ...
最新文章
- 2016cocoapods安装流程及使用
- MongoDb Windows linux平台环境及主流编程语言驱动安装同时配置mongoDb的远程连接
- C++编程优化——让你的代码飞起来
- php大数组循环嵌套的性能优化
- Keras【Deep Learning With Python】实现多元线性回归
- hbm2java和hbm2ddl的使用步骤
- WIn10+Anaconda 环境下安装 PyTorch 避坑指南
- 鲲鹏云HCIA知识总结(一)
- scikit-learn 朴素贝叶斯类库使用小结
- 怎样查看.a和so文件中的接口
- Java 使用 Timer 进行调度
- 使用Xamarin开发手机聊天程序 -- 基础篇(大量图文讲解 step by step,附源码下载)
- JAVA之运算符优先级
- mysql5.1怎么备份,MySQL 5.1升级到MySQL 5.5的步骤
- 关于解决Codeblocks中文乱码问题
- scrapy 序列化写入器 ——ItemExporter
- python零基础8分钟基础入门
- C++_类和对象_C++运算符重载_左移运算符重载_链式编程_实现直接打印对象---C++语言工作笔记056
- html5 xml文本编辑,简介XML文档的阅读与编辑
- ascii码与键盘代码的区别
热门文章
- vivado生成bit流错误:Combinatorial Loop Alert
- Android-简单单词书app
- mong命令学习记录
- 微软开发的服务器简称,AAD Connect 微软官方的描述准确吗?
- linux 显卡 压力测试软件,显卡压力测试工具 GpuTest
- C语言练习题 :猴子吃桃程序
- GradCAM神经网络可视化解释(原理和实现)
- 动物识别系统代码python_动物识别系统代码
- 测量耐力也有算法了!仅需锻炼20分钟,就能知晓自己能跑多久
- 马斯克告诉推特员工:要么继续高强度工作,要么拿遣散费走人;微信新增删除声音锁功能;Deno 1.28 发布|极客头条