%Reads in a color image, blurs it with circular PSF of radius R pixels (diameter 2R+1)
%adds noise, then deblurs it with a Wiener filter with Tikhonov parameter lambda.
%Why do we need Tikhonov regularization? Try setting lambda=0.
clear;X=imread('toucan.jpg');X=double(X);[N,M,D]=size(X);%Know D=3.
lambda=0.1;sigma=10;R=15;R1=2*R+1;
N1=N+2*R;M1=M+2*R;%blurred image is N1XM1.
%Create (roughly) circular PSF of radius R centered at MATLAB coordinates (R+1,R+1):
H(R1,R1)=0;for n=1:R1;for m=1:R1;if((n-R-1)*(n-R-1)+(m-R-1)*(m-R-1)<R*R+1);H(n,m)=1;end;end;end;
Y(:,:,1)=conv2(X(:,:,1),H);Y(:,:,2)=conv2(X(:,:,2),H);Y(:,:,3)=conv2(X(:,:,3),H);%Blur image.
Y=Y+sigma*randn(N1,M1,3);%Add 0-mean 2-D white Gaussian noise to blurred image.
FY=fft2(Y,N1,M1);FH=fft2(H,N1,M1);%zero-pad DFT of PSF to size of blurred image Y.
W=conj(FH)./(FH.*conj(FH)+lambda^2);%2-D Wiener filter wavenumber response.
FXHAT(:,:,1)=FY(:,:,1).*W;FXHAT(:,:,2)=FY(:,:,2).*W;FXHAT(:,:,3)=FY(:,:,3).*W;XHAT=real(ifft2(FXHAT));
XHAT=(XHAT-min(min(min(XHAT))))/(max(max(max(XHAT)))-min(min(min(XHAT))));%Scale for imagesc.
Y=(Y-min(min(min(Y))))/(max(max(max(Y)))-min(min(min(Y))));%Scale for imagesc.
figure,imagesc(Y),axis off,title('Blurred image')
figure,imagesc(XHAT),axis off,title('Reconstructed image')%This is zero-padded.