% session 2: Radon Transform and Filtered-Backprojection Algorithm
% values for the 'filter' variable:
% 'ram-lak', 'shepp-logan', 'cosine', 'hann', 'hamming'

clc; clear; close all;

n = 64;

% im = myMap('disk', n, 0.01, 0.2, 0.2);
% im = myMap('mickey', n, 0, 0, 0);
im = myMap('nicebumps', n, .1, 3, .4);

figure(1); dispmat(im,'image');

method = '4pt';

%---- test image rotation method, display the first 12 rotations
% nth = 12;
% thetas = 2*pi*(1/2/nth:1/nth:1-1/2/nth);
% imRots = rotateImages(im, thetas, method);
% figure(2);
% for j=1:12,
%     subplot(3,4,j); dispmat(imRots(:,:,j), ['\theta=',num2str(thetas(j))]);
% end

%---- compute the radon transform and display it
nth = 256; 
thetas = 2*pi*(1/2/nth:1/nth:1-1/2/nth);
RT = myRadon(im, thetas, method);

figure(3); dispmat(RT,'RT');

%---- compute the inverse Radon Transform and display it
filter = 'ram-lak'; d = .9;
[imRc,H] = myIRadon(RT, filter, thetas, d, n);
figure(4); 
subplot(1,3,1); dispmat(im,'inital');
subplot(1,3,2); dispmat(imRc,'reconstructed');
subplot(1,3,3); dispmat(abs(im-imRc),'error');

%---- study the impact of the parameter d and the filter
%---- d smaller: loss of high-frequencies, ringing effect 
%---- test case: d=.2 on the phantom, artifacts look like the lower bumps
% figure(5); count = 1;
% filter = 'hamming';
% for d = .1:.1:.3,
%     [imRc,H] = myIRadon(RT, filter, thetas, d, n);
%     figure(5); subplot(3,3,count); dispmat(abs(imRc-im),['d=',num2str(d)]);
%     figure(6); 
%     subplot(1,2,1); plot(fftshift(H)); hold on
%     subplot(1,2,2); plot(ifft(H)); hold on
%     count = count+1;
% end
% figure(6); title('filters');

%---- example of white noise
% noise = rand(n);
% figure(7); 
% subplot(1,2,1); dispmat(noise,'noise');
% subplot(1,2,2); dispmat(fftshift(abs(fft2(noise))),'magnitude of noise in Fourier');

%---- add some white noise and study the robustness
% SNR = .025; % signal to noise ratio
% noise = max(max(RT))/SNR*rand(size(RT));
% noisyRT = noise + RT;
% 
% filter = 'shepp-logan';
% imRc2 = myIRadon(noisyRT, filter, thetas, 1., n);
% 
% figure(8);
% ks = 1:size(RT,1);
% plot(ks,abs(fft(RT(:,1))),'b',ks,abs(fft(noise(:,1))),'r');
% legend('signal','noise');
% 
% figure(9);
% subplot(1,2,1); dispmat(noisyRT,['noisy data, SNR=',num2str(SNR)]);
% subplot(1,2,2); dispmat(imRc2,'reconstructed');

% figure(10); count = 1;
% for d = .05:.1:.955,
%     imRc2 = myIRadon(noisyRT, filter, thetas, d, n);
%     subplot(3,3,count); dispmat(imRc2);
%     count = count+1;
% end
