clc; clear; close all;

%- compute the hyperbolic Radon Transform and the Euclidean Radon
%- transform. Intertwine geometries

R = 1.2;

n = 256;

% im = myMap('nicebumps',n,0.1,3,0.5);
% im = myMap('nicedisks',n,0.1,10,0.75);
% im = myMap('mickey',n,0.1,3,0.5) ...
%     + myMap('nicedisks',n,0.1,10,0.75);
% im = myMap('disk',n,0.1,0.5,0.5) + 2*myMap('disk',n,0.1,0.5,-0.5) + 3*myMap('disk',n,0.1,0.1,-0.1);
% im = myMap('spiral',n,0.03, 0., 0.);
im = myMap('disk',n,0.2, 0.4, 0.4)- myMap('disk',n,0.1, 0.4, 0.4) +  myMap('disk',n,0.1, -0.5, -0.5);
% im = myMap('yinyang',n,0.03, 0., 0.);


nth = n; dth = 2*pi/nth; % number of thetas and theta-increment
thetas = dth/2:dth:2*pi-dth/2;
ns = n; ds = 2/ns; % number of s and s-increment
ss = -1+ds/2:ds:1-ds/2;

%%---- hyperbolic Radon transform
tic
hRT = hyperRadon(im,ss,thetas,100,R); toc

%%---- regular radon data (cropped to square size)
tic
RT = myRadon(im,-thetas,'4pt');
RT = RT(end/2-n/2+1:end/2+n/2,:); toc

figure(3);
subplot(1,4,1); dispmat(im, 'image');
subplot(1,4,2); dispmat(hRT,['hyperbolic RT, R=',num2str(R)]);
subplot(1,4,3); dispmat(RT,'RT');
subplot(1,4,4); dispmat(abs(hRT-RT),'difference');

%%---- hyperbolic backprojection
tic
imBP = hyperbackproj(hRT, ss, thetas, n, R); toc

%%---- pull out the singularities
imSing = laplacian(imBP).*myMap('disk',n,0.95,0,0);

figure(4);
subplot(1,3,1); dispmat(im,'f');
subplot(1,3,2); dispmat(imBP,'R_h^t R_h f');
subplot(1,3,3); dispmat(imSing,'\Delta R_h^t R_h f')

% %%--- mix regular Radon and hyperbolic backprojection
% imBP2 = hyperbackproj(RT, ss, thetas, n,R);
% 
% imSing2 = laplacian(imBP2).*myMap('disk',n,0.95,0,0);
% 
% figure(5);
% subplot(1,3,1); dispmat(im,'f');
% subplot(1,3,2); dispmat(imBP2,'R^t R_h f');
% subplot(1,3,3); dispmat(imSing2,'\Delta R^t R_h f')
% 
% % mix hyperbolic Radon and regular backprojection
% 
% % zero-pad the hRT data like regular radon data
% zp = ceil( (sqrt(2)-1)*n/2 );
% hRT = [zeros(zp,nth); hRT; zeros(zp,nth)];
% 
% imBP3 = myLambdaIRadon(hRT, -thetas, n, 0);
% imSing3 = laplacian(imBP3).*myMap('disk',n,0.95,0,0);
% 
% figure(6);
% subplot(1,3,1); dispmat(im,'f');
% subplot(1,3,2); dispmat(imBP3,'R_h^t R f');
% subplot(1,3,3); dispmat(imSing3,'\Delta R_h^t R f')

