%-----------------------------------------------------------%
%------------------------圖像恢復---------------------------%
%-----------------------------------------------------------%
f=double(imread('Cameraman.bmp'));%輸入實驗圖像 灰度圖像
%-----------------------------------------------------------%
%------------------------圖像模糊---------------------------%
%-----------------------------------------------------------%
[A B]=size(f);%f是A×B尺寸
%沖激函數模板 使用9*9模板 即h是9×9尺寸，設為C×D
for i=1:9
    for j=1:9
          h(i,j)=exp(sqrt((i-5)^2+(j-5)^2)/240);
    end
end
h=h/sum(sum(h));
[C D]=size(h);
%對f和h進行延拓，周期延拓為M=A+C-1，N=B+D-1
M=A+C-1;
N=B+D-1;
fe=f;
fe(M,:)=0;
fe(:,N)=0;
he=h;
he(M,:)=0;
he(:,N)=0;
%將f和h的延拓結果進行卷積得到模糊圖像
g=conv2(fe,he);
ge=g(1:M,1:N);
%添加高斯隨機噪聲
ne=randn([M,N]);
ge_8=ge+sqrt(8)*ne;%均值為0，方差為8的高斯隨機噪聲
ge_16=ge+sqrt(16)*ne;%均值為0，方差為16的高斯隨機噪聲
ge_32=ge+sqrt(32)*ne;%均值為0，方差為32的高斯隨機噪聲

figure(1)
subplot(2,2,1);
imshow(f,[]);
title('Cameraman.bmp');
subplot(2,2,2);
imshow(ge_8((M-A)/2+1:(M+A)/2,(N-B)/2+1:(N+B)/2),[]);
title('模糊并帶有方差為8高斯隨機噪聲的圖像');
subplot(2,2,3);
imshow(ge_16((M-A)/2+1:(M+A)/2,(N-B)/2+1:(N+B)/2),[]);
title('模糊并帶有方差為16高斯隨機噪聲的圖像');
subplot(2,2,4);
imshow(ge_32((M-A)/2+1:(M+A)/2,(N-B)/2+1:(N+B)/2),[]);
title('模糊并帶有方差為32高斯隨機噪聲的圖像');
%-----------------------------------------------------------%
%-------------------------逆濾波----------------------------%
%-----------------------------------------------------------%
%取頻域
G_8=fft2(ge_8);
G_16=fft2(ge_16);
G_32=fft2(ge_32);

H=fft2(he);

%為了防止H(u,v)在UV平面上取0或很小且消除振鈴效應
%如果  H(u,v)<=d,   W(u,v)=k
%      H(u,v)>d,    W(u,v)=1/H(u,v)
% 其中k和d為小于1的常數，且選的較小為好
d=0.1;
k=0.1;
for u=1:M
    for v=1:N
        if abs(H(u,v))<=d
            W(u,v)=k;
        else
            W(u,v)=1/H(u,v);
        end
    end
end

f_8_1=abs(ifft2(G_8.*W));
f_16_1=abs(ifft2(G_16.*W));
f_32_1=abs(ifft2(G_32.*W));

figure(2)
subplot(2,2,1);
imshow(f,[]);
title('Cameraman.bmp');
subplot(2,2,2);
imshow(f_8_1(1:A,1:B),[]);
title('逆濾波恢復結果(噪聲方差8)');
subplot(2,2,3);
imshow(f_16_1(1:A,1:B),[]);
title('逆濾波恢復結果(噪聲方差16)');
subplot(2,2,4);
imshow(f_32_1(1:A,1:B),[]);
title('逆濾波恢復結果(噪聲方差32)');
%-----------------------------------------------------------%
%------------------------維納濾波---------------------------%
%-----------------------------------------------------------%
F=fft2(fe);

%高斯隨機噪聲的FFT
N_8=fft2(sqrt(8)*ne); 
N_16=fft2(sqrt(16)*ne);
N_32=fft2(sqrt(32)*ne);

F_8=conj(H)./(abs(H).^2+abs(N_8).^2./abs(F).^2).*G_8;
F_16=conj(H)./(abs(H).^2+abs(N_16).^2./abs(F).^2).*G_16;
F_32=conj(H)./(abs(H).^2+abs(N_32).^2./abs(F).^2).*G_32;

f_8_2=abs(ifft2(F_8));
f_16_2=abs(ifft2(F_16));
f_32_2=abs(ifft2(F_32));

figure(3)
subplot(2,2,1);
imshow(f,[]);
title('Cameraman.bmp');
subplot(2,2,2);
imshow(f_8_2(1:A,1:B),[]);
title('維納濾波恢復結果(噪聲方差8)');
subplot(2,2,3);
imshow(f_16_2(1:A,1:B),[]);
title('維納濾波恢復結果(噪聲方差16)');
subplot(2,2,4);
imshow(f_32_2(1:A,1:B),[]);
title('維納濾波恢復結果(噪聲方差32)');