% session 4: Local tomography and limited data problems: recovering
% singularities..

clc; clear; close all;

n = 128;

% im = myMap('disk', n, 0.01, 0.2, 0.2);
im = myMap('disk', n, 0.2, .5, 0.) + myMap('disk', n, 0.1, .0, 0.5);
% im = myMap('line', n, 0.1, .0, 0.).*myMap('disk', n, 0.5, .0, 0.);
% im = myMap('nicebumps', n, .1, 3, .4);
% im = myMap('mickey',n,0.,0.,0.) + myMap('gaussian',n,0.1,0,0);

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

method = '4pt';

%--------- LAMBDA TOMOGRAPHY

%---- compute the radon transform and display it
% nth = 2*n; 
% h = 2/n;
% thetas = pi*(1/2/nth:1/nth:1-1/2/nth);
% RT = myRadon(im, thetas, method);
% 
% figure(3); dispmat(RT,'RT');
% 
% % lambda tomography
% imRc = myLambdaIRadon(RT, thetas, n,0);
% imRc2 = myLambdaIRadon(RT, thetas, n,1);
% imRc3 = myIRadon(RT, 'shepp-logan', thetas, 1., n);
% 
% figure(4); 
% subplot(1,3,1);  dispmat(imRc,'(-\Delta)^{-1/2} f');
% subplot(1,3,2);  dispmat(imRc2,'(-\Delta)^{1/2} f');
% subplot(1,3,3);  dispmat(imRc3,'inverse radon');

%--------- LIMITED ANGLE

% %---- compute the radon transform and display it
% nth = n; 
% thetas = pi*(1/2/nth:1/nth:1-1/2/nth)/2;
% 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');
% 
% figure(5);
% xs = -1+1/n:2/n:1-1/n;
% plot(xs,im(n/2,:),'b',xs,imRc(n/2,:),'r'); legend('f','f_r_c');
% title('cut plot at "y=0"')

%--------- EXTERIOR PROBLEM

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

% radius up to which we erase the Radon Transform
a = .2;

nr = size(RT,1);
ts = 1:h:(nr*h);
ts = ts + h/2-ts(nr/2+1);

% erase the central part of RT - choose how to erase
eraser = ones(size(RT,1),1);
erase_type = 1;
switch erase_type
    case 1 
        % erase linearly up to |t|=a so the limited data is still continuous
        eraser((abs(ts))<a) = abs(ts((abs(ts))<a))/a;
    case 2
        % erase linearly up to |t|=a so the limited data is still continuous
        eraser((abs(ts))<a) = (abs(ts((abs(ts))<a))/a).^2;
    otherwise
        % set the middle part to zero
        eraser((abs(ts))<a) = 0.;
end
   
RTl = RT.*eraser(:,ones(1,size(RT,2)));

figure(3); 
subplot(1,2,1); dispmat(RT,'RT');
subplot(1,2,2); dispmat(RTl,'limited RT');

% the question remains as to how to fill in the missing part in such a way
% that the reconstructions are the least affected by this missing data, or
% satisfy some smoothness conditions. 

%---- compute the inverse Radon Transform and display it
filter = 'ram-lak'; d = .9;
[imRc,H] = myIRadon(RTl, filter, thetas, d, n);

% im = im.*(myMap('disk',n,1,0,0)-myMap('disk',n,a,0,0));
imRc = imRc.*(myMap('disk',n,1,0,0)-myMap('disk',n,a,0,0));

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');

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