Algorithms

iradon.m

%% Direct reconstruction

I = iradon(proj, 1);

subplot(211)
imshow(I)
title('Reconstructed `iradon`');

%% Manually implement

% Construct the Fourier filter
m = size(proj, 1);
f = [1: (m - 1) / 2 + 1, (- m / 2 : -1)] / m;  % digital frequency
omega = 2 * pi * f;                            % angular frequency
temp = 2 * abs(f);                             % shepp-logan filter  
temp(2: end) = temp(2: end) .* sin(omega(2: end)) ./ omega(2: end);      
fourier_filter = temp;

% Apply filter in Fourier domain
projection = fft(proj) .* fourier_filter';
radon_filtered = real(ifft(projection));

% Reconstruct image by interpolation
n = size(orig, 1);
reconstructed = zeros(n);
[X, Y] = meshgrid(1: n);
xpr = X - n / 2;
ypr = Y - n / 2;

for i = 1: numel(theta)
    t = ypr * cos(deg2rad(theta(i))) - xpr * sin(deg2rad(theta(i)));
    x = (1: m) - ceil(m / 2);
    backprojected = interp1(x, radon_filtered(:, i), t);
    reconstructed = reconstructed + backprojected;
end

%% Show results

res = reconstructed * pi / (2 * numel(theta));
subplot(212)
imshow(res')
title('Reconstructed manually');