% test Radon transform using a temperature field
close all
clear all
clc
x=0:0.1:6;
y=0:0.1:6;
[xx,yy]=meshgrid(x,y);
t=200*exp(-(sqrt((xx-3.2).^2+(yy-3.2).^2)-1.5).^2/0.2);
% surf(xx,yy,t)
figure (1);
mesh(xx,yy,t)
imagesc(t); colormap(hot); colorbar
% refraction index
n = 1 + 0.292015*10^(-3)./(1+0.368184*10^(-2).*t);
figure (2);
n = n-1;
mesh(xx,yy,n)
%% radon transform for projection
theta = 0:3:179;
[R,xp] = radon(n,theta);
% figure; imagesc(theta,xp,R);
% title('R_{\theta} (X\prime)');
% xlabel('\theta (degrees)');
% ylabel('X\prime');
% set(gca,'XTick',0:20:180);
% colormap(hot);
% colorbar
I1 = iradon(R,theta);
figure;
subplot(1,2,1), imagesc(n), title('Original')
colormap(hot);
colorbar
I1 = I1(1:end-1,1:end-1);
subplot(1,2,2), imagesc(I1), title('Filtered backprojection')
colormap(hot);
colorbar
figure; mesh(xx,yy,I1)