function [image,XX,YY] = ctom( imSize, log_2_n, numPhi )
% Function to compute a tomograph of an object defined
%   by the function projection( phi, xpts )
%   P.N. Saeta, 2005-03-28
%
% Inputs
%   imSize = number of pixels in square image to generate
%   log_2_n = base-2 logarithm of pts in FFT
%             for 128-point transform, log_2_n = 7
%   numPhi = number of incident beam angles, which will
%     be equally spaced
% 
% Outputs
%   image = the 2-dimensional image made from the numPhi
%           projections
%   XX,YY = the x and y locations of the grid points of the image
   
projSize = 2^log_2_n;    % size of the projection
halfN = projSize / 2;    % half thereof
image = zeros( imSize ); % 2-dimensional image
xrange = 1.1;
xscale = inspace(-xrange, xrange, imSize ); % x coordinate of pixels

[XX,YY] = meshgrid( xscale, xscale );
dk = 1/(2 * xrange);
kscale = inspace( -0.5*dk*imSize, 0.5*dk*imSize, projSize );
projscale = inspace( -xrange, xrange, projSize );

% Now we need to loop over the projections

for i=1:numPhi
	phi = pi * i / numPhi;					% equally spaced angles
	proj = projection( phi, projscale );	% compute projection
	                                    	% using external function
	ft = fft( proj );						% use FFT to compute
	ft = ft .* abs(kscale);					%   modified projection
	rproj = real(ifft( ft ));				% revert to configuration
	                         				%   space
	% Interpolate and accumulate in image
	xiValue = cos(phi) * XX + sin(phi) * YY;
	contr = interp1( projscale, rproj, xiValue );
	image = image + contr;
    
    % display the image as it evolves
    imagesc( image );                       % draw and update image
    drawnow
end
	
function x = inspace( from, to, steps )
% Function inspace() produces a vector of equally spaced
%   points starting at 'from' and with step size
%   (to-from)/steps. Hence, the final point is one step
%   shy of 'to'.
  x = linspace( from, to - (to-from)/steps, steps );