function [L H dinfo pinfo] = DLifting(x, lsi)
% LIFTING wavelet decomposition with directional lifting
%
% USAGE:
%    [L H dinfo pinfo] = DLifting(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:
%    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.

% Lazy wavelet 
y = Split(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 = 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  = Split(x)
% Spliting along conlumns
%   Output:      y - Result of lazy wavelet
%             y{1} - Lower band
%             y{2} - Higher band
%   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
%                                           
M = size(x, 1);
M2 = floor(M/2);
lindex = 1;
hindex = 1 + mod(lindex,2);
y{1} = x(lindex:2:M, :);
y{2} = x(hindex:2:M, :);


function [h dinfo pinfo] = Predict(x, f, fd)
% primal lift: predicting along conlumns
%  direction indices and corresponding degrees
%   didx:     -3   -2   -1   0   1   2   3 
%   degree:   -18  -27  -45  90 45  27  18 
s = x{1};
d = x{2};
[r c]= size(d);
dset = -3:3;
% dset = -2:2;
% dset = -1:1;
% dset = [0, 1, -1, 2, -2, 3, -3];
ldset = length(dset);
h_all = zeros(r, c, ldset);
h = zeros(r, c);
% global bsize;
% bw = bsize;
bw = 1;
bh = bw;
br = r/bh;
bc = c/bw;
hsad = zeros(br, bc, ldset);
minsad = zeros(br,bc);
dinfo = zeros(br, bc);
pinfo = zeros(br, bc);
rextn = 3;
cextn = 1;

% 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 
% Predict for all spatial directions
for i = 1:length(dset);
    bos1 = rextn+1 + dset(i);
    bos2 = rextn+1 - dset(i);
    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
% 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;
        curblock = h_all((roff+1):(roff+bh), (coff+1):(coff+bw), 1:ldset);
       
% object function 1: SAD        
        hsad(i, j, 1:ldset) = sum(sum(abs(curblock(:,:,:)),1),2);
%         hsad(i, j, 1:ldset) = 10*fix(hsad(i, j, 1:ldset)/10);
% % object function 1: nnz
%         curblock = fix(curblock/20);
%         for m = 1:length(dset)
%         hsad(i, j, m) = nnz(curblock(:,:,m));
%         end
        [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
%  [minsad, dinfo] = min(hsad,[], 3);
 

function [l] = Update(x, f, fd, dinfo, pinfo)
% dual lift: updating along conlumns
s = x{1};
d = x{2};
dset = -3:3;
% dset = -2:2;
%  dset = -1:1;
% dset = [0, 1, -1, 2, -2, 3, -3];
ldset = length(dset);
[r c] = size(s);
% global bsize;
% bw = bsize;
% disp(['block size:  ', sprintf('%g', bw)]);
bw = 1;
bh = bw;
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
    bos1 = rextn+1 + dset(i);
    bos2 = rextn+1 - dset(i);
    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
