% Movie of reconstructing a tomogram through the Fourier slice theorem
% rather than filtered back projection (i.e., the iradon function in
% matlab). The core of the program draws heavily from the code found at
% https://stackoverflow.com/questions/61260808/implementing-a-filtered-
% backprojection-algorithm-using-the-central-slice-theorem.
% Acknowledgements go to 2fly2try, Chris Luengo, and
% Person.Woman.Man.Camera.TV (LOL)

clear all; close all
vid = VideoWriter('FourierSliceTheorem.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 15;
open(vid);
figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);

xyR = 9/16;
pos1 = [0.15*xyR 0.8 0.4*xyR 0.2];    % Subplot position for incident x-rays
pos2 = [0.15*xyR 0.4 0.4*xyR 0.4];    % Subplot position for rotating phantom
pos3 = [0.15*xyR 0.36 0.4*xyR 0.03];  % Subplot position for sinogram row
pos4 = [0.15*xyR 0.29 0.4*xyR 0.06];  % Subplot position for FT arrow
pos5 = [0.15*xyR 0.0 0.4*xyR 0.28];   % Subplot position for FT plot
pos6 = [0.64*xyR 0.05 0.9*xyR 0.9];   % Subplot position for polar FT plot 

myGold = [0.9 0.77 0];
myBlue = [0.4 0.44 0.73];

xRayArrowX = [-1  1 1 0 -1 -1];
xRayArrowY = [1 1 0.5 0 0.5 1];

% Rotation angles 0, 1, 2... 179
angs = 0:179;

initIm = imread('fstPhantom.png'); % Taken from from the public domain
% wikipedia source https://commons.wikimedia.org/wiki/Scrollable_computed_
% tomography_images_of_a_normal_brain_(case_2)

sino = radon(initIm,angs); % Create a sinogram using radon transform

% This step is necessary, so that frequency information is preserved: pad
% the sinogram with zeros so that this is ensured
sino = padarray(sino,floor(size(sino,1)/2),'both');

% Rotate the sinogram 90 degrees to be compatible with code's definition of
% view angle and radial positions
% dimension 1 -> view angle
% dimension 2 -> radial position
sino = sino';

% diagDim = length of an individual row of the sinogram
diagDim = size(sino,2);
ftTot = zeros(180,3,diagDim);   % For each angle, the x, y, and intensity 
                                % values in Fourier space. Initially zeros 
c = colormap(jet(diagDim));

% Initialize 2D Fourier transform of our final image as a matrix of zeros
fimg = zeros(diagDim);

% 1D V-shaped ramp filter with zero value at zero frequency. This
% suppresses low-frequency components and passes high-frequency components,
% thus reducing blurring due to the overlapping of the Fourier-transformed
% images in the low frequency region
rampFilter_1d  = abs(linspace(-1, 1, diagDim))';

