function lamest = backproject(sinodat, Mtheta, Mr, filt)
%Usage: lamest = backproject(sinodat, Mtheta, Mr, filt)
%Given the sinogram "sinodat" computed as sino(theta,r) with 0 <= theta < pi, 
%-1 < r < 1, backproject on a polar gridof size Mtheta by Mr. If filt = 0 
%backprojection is unfiltered, otherwise a standard ramp filter is applied.

[m,n] = size(sinodat);

%if specfied, apply a ramp filter using the fft command.
if(filt)
    for j = 1:m
        q = fft(sinodat(j,:));
        q = q.*(1-1*abs((1:n)-n/2-1)/(n/2));
        sinodat(j,:) = real(ifft(q));
    end
end

%Set up a circular grid for computing attenuation estimates.
rgrid = ((0:Mr-1)+0.5)/Mr;
tgrid = 2*pi*(0:Mtheta-1)/Mtheta;

tvals = pi*(0:m-1)/m; %Angular values in sinodat

%And set up an array to hold the recontruction
lamest = zeros(Mtheta,Mr);

%A little throw away code to use fewer than the specified number of rays.
ms = m; %Setting ms = m uses all the data.
m2 = round(m/ms);

%Loop over grid points
for j = 1:Mr
    for k = 1:Mtheta
        s = 0.0;
        for q = 1:m2:m %Loop over angular data in sinodat
            r = rgrid(j)*cos(tgrid(k)-tvals(q));
            rk = ceil(n/2*(r+1));
            s = s + sinodat(q,rk);
        end
        lamest(k,j) = s/m; %Straight back projection
    end
end
lamest = lamest/2/pi;
