l6_2.m

DigitalImageProcessing_base_on_Matlab 基于Matlab的数字图像处理

%程序L6_2.m: 迭代盲目去卷积法%
clear
%制作模糊的图像%
[h0,map]=bmpread('L6_2a.bmp');       %读取图像h0
[f0,map]=bmpread('L6_2b.bmp');       %读取图像f0
[y1,y2]=size(h0);                    %量测h0大小为y1×y2
[u1,u2]=size(f0);                    %量测f0大小为u1×u2
y=y1+y2;
u=u1+u2;
he=zeros(y,y);          %将(y1+y2)×(y1+y2) 的零矩阵令为he
fe=zeros(u,u);          %将(u1+u2)×(u1+u2) 的零矩阵令为fe    

%将h0扩张为(y1+y2)×(y1+y2) 的矩阵%
he([y/4:(y/4)+y1-1],[y/4:(y/4)+y2-1])=h0;
%将f0扩张为(u1+u2)×(u1+u2) 的矩阵%
fe([u/4:(u/4)+u1-1],[u/4:(u/4)+u2-1])=f0;
h1=he;
f1=fe;
F=fft2(f1);
H=fft2(h1);
G=F.*H;                 %将f与h做卷积得到模糊图像g
g1=ifft2(G);
%显示PSF的图像%
figure,imshow(real(h1),[min(real(h1(:))) max(real(h1(:)))]);
%显示原始图像%
figure,imshow(real(f1),[min(real(f1(:))) max(real(f1(:)))]);
%显示模糊的图像%
figure,imshow(real(g1),[min(real(g1(:))) max(real(g1(:)))]);

%开始执行迭代盲目去卷积法 %
f0_=rand(64,64);           %猜测的初始图像  
%figure,imshow(real(f0_),[min(real(f0_(:))) max(real(f0_(:)))]);
F0_=fft2(f0_);
%加上傅立叶域的条件限制%
H0=G.*(F0_.^(-1));
h0=ifft2(H0);
%加上模糊的条件限制，将猜测图像的区域范围之外设为零%
h0([1:y/4],:)=0;h0([(y/4)+y1-1:y],:)=0;
h0(:,[1:y/4])=0;h0(:,[(y/4)+y1-1:y])=0;
h0_=h0;
H0_=fft2(h0_);
u=0.9;                    %设定参数    

%进行循环迭代的步骤%
for i=1:1500              %设定迭代的次数
   i
   for k=1:64
      for l=1:64
         %加上傅立叶域的条件限制%
         if abs(G(k,l))<=0.05;     %设定noise level为0.05
            F1(k,l)=F0_(k,l);
         elseif abs(H0_(k,l))>=abs(G(k,l));
            F1(k,l)=(1- u).*F0_(k,l)+...
               u.*(G(k,l).*(H0_(k,l).^(-1)));
         else 
            F1(k,l)=(((1-u).*(F0_(k,l).^(-1)))+...
               u.*(H0_(k,l).*(G(k,l).^(-1)))).^(-1);
         end
      end
   end 
   f1=ifft2(F1);
   
   %加入图像域的条件限制，将负值的像素值设为零%
   x1=real(f1);x2=imag(f1);
   xx=find(x1<0);
   x1(xx(:))=0;
   f0__=x1;
   f0__([1:y/4],:)=0;f0__([(y/4)+y1-1:y],:)=0;
   f0__(:,[1:y/4])=0;f0__(:,[(y/4)+y1-1:y])=0;
   f0_=f0__;
   F0_=fft2(f0_);
   
   %傅立叶域的条件限制%
   H0=G.*(F0_.^(-1));
   h0=ifft2(H0);
   
   %模糊域的条件限制%
   h0([1:y/4],:)=0;h0([(y/4)+y1-1:y],:)=0;
   h0(:,[1:y/4])=0;h0(:,[(y/4)+y1-1:y])=0;
   w1=real(h0);w2=imag(h0);
   ww=find(w1<0);
   w1(ww(:))=0;
   h0_=w1;
   H0_=fft2(h0_);
end

%显示恢复后的图像%
figure,imshow(real(f0_),[min(real(f0_(:))) max(real(f0_(:)))]);
figure,imshow(real(h0_),[min(real(h0_(:))) max(real(h0_(:)))]);