动态模糊图像复原MATLAB程序

1、退化程序clc;clearall;closeall;I=imread('F:\mmw\B1\图2.jpg');%读图figure;subplot(4,3,1);imshow(I);title('原图像');LEN=30;%运动长度30THETA=30;%运动角度30%LEN=60;%THETA=60;n=2;fori=1:3forj=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('模糊化');holdonn=n+1;endend%imwrite(Blurred,'模糊∠60长60.png');%保存图2、运动角度的求解%求解模糊运动角度matlab代码closeall;clc;clearall;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;clearallim=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);fori=1:190forj=1:162fork=1:3f11(i,j,k)=im(i,j,k);endendendfori=1:130forj=501:630fork=1:3;f22(i,j-500,k)=im(i,j,k);endendfori=721:870forj=11:170fork=1:3f33(i-720,j-10,k)=im(i,j,k);endendendfori=761:830forj=561:630fork=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;clearall;closeall;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));[KDePSF]=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、截取部分图片进行对比程序clearall;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截取参考');holdon;subplot(1,3,2);imshow(b);title('eyechart1截取部分');holdon;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，计算form=1:2imsrcnamehead='area_模糊∠60长60';%源图象文件名头imsrcnameext='png';%源图象文件名扩展ifm==1%以area_eyechart1.bmp为测试图象imdstname=strcat('area_图2','.',imsrcnameext);%污染图象文件名，可修改elseifm==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计算--------------------------------------ifPSNRenable==1PSNR=0.0;%PSNR赋初值forj=1:Hfori=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计算--------------------------------ifKBlurenable==1Sin=0.0;%Sin赋初值Sout=0.0;forj=2:H-1fori=2:W-1t=doubleimsrc(j-1,i+1)+doubleimsrc(j+1,i-1)-doubleimsrc(j-1,i-1)-doubleimsrc(j+1,i+1);ift0t=-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);ift0t=-t;endendendSout=Sout+t;%污染图象邻域边缘能量计算KBlur=Sout/Sin;end%-------------------------------KBlur计算完毕------------------------------%-------------------------------质量指数Q计算------------------------------ifQenable==1Q=0.0;%Q赋初值Qnum=0;%图象以7X7块大小计算每块的Q，逐象素的移动块窗口，这里Qnum为块数量的计数forj=4:H-3fori=4:W-3midsrc=