# Showing Fourier slice theorem and Radon transform relation in MATLAB

I wrote some code to demonstrate the Fourier slice theorem and it's relation to the Radon transform. However the sampled FFT from the 2D FFT and the 1D FFT of the projection at the same angle don't match too closely.

The issue might come from two places:

Maybe the issue is undersampling and I should use zero padding for both 1 and 2D FFTs, how can I do using ifftshift(fft(fftshift(f)))? Is the way I'm using fftshift and ifftshift correct for images?

Maybe I should use another more advanced interpolation method for sampling the 2D FFT along a line? I have read it might be challenging getting a good interpolation.

Thank you,

Alex

The code is shown below:

nx     = 80;
xRange = -nx/2:nx/2-1;

dx    = 10;
angle = 20.0*pi/180.0;
sigma = 4.0;

A = makeMatrixWithLine(nx, dx,angle, sigma);

figure;
levels = 100;
contourf(xRange,xRange, A, levels, 'edgecolor','none');
colormap(flipud(bone));
xlim([-40 40])
ylim([-40 40])
pbaspect([1 1 1])
xticks(-40:10:40)
yticks(-40:10:40)
caxis([0 1])
c = colorbar;
c.Location = 'eastoutside';
c.Label.String = 'Signal magnitude';

FT_A = ifftshift(fft2(fftshift(A)));

deltak = 2*pi/nx;
kx     = (-nx/2:nx/2-1)*deltak;

figure;
contourf(kx,kx, real(FT_A), levels, 'edgecolor','none');
colormap(bone);
caxis([-400 400])
c = colorbar;
c.Location = 'eastoutside';
c.Label.String = 'Real component of FT of signal';

hold on
range = [-3.0 3.0];
plot(range, range*tan(angle), 'LineWidth',0.5, 'Color','red')

theta = 0:360;
figure
contourf(xp,theta, R', levels, 'edgecolor','none');
colormap(bone);

hold on
plot([-50.0 50.0],[angle*180/pi angle*180/pi], 'LineWidth',0.5, 'Color','red')
xlim([-50.0 50.0])

%Interpolate values
pointsInKx = xp(1:end-1)*deltak;

projectionAtAngle = R(:, round(angle*180/pi));
FTfromProjection  = ifftshift(fft(fftshift(projectionAtAngle(1:end-1))));

XI = pointsInKx*cos(angle);
YI = pointsInKx*sin(angle);
FTatAngle = interp2(kx,kx, FT_A, XI,YI);

figure
plot(pointsInKx, real(FTatAngle))
hold on
plot(pointsInKx, real(FTfromProjection))
ylim([-400 400])
legend('Projection from 2D FFT','1D FFT from Radon\newlinetransform at 20 degrees')

function A = makeMatrixWithLine(nx, dx,angle, sigma)
x     = -nx/2:nx/2-1;
[X,Y] =  meshgrid(x,x);

XY = [X(:) Y(:)];
R     = [cos(-angle) -sin(-angle);
sin(-angle)  cos(-angle)];
rotXY = XY*R';

Xqr = reshape(rotXY(:,1)-dx, size(X,1), []);
dist = abs(Xqr);
A    = exp(-dist.^2/sigma);
end