%This illustrates motion blur deconvolution on a color Christmas tree image.
%It also shows how to determine the length of the motion blur from the spectrum
%of the blurred image, by examing spacing between zero crossings.
%This may take a few seconds to run.
clear;load xmastree.mat;N=82;lamb=0.01;LH=1080+(N+1)-1;
H=ones(1,N+1)/(N+1);FH=fft2(H,691,LH);%PSF for motion blur.
%Motion-blur separately each RGB component of the Christmas tree image X:
X1=X(:,:,1);FX1=fft2(X1,691,LH);FY1=FX1.*FH;Y1=real(ifft2(FY1));Y(:,:,1)=Y1;%red
X2=X(:,:,2);FX2=fft2(X2,691,LH);FY2=FX2.*FH;Y2=real(ifft2(FY2));Y(:,:,2)=Y2;%green
X3=X(:,:,3);FX3=fft2(X3,691,LH);FY3=FX3.*FH;Y3=real(ifft2(FY3));Y(:,:,3)=Y3;%blue
%Now linearly scale so 0<Y<1 (required by imagesc for color images):
Y=(Y-min(min(min(Y))))/(max(max(max(Y)))-min(min(min(Y))));
figure,imagesc(Y),axis off,title('Motion-blurred Christmas tree image')
%Now plot 1-D slice of 2-D spectrum of motion-blurred image:
figure,plot(log10(abs(FY1(1,1:100)))),title('1-D spectrum log(G_R[k_1,0]+1)'
%Now perform motion-blur deconvolution using a Wiener filter on each RGB component:
FZ1=FY1.*conj(FH)./(abs(FH).*abs(FH)+lamb*lamb);Z1=real(ifft2(FZ1));Z(:,:,1)=Z1;
FZ2=FY2.*conj(FH)./(abs(FH).*abs(FH)+lamb*lamb);Z2=real(ifft2(FZ2));Z(:,:,2)=Z2;
FZ3=FY3.*conj(FH)./(abs(FH).*abs(FH)+lamb*lamb);Z3=real(ifft2(FZ3));Z(:,:,3)=Z3;
%Now linearly scale so 0<Z<1 (required by imagesc for color images):
Z=(Z-min(min(min(Z))))/(max(max(max(Z)))-min(min(min(Z))));
figure,imagesc(Z),axis off,title('Reconstructed Christmas tree image')


