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

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

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))),[]);

frame = getframe(gcf);

writeVideo(vid,frame);

end

end

end

% Output the movie as an mpg file

close(vid);