clear all; close all;
load('p2_proj.mat');

r = linspace(-1,1,512);

subplot(121);
imagesc(1:180,r,P);
colormap bone;
xlabel('Projection angle (degrees)');
ylabel('Projection');
title(source);

[X Y] = meshgrid(r);
F = zeros(size(X));

rf = abs(linspace(-1,1,512));
rf = fftshift(rf);

for i = 1:180
    % Filtered backprojection
    p = P(:,i)';
    pf = rf.*fft(p);
    p = real(ifft(pf));
    
    R = X*cosd(i) + Y*sind(i);
    f = interp1(r,p,R,'linear',0);
    F = F + f;
end

subplot(122);
surf(X,Y,F);
shading interp;
colormap bone;
view([0 90]);
xlabel('x (normalized)');
ylabel('y (normalized)');
title('Simple backprojection');