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

% 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);

nth = length(thetas);
ns = length(ss); 

ds = ss(2)-ss(1); % s increment

for j=1:n,
    for k=1:n,
        if x(j,k)*x(j,k)+y(j,k)*y(j,k)<1,
            %%--- non-vectorized
            
%             for t=1:nth,
%                 % find s
%                 s = hyperS(x(j,k),y(j,k),thetas(t), R);
%                 % 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
%             

            %%--- vectorized
            % find s
            s = max(min(hyperS(x(j,k),y(j,k),thetas, R),ss(end)),ss(1));
            % add to the reconstruction the contribution of that point
            s1 = (s-ss(1))/ds + 1;
            fs1 = max(min(floor(s1),ns-1),1);
            for t = 1:nth,
                temp = (fs1(t)+1-s1(t))*hRT(fs1(t),t) + (s1(t)-fs1(t))*hRT(fs1(t)+1,t);
                imBP(j,k) = imBP(j,k) + temp;
            end
            
        end
    end
end

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

end