function [imRc,H] = myIRadon(RT, filter, thetas, d, n)
%MYIRADON my Inverse Radon Transform
%   very inspired from the MatLab function iradon, except that the size of
%   input RT has is fixed wrt n, and that the reconstructed image is of
%   even size n and the reconstruction grid [-1,1]^2 is such that the 
%   center of the image is at the center of the four central pixels.
%   For now, the interpolation is linear only.
%   INPUT ARGS: RT     -  Radon data to invert for
%               filter - type of filter: 'ram-lak', 'shepp-logan', 
%               'cosine', 'hann', 'hamming'
%               thetas - angles corresponding to the projections
%               d      - in [0,1], cropping normalized frequency.
%               n      - size of the reconstructed square image
%   OUTPUT ARGS: imRc  - reconstructed image
%                H     - filter used in frequency domain
%   See also MYRADON, ROTATEIMAGE

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

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

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

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

y = (1-1/n:-2/n:-1+1/n).';
y = repmat(y, 1, n);

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

nr = size(RT,1); 

%--- Backprojection  - vectorized in (x,y), looping over theta
for i=1:length(thetas)
    proj = RT(:,i);
    t = x.*costhetas(i) + y.*sinthetas(i);
    t = (t-(-nr/n+1/n))*n/2+1;
    a = floor(t);
    
    imRc = imRc + (t-a).*proj(a+1) + (a+1-t).*proj(a); 
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,H] = filterProjections(p_in, filter, d)

p = p_in;

% Design the filter
len = size(p,1);   
H = designFilter(filter, len, d);

if strcmpi(filter, 'none')
    return;
end

p(length(H),1)=0;  % Zero pad projections 

p = fft(p);    % p holds fft of projections

for i = 1:size(p,2)
   p(:,i) = p(:,i).*H; % frequency domain filtering
end

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

end
%----------------------------------------------------------------------

%======================================================================
function filt = designFilter(filter, len, d)
% Returns the Fourier Transform of the filter which will be 
% used to filter the projections
%
% INPUT ARGS:   filter - either the string specifying the filter 
%               len    - the length of the projections
%               d      - the fraction of frequencies below the nyquist
%                        which we want to pass
%
% OUTPUT ARGS:  filt   - the filter to use on the projections


order = max(64,2^nextpow2(2*len));

if strcmpi(filter, 'none')
    filt = ones(1, order);
    return;
end

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

filt = 2*( 0:(order/2) )./order;
w = 2*pi*(0:size(filt,2)-1)/order;   % frequency axis up to Nyquist 

switch filter
    case 'ram-lak'
        % Do nothing
    case 'shepp-logan'
        % be careful not to divide by 0:
        filt(2:end) = filt(2:end) .* (sin(w(2:end)/(2*d))./(w(2:end)/(2*d)));
    case 'cosine'
        filt(2:end) = filt(2:end) .* cos(w(2:end)/(2*d));
    case 'hamming'  
        filt(2:end) = filt(2:end) .* (.54 + .46 * cos(w(2:end)/d));
    case 'hann'
        filt(2:end) = filt(2:end) .*(1+cos(w(2:end)./d)) / 2;
    otherwise
        eid = sprintf('Images:%s:invalidFilter',mfilename);
        msg = 'Invalid filter selected.';
        error(eid,'%s',msg);
end

filt(w>pi*d) = 0;                      % Crop the frequency response
filt = [filt' ; filt(end-1:-1:2)'];    % Symmetry of the filter

end
%----------------------------------------------------------------------