% DIP experiment 6
% 逆濾波,維納濾波
clc
clear all;

A = imread('Cameraman.bmp');
figure(1),subplot(221);imshow(A,[]);title('原圖像');
figure(2),subplot(221);imshow(A,[]);title('原圖像');
figure(3),subplot(221);imshow(A,[]);title('原圖像');
A = double(A);
[M,N] = size(A);

%產生模糊沖擊函數的模板,其中n為模板的大小
n = 7;
h = BlurTemplate(n);
sumh = sum(sum(h));
h = h/sumh;

%為了消除混疊效應引起的誤差，需要對模板H及原圖像的尺寸進行延拓
%對模板尺寸進行延拓
expandh = zeros(M+n-1,N+n-1);
expandh(1:n, 1:n) = h; 
%對圖像尺寸進行延拓
expandA = zeros(M+n-1,N+n-1);
expandA(1:M, 1:N) = A; 


%在頻域對圖像進行模糊
%對延拓后的expandH進行離散傅立葉變換，為了提高運算效率，故調用matlab中自帶的fft2函數
FH = fft2(expandh);
%對延拓后的圖像進行離散傅立葉變換
FA = fft2(expandA);
%退化操作
FbA = FA .* FH;
%進行傅立葉反變換得到退化后的模糊圖像
bA = abs(ifft2(FbA));
% figure(4),imshow(bA(1+ceil(n/2):M+floor(n/2),1+ceil(n/2):N+floor(n/2)),[]);
% title('與沖擊函數卷積模糊后的圖像')


%對模糊后的圖像疊加高斯噪聲
nbA1 = bA + sqrt(8)*randn([M+n-1,N+n-1]);%疊加均值為0，方差為8的高斯隨機噪聲
nbA2 = bA + sqrt(16)*randn([M+n-1,N+n-1]);%疊加均值為0，方差為16的高斯隨機噪聲
nbA3 = bA + sqrt(32)*randn([M+n-1,N+n-1]);%疊加均值為0，方差為32的高斯隨機噪聲

%畫出模糊并添加高斯噪聲后的的圖像
figure(1),subplot(222);imshow(nbA1(1+floor(n/2):M+floor(n/2),1+floor(n/2):N+floor(n/2)),[]);title('模糊并添加均值為8高斯噪聲后的圖像');
figure(2),subplot(222);imshow(nbA2(1+floor(n/2):M+floor(n/2),1+floor(n/2):N+floor(n/2)),[]);title('模糊并添加均值為16高斯噪聲后的圖像');
figure(3),subplot(222);imshow(nbA3(1+floor(n/2):M+floor(n/2),1+floor(n/2):N+floor(n/2)),[]);title('模糊并添加均值為32高斯噪聲后的圖像');

%對模糊并添加噪聲后的圖像求離散傅立葉變換
FnbA1 = fft2(nbA1);
FnbA2 = fft2(nbA2);
FnbA3 = fft2(nbA3);

%********************************逆濾波********************************

%為了防止H(u,v)在UV平面上取0或很小且消除振鈴效應，取恢復轉移函數FrH(u,v)為
% 如果  H(u,v)<d,   FrH(u,v)=k
%       H(u,v)>=d,  FrH(u,v)=1/H(u,v)
% 其中d為小于1的常數，且選的較小為好
FrH = zeros(M+n-1,N+n-1);
for i = 1:M+n-1
    for j = 1:N+n-1
        if abs(FH(i,j)) < 0.1
            FrH(i,j) = 1;
        else
            FrH(i,j) = 1/FH(i,j);
        end
    end
end
%逆濾波過程
FrnbA1 = FnbA1 .* FrH;
FrnbA2 = FnbA2 .* FrH;
FrnbA3 = FnbA3 .* FrH;
%得到逆濾波復原后的圖像
rnbA1 = abs(ifft2(FrnbA1));
rnbA2 = abs(ifft2(FrnbA2));
rnbA3 = abs(ifft2(FrnbA3));
%畫出逆濾波復原后的圖像
figure(1),subplot(223);imshow(rnbA1(1:M,1:N),[]);title('逆濾波復原的圖像');
figure(2),subplot(223);imshow(rnbA2(1:M,1:N),[]);title('逆濾波復原的圖像');
figure(3),subplot(223);imshow(rnbA3(1:M,1:N),[]);title('逆濾波復原的圖像');


%***********************************維納濾波**************************************

%為求得噪聲的功率譜對噪聲進行離散傅立葉變換
Fn1 = fft2(nbA1 - bA);
Fn2 = fft2(nbA2 - bA);
Fn3 = fft2(nbA3 - bA);
%求得維納濾波器轉移函數
W1 = conj(FH) ./ abs( FH.*conj(FH) + (Fn1.*conj(Fn1)) ./ (FA.*conj(FA)) );
W2 = conj(FH) ./ abs( FH.*conj(FH) + (Fn2.*conj(Fn2)) ./ (FA.*conj(FA)) );
W3 = conj(FH) ./ abs( FH.*conj(FH) + (Fn3.*conj(Fn3)) ./ (FA.*conj(FA)) );
%維納濾波過程
FwrnbA1 =  FnbA1 .* W1;
FwrnbA2 =  FnbA2 .* W2;
FwrnbA3 =  FnbA3 .* W3;
%求得維納濾波復原后的圖像
wrnbA1 = abs(ifft2(FwrnbA1));
wrnbA2 = abs(ifft2(FwrnbA2));
wrnbA3 = abs(ifft2(FwrnbA3));
%畫出維納濾波復原后的圖像
figure(1),subplot(224);imshow(wrnbA1(1:M,1:N),[]);title('維納濾波復原的圖像');
figure(2),subplot(224);imshow(wrnbA2(1:M,1:N),[]);title('維納濾波復原的圖像');
figure(3),subplot(224);imshow(wrnbA3(1:M,1:N),[]);title('維納濾波復原的圖像');