function RT = atRt(a,f,thetas,rotMethod)
%ATRT Attenuated Radon Transform
% computes the Radon Transform of f weighted by the attenuation term a
% arguments:
% - a: attenuation map
% - f: source map
% - thetas: angles of computation
% - rotMethod: image rotation method, '4pt' or 'spectral'

% Author: Francois Monard
% date: 06/24/2011

n = size(a,1); h = 2/n;
nth = size(thetas,2);
RT = zeros(n, nth);

for th=1:nth,
    % compute Da(.,theta) over the rotated grid
    Darot = fliplr( exp(-h*cumsum( fliplr( rotateImage(a,thetas(th)+pi/2,rotMethod) ),2) ) );
    % compute R_a f(., theta)
    RT(:,th) = h*sum(Darot.*rotateImage(f,thetas(th)+pi/2,rotMethod),2);
end

% zero-padding
zp = ceil( (sqrt(2)-1)*n/2 );
RT = [zeros(zp,nth); RT; zeros(zp,nth)];

end