function [fout,H] = hilbTrans(fin, filter, d)
%HILBTRANS regularized Hilbert transform
%    computes the filtered hilbert transform of the columns of fin. 
%    The filters available are those from the filtered back-projection
%    function.
%    INPUT ARGS:  fin    - the columns to be filtered
%                 filter - possible values : 'ram-lak', 'shepp-logan',
%                 'cosine', 'hann', 'hamming'
%                 d      - the cropping normalized frequency in [0,1]
%
%    OUTPUT ARGS: fout   - the filtered slices
%                 H      - the DFT of the filter applied to each slice
%    See also MYIRADON

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

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

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

fin = [fin; zeros(length(H)-len, size(fin,2))]; % Zero pad

% filter columns of fin in Fourier domain
fout = zeros(size(fin));
for j=1:size(fin,2),
    fout(:,j) = real(ifft(H.*fft(fin(:,j))));
end
fout = fout(1:len,:); % Truncate to the initial size

end

%======================================================================
function filt = designFilter(len, filter, d)
% Returns the Fourier Transform of the filter for the regularized Hilbert
% Transform
%
% 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


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

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

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

filt = ones(1,order/2); 
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)']/i;    % anti-symmetrize the filter

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