function [ field ] = inverseWaveletTransform( wcs, maxlevel )
%inverseWaveletTransform Reconstructs a field from the wavelet transform in wcs
%   field = inverseWaveletTransform( wcs, maxlevel ) recreates a field from the
%   wavelet coefficients in the Wavelet Coefficient Structure wcs, which
%   also contains the quadrature mirror filter that was used to calculate
%   the wavelet coefficients in the first place. Only wavelet coefficients
%   at scales less than or equal to maxlevel are used; wavelet coeffs at
%   level > maxlevel are set to zero.
%
%   field = inverseWaveletTransform( wcs ) does the same as above, but uses
%           wavelet coefficients at all scales 
%
%   A WCS can hold the wavelet transform for a padded dataset that is
%   square and power of 2 in each direction. If the original dataset was
%   not so, the field returned is the size of the original dataset, not the
%   padded one.
%
% INPUTS
%   wcs : a WCS, filled up by computeWaveletTransform(). It has the wavelet
%         coefficients of a dataset subjected to a Wavelet Tranform and the
%         QMF used for it
%
%   maxlevel, 0 <= maxlevel <= wcs.nlevels, which states that we should use
%             wavelet coefficients at scales up to maxlevel, when doing the
%             transform. If this is omitted, all levels are used
%             
%
% OUTPUTS 
%   field : a 2D array, or the same size as the original dataset containing
%           the inverse wavelet transform. If we haven't fiddled with the
%           wavelet coefficients, field == orig_data.
%
% LIBRARY REQUIREMENTS :  Needs WaveLab850 to run
%
%  Jaideep Ray, SNL, CA, 2/9/2011.
% =========================================================================

% ---- Extract the wavelet coeffs from the WCS struct.
wc = zeros( wcs.paddedDataSize ) ;
wc(1, 1)  = wcs.baseVal ;

if (maxlevel > 0)
    limit = wcs.nlevels ; % by default, we use all levels of wavelet coeffs
    if (nargin == 2)
        limit  = min( [maxlevel, wcs.nlevels] ) ;
    end
    for l = 1 : limit
        dyadNumber = l - 1 ; % counted from zero
        lb = min( dyad( dyadNumber ) ) ;
        ub = max( dyad( dyadNumber ) ) ;
        wc(1:ub-lb+1, lb:ub) = wcs.coeffs{l, 1} ;
        wc(lb:ub, 1:ub-lb+1) = wcs.coeffs{l, 2} ;
        wc(lb:ub, lb:ub)     = wcs.coeffs{l, 3} ;
    end
end

% ---- Do the actual inverse transform
D = IWT2_PO(wc, 0, wcs.qmf) ;

% ---- Cut out the pad region and return
field = D(1:wcs.origDataSize(1), 1:wcs.origDataSize(2)) ;
end

