% Script to create a software phantom, its parallel projections and its
% reconstruction using Filtered BackProjection
% Created by Edgar Garduño Angeles, 2008
%
clear all %Cleaning all variables
clc       % Clearing the main screen

IDx = 400;
IDy = 400;
L   = 400;
TI  = 0;
TO  = 179;
DA  = 1;
A = [TI:DA:TO];

[S,PI] = sinogram_shepp_logan(IDx,IDy,TI,TO,DA,L);

set(figure(1),'Position',[1 1 1000 1000])
imagesc(PI);
colormap(gray)
title('Original Phantom');

set(figure(2),'Position',[1000 1 1000 1000])
imagesc(S)
colormap(gray)
axis square
title('SinoGram')

% Filter
a = 1; % Hamming Filter, this value affects the signal frequency, increasing a increases the frequency
w =[-pi:(2*pi)/size(S,2):pi-(2*pi)/size(S,2)]; % Defines the filter range
% Components of the filter
r1 = abs(2/a*sin(a.*w./2));
r2 = sin(a.*w./2);
r3 = (a*w)./2;
r = r1*(r2/r3)^2;
f = fftshift(r); % Fourier Transform
set(figure(3),'Position',[1000 1 1000 1000])
plot(r,'b','LineWidth',2)
title('Filter')

