动态模糊图像复原MATLAB程序

1、 退化程序clc;clear all;close all;I=imread(F:mmwB1图2.jpg);%读图figure;subplot(4,3,1);imshow(I);title(原图像);LEN=30;%运动长度30THETA=30;%运动角度30% LEN=60;% THETA=60;n=2;for i=1:3 for j=1:3PSF=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 onn=n+1; endend%imwrite(Blurred,模糊60长60.png);%保存图2、 运动角度的求解%求解模糊运动角度matlab代码close all;clc;clear all;im=imread(F:mmwB1模糊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:mmwB1模糊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 allim=imread(F:mmwB1图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:190for j=1:162for k=1:3f11(i,j,k)=im(i,j,k);endendendfor i=1:130for j=501:630fork=1:3;f22(i,j-500,k)=im(i,j,k);endendfor i=721:870for j=11:170for k=1:3f33(i-720,j-10,k)=im(i,j,k);endendendfor i=761:830for j=561:630for k=1:3f33(i-760,j-560,k)=im(i,j,k);endendendfigure;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:mmwB1图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:mmwB1模糊60长60.png); %未处理质量较差图像b=a(64:120,67:126);a=imread(F:mmwmin_recover1.png);%算法处理后质量较好图象c=a(64:120,67:126);%从eyechart3中截取测试参考图象，截取部分需要进行缩放-%使之与eyechart1,eyechart2截取部分大小匹配-a=imread(F:mmwB1图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:2imsrcnamehead=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=1PSNR=0.0; %PSNR赋初值for j=1:Hfor i=1:WPSNR=PSNR+double(doubleimsrc(j,i)-doubleimdst(j,i)*(doubleimsrc(j,i)-doubleimdst(j,i);endendPSNR=PSNR/W/H;PSNR=10*log10(255*255/PSNR);%-PSNR计算完毕-end%-模糊系数KBlur计算-if KBlurenable=1Sin=0.0; %Sin赋初值Sout=0.0;for j=2:H-1for i=2:W-1t=doubleimsrc(j-1,i+1)+doubleimsrc(j+1,i-1)-doubleimsrc(j-1,i-1)-doubleimsrc(j+1,i+1);if t0t=-t;endSin=Sin+t; %源图象邻域边缘能量计算t=doubleimdst(j-1,i+1)+doubleimdst(j+1,i-1)-doubleimdst(j-1,i-1)-doubleimdst(j+1,i+1);if t0t=-t;endendendSout=Sout+t; %污染图象邻域边缘能量计算 KBlur=Sout/Sin;end%-KBlur计算完毕-%-质量指数Q计算-if Qenable=1Q=0.0; %Q赋初值Qnum=0; %图象以7X7块大小计算每块的Q，逐象素的移动块窗口，这里Qnum为块数量的计数for j=4:H-3for i=4:W-3midsrc=0.0;middst=0.0;varsrc=0.0;vardst=0.0; %源图象和污染图象块内的平均值和方差赋初值varsrcdst=0.0;%源图象和污染图象块内的协方差赋初值for n=-3:3for m=-3:3 midsrc=midsrc+doubleimsrc(j+n,i+m);middst=middst+doubleimdst(j+n,i+m);endendmidsrc=midsrc/49;middst=middst/49;%源图象和污染图象块内的平均值计算-for n=-3:3for m=-3:3varsrc=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);endendvarsrc=va