function x = iQDLifting_2(L, H, dinfo, pinfo, lsi)
% LIFTING wavelet decomposition with directional lifting on
% qincux lattice
% 
% USAGE:
%    x = iQDLifting(L, H, dinfo, pinfo, lsi)
% INPUT:   
%    L    - Lower band
%    H    - Higher band
%  dinfo  - directional information 
%  pinfo  - partition information
%    lsi  - Lifting lattice
% OUTPUT:
%    x    - reconstructed signal
% 
% REFERENCE:
%    W. Ding, F. Wu, and S. Li, ¡°Lifting-based wavelet transform
%    with directionally spatial prediction,¡± in Proc. Picture 
%    Coding Symposium 2004, San Francisco, CA, USA, Dec. 2004.

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

N = size(lsi, 1);
y = Normalize(y, lsi{N,1}, lsi{N,2});
for i = (N-1):-1: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} = Predict_2(y, f, dinfo, pinfo, fd);
    elseif lsi{i,1}=='u'
         y{1} = Update(y, f, dinfo, pinfo, fd);
    end
end
x = QMerge_2(y{1}, y{2}); 


function [h] = Predict_2(x, f, dinfo, pinfo, fd)
% primal lift: predicting along conlumns
%  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);
dset1 = [0 0; 0   1; 0  0; 1  0; 0 0;   1 -1; 0 -1; 0 -2];
dset2 = [0 0; -1 -2; -1 -1;-2 -1; 0 0; -2  0; -1 0; -1 1];
ldset = length(dset1);

h_all = zeros(r, c, ldset);
h = zeros(r, c);
bw = 4;
bh = 4;
br = r/bh;
bc = c/bw;
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~=1&i~=5)
    h_all(1:r, :, i) = d(1:r, :) + 0.5*extns( (1:r)+cextn + dset1(i,1), (1:c)+rextn + dset1(i,2) ) ...
                                 + 0.5*extns( (1:r)+cextn + dset2(i,1), (1:c)+rextn + dset2(i,2) );
%     end
end

% generate half pixel value for direction 1
halfh1 = zeros(r, 1+c);
i = 1:r;
j = 1:(1+c);
        half1(i, j) = (-extns(cextn + i - 2, rextn + j -1) + 3*extns(cextn + i - 1, rextn + j -1) + ...
            3*extns(cextn + i , rextn + j - 1) - extns(cextn + i + 1, rextn + j - 1) )/4;
        
% predict for direction 1
    h_all(1:r, 1:c, 1) = d(1:r, 1:c) + 0.5*half1(1:r, 1:c) + 0.5*half1(1:r, (1:c)+1);
    
% generate half pixel value for direction 5
halfh5 = zeros(1+r, c);
i = 1:(1+r);
j = 1:c;
        half5(i, j) = (-extns(cextn + i - 1, rextn + j - 2) + 3*extns(cextn + i -1, rextn + j -1) + ...
            3*extns(cextn + i - 1, rextn + j ) - extns(cextn + i - 1, rextn + j + 1) )/4;

% predict for direction 5
    h_all(1:r, 1:c, 5) = d(1:r, 1:c) + 0.5*half5(1:r, 1:c) + 0.5*half5((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;
        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, dinfo, pinfo, fd)
% dual lift: updating along conlumns
s = x{1};
d = x{2};
[r c]= size(d);

dset1 = [0 0;  1  2; 1  1; 2  1; 0 0;  2  0; 1 0; 1 -1];
dset2 = [0 0;  0 -1; 0  0;-1  0; 0 0; -1  1; 0 1; 0  2];
ldset = length(dset1);

l_all = zeros(r, c, ldset);
l = zeros(r, c);
bw = 4;
bh = 4;
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 directiona 
% 1 and 5 since they lie in half pixel position
for i = 1:ldset;
%     if (i~=1&i~=5)
    l_all(1:r, :, i) = s(1:r, :) - 0.25*extnd( (1:r)+cextn + dset1(i,1), (1:c)+rextn + dset1(i,2) ) ...
                                 - 0.25*extnd( (1:r)+cextn + dset2(i,1), (1:c)+rextn + dset2(i,2) );
%     end
end
% generate half pixel value for direction 1
halfh1 = zeros(r, 1+c);
i = 1:r;
j = 1:(1+c);
        half1(i, j) = (-extnd(cextn + i - 1, rextn + j) + 3*extnd(cextn + i, rextn + j) + ...
            3*extnd(cextn + i + 1, rextn + j) - extnd(cextn + i + 2, rextn + j) )/4;
    end
end
% update for direction 1
    l_all(1:r, 1:c, 1) = s(1:r, 1:c) - 0.25*half1(1:r, 1:c) - 0.25*half1(1:r, (1:c)+1);
    
% update half pixel value for direction 5
halfh5 = zeros(1+r, c);
for i = 1:(1+r)
    for j = 1:c
        half5(i, j) = (-extnd(cextn + i, rextn + j - 1) + 3*extnd(cextn + i, rextn + j) + ...
            3*extnd(cextn + i, rextn + j + 1) - extnd(cextn + i, rextn + j + 2) )/4;
    end
end
% update for direction 5
    l_all(1:r, 1:c, 5) = s(1:r, 1:c) - 0.25*half5(1:r, 1:c) - 0.25*half5((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


function y = Normalize(x, Lscale, Hscale)

y{1} = (1/Lscale)*x{1};
y{2} = (1/Hscale)*x{2};

function x  = QMerge_2(L, H)
% Spliting along conlumns
%    Input:     L - Lower band
%               H - Lower band
%   Output:      x - recontructed signal
%   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                                     
x = zeros(size([H,L]));
x(:,1:2:end) = H(:,:);
x(:,2:2:end) = L(:,:);

 
