function imRc = UWIRadon(RT, thetas, d, n)

%--- filter each projection
RT = filterProjections(RT, d);

%--- allocate memory for the image
imRc = zeros(n); 

% define the reconstruction grid
x = -1+1/n:2/n:1-1/n;
[x,y] = meshgrid(x,x);

costhetas = cos(thetas); sinthetas = sin(thetas);

nr = size(RT,1); 
%--- Backprojection  - vectorized in (x,y), looping over theta
for i=1:length(thetas),
    for j=1:n,
        for k=1:n,
            xp = x(j,k);
            yp = y(j,k);
            
            t = xp.*costhetas(i) + yp.*sinthetas(i);
            t = (t-(-nr/n+1/n))*n/2+1;
            a = floor(t);
            
            %% a designates between which two indices t the value is
            %% located
            % update imRc(j,k)
            imRc(j,k) = imRc(j,k) + %%% COMPLETE
        end
    end
end

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

%--- take out the part outside the unit disk
imRc = (n/2)*imRc.*double(x.^2+y.^2<=1);

end

%======================================================================
function p = filterProjections(p_in, d)

p = p_in;

%---- Design the filter
len = size(p,1);   

% make its length a power of 2 for computing purposes
order = max(64,2^nextpow2(2*len));

%---- First create a ramp filter - go up to the next highest power of 2

%% COMPLETE
% create a ramp filter of half-size, crop the "normalized frequencies" that
% are greater to d, and then symmetrize the filter

filt = ;

%---- filter the projections
p(length(filt),1)=0;  % Zero pad projections 

% frequency domain filtering on each column

for j=1:size(p,2),
    % filter the j-th column by applying the filter in Fourier
    %%% COMPLETE
end

%---- Truncate the filtered projections
p(len+1:end,:) = [];   

end
