function y = nsctdec(x, levels, dfilt, pfilt )
% NSSCDEC   Nonsubsampled Contourlet Transform Decomposition
%
%	y = nsctdec(x, levels, [dfilt, pfilt] )
%
% Input:
%   x:      
%       a matrix, input image
%   levels:  
%       vector of numbers of directional filter bank decomposition levels 
%       at each pyramidal level (from coarse to fine scale).
%       If the number of level is 0, a critically sampled 2-D wavelet 
%       decomposition step is performed.
%       The support for the wavelet decomposition has not been verified!!!!!!
%   dfilt:  
%       a string, filter name for the directional decomposition step.
%       It is optional with default value 'dmaxflat7'. See dfilters.m for all
%       available filters.
%   pfilt:  
%       a string, filter name for the pyramidal decomposition step.
%       It is optional with default value 'maxflat'. See atrousfilters.m for 
%       all available filters. 
%
% Output:
%   y:  a cell vector of length length(nlevs) + 1, where except y{1} is 
%       the lowpass subband, each cell corresponds to one pyramidal
%       level and is a cell vector that contains bandpass directional
%       subbands from the DFB at that level.
%
% Index convention:
%   Suppose that nlevs = [l_J,...,l_2, l_1], and l_j >= 2.
%   Then for j = 1,...,J and k = 1,...,2^l_j 
%       y{J+2-j}{k}(n_1, n_2)
%   is a contourlet coefficient at scale 2^j, direction k, and position
%       (n_1 * 2^(j+l_j-2), n_2 * 2^j) for k <= 2^(l_j-1), 
%       (n_1 * 2^j, n_2 * 2^(j+l_j-2)) for k > 2^(l_j-1).
%   As k increases from 1 to 2^l_j, direction k rotates clockwise from
%   the angle 135 degree with uniform increment in cotan, from -1 to 1 for
%   k <= 2^(l_j-1), and then uniform decrement in tan, from 1 to -1 for 
%   k > 2^(l_j-1).
%
% See also:	ATROUSFILTERS, DFILTERS, NSCTREC, NSFBDEC, NSDFBDEC.

%  History:
%     02/17/04 Created by Jianping Zhou.
%     08/07/04 Modified by Jianping Zhou. Incorporte the fast implementation 
%              of convolution algorithm by Jason Laska. 
%     08/30/04 Modified by Arthur. L. Cunha, added more pyramid filters.
%     10/17/04 Modified by Arthur. L. Cunha, replaced periodic with symmetric extension
%     01/24/05 Modified by Jianping Zhou, changed function names and corrected a bug. 
%     10/31/05 Modified by Arthur L Cunha, corrected a bug in the pyramid decomposition (a wrong index)

% Check input
if ~isnumeric( levels )
    error('The decomposition levels shall be integers');
end
if isnumeric( levels )
    if round( levels ) ~= levels
        error('The decomposition levels shall be integers');
    end
end
        
if ~exist('dfilt', 'var')
    dfilt = 'dmaxflat7' ;
end;

if ~exist('pfilt', 'var')
    pfilt = 'maxflat' ; 
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get fan filters, parallelogram filters, and pyramid filters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get the directional filters for the critically sampled DFB.
filters = cell(4) ;
[h1, h2] = dfilters(dfilt, 'd');
% A scale is required for the nonsubsampled case.
h1 = h1./sqrt(2) ;
h2 = h2./sqrt(2) ;

% Generate the first-level fan filters by modulations.
filters{1} = modulate2(h1, 'c');
filters{2} = modulate2(h2, 'c'); 

% Obtain the parallelogram filters from the diamond filters
[filters{3}, filters{4}] = parafilters( h1, h2 ) ;

[h1, h2, g1, g2] = atrousfilters(pfilt); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nonsubsampled Contourlet transform with tree structure filter banks
% Nonsubsampled pyramids make multiresolution decomposition.
% Nonsubsampled directional filter banks make directional decomposition.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Numer of levels:
clevels = length( levels ) ;
nIndex = clevels + 1 ;
% Initialize the output
y = cell(1, nIndex) ;

% Nonsubsampled pyramid decomposition
for i= 1 : clevels    
    
    % Nonsubsampled Contourlet transform
    % Nonsubsampled pyramids decomposition
    [xlo, xhi] = nsfbdec(x, h1, h2, i-1) ;
        
    if levels(nIndex-1) > 0        
        % Nonsubsampled DFB decomposition on the bandpass image
        xhi_dir = nsdfbdec(xhi, filters, levels(nIndex-1));        
        y{nIndex}=xhi_dir ;
    else
        % Copy the result directly:
        y{nIndex}=xhi ;
    end
    
    % Update the index for the Nonsubsampled Pyramdis
    nIndex = nIndex - 1 ;
    
    % Prepare for next iteration
    x = xlo ;
    
end

% The lowpass output
y{1}=x;