function  Y = mylwt2(X, ls, nivaa)
% mylwt2          My variant of lwt2(..), use as lwt2 with Lifting
% Here mirror extension is done, while matlab seems to do zero-padding.
% Some alternative lifting schemes:
% lsdb1 = {'d',-1,0;'p',1/2,0;sqrt(2),1/sqrt(2),[]};  % db1 wavelet (Haar)
% lshaar = {'d',-1,0;'p',1/2,0;2,1,[]};  % same with "integer" results
% ls53 = {'d',-[1,1]/2,1;'p',[1,1]/4,0;1,1,[]};  % 5/3 wavelet
% 
% use: 
%   Y = mylwt2(A, ls, level);   
%----------------------------------------------------------------------
% A     the image (matrix of doubles), size MxN, both M and N should have 
%       2^level (typically 8, 16 or 32) as a factor. If not myimadjust.m
%       may be used to extend the image.
% Y     the coefficients, same size as X
% ls    cell array or a name as in liftwave (Matlab Wavelet Toolbox)
%       the name may also be 'j97' to get the wavelet as used in JPEG-
%       2000, M.D. Adams and R.Ward: Wavelet Transforms in the JPEG-2000 
%       Standard. (the obsolete name 'db97' is the same)
%       Note that 'jp2' is not the same as '9.7' in liftwave.
%       Also, the name 'm97' (better name could be 'e97' or 'ee97') let
%       the energi levels of low-pass and high-pass be more equal to
%       each other, (JPEG favours low-pass as it gives better visual
%       reconstruction when compressing)
% level number of levels for the wavelet, typically 3, 4 or 5
%----------------------------------------------------------------------
% Examples:
% A = double(imread('lena.bmp')) - 128;
% Y = mylwt2(A, 'j97', 3);
% Ar = myilwt2(Y, 'j97', 3);  
% if (norm(A-Ar,'fro') < 1e-10); disp(' O K '); else disp('NOT OK'); end;

%----------------------------------------------------------------------
% Copyright (c) 2009.  Karl Skretting.  All rights reserved.
% University of Stavanger (Stavanger University), Signal Processing Group
% Mail:  karl.skretting@uis.no   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  04.06.2009  KS: m-file made 
% Ver. 1.1  08.03.2010  KS: made the function 'cleaner'
%----------------------------------------------------------------------

% image A is denoted X below, and level is denoted nivaa
[M,N] = size(X);
K = 2^nivaa;
Y = X;       % do transform in-place in Y

% check arguments
if ~isnumeric(X)
    disp('mylwt2: The supplied image A is not numeric.');
    return
end
if ~((mod(M,K)==0) && (mod(N,K)==0))
    disp(['mylwt2: Image size has not ',int2str(K),' as factor.']);
    return
end

if ischar(ls) 
    % 
    i = strfind(ls,'79');    % allow for '97' to be flipped in ls
    if (numel(i) == 1); ls(i:(i+1)) = '97'; end;
    i = strfind(ls,'7.9');    % allow for '9.7' to be flipped in ls
    if (numel(i) == 1); ls(i:(i+2)) = '9.7'; end;
    %
    if (strcmpi(ls,'j97') || strcmpi(ls,'db97') || strcmpi(ls,'m97'))
        p1 = -1.586134342059924;
        u1 = -0.052980118572961;
        p2 = 0.882911075530934;
        u2 = 0.443506852043971;
        if strcmpi(ls,'m97')
            sc = 1.1496;       % more equal energy in subbands
        else
            sc = 1.230174104914001;  % as in JPEG-2000
        end
        ls = {'d', [p1,p1], 1; ...      % z^1, z^0
              'p', [u1,u1], 0; ...      % z^0, z^{-1}
              'd', [p2,p2], 1; ...
              'p', [u2,u2], 0; ...
               sc,    1/sc, [] };  
    else
        ls = liftwave(ls);   % if string is invalid an error occur
    end
end

if ~(iscell(ls) && (size(ls,2) == 3)) 
    disp('mylwt2: No valid lifting scheme in second argument, see help.');
    return
end
sc1 = ls{end,1};
sc2 = ls{end,2};

% do transform inplace in Y
for k = 1:nivaa
    %
    if (rem(M,2^k) == 0) % the columns
        all = 1:(2^(k-1)):N;   % all columns at this level
        lpI = 1:(2^k):M;       % low-pass rows (elements in columns)
        hpI = lpI+(2^(k-1));
        for i=1:(size(ls,1)-1)
            if (ls{i,1} == 'd')  % dual, update hp
                for j = 1:numel(ls{i,2})
                    offset = 1-j+ls{i,3};
                    if (offset == 0)
                        I = lpI;
                    elseif offset > 0
                        I = [lpI((1+offset):end),lpI(end:(-1):(end+1-offset))];
                    elseif offset < 0
                        I = [lpI((-offset):(-1):1),lpI(1:(end+offset))];
                    end
                    %
                    Y(hpI,all) = Y(hpI,all) + Y(I,all)*ls{i,2}(j);
                end
            end
            if (ls{i,1} == 'p')  % primal, update lp
                for j = 1:numel(ls{i,2})
                    offset = 1-j+ls{i,3};
                    if (offset == 0)
                        I = hpI;
                    elseif offset > 0
                        I = [hpI((1+offset):end),hpI(end:(-1):(end+1-offset))];
                    elseif offset < 0
                        I = [hpI((-offset):(-1):1),hpI(1:(end+offset))];
                    end
                    %
                    Y(lpI,all) = Y(lpI,all) + Y(I,all)*ls{i,2}(j);
                end
            end
        end
        Y(lpI,all) = Y(lpI,all)*sc1;  % scaling low-pass
        Y(hpI,all) = Y(hpI,all)*sc2;  % scaling high-pass
    end
    %
    if (rem(N,2^k) == 0) % the rows
        all = 1:(2^(k-1)):M;  % all rows at this level
        lpI = 1:(2^k):N;      % low-pass columns (elements in rows)
        hpI = lpI+(2^(k-1));
        for i=1:(size(ls,1)-1)
            if (ls{i,1} == 'd')  % dual, update hp
                for j = 1:numel(ls{i,2})
                    offset = 1-j+ls{i,3};
                    if (offset == 0)
                        I = lpI;
                    elseif offset > 0
                        I = [lpI((1+offset):end),lpI(end:(-1):(end+1-offset))];
                    elseif offset < 0
                        I = [lpI((-offset):(-1):1),lpI(1:(end+offset))];
                    end
                    %
                    Y(all,hpI) = Y(all,hpI) + Y(all,I)*ls{i,2}(j);
                end
            end
            if (ls{i,1} == 'p')  % primal, update lp
                for j = 1:numel(ls{i,2})
                    offset = 1-j+ls{i,3};
                    if (offset == 0)
                        I = hpI;
                    elseif offset > 0
                        I = [hpI((1+offset):end),hpI(end:(-1):(end+1-offset))];
                    elseif offset < 0
                        I = [hpI((-offset):(-1):1),hpI(1:(end+offset))];
                    end
                    %
                    Y(all,lpI) = Y(all,lpI) + Y(all,I)*ls{i,2}(j);
                end
            end
        end
        Y(all,lpI) = Y(all,lpI)*sc1;
        Y(all,hpI) = Y(all,hpI)*sc2;
    end
end

return;