pdtdfbdec.m

PDTDFB toolbox The filter bank is described in: The Shiftable Complex Directional Pyramid—Pa

function y = pdtdfbdec(x, nlevs, lpfname, dfname, wfname, res)
% PDTDFBDEC   Pyramidal Dual Tree Directional Filter Bank Decomposition
%
%	y = pdtdfbdec(x, nlevs, [lpfname], [dfname], [wfname], [res])
%
% Input:
%   x:      Input image
%   nlevs:  vector of numbers of directional filter bank decomposition levels
%           at each pyramidal level (from coarse to fine scale).
%           0 : Laplacian pyramid level
%           1 : Wavelet decomposition
%           n : 2^n directional PDTDFB decomposition
%   lpfname : [optional] The name of the laplacian lowpass filter that is used, 
%           default is 'nalias'
%           if lpfname is 'nalias' then the lpftype = 2
%   dfname  :[optional] The name of the diamond or fan filter that is used,
%           'meyer': default (very big filters - very slow)
%           'pkva''pkva6', 'pkva8', 'pkva12' (Contourlet toolbox, using ladder structure) 
%
%   wfname: [optional] The name of the wavelet filter that is used, default is '9-7'
%   res :   [optional] Residual high passband, default 1
%           0 is no, 1 is yes
%   mode :
%
% 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.
%
% Call function: LPDEC, 
% See also:	PDTDFBREC, PDFBDEC, TDFBDEC,
%
% Note : PDTDFB data structure y{resolution}{1}{1-2^n} : primal branch
%                              y{resolution}{2}{1-2^n} : dual branch

% running time for 1024 image size [3]: 31 sec 
%                                  [5]: 33 sec
%                             [4 7 7 5]  : 41 sec
%

if ~exist('res','var')
    res = 0 ; % default implementation have a highpass residual band
end

if ~exist('lpfname','var')
    lpfname = 'nalias' ; % default implementation by meyer type FB
end

if ~exist('dfname','var')
    dfname = 'meyer' ; % default implementation by meyer type FB
    % Note: meyer diamond is different from meyer type low pass
end

if ~exist('wfname','var')
    wfname = '9-7' ; % default implementation by meyer type FB
end

if ~exist('mode','var')
    mode = 'pdtdfb' ; % default implementation is PDTDFB decomposition
end

if length(nlevs) == 0
    y = {x};

