动态模糊图像复原MATLAB程序

1、退化程序

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');%保存图

2、运动角度的求解

%求解模糊运动角度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);

%模糊图像

B1=fft2(double(B));

B2=mat2gray(log(abs(B1)));

subplot(122);imshow(B2);

%模糊图像的频谱图

C=sum(B2,1);

%对频谱图求列和

[m,n]=size(C); x=0:1:n-1; y=C;

figure,plot(x,y);title('频谱列和曲线图1') %绘制频谱列和曲线图

3、运动长度算法

%求解模糊运动长度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('模糊运动长度');

4、噪声分析

%噪声分析

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');

5、去噪还原

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');

6、截取部分图片进行对比程序

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截取部分');

7、模糊系数、质量指数、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 i=2:W-1

t=doubleimsrc(j-1,i+1)+doubleimsrc(j+1,i-1)-doubleimsrc(j-1,i-1)-doub leimsrc(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)-doub leimdst(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

varsrc=varsrc+(doubleimsrc(j+n,i+m)-midsrc)*(doubleimsrc(j+n,i+m)-mid src);

vardst=vardst+(doubleimdst(j+n,i+m)-middst)*(doubleimdst(j+n,i+m)-mid dst);

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