function  A = myilwt2(Y, ls, nivaa)
% myilwt2         My variant of ilwt2(..), use as ilwt2 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: 
%   Ar = myilwt2(Y, ls, level);   
%----------------------------------------------------------------------
% Y     the coefficients, size MxN, both M and N should have 
%       2^level (typically 8, 16 or 32) as a factor.
% A     the reconstructed image (matrix of doubles), same size as Y.
% 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'
%----------------------------------------------------------------------

[M,N] = size(Y);
K = 2^nivaa;

% check arguments
if ~isnumeric(Y)
    disp('myilwt2: The coefficient matrix Y is not numeric.');
    A = Y;
    return;
end
if ~((mod(M,K)==0) && (mod(N,K)==0))
    disp(['myilwt2: Size of coefficient matrix Y has not ',...
           int2str(K),' as factor.']);
    A = Y;
    return;
end

if ischar(ls) 
    % 
    i = strfind(ls,'79');    % allow for 97 to be swapped in ls
    if (numel(i) == 1); ls(i:(i+1)) = '97'; end;
    i = strfind(ls,'7.9');    % allow for 9.7 to be swapped 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('myilwt2: 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 = nivaa:(-1):1
    %
    if (rem(N,2^k) == 0) % the rows
        all = 1:(2^(k-1)):M;
        lpI = 1:(2^k):N;
        hpI = lpI+(2^(k-1));
        Y(all,lpI) = Y(all,lpI)/sc1;
        Y(all,hpI) = Y(all,hpI)/sc2;
        for i=(size(ls,1)-1):(-1):1
            if (ls{i,1} == 'd')  % dual, update hp
                for j = numel(ls{i,2}):(-1):1
                    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 = numel(ls{i,2}):(-1):1
                    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
    end
    %
    if (rem(M,2^k) == 0) % the columns
        all = 1:(2^(k-1)):N;
        lpI = 1:(2^k):M;
        hpI = lpI+(2^(k-1));
        Y(lpI,all) = Y(lpI,all)/sc1;
        Y(hpI,all) = Y(hpI,all)/sc2;
        for i=(size(ls,1)-1):(-1):1
            if (ls{i,1} == 'd')  % dual, update hp
                for j = numel(ls{i,2}):(-1):1
                    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 = numel(ls{i,2}):(-1):1
                    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
    end
end

A = Y;

return;