% session on the reconstruction of the attenuated Radon transform

clc; clear; close all;

n = 128;

% a = zeros(n,n);
a = myMap('disk',n,1,0,0) + myMap('nicebumps',n,0.1,3,0.4);
% a = myMap('gaussian',n,0.1,.5,0.);
f = myMap('mickey',n,1,0,0) + 5*myMap('disk',n,0.05,.5,0);
% f = myMap('nicebumps',n,0.1,3,0.4);

filter = 'cosine';

figure(1);
subplot(1,2,1); dispmat(a,'attenuation');
subplot(1,2,2); dispmat(f,'source');

nth = 2*n;
thetas = 2*pi*(1/2/nth:1/nth:1-1/2/nth);

% attenuated Radon Transform and comparison with the regular RT
rotMethod = '4pt';
atRT = atRt(a,f,thetas,rotMethod);
RT = myRadon(f, thetas, rotMethod);

figure(3); 
subplot(1,3,1); dispmat(atRT,'atRt');
subplot(1,3,2); dispmat(RT,'RT');
subplot(1,3,3); dispmat(RT-atRT,'RT-atRT');

% inversion of the atRt
d = .9; % regularization parameter
fRc = iAtRt(a, atRT, thetas, rotMethod, filter, d);

figure(4); 
subplot(1,3,1); dispmat(fRc,'fRc');
subplot(1,3,2); dispmat(f,'f');
subplot(1,3,3); dispmat(abs(fRc-f),'error');

figure(5);
xs = -1+1/n:2/n:1-1/n;
plot(xs,f(n/2,:),'b',xs,fRc(n/2,:),'r'); legend('f','f_r_c');
title('cut plot at "y=0"')

%---- impact of the regularization parameter

% xs = -1+1/n:2/n:1-1/n;
% figure(6);
% plot(xs,f(n/2,:)); hold on;
% 
% for filter={'ram-lak','shepp-logan','cosine','hamming','hann'},%d=0.1:0.1:0.9,
%     fRc = iAtRt(a, atRT, thetas, rotMethod, char(filter), d);
%     plot(xs,fRc(n/2,:)); hold on;
% end
%     
