% EX_DECONV.M

figure(1)
clf

figure(2)
clf

load lenna.mat % (15,40,2)
load bible.mat % (25,40,1)
load paris.mat % (11,40,2) (51,40,4) (11,20,5/6) 
load annak1.mat % (13,50,1) (13,40,3) (13,40,4) (9,20,5)
load annak2.mat % (25,40,1)
load galaxie.mat % (13,50,1)

signal = 6 ; % 1 = lenna; 2 = bible; 3 = Paris; 4 = anna k1; 5 = anna k2 ; 6 = galaxie
d = 13 ;     % blur size
SNR = 50 ;   % SNR
choice = 1 ;
; % 1 = Gaussian; 2 = square; 3 = vert. echo; 
             % 4 = horiz. echo; 5 = horiz. move; 6 = vert. move

correc= 10.^(-SNR/20);

if signal == 1
	
	A = i1;
	
elseif signal == 2

	A = KK;
	
elseif signal == 3

	A = PP;
	
elseif signal == 4

	A = AKK1;
	
elseif signal == 5

	A = AKK2;
    
elseif signal == 6
    
    A = GG ;
	
end

[N,M] = size(A);
n = log2(N);

P = 10^(-SNR/10); 

if choice == 1;
	
%	H = window(d,'Gauss')*window(d,'Gauss')';
    H = window(@gausswin,d)*window(@gausswin,d)';

elseif choice == 2

	H = ones(d);
	
elseif choice == 3

	H = ones(d);
	H(2:d-1,:) = zeros(d-2,d);
	
elseif choice == 4

	H = ones(d);
	H(:,2:d-1) = zeros(d,d-2);
	
elseif choice == 5

	H = zeros(d);
    H((d+1)/2,:) = ones(1,d);
	
elseif choice == 6

	H = zeros(d);
	H(:,(d+1)/2) = ones(d,1);
	
end

figure(1)

subplot(221)
imagesc(A)
colormap(gray)
axis('square')
set(gca,'XTick',[])
set(gca,'YTick',[])
title(['original ',int2str(N),' x ',int2str(N)])

pause

figure(2)

Af = fft2(A,2^(n+1),2^(n+1));

subplot(331)
imagesc(fftshift(log(abs(Af))))
axis('square')
set(gca,'XTick',[])
set(gca,'YTick',[])
title('original')

pause

B = conv2(A,H);
EB = sum(sum(B.^2));
F = randn(size(B));
EF = sum(sum(F.^2));
B = B + P*F/sqrt(EF)*sqrt(EB);

figure(1)

subplot(222)
imagesc(B(1+(d-1)/2:N+(d-1)/2,1+(d-1)/2:N+(d-1)/2))
colormap(gray)
axis('square')
set(gca,'XTick',[])
set(gca,'YTick',[])

if choice == 1
	
	title(['with Gaussian blur ',int2str(d),' x ',int2str(d),' ; SNR = ',int2str(fix(SNR)),' dB'])
	
elseif choice == 2

	title(['with square blur ',int2str(d),' x ',int2str(d),' ; SNR = ',int2str(fix(SNR)),' dB'])
	
elseif choice == 3

	title(['with vertical echo ',int2str(d),' x ',int2str(1),' ; SNR = ',int2str(fix(SNR)),' dB'])
	
elseif choice == 4

	title(['with horizontal echo ',int2str(1),' x ',int2str(d),' ; SNR = ',int2str(fix(SNR)),' dB'])
	
end

pause

Bf = fft2(B,2^(n+1),2^(n+1));
Hf = fft2(H,2^(n+1),2^(n+1));
% AAf = Bf./(Hf+10*P);
AAf = Bf./Hf;
AA = ifft2(AAf);


figure(2)

subplot(332)
imagesc(fftshift(log(max(eps,abs(Hf)))))
axis('square')
set(gca,'XTick',[])
set(gca,'YTick',[])
title('blur')

subplot(333)
imagesc(fftshift(log(max(eps,abs(Bf)))))
axis('square')
set(gca,'XTick',[])
set(gca,'YTick',[])
title('blurred image')

pause

subplot(335)
imagesc(fftshift(-log(max(eps,abs(Hf)))))
axis('square')
set(gca,'XTick',[])
set(gca,'YTick',[])
title('inverse filter')

pause

subplot(334)
imagesc(fftshift(log(max(eps,abs(AAf)))))
axis('square')
set(gca,'XTick',[])
set(gca,'YTick',[])
title('inverse filtered')

pause

dd = 3;
reg = window(@gausswin,dd)*window(@gausswin,dd)';
Hff = conv2(Hf,reg);
Hf = Hff(1+(dd-1)/2:2^(n+1)+(dd-1)/2,1+(dd-1)/2:2^(n+1)+(dd-1)/2);

figure(1)

subplot(224)
imagesc(real(AA(1:N,1:N)))
colormap(gray)
axis('square')
set(gca,'XTick',[])
set(gca,'YTick',[])
title('inverse filtered')

% -10*log10(sum(sum((A-real(AA(1:N,1:N))).^2))/sum(sum(A.^2)))

pause

AAAf = Bf./(Hf+correc./conj(Hf));
%AAAf = Bf./(Hf+correc);
AAA = ifft2(AAAf);

figure(2)

subplot(338)
imagesc(fftshift(-log(max(eps,abs(Hf+correc)))))
axis('square')
set(gca,'XTick',[])
set(gca,'YTick',[])
title('regularized')

pause

subplot(337)
imagesc(fftshift(log(max(eps,abs(AAAf)))))
axis('square')
set(gca,'XTick',[])
set(gca,'YTick',[])
title('restored')

pause

figure(1)

subplot(223)
imagesc(real(AAA(1:N,1:N)))
colormap(gray)
axis('square')
set(gca,'XTick',[])
set(gca,'YTick',[])
title('restored')

% -10*log10(sum(sum((A-real(AA(1:N,1:N))).^2))/sum(sum(A.^2)))