%CODE TO COMPUTE PROJECTIONS OF A SHEPP-LOGAN PHANTOM

% THIS CODE CALLS THE ROUTINE projection2 
% INPUT IS THE MATRIX image
% OUTPUT IS THE MATRIX PR WHOSE COLUMNS ARE PROJECTIONS AT EACH ANGLE 
% Create a phantom of size K*K
iptsetpref('ImshowAxesVisible','on')
K=128;
 image = phantom(K);
         figure(1), imshow(image)
         title('Shepp-Logan phantom')
% Set the number of rotation angles M 
M=128;
theta=1:180/M:180;
PR = projection2(image,theta);
L=ceil(sqrt(2) * K);t=1:L;
figure(2), imshow(PR,[],'Xdata',90-theta,'Ydata',t-ceil(L/2),'InitialMagnification','fit')
 xlabel('\theta (degrees)')
 ylabel('t')

% THIS CODE CALLS THE ROUTINE radon from MATLAB 
% INPUT IS THE MATRIX image
% OUTPUT IS THE MATRIX [R,xp] 
% Create a phantom of size K*K
iptsetpref('ImshowAxesVisible','on')
K=128;
 image = phantom(K);
         figure(1), imshow(image)
         title('Shepp-Logan phantom')
% Set the number of rotation angles M 
M=128;
theta=1:180/M:180;
[R,xp] = radon(image,theta);
figure(3), imshow(R,[],'Xdata',theta,'Ydata',xp,'InitialMagnification','fit')
 xlabel('\theta (degrees)')
 ylabel('t')

%% This MATLAB function takes an image matrix and vector of angles and 
%% computes the 1D projection (Radon transform) at each of the angles.  
%% Output is a matrix whose columns are the projections at each angle.

function PR = projection2(IMG,THETA)

% pad the image with zeros so we don't lose anything when we rotate.
[iLength, iWidth] = size(IMG);
iDiag = sqrt(iLength^2 + iWidth^2);
LengthPad = ceil(iDiag - iLength) + 2;
WidthPad = ceil(iDiag - iWidth) + 2;
padIMG = zeros(iLength+LengthPad, iWidth+WidthPad);
padIMG(ceil(LengthPad/2):(ceil(LengthPad/2)+iLength-1), ...
       ceil(WidthPad/2):(ceil(WidthPad/2)+iWidth-1)) = IMG;

% loop over the number of angles, rotate 90-theta (because we can easily sum
% if we look at data from the top), and then add up.
n = length(THETA);
PR = zeros(size(padIMG,2), n);
for i = 1:n
   tmpimg = imrotate(padIMG, 90-THETA(i), 'bilinear', 'crop');
   PR(:,i) = (sum(tmpimg))';
end