function ut = mypred(inn, varargin)
% mypred          Prediction in matrix (image) into cell array, or inverse.
%
% The sequences will be suited for entropy coding, see estimateBits.m
% If first argument is a cell array reshape is from cell array to matrix
% else first argument is matrix and reshape is from matrix to cell array. 
%
% Used on images it works similar to CALIC but slower, a little bit larger
% neighborhood but simpler rules. 5 predictions may be used: p1 (left), p2
% (above), p3 (left above), p4 (right above) and (3*p1+2*p4)/5. Which to
% use is the one that was best for the four known neighbors. The sequences
% returned is divided based on estimated standard deviation.
% The more CALIC like implementation, calic.m, often works better on images!
% 
% use: 
%  xC = mypred(X, argName, argVal, ...);   
%  Xr = mypred(xC);      % decoding only need xC  
%------------------------------------------------------------------------
%  xC     cell array of integer sequences 
%  X      MxN matrix of integers, (an image)
%         Options:
%  'nofS' is number of sequences to split X into. Since first sequence in 
%         xC is side information, the number of sequences in xC will be 
%         one more, (nofS+1). Default: nofS = 3.
%  'v' or 'verbose' to indicate verboseness, default 0
%------------------------------------------------------------------------
% examples:
% A = double(imread('lena.bmp')) - 128; 
% xC = mypred(A, 'nofS',5, 'v',1);
% [b,r] = estimateBits(xC, 'Huff06', 2);
% Ar = mypred(xC );
%
% r1 = imageapprox(A, 'd','ex312Jan191930.mat', 'tPSNR',38, 'v',1, 'del',-10);
% dc = reshape(r1.Zdc,64,64);
% xC = mypred(dc, 'v',1);
% dcr = mypred(xC, 'v',1); if (norm(dc-dcr)==0); disp(' -- O K -- '); end;

%% default options
verbose = 0;
M = 0; 
N = 0;
w = [3,3,2,2]';   % weights
nofS = 3;
vLimRange = 4000;

%% get the options
nofOptions = nargin-1;
optionNumber = 1;
fieldNumber = 1;
while (optionNumber <= nofOptions)
    if isstruct(varargin{optionNumber})
        sOptions = varargin{optionNumber}; 
        sNames = fieldnames(sOptions);
        opName = sNames{fieldNumber};
        opVal = sOptions.(opName);
        % next option is next field or next (pair of) arguments
        fieldNumber = fieldNumber + 1;  % next field
        if (fieldNumber > numel(sNames)) 
            fieldNumber = 1;
            optionNumber = optionNumber + 1;  % next pair of options
        end
    elseif iscell(varargin{optionNumber})
        sOptions = varargin{optionNumber}; 
        opName = sOptions{fieldNumber};
        opVal = sOptions{fieldNumber+1};
        % next option is next pair in cell or next (pair of) arguments
        fieldNumber = fieldNumber + 2;  % next pair in cell
        if (fieldNumber > numel(sOptions)) 
            fieldNumber = 1;
            optionNumber = optionNumber + 1;  % next pair of options
        end
    else
        opName = varargin{optionNumber};
        opVal = varargin{optionNumber+1};
        optionNumber = optionNumber + 2;  % next pair of options
    end
    % interpret opName and opVal
    if isnumeric(inn) && strcmpi(opName,'vLim') 
        if isnumeric(opVal)
            vLim = opVal;
            nofS = numel(vLim)+1;
            disp(['mypred: Length of ''vLim'' is used to set ''nofS'',',...
                ' but values are ignored.']);
        else
            disp('mypred: illegal option ''vLim''. Option ignored.');
        end
    end
    if strcmpi(opName,'nofS') 
        if isnumeric(opVal)
            nofS = max(1,floor(abs(opVal(1)))); 
        else
            disp('mypred: illegal option ''nofS''. Option ignored.');
        end
    end
    if strcmpi(opName,'verbose') || strcmpi(opName,'v')
        if (islogical(opVal) && opVal); verbose = 1; end;
        if isnumeric(opVal); verbose = opVal(1); end;
    end
