%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Yue M. Lu and Minh N. Do
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%	PrySDrec_onestep.m
%	
%   First Created: 10-11-05
%	Last Revision: 07-13-09
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function rec = PrySDrec_onestep(Lp, Hp, w, tbw, D, smooth_func)

%   N-dimensional multiscale pyramid reconstruction - One Step
%
%   INPUT:
%
%     Lp: the lowpass subband in the Fourier domain with the symmetric part
%         removed along the last dimension.
%
%     Hp: the highpass subband in the Fourier domain with the symmetric part
%         removed along the last dimension.
%  
%     w: The ideal cutoff frequency of the lowpass subband. w is specified as a
%        fraction of PI, e.g. 0.5, 1/3, ....
%
%     tbw: Width of the transition band, specified as a fraction of PI.
%
%     D: downsampling factor. For example, 1, 1.5, 2, ......
%        Note: D should divide the size of the FFT array.
%
%     smooth_func: function handle to generate the smooth passband theta
%                  function
%
%   OUTPUT:
%
%     rec: the reconstructed signal in the Fourier domain. Note: Since the 
%          signal is real-vaued in the spatial domain, we remove the 
%          symmetric part of rec along the last dimension to save 
%          computation time. This is important when rec is a large 3-D 
%          array. See ccsym.m for details.
%
%   See also:
%
%   PyrSDdec_onestep.m, PyrNDDec_mm.m, PyrNDRec_mm.m


%% The dimension of the problem
N = ndims(Hp);
szX = size(Hp);

%% the original full size
szF = szX;
szF(end) = (szX(end) - 1) * 2;

%% Passband index arrays
pbd_array = cell(N, 1);

%% Lowpass filter (nonzero part)
szLf = 2 * ceil(szF / 2 * (w + tbw)) + 1;
szLf(end) = (szLf(end) + 1) / 2;
Lf = ones(szLf);

szR = ones(1, N); %% for resizing

for n = 1 : N
    %% current Fourier domain resolution
    nr = szF(n);
    
    szR(n) = szLf(n);
    
    %% passband
    pbd = [0 : ceil(nr / 2 * (w + tbw))] .';
    
    %% passband value
    pbd_value = realsqrt(smooth_func((pbd(end) - pbd) ./ (ceil(nr / 2 * (w + tbw)) ...
        - floor(nr / 2 * (w - tbw)))));
    
    %% See if we need to consider the symmetric part
    if n ~= N
        pbd = [pbd ; [nr - pbd(end) : nr - 1].'];
        pbd_value = [pbd_value ; flipud(pbd_value(2:end))];
    end
    
    pbd_array{n} = pbd + 1;
    pbd_value = reshape(pbd_value, szR);
    Lf = Lf .* repmat(pbd_value, szLf ./ szR);
    
    szR(n) = 1;
end

%% Get the reconstruction
rec = Hp;
clear Hp;
rec(pbd_array{:}) = realsqrt(1 - realpow(Lf, 2)) .* rec(pbd_array{:});

%% Get the lowpass subband
szLp = szF ./ D;
if any(fix(szLp ~= szLp))
    error('The downsampling factor must be able to divide the size of the FFT matrix!');
end
szLp(end) = szLp(end) / 2 + 1;

pbd_array_sml = pbd_array;
for n = 1 : N - 1
    pbd = pbd_array{n};
    pbd((length(pbd) + 3) / 2 : end) = pbd((length(pbd) + 3) / 2 : end) ...
        + szLp(n) - szX(n);
    pbd_array_sml{n} = pbd;
end

rec(pbd_array{:}) = rec(pbd_array{:}) + Lp(pbd_array_sml{:}) .* (Lf .* D^(N/2));



