function hRT = hyperRadon(f, ss, thetas, npts)

n = size(f,1);
xs = -1+1/n:2/n:1-1/n;

nss = length(ss);
nth = length(thetas);
hRT = zeros(nss,nth);

%--- version that uses the slow matlab interp2 function
% xx = repmat(xs, n, 1);
% yy = repmat(xs', 1, n);
% for j=1:nss,
%     for k=1:nth,
%         [x,y,r0,alpha0] = geodesic(ss(j),thetas(k),npts);
%         vals = interp2(xx,yy,f,x(x.^2+y.^2<.81),y(x.^2+y.^2<.81));
%         
%         if any(isnan(vals)),
%             x( isnan(vals) )
%             y( isnan(vals) )
%             error('not a number !?!?!?!?!?')
%         end
%         
%         hRT(j,k) = (2*alpha0*r0)*sum(vals)/npts;
%     end
% end

%--- faster hardcoded version
for j=1:nss,
    for k=1:nth,
        [x,y,r0,alpha0] = geodesic(ss(j),thetas(k),npts);
        
        for p=1:length(x),
            if x(p)*x(p) + y(p)*y(p) <= .81,
                x_rs = n/2. * (x(p) + 1.-1./n);
                y_rs = n/2. * (y(p) + 1.-1./n);
                i0 = floor(y_rs);
                j0 = floor(x_rs);
                xbar = x_rs - j0;
                ybar = y_rs - i0;
                i0 = i0+1; j0=j0+1;
                hRT(j,k) = hRT(j,k) ...
                    + f(i0, j0) * (1.-xbar) * (1.-ybar) ...
                    + f(i0+1, j0) * ybar * (1.-xbar) ...
                    + f(i0, j0+1) * xbar * (1.-ybar) ...
                    + f(i0+1, j0+1) * xbar * ybar;
            end
        end
        hRT(j,k) = (2*alpha0*r0)*hRT(j,k)/npts;
    end
end

end