%This blurs a MRI head image with a disk PSF, adds noise, and deconvolves it
%using (a) deterministic Wiener filter and then (b) stochastic Wiener filter.
%The stochastic reconstruction is visibly better than the deterministic reconstruction.
%Note: this may take a few seconds to run.
clear;load MRI.mat;sigma=50;lambda=1;
%Next 2 lines:PSF is disk-shaped.
H(25,25)=0;H(1,1)=1;for I=1:25;for J=1:25;
if((I-13)*(I-13)+(J-13)*(J-13)<145);H(I,J)=1;end;end;end;
Y=conv2(X,H);M=size(Y,1);N=sigma*randn(M,M);Y=Y+N;%Blur image and add noise.
figure,imagesc(X),colormap(gray),axis off,title('Original image f[n,m]')
figure,imagesc(Y),colormap(gray),axis off,title('Noisy blurred image g[n,m]')
FH=fft2(H,M,M);FY=fft2(Y,M,M);%Zero-pad to size MXM of blurred image.
%Next 2 lines implement Eq. (9.112):
VH=ones(M,1)*[1:M/2];VH=[VH fliplr(VH)];VV=[1:M/2]'*ones(1,M);VV=[VV;flipud(VV)];
VVV=(VV.*VV+VH.*VH)*4*pi*pi/M/M;%Omega_i=2(pi)k_i/N in each direction.
%Next 2 lines are dterministic and stochastic Wiener filters:
FW1=FY.*conj(FH)./(abs(FH).*abs(FH)+lambda*lambda);W1=real(ifft2(FW1));%Deterministic
FW2=FY.*conj(FH)./(abs(FH).*abs(FH)+VVV*lambda*lambda);W2=real(ifft2(FW2));%Stochastic
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')