end

if isnumeric(inn)   
    %% start encoding
    X = double(inn);    % make sure we can calculate on X
    clear inn;
    if (((M ~= 0) && (M ~= size(X,1))) || ((N ~= 0) && (N ~= size(X,2))))
        disp('Actual size of A overrule input given by ''size'' or ''M'' option');
    end
    [M,N] = size(X);
    range = 2^ceil(log2(max(X(:)) - min(X(:))));
    if verbose
        disp(['mypred: from ',int2str(M),'x',int2str(N),...
              ' matrix/image (range = ',int2str(range),...
              ', min = ',int2str(min(X(:))),...
              ', max = ',int2str(max(X(:))),...
              ') to ',int2str(nofS),' sequences in cell array.']);
    end
    %
    if ((M < 6) || (N < 6)) || ((M*N) < 100)
        % may be extended in future versions
        if verbose
            disp('mypred: just use DPCM here.');
        end
        X(:,2:2:end) = flipud( X(:,2:2:end) );
        ut = { [M; log2(range); 1]; [X(1); diff(X(:))] };
        return;
    end
    %
    Y = zeros(M,N);     % start using Y as prediction
    V = zeros(M,N);     % est. of standard deviation for prediction error
    
    i=1;  % first row
    j=1; v = range/8; 
    V(i,j) = v; 
    Y(i,j) = 0;
    j=2; v = range/16; 
    V(i,j) = v; 
    Y(i,j) = X(i,j-1);
    j=3; v = (range/16 + abs(X(i,2) - X(i,1)))/2;
    V(i,j) = v; 
    Y(i,j) = X(i,j-1);
    j=4; v = (range/16 + sum(abs(X(i,2:3) - X(i,1:2))))/3;
    V(i,j) = v; 
    Y(i,j) = X(i,j-1);
    for j=5:N  
        v = sum(abs(X(i,(j-3):(j-1)) - X(i,(j-4):(j-2))))/3;
        V(i,j) = v; 
        Y(i,j) = X(i,j-1);
    end
    %
    i=2;  % second row
    j=1;  v = range/16; 
    V(i,j) = v; 
    Y(i,j) = X(i-1,j);
    j=2; 
    % ***** Block, for i=2,j=2, 2 lines identical in encode/decode ********
    x1=X(i,j-1); x2=X(i-1,j); x3=X(i-1,j-1);                            % *
    d1=abs(x2-x3); d2=abs(x1-x3); v=(d1+d2)/2;                          % *
    % *********************************************************************
    V(i,j) = v; 
    if (d1<=d2); Y(i,j) = x1; else Y(i,j) = x2; end
    for j=3:N
        % ***** Block, for i=2,j>2, 5 lines identical in encode/decode ****
        x1=X(i,j-1); x2=X(i-1,j); x3=X(i-1,j-1);                        % *
        x5=X(i,j-2); x6=X(i-1,j-2);                                     % *
        d1 = (abs(x2-x3)+abs(x5-x1))/2;                                 % *
        d2 = (3*abs(x1-x3)+2*abs(x6-x5))/5;                             % *
        v=(d1+d2)/2;                                                    % *
        % *****************************************************************
        V(i,j) = v; 
        if (d1<=d2); Y(i,j) = x1; else Y(i,j) = x2; end
    end
    %
    for i=3:M
        j=1;
        if (i==3) 
            v = (range/16 + abs(X(2,j) - X(1,j)))/2;
            V(i,j) = v; 
            Y(i,j) = X(i-1,j);
        end
        if (i==4)
        	v = (range/16 + sum(abs(X(2:3,j) - X(1:2,j))))/3;
            V(i,j) = v; 
            Y(i,j) = X(i-1,j);
        end
        if (i > 4)
        	v = sum(abs(X((i-1):(i-3),j) - X((i-2):(i-4),j)))/3;
            V(i,j) = v; 
            Y(i,j) = X(i-1,j);
        end
        j=2; 
        % ***** Block, for i>2,j=2, 5 lines identical in encode/decode ****
        x1=X(i,j-1); x2=X(i-1,j); x3=X(i-1,j-1);                        % *
        x8=X(i-2,j-1); x9=X(i-2,j);                                     % *  
        d1 = (3*abs(x2-x3)+2*abs(x8-x9))/5;                             % *
        d2 = (abs(x1-x3)+abs(x2-x9))/2;                                 % *
        v=(d1+d2)/2;                                                    % *
        % *****************************************************************
        V(i,j) = v; 
        if (d1<=d2); Y(i,j) = x1; else Y(i,j) = x2; end
        %
        for j = 3:N
            % *** Block, for i>2,j>2, lines identical in encode/decode ****
            x7=X(i-2,j-2); x8=X(i-2,j-1); x9=X(i-2,j);                  % *
            x6=X(i-1,j-2); x3=X(i-1,j-1); x2=X(i-1,j);                  % *
            x5=X(i,j-2);   x1=X(i,j-1);   % X(i,j)                      % *
            if (j==N); x4=x3; x10=x8;                                   % *
            else x4=X(i-1,j+1); x10=X(i-2,j+1); end;                    % *
            if ((j+1)>=N); x11=x2; else x11=X(i-2,j+2); end;            % *
            d = abs( ones(5,1)*[x1,x2,x3,x4] - ...                      % * 
                [ x5, x3,  x6, x2;  x3, x9,  x8, x10; ...               % *
                x6, x8,  x7, x9;  x2, x10, x9, x11; ...                 % *
                (3*x5+2*x2)/5, (3*x3+2*x10)/5, ...                      % *
                    (3*x6+2*x9)/5, (3*x2+2*x11)/5 ]);                   % *
            [temp,pNo] = min(d*w);                                      % *  
            v = mean(d(pNo,:));                                         % *
            % *************************************************************
            V(i,j) = v;
            switch pNo
                case 1
                    Y(i,j) = x1;
                case 2
                    Y(i,j) = x2;
                case 3
                    Y(i,j) = x3;
                case 4
                    Y(i,j) = x4;
                case 5
                    Y(i,j) = floor((3*x1+2*x2)/5+0.45);
            end
        end
    end
    
    xC = cell(nofS+1,1);
    % limits to store
    t = sort(V(:));
    vLim = ceil(t( floor((1:(nofS-1))*((M*N)/nofS)) )*vLimRange/range)';
    xC{1} = [M, log2(range), nofS, vLim];
    % limits to use
    vLim = [xC{1}(4:end)*range/vLimRange, inf];
    Sekv = zeros(M,N);
    for i = 1:M
        for j = 1:N
            Sekv(i,j) = find(V(i,j) <= vLim, 1);
        end
    end
    
    % make Y prediction error (and transpose it to get elements orderd by row)
    Y = (X-Y)';
    % sequence also transposed
    Sekv = Sekv';
    for i=2:(nofS+1)
        xC{i} = Y(Sekv==(i-1));
        % xC{i} = makePositive( Y(Sekv==(i-1)) );
        if verbose
            disp(['Seq. ',int2str(i),' has ',int2str(numel(xC{i})),' elements.',...
                '  min = ',int2str(min(xC{i})),'  max = ',int2str(max(xC{i})),...
                '  mean = ',num2str(mean(xC{i})), '  std = ',num2str(std(xC{i}))]);
        end
    end
    ut = xC;
    %
