function [y,dinfo, pinfo]= QDLifting2D(x, J, lsi)
% 2D wavelet decomposition with directional lifting scheme 
% on quincux lattice
%
% USAGE:
%    [y,dinfo, pinfo] = QDLifting2D(x, J, lsi)
% 
% INPUT:   
%    x    - input 2D signals
%    J    - Decomposition levels
%    lsi  - lifting lattice
% OUTPUT:
%    y       - Decomposition Coefficients
%    y{J}    - Lower band
%    y{i}{1} - Lower-Higher band
%    y{i}{2} - Higher-Lower band
%    y{i}{3} - Higher-Higher band
% NOTES: 
%    Spatial correlation of high-pass band after quincux
%    lifting usually is very weak. So it's not necessary
%    to be further decomposed with quincux lifting. To 
%    generate pyramid as conventional wavelet, the high-
%    pass band are decomposed with usual 5/3 wavelet.

LL = x;
for i = 1:J
[L, H, dinfo{i}{1}, pinfo{i}{1}] = QDLifting_1(LL, lsi);
[LL, LH, dinfo{i}{2}, pinfo{i}{2}] = QDLifting_2(L, lsi);
lsi53 = {'p', [-0.5 -0.5], [1], [0];'u', [0.25 0.25], [0], [0]; [sqrt(2)], [sqrt(2)/2], [], []};
[HL, HH] = Lifting1D(H', lsi53);
y{i}{1} = LH;
y{i}{2} = HL';
y{i}{3} = HH';
LL = LL;
end
y{J+1} = LL;