function [L H dinfo pinfo] = QDLifting_1(x, lsi)
% Stage 1 decomposition: decompose original 2-D signal into Lower 
% band and Higher band with directional lifting  on qincux lattice
%
% USAGE:
%    [L H dinfo pinfo] = QDLifting(x, lsi)
% INPUT:   
%    x    - Input signal could be vector or matrix. If x is 
%           a matrix, lifting will be applied conlumn by conlumn
%    lsi  - Lifting lattice
% OUTPUT:
%    L    - Lower band
%    H    - Higher band
%  dinfo  - directional information 
%  pinfo  - partition information
% 
% REFERENCE:
%    Chuo-Ling Chang, Arian Maleki, and Bernd Girod, ˇ°Adaptive
%    wavelet transform for image for image compression via
%    directional quincux lifting", in Proc. IEEE International 
%    Workshop on Multimedia Signal Processing, MMSP-05, Shanghai,
%    China, November 2005.

% Lazy wavelet 

y = QSplit(x);  

N = size(lsi, 1);
for i = 1:(N-1)
     f  = lsi{i,2}; % filter taps
     fd = lsi{i,3}; % filter delay: z(-fd)
     cd = lsi{i,4}; % component delay which corresponds non-one diagonal components: z(-cd)
    if lsi{i,1}=='p'
        [y{2}, dinfo, pinfo] = Predict(y, f, fd);
    elseif lsi{i,1}=='u'
         y{1} = Update(y, f, fd, dinfo, pinfo);
    end
end

y  =  QExtract(y);

y = Normalize(y, lsi{N,1}, lsi{N,2});

L = y{1};
H = y{2};


function y = Normalize(x, Lscale, Hscale)
y{1} = Lscale*x{1};
y{2} = Hscale*x{2};

function y  =  QSplit(x)
% Qincux sampling
%   Output:      y - Result of lazy wavelet
%             y{1} - Even field of original size with zeros in odd field
%                    positions
%             y{2} - Odd field of original size with zeros in even field
%                    positions
%   Notes: the indices of the signal is supposed to
%          be 0, 1, 2, ..., M-1 and  0, 1, 2, ..., N-1
%          for the size of MxN. This is not consistent
%          with the convention of MATLAB
[h, w] = size(x);
y{1} = zeros(size(x));
y{2} = zeros(size(x));
xe = x(:, 1:2:end);
xo = x(:, 2:2:end);
xel = downsample(xe, 2);
xol = downsample(xo, 2, 1);
le = upsample(xel, 2);
lo = upsample(xol, 2, 1);
y{1}(:,1:2:end) = le(:,:);
y{1}(:,2:2:end) = lo(:,:);
xeh = downsample(xe, 2,1);
xoh = downsample(xo, 2);
he = upsample(xeh, 2, 1);
ho = upsample(xoh, 2);
y{2}(:,1:2:end) = he(:,:);
y{2}(:,2:2:end) = ho(:,:);



function y  =  QExtract(x)
[h, w] = size(x{1});
Let = x{1}(:, 1:2:w);
Lot = x{1}(:, 2:2:w);
Le = downsample(Let, 2);
Lo = downsample(Lot, 2, 1);
L(:,1:2:w) = Le(:,:);
L(:,2:2:w) = Lo(:,:);
y{1} = L;
[h, w] = size(x{2});
Het = x{2}(:, 1:2:w);
Hot = x{2}(:, 2:2:w);
He = downsample(Het, 2, 1);
Ho = downsample(Hot, 2);
H(:,1:2:w) = He(:,:);
H(:,2:2:w) = Ho(:,:);
y{2} = H;

function [h dinfo pinfo] = Predict(x, f, fd)
% primal lift
%  direction indices and corresponding degrees
%   didx:      1    2    3    4    5    6    7     8
%   degree:    0  -27  -45  -63  -90 -117  -135  -153 
s = x{1};
d = x{2};
[r c]= size(d);
dset = [0 1; 1 2; 0 0; 2 1; 1 0; 2 -1; 0 0; 1 -2];
ldset = length(dset);

h_all = zeros(r, c, ldset);
h = zeros(r, c);
bw = 8;
bh = 8;
br = r/bh;
bc = c/bw;
hsad = zeros(br, bc, ldset);
minsad = zeros(br,bc);
dinfo = zeros(br, bc);
pinfo = zeros(br, bc);
rextn = 2;
cextn = 2;

