% Program to plot development of a sinogram, simple back projection, and 
% spatially-filtered back projection of a cartoon fly 

clear
close all

vid = VideoWriter('SimpleVfilteredBP.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 60;
open(vid);

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]); 
colormap bone

phantomData = imread('cartoonFly.png');
phantomData = double(phantomData - min(phantomData(:))); % set min of image to zero
% pad the image with zeros so we don't lose anything when we rotate.
padDims = ceil(norm(size(phantomData)) - size(phantomData)) + 2;
P       = padarray(phantomData,padDims);

pos1 = [0.23 0.5 0.27 0.45]; % Original
pos2 = [0.5 0.5 0.27 0.45];  % Sinogram
pos3 = [0.23 0.0 0.27 0.45]; % Simple back projection
pos4 = [0.5 0.0 0.27 0.45];  % Filtered back projection

%__________________________________________________________________________

% Plot original image
subplot('position',pos1)
imagesc(P);
title('Original','FontSize',20)
axis equal
axis off
% Set angular step size
freq = 2; % [1/degree]
thetas = 0:1/freq:180-1/freq;

%__________________________________________________________________________

subplot('position',pos2)
% Compute sinogram
numOfAngularProjections  = length(thetas);
numOfParallelProjections = size(P,1);
sinogram = zeros(numOfParallelProjections,numOfAngularProjections);
% Loop over the number of angles
for i = 1:length(thetas)
    % Rotate image
    tmpImage = imrotate(P,-thetas(i),'bilinear','crop');
    
    % Fill sinogram
    sinogram(:,i) = sum(tmpImage,2);
    imagesc(sinogram);
    title('Sinogram','FontSize',20)
    axis square
    axis off
    drawnow
    
    frame = getframe(gcf);
    writeVideo(vid,frame);
end

%__________________________________________________________________________

subplot('position',pos3)
% Simple backprojection

numOfParallelProjections = size(sinogram,1);
numOfAngularProjections  = length(thetas);
% Convert thetas to radians
thetas = (pi/180)*thetas;
% Set up the backprojected image
BPI = zeros(numOfParallelProjections,numOfParallelProjections);
% Find the middle index of the projections
midindex = floor(numOfParallelProjections/2) + 1;
% Set up the coords of the image
[xCoords,yCoords] = meshgrid(ceil(-numOfParallelProjections/2):ceil(numOfParallelProjections/2-1));
% Loop over each projection angle
for i = 1:numOfAngularProjections
    % Figure out which projections to add to which spots
    rotCoords = round(midindex + xCoords*sin(thetas(i)) + yCoords*cos(thetas(i)));
    % Check which coords are in bounds
    indices   = find((rotCoords > 0) & (rotCoords <= numOfParallelProjections));
    newCoords = rotCoords(indices);
    
    % Summation
    BPI(indices) = BPI(indices) + sinogram(newCoords,i)./numOfAngularProjections;
    imagesc(BPI)
    title('Simple backprojection','FontSize',20)
    axis equal
    axis off
    drawnow
    
    frame = getframe(gcf);
    writeVideo(vid,frame);    
end

%__________________________________________________________________________

subplot('position',pos4)
% Reconstruction with filtered back projection in spatial domain

BPI = zeros(numOfParallelProjections,numOfParallelProjections);
% Find the middle index of the projections
midindex = floor(numOfParallelProjections/2) + 1;
% Set up the coords of the image
[xCoords,yCoords] = meshgrid(ceil(-numOfParallelProjections/2):ceil(numOfParallelProjections/2-1));
% Set up filter: now for the spatial domain!!!
if mod(numOfParallelProjections,2) == 0
    halfFilterSize = floor(1 + numOfParallelProjections);
else
    halfFilterSize = floor(numOfParallelProjections);
end
filter = zeros(1,halfFilterSize);
filter(1:2:halfFilterSize) = -1./((1:2:halfFilterSize).^2 * pi^2);
filter = [fliplr(filter) 1/4 filter]; 

% Loop over each projection
for i = 1:numOfAngularProjections
    % Figure out which projections to add to which spots
    rotCoords = round(midindex + xCoords*sin(thetas(i)) + yCoords*cos(thetas(i)));
    % Check which coords are in bounds
    indices   = find((rotCoords > 0) & (rotCoords <= numOfParallelProjections));
    newCoords = rotCoords(indices);
    
    % Filter
    filteredProfile = conv(sinogram(:,i),filter,'same');
    % Summation
    BPI(indices) = BPI(indices) + filteredProfile(newCoords)./numOfAngularProjections;
    imagesc(BPI)
    title('Ram-Lak filtered backprojection','FontSize',20)
    axis square
    axis off
    drawnow
    
    frame = getframe(gcf);
    writeVideo(vid,frame);
end
close(vid);