else  % inn is not numeric (it should be cell array)
    %% start  decoding    
    xC = inn;     
    clear inn;
    % get side info
    M = xC{1}(1);
    range = 2^xC{1}(2);
    nofS = xC{1}(3);
    if (nofS+1) ~= numel(xC);
        error('mypred: ERROR decoding xC, (nofS+1) ~= number of sequences in xC.');
    end
    nofElements = 0;
    for i = 2:(nofS+1);
        nofElements = nofElements + numel(xC{i});
    end
    N = floor(nofElements/M);
    if (N*M) ~= nofElements
        error('mypred: ERROR decoding xC, (N*M) ~= nofElements in xC.');
    end
    %        
    if verbose
        disp(['mypred: decoding from ',int2str(nofS),...
              ' sequences in xC into ',int2str(M),'x',int2str(N),...
              ' matrix/image (range = ',int2str(range),')']);
    end
    % 
    if ((M < 6) || (N < 6)) || ((M*N) < 100)
        % may be extended in future versions
        if verbose
            disp('mypred: just use DPCM here.');
        end
        X = reshape(xC{2},M,N);
        for i = 2:(M*N)
            X(i) = X(i) + X(i-1);
        end
        X(:,2:2:end) = flipud( X(:,2:2:end) );
        ut = X;
        return;
    end
    %
    % The 'normal' decoding
    vLim = [reshape(xC{1}(4:end),1,numel(xC{1})-3)*range/vLimRange, inf];
    X = zeros(M,N); % restored matrix/image (to be filled by rows!)
    xCi = zeros(nofS+1,1);
    % for i=2:(nofS+1);      xC{i} = unmakePositive( xC{i} );    end;
    %
    i=1;  % first row
    j=1; v = range/8; 
    Sekv = find(v <= vLim, 1)+1; 
    xCi(Sekv)=xCi(Sekv)+1; X(i,j) = xC{Sekv}(xCi(Sekv));
    j=2; v = range/16; 
    Sekv = find(v <= vLim, 1)+1;
    xCi(Sekv)=xCi(Sekv)+1; X(i,j) = X(i,j-1) + xC{Sekv}(xCi(Sekv));
    j=3; v = (range/16 + abs(X(i,2) - X(i,1)))/2;
    Sekv = find(v <= vLim, 1)+1;
    xCi(Sekv)=xCi(Sekv)+1; X(i,j) = X(i,j-1) + xC{Sekv}(xCi(Sekv));
    j=4; v = (range/16 + sum(abs(X(i,2:3) - X(i,1:2))))/3;
    Sekv = find(v <= vLim, 1)+1;
    xCi(Sekv)=xCi(Sekv)+1; X(i,j) = X(i,j-1) + xC{Sekv}(xCi(Sekv));
    for j=5:N   
        v = sum(abs(X(i,(j-3):(j-1)) - X(i,(j-4):(j-2))))/3;
        Sekv = find(v <= vLim, 1)+1;
        xCi(Sekv)=xCi(Sekv)+1; X(i,j) = X(i,j-1) + xC{Sekv}(xCi(Sekv));
    end
    %
    i=2;  % second row
    j=1;  v = range/16; 
    Sekv = find(v <= vLim, 1)+1;
    xCi(Sekv)=xCi(Sekv)+1; X(i,j) = X(i-1,j) + xC{Sekv}(xCi(Sekv));
    j=2; 
    % ***** Block, for i=2,j=2, 2 lines identical in encode/decode ********
    x1=X(i,j-1); x2=X(i-1,j); x3=X(i-1,j-1);                            % *
    d1=abs(x2-x3); d2=abs(x1-x3); v=(d1+d2)/2;                          % *
    % *********************************************************************
    Sekv = find(v <= vLim, 1)+1;
    xCi(Sekv)=xCi(Sekv)+1; p = xC{Sekv}(xCi(Sekv));
    if (d1<=d2); X(i,j) = x1+p; else X(i,j) = x2+p; end
    for j=3:N
        % ***** Block, for i=2,j>2, 5 lines identical in encode/decode ****
        x1=X(i,j-1); x2=X(i-1,j); x3=X(i-1,j-1);                        % *
        x5=X(i,j-2); x6=X(i-1,j-2);                                     % *
        d1 = (abs(x2-x3)+abs(x5-x1))/2;                                 % *
        d2 = (3*abs(x1-x3)+2*abs(x6-x5))/5;                             % *
        v=(d1+d2)/2;                                                    % *
        % *****************************************************************
        Sekv = find(v <= vLim, 1)+1;
        xCi(Sekv)=xCi(Sekv)+1; p = xC{Sekv}(xCi(Sekv));
        if (d1<=d2); X(i,j) = x1+p; else X(i,j) = x2+p; end
    end
    %
    for i=3:M
        j=1;
        if (i==3) 
            v = (range/16 + abs(X(2,j) - X(1,j)))/2;
            Sekv = find(v <= vLim, 1)+1;
            xCi(Sekv)=xCi(Sekv)+1; X(i,j) = X(i-1,j) + xC{Sekv}(xCi(Sekv));
        end
        if (i==4)
        	v = (range/16 + sum(abs(X(2:3,j) - X(1:2,j))))/3;
            Sekv = find(v <= vLim, 1)+1;
            xCi(Sekv)=xCi(Sekv)+1; X(i,j) = X(i-1,j) + xC{Sekv}(xCi(Sekv));
        end
        if (i > 4)
        	v = sum(abs(X((i-1):(i-3),j) - X((i-2):(i-4),j)))/3;
            Sekv = find(v <= vLim, 1)+1;
            xCi(Sekv)=xCi(Sekv)+1; X(i,j) = X(i-1,j) + xC{Sekv}(xCi(Sekv));
        end
        j=2; 
        % ***** Block, for i>2,j=2, 5 lines identical in encode/decode ****
        x1=X(i,j-1); x2=X(i-1,j); x3=X(i-1,j-1);                        % *
        x8=X(i-2,j-1); x9=X(i-2,j);                                     % *  
        d1 = (3*abs(x2-x3)+2*abs(x8-x9))/5;                             % *
        d2 = (abs(x1-x3)+abs(x2-x9))/2;                                 % *
        v=(d1+d2)/2;                                                    % *
        % *****************************************************************
        Sekv = find(v <= vLim, 1)+1;
        xCi(Sekv)=xCi(Sekv)+1; p = xC{Sekv}(xCi(Sekv));
        if (d1<=d2); X(i,j) = x1+p; else X(i,j) = x2+p; end
        %
        for j=3:N  % column 3 and more
            %
            % *** Block, for i>2,j>2, lines identical in encode/decode ****
            x7=X(i-2,j-2); x8=X(i-2,j-1); x9=X(i-2,j);                  % *
            x6=X(i-1,j-2); x3=X(i-1,j-1); x2=X(i-1,j);                  % *
            x5=X(i,j-2);   x1=X(i,j-1);   % X(i,j)                      % *
            if (j==N); x4=x3; x10=x8;                                   % *
            else x4=X(i-1,j+1); x10=X(i-2,j+1); end;                    % *
            if ((j+1)>=N); x11=x2; else x11=X(i-2,j+2); end;            % *
            d = abs( ones(5,1)*[x1,x2,x3,x4] - ...                      % * 
                [ x5, x3,  x6, x2;  x3, x9,  x8, x10; ...               % *
                x6, x8,  x7, x9;  x2, x10, x9, x11; ...                 % *
                (3*x5+2*x2)/5, (3*x3+2*x10)/5, ...                      % *
                    (3*x6+2*x9)/5, (3*x2+2*x11)/5 ]);                   % *
            [temp,pNo] = min(d*w);                                      % *  
            v = mean(d(pNo,:));                                         % *
            % *************************************************************
            Sekv = find(v <= vLim, 1)+1;
            xCi(Sekv)=xCi(Sekv)+1; p = xC{Sekv}(xCi(Sekv));
            switch pNo
                case 1
                    X(i,j) = x1+p;
                case 2
                    X(i,j) = x2+p;
                case 3
                    X(i,j) = x3+p;
                case 4
                    X(i,j) = x4+p;
                case 5
                    X(i,j) = floor((3*x1+2*x2)/5+0.45)+p;
            end
        end
    end
    %  
    ut = X;
end

return;

% function y = makePositive(x)
% %    ... -2  -1  0  1  2 ... to
% %         4   2  0  1  3 ...
% y = x(:);
% y(y>0) = 2*y(y>0)-1;
% y(y<0) = -2*y(y<0);
% return;
% 
% function x = unmakePositive(y)
% x = y(:);
% x(mod(x,2)==0) = -x(mod(x,2)==0)/2;
% x(x>0) = (x(x>0)+1)/2;
% return;

