1
$\begingroup$

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;
[R,xp] = radon(A,-theta);
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
$\endgroup$

1 Answer 1

0
$\begingroup$

I figured this out; the first issue was solved by zero padding both the 2D and 1D FFTs; this makes the transforms approach the Fourier transform. This produces the expected Fourier transform of a shifted Gaussian, both results have a better shape without the sampling related features seen in output of the current code.

However there was still a linear phase difference in the signal even with the line angle being zero; this discounts any issues with sampling a 2D FFT with an angled line. The error came down to the Radon transform moving the peak of the projection one pixel; the line at the x axis was a distance of 10 pixels away while the Radon transform at the x axis had a peak at 11 pixels. Shifting the projection one pixel when doing the zero padding produces a much closer Fourier transform to the 2D FT sampled at a line with the corresponding angle

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.