%This motion-blurs a fractal tree image, adds noise, and deblurs the noisy blurred image
%using (a) deterministic Wiener filter and then (b) stochastic Wiener filter.
%It also plots S_x(0,\Omega_2) on a log-log scale, showing the tree is a fractal with a=2.
%The stochastic reconstruction is visibly better than the deterministic reconstruction.
%Note: this may take a few seconds to run.
clear;load tree.mat;%Fractal tree image.
figure,imagesc(XX),colormap(gray),axis off,title('Fractal tree image')
H(1,1:151)=ones(1,151);%Motion blur PSF of duration 151.
Y=conv2(XX,H);N=100*randn(960,1430);Z=Y+N;%Noisy motion-blurred.
figure,imagesc(Y),colormap(gray),axis off,title('Blurred image')
figure,imagesc(Z),colormap(gray),axis off,title('Noisy blurred image')
FX=fft2(XX,960,1430);FZ=fft2(Z,960,1430);FH=fft2(H,960,1430);%zero-pad to size of blurred image.
%Log-log plot of S_x(0,\Omega_2) vs. \Omega_2 shows the tree image is indeed a fractal with a=2:
figure,plot(log10([2:480]*2*pi/960),log10(abs(FX(2:480,1)).^2/960/1280)),...
title('log(S_x(0,\Omega_2)) vs. log(\Omega_2)')
%Next 2 lines implement Eq. (9.112).
VH=ones(960,1)*[1:715];VH=[VH fliplr(VH)];VV=[1:480]'*ones(1,1430);VV=[VV;flipud(VV)];
VVV=(VV.^2+VH.^2)*4*pi^2/960/1430;%Omega_i=2(pi)k_i/N in each direction.
FW1=FZ.*conj(FH)./(abs(FH).^2+1);W1=real(ifft2(FW1));%Deterministic Wiener filter
FW2=FZ.*conj(FH)./(abs(FH).^2+VVV);W2=real(ifft2(FW2));%Stochastic Wiener filter
figure,imagesc(W1),colormap(gray),axis off,...
title('Reconstructed using deterministic Wiener filter')
figure,imagesc(W2),colormap(gray),axis off,...
title('Reconstructed using stochastic Wiener filter')

