n = 201;     % Number of pixels in x and y direcrtions
n_phi = 128; % nunmber of angles for the Radon Transform
%---
h_d = 1; % size of the detector's pixel
%---
x = h_d*linspace(-(n-1)/2,(n-1)/2,n).'; % X coordinates of pixel in the Phantom
y =h_d* linspace(-(n-1)/2,(n-1)/2,n).'; % Y coordinates of pixel in the Phantom and X-ray detecto
phi = linspace(0, pi, n_phi);           % Angles for the Radon transform
h_phi= phi(2)-phi(1);                   % Angular resolution
%---
k = linspace(-pi/h_d,pi/h_d,n).'; % wave vectors conjugate to the Y coordinate
h_k = k(2)-k(1);                  % The resolution in hte k space
%---
M = head_phantom(n, phi(1)); % Function that creates custom Head phantom rotated by an angle phi
%----
% Display original phantom
imagesc([x(1),x(end)],[y(1),y(end)],M);
axis square;
set(gca,'Ydir','normal');
xlabel('x (px)');
ylabel('y (px)');
title('Head Phantom');
colormap(bone);
%---
% Calculate Radon transform and show the resulting sinogram
S = zeros(n,n_phi);
S(:,1) = sum(M,2)*h_d;
for i= 2:n_phi
    M = head_phantom(n, phi(i));
    S(:,i) = sum(M,2)*h_d;
end
%---
figure
imagesc([0,pi],[x(1),x(end)], S);
xlabel('rotation angle');
ylabel('detector position (px)');
title('Sinogram')
colormap(bone);
%---
% Reconstruction
%---
% Calculate Fourier transform of sinogram
FS = zeros(size(S));
%---
for l = 1:numel(phi)
    for p=1:numel(k)
        FS(p,l) = sum(S(:,l).*exp(-1i*k(p)*x))*h_d; 
    end
end
%---
% show resulting spectrum on the obtained polar grid
figure
imagesc([0,pi],[k(1),k(end)],log(abs(FS)));
xlabel('rotation angle');
ylabel('k_r (1/px)');
title('Siniogram Spectrum');
colormap('bone');
%---
% Naive inplementation of hte inverse Fourier transform in polar
% coordinates
%---
RI = zeros(n);
for i=1:numel(x)
    for j=1:numel(y)
        for l = 1:numel(phi)
            RI(j,i) = RI(j,i) + 1/(2*pi)^2*sum(FS(:,l).*abs(k).*exp(1i*k*cos(phi(l)+pi/2)*x(i)+1i*k*sin(phi(l)+pi/2)*y(j)))*h_k*h_phi;
        end
    end
end
%---
figure
imagesc([x(1),x(end)],[y(1),y(end)],abs(RI));
axis square;
set(gca,'Ydir','normal');
xlabel('x (px)');
ylabel('y (px)');
title('Reconstructed Phantom');
colormap(bone);
%---
% Optimized Matlab's implementations of the Radon transform and inverse
% Radon transform
%---
RM = radon(head_phantom(n,0),rad2deg(phi));
IR = iradon(RM,rad2deg(phi),"linear","Ram-Lak",1,n);
%---
% Display obtained reconstructed phantom
%---
figure
imagesc([x(1),x(end)],[y(1),y(end)],IR);
axis square;
set(gca,'Ydir','normal');
xlabel('x (px)');
ylabel('y (px)');
title('Optimized Reconstructed Phantom');
colormap(bone);
%---


function M = head_phantom(n,phi)
    %---
    % Rot matrix
    R = [cos(phi),sin(phi);-sin(phi),cos(phi)];
    %---
    % Elipsoid representing the head
    %---
    mu_0 = [0;0];
    Sigma_0 = [1.5^2,0;0,1^2];
    %---
    mu_0 = R*mu_0; % Rotate ellipse
    Sigma_0 = R*Sigma_0*R.';
    %---
    s = sqrt(2)/min(eig(Sigma_0)); % determine the range of the rotated head
    D_r = 1.2*max([max(mu_0+s),max(mu_0-s)]);
    %---
    [X, Y]=meshgrid(linspace(-D_r,D_r,n)); % define the computational grid
    M = zeros(n);    
    El = char_ellips(X, Y, mu_0, Sigma_0); % Function representing the characteristic ellipse
    % Fill the inside of the ellipse El
    M(El<=1) = 1; 
    El = char_ellips(X, Y, mu_0, 1.05^2*Sigma_0);
    % Fill the inside of the ellipse El
    M(El<=1) = 0.5*exp(-El(El<=1));
    %----
    % Ellipses representing some features of hte phantom
    mu_1 = [0.4;0];
    Sigma_1 = [7^2,0;0,4^2];
    %---
    mu_1 = R*mu_1; % Rotate ellipse    
    Sigma_1 = R*Sigma_1*R.';  % Rotate 
    %---
    El = char_ellips(X, Y, mu_1, Sigma_1);
    %---
    % Fill the inside of the ellipse El
    M(El<=1) = M(El<=1) - 0.2*exp(-El(El<=1));
    %---
    mu_1 = [-0.2;-0.2];
    Sigma_1 = [4^2,1;2,3^2];
    %---
    mu_1 = R*mu_1; % Rotate ellipse 
    Sigma_1 = R*Sigma_1*R.';
    %---
    El = char_ellips(X, Y, mu_1, Sigma_1);
    % Fill the inside of the ellipse El
    M(El<=1) = M(El<=1) + 0.4*exp(-El(El<=1));
end

% Function calculating the charasteristic polynomial of the ellipse.
function El = char_ellips(X,Y, mu, Sig) 
    El = 1/2*( (X-mu(1)).*Sig(1,1).*(X-mu(1)) + 2*(X-mu(1)).*Sig(1,2).*(Y-mu(2)) ...
        +(Y-mu(2)).*Sig(2,2).*(Y-mu(2)) );
end