%======================================================================
function [p_out,H] = filterProjections(p_in, filter)

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

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


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

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

    p = bsxfun(@times, p, H); % faster than for-loop

    p = ifft(p,'symmetric');  % p is the filtered projections

    p(len+1:end,:) = [];      % Truncate the filtered projections
    p_out = p;
    
elseif ndims(p)==3
    
    p_out = zeros(length(H),size(p,2),size(p,3));
    
    for nz=1:size(p,3)
        p_z = p(:,:,nz);
        
        p_z(length(H),1)=0;
        
        p_z = fft(p_z);
        
        p_z = bsxfun(@times, p_z, H);
        p_z = ifft(p_z,'symmetric');

        p_out(:,:,nz) = p_z;
    end
    
    p_out(len+1:end,:,:) = [];
    
    
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 bandlimited ramp filter (Eqn. 61 Chapter 3, Kak and
% Slaney) - go up to the next highest power of 2.

n = 0:(order/2); % 'order' is always even. 
filtImpResp = zeros(1,(order/2)+1); % 'filtImpResp' is the bandlimited ramp's impulse response (values for even n are 0)
filtImpResp(1) = 1/4; % Set the DC term 
filtImpResp(2:2:end) = -1./((pi*n(2:2:end)).^2); % Set the values for odd n
filtImpResp = [filtImpResp filtImpResp(end-1:-1:2)]; 
filt = 2*real(fft(filtImpResp)); 
filt = filt(1:(order/2)+1);

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;
    case 'superlinear'
        filt(2:end) = filt(2:end).^3;

    otherwise
        error(message('images:iradon:invalidFilter'))
end

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