function [ im ] = bp(fpm,n )
%
% Very simple backprojection routine, no rho filter
% 
%   im = bp(fpm,n )
%
% fpm -- projection data in columns, of even length
% n   -- size of the reconstructed image
% im  -- backprojection, n x n

% lp is length of the projections
% np is the number of projections
[lp np] = size(fpm);

% make sure the projections are of even length to make things easier
if floor(lp/2)*2 ~= lp,
    display('Projections must be of even length');
    return;
end

% 2D indexing array to be used below, it is just nx + i*ny
y = [-n/2:n/2-1]/(n/2); 
x = y'; 
xy = x*ones(1,n)+i*ones(n,1)*y; 


% zero out the images we are backprojecting into  
im = zeros(n,n);

% backproject each projection
for jj=1:np, 
    
    % assume that the angles are uniform over pi radians
    proj_angle = (jj-1)*pi/np;

    % find the locations in the current projection for each point
    % in the image at the current angle
    xoffset = real(xy*exp(i*proj_angle));
    ndx = round(xoffset*lp/2)+lp/2 + 1;

    % Clip the indices to the ends of the projection
    ndx = min(ndx, lp);
    ndx = max(1, ndx);
  
    % The current projection
    fp = fpm(:,jj); 
    % zero out the end samples, so they don't do anything
    fp(1) = 0; 
    fp(lp) = 0;
    
    % Backrpoject!  This line does everything
    im = im+fp(ndx); 
end

end

