backpro.m

function img=backpro(sino,theta)

    
    %theta=theta+pi;  probably wrong !
	%interp='linear';
    interp='nearest neighbor';
    
	N=size(sino,1);


    img = zeros(N);        % Allocate memory for the image.

	% Define the x & y axes for the reconstructed image so that the origin
	% (center) is in the spot which RADON would choose.
	xax = (1:N)-ceil(N/2);
	x = repmat(xax, N, 1);    % x coordinates, the y coordinates are rot90(x)
	y = rot90(x);

	costheta = cos(theta);
	sintheta = sin(theta);
	ctrIdx = ceil(N/2);     % index of the center of the projections

	% Zero pad the projections to size 1+2*ceil(N/sqrt(2)) if this
	% quantity is greater than the length of the projections
	imgDiag = 2*ceil(N/sqrt(2))+1;  % largest distance through image.
	if size(sino,1) < imgDiag
		rz = imgDiag - size(sino,1);  % how many rows of zeros
		sino = [zeros(ceil(rz/2),size(sino,2)); sino; zeros(floor(rz/2),size(sino,2))];
		ctrIdx = ctrIdx+ceil(rz/2);
	end

	% Backprojection - vectorized in (x,y), looping over theta
	if strcmp(interp, 'nearest neighbor')
		for i=1:length(theta)
			proj = sino(:,i);
			t = round(x*costheta(i) + y*sintheta(i));
			img = img + proj(t+ctrIdx);
		end
	elseif strcmp(interp, 'linear')
		for i=1:length(theta)
			proj = sino(:,i);
			t = x.*costheta(i) + y.*sintheta(i);
			a = floor(t);
			img = img + (t-a).*proj(a+1+ctrIdx) + (a+1-t).*proj(a+ctrIdx);
		end
	elseif strcmp(interp, 'spline')
		for i=1:length(theta)
			proj = sino(:,i);
			taxis = (1:size(sino,1)) - ctrIdx;
			t = x.*costheta(i) + y.*sintheta(i);
			projContrib = interp1(taxis,proj,t(:),'*spline');
			img = img + reshape(projContrib,N,N);
		end
	end