Cartoon of the Fourier slice theorem 

Matlab code 

% 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);