% tomography.m
% Reconstruct density from tomography
% JMS, Jan.2020

% Set data
clear
nx = 9;    % pixels along x axis
ny = nx;   % pixels along y axis

% Set ray angles and lengths
a = [0; atan2(1,2); pi/4; atan2(2,1); pi/2; atan2(2,-1); 3*pi/4; atan2(1,-2)];
d = [1;sqrt(5);sqrt(2);sqrt(5);1;sqrt(5);sqrt(2);sqrt(5)];  
na = numel(a);  % number of angles

% Read tomography data
b = load('tomography.dat');
nr = numel(b);               % number of rays

% Find matrix A such that b=A*rho
A = zeros(nr,nx*ny);
for i = 1:nx                      % pixel index along x
  for j = 1:ny                    % pixel index along y
    p = (i-1)*ny+j;               % combined pixel index
    r(1) = j;                     % ray index for angle 1
    r(2) = nx+ny-1-i+2*j;         % ray index for angle 2
    r(3) = 2*nx+3*ny-2-i+j;
    r(4) = 4*nx+4*ny-3-2*i+j;
    r(5) = 4*nx+5*ny-5+i;
    r(6) = 5*nx+5*ny-7+2*i+j;
    r(7) = 7*nx+6*ny-8+i+j;
    r(8) = 8*nx+7*ny-10+i+2*j;    % ray index for angle 8
    A(r,p) = d;
  end
end

% Invert tomography signal to find density
rho = A\b;
rho = reshape(rho,ny,nx);

% Plot density using image
figure(1)
image(rho*256)
colorbar
axis equal

% Plot density using surf
figure(2)
x = 1:nx;
y = -(1:ny);            % minus sign to plot matrix: y=-j
[x,y] = meshgrid(x,y);
surf(x,y,rho)

% Plot density using contour
figure(3)
contour(x,y,rho)
axis equal