else %----------------------------------------------------------------

    if res % remove the highpass component to residual image
        N = 128;
        cutoff = 2*pi;
        % cutoff =    pi;
        rev = (8*pi^2/3)/cutoff;
        xa = 0:rev/(N):rev;
        % Compute support of Fourier transform of phi.
        int1 = find((xa < 2*pi/3));
        int2 = find((xa >= 2*pi/3) & (xa < 4*pi/3));

        % Compute Fourier transform of phi.
        phihat = zeros(1,N+1);
        phihat(int1) = ones(size(int1));
        phihat(int2) = cos(pi/2*meyeraux(3/2/pi*xa(int2)-1));

        phihat = [phihat,phihat((end-1):-1:2)];
        psihat = (1-phihat).^(1/2);
        h = ifftshift(ifft(phihat));
        g = ifftshift(ifft(psihat));

        % [h,g] = wfilters('dmey','l');
        h = fitmat(h,[1,32]); h = h(2:end);
        h = h./sum(h);
        H = h'*h;
        G = fftshift(ifft2( (1 - abs(fft2(H)).^2).^(1/2) ) );
        G = fitmat(real(G),31);
        x_res = efilter2(x,G,'sym');
        x = efilter2(x,H,'sym');
    end

    % determine the kind of laplacian pyramid decomposition based on
    % lpfname decide type of pyramidal decomposition
    switch lower(lpfname)
       case ('special')
            lptype = 3; % 
            disp('Frequency implemenation for the first two level')
       case ('nalias')
            lptype = 2; % noaliasing type pyramid
            disp('Noalias frame pyramid')
        case ('fp')
            lptype = 0; % Burt and Aldeson pyramid
            disp('framming pyramid')
        otherwise
            lptype = 1; % Burt and Aldeson pyramid
            disp('Burt and Aldeson pyramid')
    end
    
    switch nlevs(end)
        case {0}
            % Get the pyramidal filters
            [h,g] = pfilters(lpfname);

            % keep the laplacian highpass
            [xlo, xhi] = lpdec(x,h,g,lptype);
            xhi_dir = xhi;
        case {1}
            % Get the pyramidal filters
            [h,g] = pfilters(lpfname);
            % wavelet decomposittion
            [h, g] = pfilters(wfname);

            [xlo, xLH, xHL, xHH] = wfb2dec(x, h, g);
            xhi_dir = {xLH, xHL, xHH};
        otherwise
            if (lptype~=3)
                % laplacian pyramid
                % Get the pyramidal filters
                [h,g] = pfilters(lpfname);
                [xlo, xhi] = lpdec(x,h,g,lptype);

                % DFB on the bandpass image
                switch dfname        % Decide the method based on the filter name
                    case {'pkva6', 'pkva8', 'pkva12', 'pkva'}
                        % Use the ladder structure (whihc is much more efficient)
                        % disp('ladder')
                        % dual tree dfb
                        tmp = tdfbdec_l(xhi, nlevs(end),'primal', dfname);
                        xhi_dir{1} = tmp;
                        tmp = tdfbdec_l(xhi, nlevs(end),'dual', dfname);
                        xhi_dir{2} = tmp;
                    otherwise
                        % General case
                        % dual tree dfb
                        tmp = tdfbdec(xhi, nlevs(end),'primal', dfname);
                        xhi_dir{1} = tmp;
                        tmp = tdfbdec(xhi, nlevs(end),'dual', dfname);
                        xhi_dir{2} = tmp;
                end

            else
                % frequency implementation

                alpha = 0.15;

                s = size(x);

                % create the grid and transform to circle grid
                S1 = -1.5*pi:pi/(s(1)/2):1.5*pi;
                S2 = -1.5*pi:pi/(s(2)/2):1.5*pi;
                [x1, x2] = meshgrid(S2,S1);
                r = [0.4 0.5 1-alpha, 1+ alpha];

                [x1, x2]=tran_sf(x1,x2);

                rd = sqrt(x1.*x1+x2.*x2);
                theta = angle(x1+sqrt(-1)*x2);

                % Low pass window
                sz = size(rd);
                cen = rd((sz(1)+1)/2,:);
                cen = abs(cen);
                fl = fun_meyer(cen,pi*[-2 -1 r(1:2)]);
                FL =  fl'*fl;

                % high pass window
                ang = pi/4*[-alpha alpha];
                ang = [-pi/4+ang(1:2), 3*pi/4+ang(1:2)];

                f3 = fun_meyer_curv(rd, theta, pi*r, ang, 's');
                f3 = periodize(f3,s, 's');

                % take out the center and square root
                FL = sqrt(fftshift(FL(s(1)/4+1:s(1)/4+s(1),s(2)/4+1:s(2)/4+s(2))));
                f3 = sqrt(fftshift(f3(s(1)/4+1:s(1)/4+s(1),s(2)/4+1:s(2)/4+s(2))));

                % actual transform by FFT
                imf = fft2(x);
                sz = s/2;
                dm = [2 2];

                imf2 = (1/prod(dm))*periodize(imf.*FL, sz);

                xlo = ifft2(imf2(1:sz(1),1:sz(2)));
                xhi = ifft2(imf.*f3);

                %now doing normal DFB on both branches
                tmp = tdfbdec_l(real(xhi), nlevs(end),'primal', dfname);
                xhi_dir{1} = tmp;
                tmp = tdfbdec_l(imag(xhi), nlevs(end),'primal', dfname);
                xhi_dir{2} = tmp;
            end

    end

    % Recursive call on the low band
    res2 = 0; % from the second level , there is no residual band
    ylo = pdtdfbdec(xlo, nlevs(1:end-1), lpfname, dfname, wfname, res2);

    % Add bandpass directional subbands to the final output
    y = {ylo{:}, xhi_dir};

    if res % add back the residual band to the end of xhi_dir
        y{end+1} = x_res;
    end
end