(see also plotting using 3D bars)
% simple image
P = zeros(11,11);
P(6, 6) = 1;
imshow(P,'InitialMagnification', 2000), title('original sample figure')
surf(P) colormap(hsv)
% sinogram
theta = 1:1:179;
[R,xp] = radon(P, theta);
imshow(R,[],'Xdata',theta,'Ydata',xp,...
'InitialMagnification','fit')
xlabel('\theta (degrees)')
ylabel('x''')
colormap(bone), colorbar
mesh(R) colormap(hsv)
I1 = iradon(R,theta,'linear','none');
I1 = I1 ./max(max(I1));
surf(I1)
colormap(hsv)
title('Unfiltered Backprojection')
I = iradon(R, theta);
I = I ./max(max(I));
surf(I)
colormap(hsv)
title('Filtered Backprojection')
I1 = iradon(R,theta,'linear','none');
I1 = I1 ./max(max(I1));
I = iradon(R, theta);
I = I ./max(max(I));
diff = I1 - I;
surf(diff)
colormap(hsv)
title('Difference: Unfiltered - Filtered') % <==================== difference
% simple image
P = zeros(11,11);
P(6, 6) = 1;
theta = 1:1:179;
[R,xp] = radon(P, theta);
I1 = iradon(R,theta,'linear','none');
I = iradon(R, theta);
% fft at center of figure
y = I(:,6);
[L,W] = size(y);
Fs = L;
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2);
% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2)))
title('Single-Sided Amplitude Spectrum - Filtered')
xlabel('Frequency (Hz)')
% fft at center of figure
y = I1(:,6);
[L,W] = size(y);
Fs = L;
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2);
% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2)))
title('Single-Sided Amplitude Spectrum - Unfiltered')
xlabel('Frequency (Hz)')
http://amath.colorado.edu/courses/4380/2004fall/FB.ppt#279,10,Why do we Apply a Filter ?
The Fourier Slice Theorem (ref)
The Fourier Slice Theorem states that the Fourier transform of a projection of a function f(x,y), seen from an angle θ, equals the slice of the Fourier transform of f(x,y), F(f(x,y)) = F(ωx, ωy), under that angle θ.
Under the projection we understand the function P(θ, t) which results from integrating along all parallel lines perpendicular to a line (x, y) = t * (cos(θ), sin(θ)) through the origin of the coordinate axes and forming the angle θ with the x-axis.
It is clear that if projections from all angles 0 <= θ < π are given, their Fourier transforms will completely cover the Fourier transform of the function f(x,y). Thus, the function f(x,y) can be determined by its projections when assembling the Fourier transforms of the projections correctly and then back-transform to the space domain. A standard application is the mapping of internal body tissue, where X-ray images as projections of, for instance, a tumor, are used to reconstruct a three dimensional image of the tumor itself.
http://en.wikipedia.org/wiki/Filtered_back_projection#Filtered_back-projection
Study of Radon Transformation and Application of its Inverse to NMR - http://cgi.di.uoa.gr/~erga/mobio/05/proj2/Flaounas_Venturas_Final.doc local
Discrete Radon transform has an exact, fast inverse and generalizes to operations other than sums along lines - http://www.nr.com/whp/radonpreprint.pdf local