基于MATLAB的图像模糊复原系统
基于MATLAB的图像模糊复原系统
源码如下：
- 退化程序
clc;
clear all;
close all;
I=imread('F:\mmw\B1\图2.jpg');%读图
figure;
subplot(4,3,1);
imshow(I);
title('原图像');
LEN=30;%运动长度30
THETA=30;%运动角度30
% LEN=60;
% THETA=60;
n=2;
for i=1:3
for j=1:3
PSF=fspecial('motion',LEN*i,THETA*j);%退化并研究运动角度和长度对图片模糊程度的影响
PSF=fspecial('motion',LEN,THETA);
Blurred=imfilter(I,PSF,'circular','conv');
subplot(4,3,n);
imshow(uint8(Blurred));
title('模糊化');
hold on
n=n+1;
end
end
%imwrite(Blurred,'模糊∠60长60.png');%保存图
- 运动角度的求解
%求解模糊运动角度matlab代码
close all;clc;clear all;
im=imread('F:\mmw\B1\模糊∠60长60.png');
img_gray=rgb2gray(im);%灰度化
img_fft=fftshift(fft2(img_gray));
N=abs(img_fft);
P=(N-min(min(N)))/(max(max(N))-min(min(N)))*225;
figure;
imshow(P);
title('频谱图（运动角度与光斑方向垂直）');
len=35;
theta=0;
PSF=fspecial('motion',len,theta);
B=imfilter(img_gray,PSF,'circular','conv');
subplot(121);imshow(B);
%模糊图像的频谱图
C=sum(B2,1);
%对频谱图求列和
[m,n]=size(C); x=0:1:n-1; y=C;
figure,plot(x,y);title('频谱列和曲线图1')
%绘制频谱列和曲线图
- 运动长度算法
%求解模糊运动长度matlab代码：
im=imread('F:\mmw\B1\模糊∠60长60.png');
img_gray=rgb2gray(im);%灰度化
h=fspecial('sobel');%sobel边缘检测
img_double=double(img_gray);
J=conv2(img_double,h,'same');
IP=abs(fft2(J));
S=fftshift(real(ifft2(IP)));
figure;plot(S);
title('模糊运动长度');
- 噪声分析
%噪声分析
clc;
clear all
im=imread('F:\mmw\B1\图1.png');
[m,n,h]=size(im);
f11=ones(192,162,3);
f22=ones(130,130,3);
f33=ones(100,100,3);
f44=ones(70,70,3);
for i=1:190
for j=1:162
for k=1:3
f11(i,j,k)=im(i,j,k);
end
end
end
for i=1:130
for j=501:630
fork=1:3;
f22(i,j-500,k)=im(i,j,k);
end
end
for i=721:870
for j=11:170
for k=1:3
f33(i-720,j-10,k)=im(i,j,k);
end
end
end
for i=761:830
for j=561:630
for k=1:3
f33(i-760,j-560,k)=im(i,j,k);
end
end
end
figure;
subplot(221),hist(f11,100);
subplot(222),hist(f22,100);
subplot(223),hist(f33,100);
subplot(224),hist(f44,100);
title('噪声分析2');
- 去噪还原
clc;
clear all;
close all;
I=imread('F:\mmw\B1\图1.png');%读图
Len=60;
Theta=60;
PSF=fspecial('motion',Len,Theta); %模糊化
BlurredA=imfilter(I,PSF,'circular','conv');
wnr1=deconvwnr(BlurredA,PSF);%维纳滤波
BlurredD=imfilter(I,PSF,'circ','conv');
INITPSF=ones(size(PSF));
[K DePSF]=deconvblind(BlurredD,INITPSF,30);%盲去卷积法
BlurredB=imfilter(I,PSF,'conv');
v=0.02;
Blurred_I_Noisy=imnoise(BlurredB,'gaussian',0,v);
NP=v*prod(size(I));
J=deconvreg(Blurred_I_Noisy,PSF,NP);%最小二乘法
BlurredC=imfilter(I,PSF,'symmetric','conv');
v=0.002;
BlurredNoisy=imnoise(BlurredC,'gaussian',0,v);
Luc=deconvlucy(BlurredNoisy,PSF,5);%L_Rl滤波
subplot(221);imshow(I);title('原图');
subplot(222);imshow(BlurredA);title('模糊化');
%subplot(233);imshow(wnr1);title('维纳滤波');
subplot(223);imshow(J);title('最小二乘法');imwrite(J,'min_recover1.png');
subplot(224);imshow(Luc);title('L_R法');imwrite(Luc,'LR_recover1.png');
- 截取部分图片进行对比程序
clear all;clc;
a=imread('F:\mmw\B1\模糊∠60长60.png');
%未处理质量较差图像
b=a([64:120],[67:126]);
a=imread('F:\mmw\min_recover1.png');
%算法处理后质量较好图象
c=a([64:120],[67:126]);
%%从eyechart3中截取测试参考图象，截取部分需要进行缩放------------------
%%使之与eyechart1,eyechart2截取部分大小匹配-----------------------
a=imread('F:\mmw\B1\图2.jpg');
%高清晰参考图象
d=a([64:120],[67:126]);
e=imresize(d,[length(b(:,1)),length(b(1,:))],'bicubic');
%调整
imwrite(b,'area_模糊∠60长60.png');
imwrite(c,'area_最小二乘法复原图.png');
imwrite(e,'area_图2.png');
subplot(1,3,1);
imshow(e);
title('模糊∠60长60截取参考');
hold on;
subplot(1,3,2);
imshow(b);
title('eyechart1截取部分');
hold on;
subplot(1,3,3);
imshow(c);
title('eyechart2截取部分');
- 模糊系数、质量指数、PSNR的计算（评价）
clc;
clear;
PSNRenable=1; %PSNR计算使能，为0不计算，为1，计算
KBlurenable=1; %模糊系数KBlur计算使能，为0不计算，为1，计算
Qenable=1; %质量指数Q计算使能，为0不计算，为1，计算
for m=1:2
imsrcnamehead='area_模糊∠60长60'; %源图象文件名头
imsrcnameext='png'; %源图象文件名扩展
if m==1 %以area_eyechart1.bmp为测试图象
imdstname=strcat('area_图2','.',imsrcnameext);%污染图象文件名，可修改
elseif m==2%以area_eyechart2.bmp为测试图象
imdstname=strcat('area_最小二乘法复原图','.',imsrcnameext);%污染图象文件名，可修改
end
%--------------------------------------------------------------------------
iminfo=imfinfo(strcat(imsrcnamehead,'.',imsrcnameext));%源图象信息读取
imsrc=imread(strcat(imsrcnamehead,'.',imsrcnameext)); %源图象读取
imdst=imread(imdstname,imsrcnameext); %污染图象读取
doubleimsrc=double(imsrc); %转换为浮点类型
doubleimdst=double(imdst); %转换为浮点类型
%----------------------------------------------------源图象和污染图象读取
W=iminfo.Width; %图象宽
H=iminfo.Height; %图象高
%----------------------------PSNR计算--------------------------------------
if PSNRenable==1
PSNR=0.0; %PSNR赋初值
for j=1:H
for i=1:W
PSNR=PSNR+double((doubleimsrc(j,i)-doubleimdst(j,i))*(doubleimsrc(j,i)-doubleimdst(j,i)));
end
end
PSNR=PSNR/W/H;
PSNR=10*log10(255*255/PSNR);
%---------------------------PSNR计算完毕-----------------------------------
end
%-------------------------模糊系数KBlur计算--------------------------------
if KBlurenable==1
Sin=0.0; %Sin赋初值
Sout=0.0;
for j=2:H-1
for i=2:W-1
t=doubleimsrc(j-1,i+1)+doubleimsrc(j+1,i-1)-doubleimsrc(j-1,i-1)-doubleimsrc(j+1,i+1);
if t<0
t=-t;
end
Sin=Sin+t; %源图象邻域边缘能量计算
t=doubleimdst(j-1,i+1)+doubleimdst(j+1,i-1)-doubleimdst(j-1,i-1)-doubleimdst(j+1,i+1);
if t<0
t=-t;
end
end
end
Sout=Sout+t; %污染图象邻域边缘能量计算
KBlur=Sout/Sin;
end
%-------------------------------KBlur计算完毕------------------------------
%-------------------------------质量指数Q计算------------------------------
if Qenable==1
Q=0.0; %Q赋初值
Qnum=0; %图象以7X7块大小计算每块的Q，逐象素的移动块窗口，这里Qnum为块数量的计数
for j=4:H-3
for i=4:W-3
midsrc=0.0;
middst=0.0;
varsrc=0.0;
vardst=0.0; %源图象和污染图象块内的平均值和方差赋初值
varsrcdst=0.0;%源图象和污染图象块内的协方差赋初值
for n=-3:3
for m=-3:3
midsrc=midsrc+doubleimsrc(j+n,i+m);
middst=middst+doubleimdst(j+n,i+m);
end
end
midsrc=midsrc/49;
middst=middst/49;
%%源图象和污染图象块内的平均值计算----------------------------------
for n=-3:3
for m=-3:3
varsrc=varsrc+(doubleimsrc(j+n,i+m)-midsrc)*(doubleimsrc(j+n,i+m)-midsrc);
vardst=vardst+(doubleimdst(j+n,i+m)-middst)*(doubleimdst(j+n,i+m)-middst);
varsrcdst=varsrcdst+(doubleimsrc(j+n,i+m)-midsrc)*(doubleimdst(j+n,i+m)-middst);
end
end
varsrc=varsrc/48;
vardst=vardst/48;
varsrcdst=varsrcdst/48;
if((varsrc+vardst)*(midsrc*midsrc+middst*middst))~=0%分母不为零的块才计算质量指数Q
Q=Q+4*varsrcdst*midsrc*middst/((varsrc+vardst)*(midsrc*midsrc+middst*middst));
%%源图象和污染图象块内Q计算完毕------------------------------------
Qnum=Qnum+1; %块计数加1
end
end
end
Q=Q/Qnum;
end
end
Q
PSNR
KBlur