function [ wcs ] = createEmptyWCS( data, waveletType, waveletOrder )
%createEmptyWCS Create an empty wavelet coefficient structure (WCS)
%   wcs = createEmptyWCS( data, waveletType, waveletOrder ) creates an
%   empty wavelet coefficient structure (a struct for people who know C).
%   The details of this struct are described below in the output section.
%
% INPUTS :
%    data : an array, NR x NC in size. The array can be empty since we are
%           not going to fill up the wavelet coefficient matrix.
%
%    waveletType, waveletOrder : specifies the quadrature mirror filter to
%           use. The allowed combinations of (Type, Order) are:
%
%         waveletType      | waveletOrder
%
%         'Haar'           |  2
%         'Beylkin'        |  1, 2, 3, 4, 5
%         'Daubechies'     |  4, 6, 8, 10, .... 20
%         'Symmlets'       |  4, 5, 6, ...10
%         'Vaidyanathan'   |  0
%         'Battle'         |  1, 3, 5
%
% CAVEATS :
% 
%    For the time being, we will insist that NR = NC = power-of-two.
%
% OUTPUT:
%   wcs, a variable of type Wavelet Coefficient Structure, containing the
%        following fields:
%
%   wcs.waveletType : Wavelet type used for decomposing the data. This WCS 
%                     should contain wavelet coefficients which were
%                     obtained by projecting the data on this type of
%                     wavelet basis set.
%   wcs.waveletOrder: Contains the order of the wavelet basis set. You need
%                     to know the wavelet type and order to fully specify a
%                     wavelet basis set.
%   wcs.qmf         : The quadrature mirror filter that should be used to
%                     decompose data onto a wavelet basis set. This process
%                     is also called doing a WAVELET TRANSFORM. In reality
%                     qmf should be an object, with wavelet type and order
%                     as its attributes. Unfortunately, it's a 1D array
%                     with the profile of a wavelet in it, and is used for
%                     convolving the data when doing a wavelet transform.
%   wcs.nlevels     : number of levels L in the wavelet structure. L = max(
%                     [ nextpow2(NR), nextpow2(NC) ] )
%   wcs. baseVal    : A residual value (roughly, a representative value)
%                     for the entire data. The wavelets add structure 
%                     around this value.
%   wcs.coeffs      : A cell array, (nlevels, 3) in size. The value of each
%                     cell coeffs{l, i}, is a 2D array, 2^{l-1} x 2^{l-1}
%                     in size. coeffs{l, 1}, a 2D array, contains wavelet
%                     coefficients at scale l (higher l means finer spatial
%                     scales) decribing variations in the x-direction,
%                     coeffs{l, 2} does so for variations in the
%                     y-direction and coeffs{l, 3} contains variations in
%                     x-y direction.
%   wcs.origDataSize: (NR, NC)
%   wcs.paddedDataSize: The data size used when we generated the qmf. qmf
%                     requires that the data set be power-of-two and we
%                     will have to pad it some day.
%
%  LIBRARY REQUIREMENTS :  Needs WaveLab850 to run
%
%  Jaideep Ray, SNL, CA, 2/9/2011.
% =========================================================================

% ---- Fill up the description fields
wcs.waveletType  = waveletType ;
wcs.waveletOrder = waveletOrder ;
wcs.qmf          = MakeONFilter(wcs.waveletType, wcs.waveletOrder) ;
wcs.origDataSize = size(data) ;
fprintf('createEmptyWCS() : Original size of the data set is %d X %d\n', ...
        wcs.origDataSize(1), wcs.origDataSize(2)) ;
    
% ---- If the data set is not power of two, get the padded dataset size.
% Also, ensure it is square
dyadicSize1 = ceil( log(wcs.origDataSize(1)) / log(2) ) ;
dyadicSize2 = ceil( log(wcs.origDataSize(2)) / log(2) ) ;
dyadicSize = max( [dyadicSize1 dyadicSize2] ) ;
if ( (2^dyadicSize1 ~= wcs.origDataSize(1)) || ...
     (2^dyadicSize2 ~= wcs.origDataSize(2)) )
    wcs.paddedDataSize = [2^dyadicSize 2^dyadicSize] ;
else
    wcs.paddedDataSize = wcs.origDataSize ;
end
fprintf('createEmptyWCS() : The padded size is %d X %d, the dyadic size is %d X %d\n', ...
        wcs.paddedDataSize(1), wcs.paddedDataSize(2), dyadicSize, ...
        dyadicSize) ;
    
% ---- Calculate the number of levels I will have. 
wcs.nlevels = dyadicSize ;

% ---- Allocate the cell array and zero it out at all levels.
wcs.baseVal = 0.0 ;
wcs.coeffs = cell(wcs.nlevels, 3) ;
for l = 1:wcs.nlevels
    wcs.coeffs{l, 1} = zeros(2^(l-1), 2^(l-1)) ;
    wcs.coeffs{l, 2} = zeros(2^(l-1), 2^(l-1)) ;
    wcs.coeffs{l, 3} = zeros(2^(l-1), 2^(l-1)) ;
end

end

