%
% Optical Imaging and Spectroscopy
%
% David J. Brady
% Duke University
% www.opticalimaging.org
%
% Figure 8.10 Tikhonov and tSVD reconstruction for shift codes
%
%
figure(1);
set(gcf,'color','white');
colormap gray;
load flowers;
dog=dog/max(max(dog));
sigma=0.01;
lambda=0.3;
%
% Tikhonov recon
%
N=4;
Hl=toeplitz([1/N zeros(1,255)],[ones(1, N)/N,zeros(1,256-N)]);
Hr=Hl;
dData=Hl*dog*Hr'+sigma*randn(size(dog));
dNoise=pinv(Hl'*Hl+lambda^2*eye(256))*Hl'*dData*Hr*pinv(Hr*Hr'+lambda^2*eye(256));
subplot(3,2,1);imagesc(dNoise(100:150,100:150));axis off; axis image;
title(['Tikhonov mse=' num2str(mean(mean((dNoise-dog).^2)),'%5.2e')]);
%
%
Hl=toeplitz([1, zeros(1,255)],[1 1 0 0 0 1 0 1,zeros(1,256-8)])/4;
Hr=Hl;
dData=Hl*dog*Hr'+sigma*randn(size(dog));
dNoise=pinv(Hl'*Hl+lambda^2*eye(256))*Hl'*dData*Hr*pinv(Hr*Hr'+lambda^2*eye(256));
subplot(3,2,3);imagesc(dNoise(100:150,100:150));axis off; axis image;
title(['1100101 mse=' num2str(mean(mean((dNoise-dog).^2)),'%5.2e')]);
%
%
Hl=toeplitz([1, zeros(1,255)],[1 1 1 0 0 1 0 0,zeros(1,256-8)])/4;
Hr=Hl;
dData=Hl*dog*Hr'+sigma*randn(size(dog));
dNoise=pinv(Hl'*Hl+lambda^2*eye(256))*Hl'*dData*Hr*pinv(Hr*Hr'+lambda^2*eye(256));
subplot(3,2,5);imagesc(dNoise(100:150,100:150));axis off; axis image;
title(['11100100 mse='  num2str(mean(mean((dNoise-dog).^2)),'%5.2e')]);
%
% tsvd
%
rl=125;
rr=125;
N=4;
Hl=toeplitz([1/N zeros(1,255)],[ones(1, N)/N,zeros(1,256-N)]);
Hr=Hl;
dData=Hl*dog*Hr'+sigma*randn(size(dog));
[al bl cl]=svd(Hl);[ar br cr]=svd(Hr);
hlpi=cl(:,1:rl)*pinv(bl(1:rl,1:rl))*al(:,1:rl)';
dNoise=hlpi*dData*hlpi';
subplot(3,2,2);imagesc(dNoise(100:150,100:150));axis off; axis image;
title(['TSVD mse=' num2str(mean(mean((dNoise-dog).^2)),'%5.2e')]);
%
%
Hl=toeplitz([1, zeros(1,255)],[1 1 0 0 0 1 0 1,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);
hlpi=cl(:,1:rl)*pinv(bl(1:rl,1:rl))*al(:,1:rl)';
dNoise=hlpi*dData*hlpi';
subplot(3,2,4);imagesc(dNoise(100:150,100:150));axis off; axis image;
title(['11000101 mse=' num2str(mean(mean((dNoise-dog).^2)),'%5.2e')]);
%
%
Hl=toeplitz([1, zeros(1,255)],[1 1 1 0 0 1 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);
hlpi=cl(:,1:rl)*pinv(bl(1:rl,1:rl))*al(:,1:rl)';
dNoise=hlpi*dData*hlpi';
subplot(3,2,6);imagesc(dNoise(100:150,100:150));axis off; axis image;
title(['11100100 mse='  num2str(mean(mean((dNoise-dog).^2)),'%5.2e')]);