%
% Optical Imaging and Spectroscopy
%
% David J. Brady
% Duke University
% www.opticalimaging.org
%
% Figure 8.14
% 
%
%
figure(2);
set(gcf,'color','white');
colormap gray;
%dog=double(imread('cameraman.tif'));
%dog=double(imread('nycN.jpg'));
load flowers;
dog=dog/max(max(dog));
imagesc(dog);axis image; axis off;
title('original image');
sigma=0.01;
lambda=0.3;
N=4;
N=4;
Hl=toeplitz([1, zeros(1,255)],[1 1 1 1 0 0 0 0,zeros(1,256-8)])/4;
Hr=Hl;
dData=Hl*dog*Hr'+sigma*randn(size(dog));
[al bl cl]=svd(Hl);[ar br cr]=svd(Hr);
rl=100;
rr=100;
hlpi=cl(:,1:rl)*pinv(bl(1:rl,1:rl))*al(:,1:rl)';
hrpi=cr(:,1:rr)*pinv(br(1:rr,1:rr))*ar(:,1:rr)';
dNoise=hlpi*dData*hrpi';
%dNoise=pinv(Hl'*Hl+lambda^2*eye(256))*Hl'*dData*Hr*pinv(Hr*Hr'+lambda^2*eye(256));
fuLG=nullSpaceSmoothing(null(hlpi),dNoise);
fuLG=nullSpaceSmoothing(null(hrpi),fuLG')';
subplot(1,2,1);imagesc(dNoise(100:150,100:150),[0 1]);axis image; axis off;
title(['Tikhonov mse=' num2str(mean(mean((dNoise-dog).^2)),'%5.2e')]);
subplot(1,2,2);imagesc(fuLG(100:150,100:150),[0 1]);axis image; axis off;
title(['LG mse=' num2str(mean(mean((fuLG-dog).^2)),'%5.2e')]);