%
% Method to create the Shepp-Logan Phantom
%
function [S,P] = sinogram_shepp_logan(N,M,AI,AO,DA,L)
LH = [-1 1];
LV = [-1 1];
D = [(LH(2) - LH(1))/N (LV(2) - LV(1))/M];
O = [(LH(2) + LH(1))/2 (LV(2) + LV(1))/2];
C = [abs(LH(1)+O(1))/D(1) abs(LV(1)+O(2))/D(2)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% List of Ellipses
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s.x=0; % Center coordinate x
s.y=0; % Center coordinate y
s.a=0; % Minor axis
s.b=0; % Major axis
s.d=0; % Density
s.t=0; % Angle
% Number of Ellipses
Ne=10;
E=repmat(s,1,Ne);
% First Ellipse
E(1).x= 0;
E(1).y= 0;
E(1).a= 0.69;
E(1).b= 0.92;
E(1).d= 2.0;
E(1).t= 0.0;
% Second Ellipse
E(2).x=  0.0;
E(2).y= -0.0184;
E(2).a=  0.6624;
E(2).b=  0.874;
E(2).d= -0.98;
E(2).t=  0;
% Third Ellipse
E(3).x= 0.22;
E(3).y= 0.0;
E(3).a= 0.11;
E(3).b= 0.31;
E(3).d=-0.4; % -0.02
E(3).t= -18;
% Fourth Ellipse
E(4).x= -0.22;
E(4).y=  0.0;
E(4).a= 0.16;
E(4).b= 0.41;
E(4).d=-0.4; % -0.02
E(4).t= 18;
% Fifth Ellipse
E(5).x= 0.0;
E(5).y= 0.35;
E(5).a= 0.21;
E(5).b= 0.25;
E(5).d= 0.2; % 0.01
E(5).t= 0.0;
% Sixth Ellipse
E(6).x= 0.0;
E(6).y= 0.1;
E(6).a= 0.046;
E(6).b= 0.046;
E(6).d= 0.2; % 0.01
E(6).t= 0.0;
% Seventh Ellipse
E(7).x= 0.0;
E(7).y=-0.1;
E(7).a= 0.046;
E(7).b= 0.046;
E(7).d= 0.2; % 0.01
E(7).t= 0.0;
% Eighth Ellipse
E(8).x=-0.08;
E(8).y=-0.605;
E(8).a= 0.046;
E(8).b= 0.023;
E(8).d= 0.2; % 0.01
E(8).t= 0.0;
% Ninth Ellipse
E(9).x= 0.0;
E(9).y=-0.605;
E(9).a= 0.023;
E(9).b= 0.023;
E(9).d= 0.2; % 0.01
E(9).t= 0.0;
% Tenth Ellipse
E(10).x= 0.06;
E(10).y=-0.605;
E(10).a= 0.023;
E(10).b= 0.046;
E(10).d= 0.2; % 0.01
E(10).t= 0.0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Creating Phantom on the Screen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
P=zeros(M,N);
for k=1:Ne
    radius = max([E(k).a E(k).b]);
    My=max([(fix((O(2) + E(k).y + radius)/D(2) + 0.5) + C(2)) M]);
    my=max([(fix((O(2) + E(k).y - radius)/D(2) + 0.5) + C(2)) 1]);
    for j=my:My
        y=(j-C(2))*D(2);
        a=1;
        b=-2*E(k).x;
        c=E(k).x*E(k).x - E(k).a*E(k).a*(1 - (y-E(k).y)*(y-E(k).y)/(E(k).b*E(k).b));
        d=b*b - 4*a*c;
        if d>=0
            x1=(-b - sqrt(d))/(2*a);
            x2=(-b + sqrt(d))/(2*a);
            Mx=floor(((O(1) + x2)/D(1) + 1) + C(1));
            mx=floor(((O(1) + x1)/D(1) + 1) + C(1));
            for i=mx:Mx
                P(M - j,i) = P(M - j,i) + E(k).d;
            end
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Creating Projections
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dt=DA;
Theta=AI:Dt:AO;
Nl=L;
h=sqrt((LH(2)-LH(1))*(LH(2)-LH(1)) + (LV(2)-LV(1))*(LV(2)-LV(1)));
Dl=h/(Nl+1);
MX=h/2;           % Half of the Projection
Cp=fix(MX/Dl+1);  % Center of the Projection in Pixels
S = zeros(size(Theta,2),Nl);

for t=1:size(Theta,2)
    r=Theta(t)*unitsratio('rad','deg');
    o=[cos(r) sin(r)];
    x=[-sin(r) cos(r)];
    R = zeros(Nl,1);
    for k=1:Ne
        if E(k).b > E(k).a
            p=O + [E(k).x E(k).y] + E(k).b*x;
            d2=ceil(dot(x,p)/Dl + Cp);
            p=O + [E(k).x E(k).y] - E(k).b*x;
            d1=floor(dot(x,p)/Dl + Cp);
        else
            p=O + [E(k).x E(k).y] + E(k).a*x;
            d2=ceil(dot(x,p)/Dl + Cp);
            p=O + [E(k).x E(k).y] - E(k).a*x;
            d1=floor(dot(x,p)/Dl + Cp);
        end
        for l=d1:d2
            %Checking whether this line intersects the k Ellipse
            [e,T]=IntersectsEllipse((l-Cp)*Dl*x,o,[E(k).x E(k).y],E(k).a,E(k).b);
            if e == true
                R(l) = R(l) + E(k).d*abs(T(2)-T(1));
            end
        end
    end
    S(t,:) = R;
end
end

% Subroutine to check whether a line intersects an ellipse.
% The ellipse is centered at point c (a vector) with radii w and h
% The line is determined by X and o by p=X+t*o
% The function returns three values, e to inform whether there was
% anintersection and a vector T with values t1 and t2 of intersection
% Created by Edgar Garduño Angeles, 2008
%
function [e,T]=IntersectsEllipse(X,o,c,w,h)

a=(o(1)*o(1))/(w*w) + (o(2)*o(2))/(h*h);
b=(2*X(1)*o(1) - 2*c(1)*o(1))/(w*w) + (2*X(2)*o(2) - 2*c(2)*o(2))/(h*h);
c=(X(1)*X(1) - 2*c(1)*X(1) + c(1)*c(1))/(w*w) + (X(2)*X(2) - 2*c(2)*X(2) + c(2)*c(2))/(h*h) - 1;

d=b*b - 4*a*c;
%fprintf(1,'Disc:%f\n',d);
if d<0
    e=false;
    T(1)=0;
    T(2)=0;
else
    e=true;
    T(2)=(-b+sqrt(d))/(2*a);
    T(1)=(-b-sqrt(d))/(2*a);
    if T(1)>T(2)
        T=[T(2) T(1)];
    end
    %fprintf(1,'%f,%f,Disc:%f\n',t1,t2,d);
end
end