function [ wcs ] = computeWaveletTransform( data, wcs_in )
%computeWaveletTransform Take data in and compute its wavelet transform
%   wcs = computeWaveletTransform( data, wcs_in ) takes in a 2D array
%   (data) and computes its wavelet transform. The wavelet coefficients are
%   stored in the WCS struct that is returned. The QMF to be used for
%   computing the transform is also stored in the WCS. The incoming WCS,
%   wcs_in, contains info about the dataset, the QMf that you want to use
%   etc. In order to call this function, you should have called
%   createEmptyWCS beforehand and initialized the incoming WCS. The
%   outgoing WCS is just the incoming one, with the coeffs field filled in.
%
% INPUTS
%    data   : a 2D array whose wavelet transform has to be taken.
%
%    wcs_in : a Wavelet Coefficient Structure that has been initialized by
%             a call to createEmptyWCS(). Make sure that the size of the
%             dataset and that used to create wcs_in (using
%             createEmptyWCS()) are the same!
%
% OUTPUT :
%    wcs    : The output WCS. It is the same as the incoming WCS, wcs_in,
%             but now has wcs.coeffs and wcs.baseVal filled in with wavelet 
%             coefficients, as a result of the  wavelet transform.
% 
%  LIBRARY REQUIREMENTS :  Needs WaveLab850 to run
%
%  Jaideep Ray, SNL, CA, 2/9/2011.
% =========================================================================

% ---- Make sure that the incoming dataset has the same size that the WCS
% was created for.
diff = size(data) - wcs_in.origDataSize ;
if ( norm(diff) ~= 0 )
    error('computeWaveletTransform() : the incoming dataset and WCS should have same size') ;
end

% ---- At this point, if not dyadic, pad the dataset to be dyadic. It
% should always be square, and cws_in.paddedDataSize should indicate a
% padded, dyadic, square array. Let the pad be at one end of the array
D = zeros(wcs_in.paddedDataSize) ;
D(1:size(data,1), 1:size(data, 2)) = data ;

% ---- Take the 2D wavelet transform
coarseLevel = 0 ;
wc = FWT2_PO(D, coarseLevel, wcs_in.qmf) ;

% ---- Unroll the 2D wc array into the outgoing wcs
wcs = wcs_in ;
wcs.baseVal = wc(1, 1) ;
for l = 1 : wcs.nlevels
    dyadNumber = l - 1 ; % counted from zero
    lb = min( dyad( dyadNumber ) ) ;
    ub = max( dyad( dyadNumber ) ) ;
    wcs.coeffs{l, 1} = wc(1:ub-lb+1, lb:ub) ;
    wcs.coeffs{l, 2} = wc(lb:ub, 1:ub-lb+1) ;
    wcs.coeffs{l, 3} = wc(lb:ub, lb:ub) ;
end
end