% symmetric extension
extns  = extn2D(s',1:(rextn+1), 'symodd')'; % odd symmtric extension with rextn samples for each row
extns  = extn2D(extns, 1:(cextn+1), 'symodd'); % odd symmtric extension with cextn samples for each row 
% generate 
% Predict for all spatial directions except directiona 
% 3 and 7 since they lie in half pixel position

for i = 1:ldset;
%     if (i~=3&i~=7)
    h_all(:, :, i) = d(:, :) - 0.5*extns( (1:r)+cextn + dset(i,1), (1:c)+rextn + dset(i,2) ) ...
                                 - 0.5*extns( (1:r)+cextn - dset(i,1), (1:c)+rextn - dset(i,2) );
%     end
end

% generate half pixel value for direction 3
halfh3 = zeros(1+r, 1+c);
half3(1:(1+r), 1:(1+c)) = (-extns(cextn + 1 + (1:(1+r)), rextn -2 + (1:(1+c))) + 3*extns(cextn + (1:(1+r)), rextn -1 + (1:(1+c))) + ...
            3*extns(cextn + (1:(1+r)) - 1, rextn + (1:(1+c))) - extns(cextn + (1:(1+r)) - 2, rextn + (1:(1+c)) + 1) )/4;


% predict for direction 3
    h_all(1:r, 1:c, 3) = d(1:r, 1:c) - 0.5*half3(1:r, 1:c) - 0.5*half3((1:r)+1, (1:c)+1);
    
% generate half pixel value for direction 7
halfh7 = zeros(1+r, 1+c);
half7(1:(1+r), 1:(1+c)) = (-extns(cextn + (1:(1+r)) - 2 , rextn + (1:(1+c)) - 1 -1) + 3*extns(cextn + (1:(1+r)) -1, rextn + (1:(1+c)) -1) + ...
            3*extns(cextn + (1:(1+r)) , rextn + (1:(1+c)) + 1 - 1) - extns(cextn + (1:(1+r)) + 1, rextn + (1:(1+c)) + 2 -1) )/4;

% predict for direction 7
    h_all(1:r, 1:c, 7) = d(1:r, 1:c) - 0.5*half7((1:r), (1:c)+1) - 0.5*half7((1:r)+1, (1:c));

% finding the best direction which results in 
% smallest residue in SAD

for i = 1:br
for j = 1:bc
roff = (i-1)*bh;
coff = (j-1)*bw;
hsad(i, j, 1:ldset) = sum(sum(abs(h_all((roff+1):(roff+bh), (coff+1):(coff+bw), 1:ldset)),1),2);
[minsad(i,j) dinfo(i,j)] = min(hsad(i, j, :));
h((roff+1):(roff+bh), (coff+1):(coff+bw)) = h_all((roff+1):(roff+bh), (coff+1):(coff+bw), dinfo(i,j));
end
end


 

function [l] = Update(x, f, fd, dinfo, pinfo)
% dual lift
s = x{1};
d = x{2};
[r c]= size(d);
dset = [0 1; 1 2; 0 0; 2 1; 1 0; 2 -1; 0 0; 1 -2];
ldset = length(dset);

l_all = zeros(r, c, ldset);
l = zeros(r, c);
bw = 8;
bh = 8;
br = r/bh;
bc = c/bw;
rextn = 2;
cextn = 2;

% symmetric extension
extnd  = extn2D(d',1:(rextn+1), 'symodd')'; % odd symmtric extension with rextn samples for each row
extnd  = extn2D(extnd, 1:(cextn+1), 'symodd'); % odd symmtric extension with cextn samples for each row 

% Update for all spatial directions except direction 
% 3 and 7 since they lie in half pixel position

for i = 1:ldset;
%     if (i~=3&i~=7)
    l_all(:, :, i) = s(:, :) + 0.25*extnd( (1:r)+cextn + dset(i,1), (1:c)+rextn + dset(i,2) ) ...
                                 + 0.25*extnd( (1:r)+cextn - dset(i,1), (1:c)+rextn - dset(i,2) );
%     end
 end

% generate half pixel value for direction 3
half3 = zeros(1+r, 1+c);
i = 1:(1+r);
j = 1:(1+c);
        half3(i, j) = (-extnd(cextn + i + 1, rextn + j -2) + 3*extnd(cextn + i, rextn + j -1) + ...
            3*extnd(cextn + i - 1, rextn + j) - extnd(cextn + i - 2, rextn + j + 1) )/4;

% predict for direction 3
    l_all(1:r, 1:c, 3) = s(1:r, 1:c) + 0.25*half3(1:r, 1:c) + 0.25*half3((1:r)+1, (1:c)+1);
 
% generate half pixel value for direction 7
half7 = zeros(1+r, 1+c);
i = 1:(1+r);
j = 1:(1+c);
        half7(i, j) = (-extnd(cextn + i - 2 , rextn + j - 1 -1) + 3*extnd(cextn + i -1, rextn + j -1) + ...
            3*extnd(cextn + i , rextn + j + 1 - 1) - extnd(cextn + i + 1, rextn + j + 2 -1) )/4;
        
% predict for direction 7
    l_all(1:r, 1:c, 7) = s(1:r, 1:c) + 0.25*half7((1:r), (1:c)+1) + 0.25*half7((1:r)+1, (1:c));

% Pick the right one with dinfo

for i = 1:br
    for j = 1:bc
        roff = (i-1)*bh;
        coff = (j-1)*bw;
        l((roff+1):(roff+bh), (coff+1):(coff+bw)) = l_all((roff+1):(roff+bh), (coff+1):(coff+bw), dinfo(i,j));
    end
end
