function imBP = hyperbackproj(hRT, ss, thetas, n)

% initialize backprojected image
imBP = zeros(n);

% define the reconstruction grid
xs = -1+1/n:2/n:1-1/n;
x = repmat(xs, n, 1);
y = repmat(xs', 1, n);

ds = ss(2)-ss(1);

for j=1:n,
    for k=1:n,
        if x(j,k)*x(j,k)+y(j,k)*y(j,k)<=.81,
            for t=1:length(thetas),
                % find s
                s = hyperS(x(j,k),y(j,k),thetas(t));
                % add to the reconstruction the contribution of that point
                if s>=ss(1) && s<=ss(end)
                    s1 = (s-ss(1))/ds + 1;
                    fs1 = floor(s1);
                    temp = (fs1+1-s1)*hRT(fs1,t) + (s1-fs1)*hRT(fs1+1,t);
                    imBP(j,k) = imBP(j,k) + temp;
                end
%                 imBP(j,k) = imBP(j,k) + interp1(ss,hRT(:,t),s);
            end
        end
    end
end

imBP = imBP*2*pi/length(thetas);

end