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