%程序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擴(kuò)張為(y1+y2)×(y1+y2) 的矩陣%
he([y/4:(y/4)+y1-1],[y/4:(y/4)+y2-1])=h0;
%將f0擴(kuò)張為(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(:)))]);

%開始執(zhí)行迭代盲目去卷積法 %
f0_=rand(64,64);           %猜測的初始圖像  
%figure,imshow(real(f0_),[min(real(f0_(:))) max(real(f0_(:)))]);
F0_=fft2(f0_);
%加上傅立葉域的條件限制%
H0=G.*(F0_.^(-1));
h0=ifft2(H0);
%加上模糊的條件限制，將猜測圖像的區(qū)域范圍之外設(shè)為零%
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;                    %設(shè)定參數(shù)    

%進(jìn)行循環(huán)迭代的步驟%
for i=1:1500              %設(shè)定迭代的次數(shù)
   i
   for k=1:64
      for l=1:64
         %加上傅立葉域的條件限制%
         if abs(G(k,l))<=0.05;     %設(shè)定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);
   
   %加入圖像域的條件限制，將負(fù)值的像素值設(shè)為零%
   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

%顯示恢復(fù)后的圖像%
figure,imshow(real(f0_),[min(real(f0_(:))) max(real(f0_(:)))]);
figure,imshow(real(h0_),[min(real(h0_(:))) max(real(h0_(:)))]);