rowIndex = 1;
for nn = angs
    % ******************** Maths part **********************
    imContrib = zeros(diagDim);
    
    % Get the current row of the sinogram - use rowIndex
    curRow = sino(rowIndex,:);
    
    % Take the 1D-FT of the current row - be careful; one needs to perform
    % first ifftshift and fftshift, as Matlab tends to place zero-frequency
    % components of a spectrum at the edges
    fourierCurRow = fftshift(fft(ifftshift(curRow)));
    
    % Take the Fourier-transformed sino row and place it at the center of
    % the next image contribution. Multiply by the ramp filter in the
    % Fourier domain. imContrib is a temporary parking place for the
    % contribution currently being worked on
    imContrib(round(diagDim/2),:) = fourierCurRow;
    imContrib(round(diagDim/2),:) = imContrib(round(diagDim/2), :)'...
        .* rampFilter_1d;
    
    % Rotate the current image contribution to be at the correct angle in
    % the 2D Fourier-domain image.
    imContrib = imrotate(imContrib, nn, 'crop');
    
    % Add the current image contribution to the running representation of
    % the image in Fourier space
    fimg = fimg + imContrib;
    
    % **************** Graphics and video creating part ******************

    % Top top left: draw x-ray arrow pointing downwards     
    subplot('position',pos1);
    newplot
    fill(xRayArrowX,xRayArrowY,myGold,'LineStyle','none');
    text(0, 0.6,'x-rays','FontSize',28','HorizontalAlignment','center');
    axis off
    axis equal
    
    % Top left: draw rotating original object being irradiated from above 
    subplot('position',pos2);
    newplot
    hold on
    rotIm = imrotate(initIm,-nn,'crop');
    imshow(rotIm,[]);
    axis equal
    
    % Middle left: draw current sinogram row    
    subplot('position',pos3);
    newplot
    hold on
    sst = round(diagDim/2)-size(initIm,1)/2;
    sen = round(diagDim/2)+size(initIm,1)/2;
    imshow(curRow(sst:sen),[],'XData', [-10 10], 'YData', [0 0.2]);
    axis equal
    
    % Bottom left: draw arrow indicating transform to FT    
    subplot('position',pos4);
    newplot
    xRayArrowX = [-1  1 1 0 -1 -1];
    xRayArrowY = [1 1 0.5 0 0.5 1];
    fill(xRayArrowX,xRayArrowY,myBlue,'LineStyle','none');
    text(0,0.6,'FT(Tr)','Color','w','FontSize',20',...
        'HorizontalAlignment','center','VerticalAlignment','middle');    
    axis off
    axis equal
    
    % Bottom bottom left: FT plot (log scale)    
    subplot('position',pos5);
    newplot
    semilogy(abs(fourierCurRow),'color','r', 'LineWidth', 2.0)
    axis off
    axis square
    
    % Right: buildup of Fourier components in Fourier space
    subplot('position',pos6);
    newplot
    FTx = -floor(diagDim/2):floor(diagDim/2);
    ftIm = scatter3(FTx*cosd(-nn),FTx*sind(-nn),log(abs(fourierCurRow)),2,...
        log(abs(fourierCurRow)),'filled');
    ftTot(nn+1,1,:) = FTx*cosd(-nn);
    ftTot(nn+1,2,:) = FTx*sind(-nn);
    ftTot(nn+1,3,:) = log(abs(fourierCurRow));
    xxx = reshape(ftTot(:,1,:),180*diagDim,1);
    yyy = reshape(ftTot(:,2,:),180*diagDim,1);
    zzz = reshape(ftTot(:,3,:),180*diagDim,1);
    
    hold on
    caxis([min(min(log(abs(fimg)))),max(max(log(abs(fimg))))]);
    axis off
    axis equal
    view(0,90); 
    
    % Move to next row in sinogram
    rowIndex = rowIndex + 1;

    frame = getframe(gcf);
    writeVideo(vid,frame);
    
    if (nn == 179)
         for fade = 1:-0.05:0
            hold off
            ftIm = scatter3(xxx,yyy,zzz,2,zzz,'filled');
            caxis([min(min(log(abs(fimg)))),max(max(log(abs(fimg))))]);
            axis off
            axis equal
            view(0,90);           
            set(ftIm, 'MarkerEdgeAlpha',fade, 'MarkerFaceAlpha',fade);
            frame = getframe(gcf);
            writeVideo(vid,frame);
        end
        hold off
        % Finally, take the inverse 2D Fourier transform of the image.
        % Remember to implement the fftshift or ifftshift here
        for fade = 0:0.05:1
            rcon = fftshift(ifft2(ifftshift(fimg)));
            finalIm = imshow(abs(rcon(round(0.25*diagDim):round(0.75*diagDim),...
                round(0.25*diagDim):round(0.75*diagDim))),[]);
            alpha(finalIm,fade);
            frame = getframe(gcf);
            writeVideo(vid,frame);
        end
    end
end
% Output the movie as an mpg file
close(vid);