function x = iDLifting_comb(L, H, dinfo, pinfo, lsi)
% LIFTING wavelet decomposition with directional lifting
%
% USAGE:
%    x = iDLifting(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(y, f, dinfo, pinfo, fd);
    elseif lsi{i,1}=='u'
        y{1} = Update(y, f, dinfo, pinfo, fd);
    end
end
x = Merge(y{1}, y{2}); 


function [h] = Predict(x, f, dinfo, pinfo, fd)
% primal lift: predicting along conlumns
s = x{1};
d = x{2};
[r,c] = size(d);
h = zeros(size(d));
dset = -5:5;
ldset = length(dset);
h_all = zeros(r, c, ldset);
bw = 16;
bh = 16;
br = r/bh;
bc = c/bw;

rextn = 3;
cextn = 1;
extns  = extn2D(s',1:(rextn+1), 'symodd')'; % odd symmtric extension with 3 samples for each row
extns  = extn2D(extns, 1:(cextn+1), 'symodd'); % odd symmtric extension with 1 samples for each row

% Predict for all spatial directions
for i = 1:length(dset);
    if abs(dset(i))<4
    bos1 = rextn+1 + dset(i);
    bos2 = rextn+1 - dset(i);
    else
        if dset(i)<0
    bos1 = 8 + dset(i);
    bos2 = bos1 + 1;
        else
            bos1 = dset(i);
            bos2 = bos1 - 1;
        end 
    end
    h_all(1:r, :, i) = d(1:r, :) + 0.5*extns((1:r)+cextn, bos1:(bos1+c-1)) + 0.5*extns((1:r)+cextn + 1, bos2:(bos2+c-1));
end

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};
dset = -5:5;
ldset = length(dset);
[r c] = size(s);
bw = 16;
bh = 16;
br = r/bh;
bc = c/bw;
l = zeros(r, c);
l_all  = zeros(r, c, ldset);

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

% Update for all spatial directions
for i = 1:ldset
    if abs(dset(i))<4
    bos1 = rextn+1 + dset(i);
    bos2 = rextn+1 - dset(i);
    else
        if dset(i)<0
    bos1 = 8 + dset(i);
    bos2 = bos1 + 1;
        else
            bos1 = dset(i);
            bos2 = bos1 - 1;
        end 
    end
    l_all(1:r, :, i) = s(1:r, :) - 0.25*extnd((1:r)+cextn - 1, bos1:(bos1+c-1)) - 0.25*extnd((1:r)+cextn , bos2:(bos2+c-1));
end

% 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  = Merge(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([L;H]));
M = size(x,1);

lindex = 1;
hindex = 1 + mod(lindex,2);
 x(lindex:2:M, :) = L;
 x(hindex:2:M, :) = H;