function  ut = myreshape(inn, varargin)
% myreshape       Reshape between matrix and cell array of sequences.
% The sequences should be fitted for entropy coding, see estimateBits.m
% If first argument is a cell array reshape is from cell array to matrix
% else first argument should be matrix and reshape is from matrix to cell
% array.
% see also: Mat2Vec in Comp toolbox (made by Karl Skretting)
% 
% use: 
%  xC = myreshape(W, argName, argVal, ...);   
%  W  = myreshape(xC, argName, argVal, ...);       
%------------------------------------------------------------------------
%  xC     cell array of integer sequences 
%         size of W and method are stored as xC{1}, i.e [K,L,met]
%         more options may also be stored in xC{1}
%  W      KxL matrix of integers, i.e. quantized transform coeffients
%         Options:
%  'm' or 'method' gives the method to use, default 1
%         1 : my own variant of End-Of-Block (EOB) coding into several
%             sequences. Suited for quantized transform (or wavelet)
%             coefficients with 3 levels, i.e. size(W,1)=64.
%             Additional arguments may be: 'sortRows', 'sDC'
%             (Note, algorithm is slow for large W)
%         2 : index + value sequences, suited for sparse matrices with no
%             DC component and size(W,1) might as well be large (ex. 440).
%         3 : is like method 1, with row sorted to match the order in 2D
%             wavelet, i.e. myim2col (with wavelet transform mylwt2) or as
%             lwt2 (inplace) followd by im2col. K = 16, 64, 256, 1024
%             This is usually only better than method 1 for K = 256 or 1024
%  'sortRows' for method 1. true (default) or false. 
%  'sDC' is size of DC image, [Mdc, Ndc]. If this is correct given such
%         that (Mdc*Ndc == L), prediction of DC component (first row) will
%         be by mypred.m, otherwise DPCM will be done.
%  'largeLimit' for method 2. A number (ex. 100-1000), default 440.
%  'v' or 'verbose' to indicate verboseness, default 0
%------------------------------------------------------------------------
% examples:
% A = double(imread('lena.bmp')) - 128; peak = 255;
% X = myim2col(A, 't','m79', 's',8);
% del = 10; W = UniQuant(X, del, del, 2000); 
% xCold = Mat2Vec(W, 5); [bold , rold] = estimateBits(xCold, 'Huff06', 2);
% xCnew = myreshape(W);  [bnew , rnew] = estimateBits(xCnew, 'Huff06', 2);
% Wr = myreshape(xCnew); if (norm(W-Wr) == 0); disp(' -- O K --'); end;
% xCv11 = myreshape(W, 'm',1, 'sDC',size(A)/8);   % version 1.1 
% [bv11 , rv11] = estimateBits(xCv11, 'Huff06', 2);
% Wr = myreshape(xCv11); if (norm(W-Wr) == 0); disp(' -- O K --'); end;
% Ar = mycol2im(UniQuant(Wr,del,del), 't','m79', 'imsize',size(A), 's',8);
% R = A - Ar; PSNR = 10*log10(((numel(A)*peak^2)/sum(sum(R.*R))));


%------------------------------------------------------------------------
% Copyright (c) 2010.  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  12.02.2010  KS: m-file made
% Ver. 1.1  22.02.2010  KS: Added 'sDC', may do better prediction in DC
% Ver. 1.2  04.03.2010  KS: Added method 3
%------------------------------------------------------------------------


%% default options
method = 1;
sortrows = true;
largeLimit = 440;
verbose = 0;
Mdc = 0;
Ndc = 0;
nofSdc = 3;   % may be used as nofS in mypred
m3variant = 1;

%% 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 strcmpi(opName,'method') || strcmpi(opName,'m')
        if isnumeric(opVal)
            method = max(1, floor(opVal(1)));
        else
            disp('myreshape: illegal option method. We ignore it.');
        end
    end
    if strcmpi(opName,'sortRows') 
        if islogical(opVal); sortrows = opVal(1); end;
        if isnumeric(opVal); sortrows = (opVal(1) ~= 0); end;
    end
    if strcmpi(opName,'sDC') 
        if isnumeric(opVal)
            if (numel(opVal) >= 2)
                Mdc = opVal(1);
                Ndc = opVal(2);
            else
                Mdc = opVal(1);
            end
        end
    end
    if strcmpi(opName,'largeLimit') 
        if isnumeric(opVal)
            largeLimit = max(5, floor(opVal(1)));
        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

%% start of function
if isnumeric(inn)   % encoding
    W = double(inn);    % make sure we can calculate on W
    clear inn;
    [K,L] = size(W);
    if verbose
        disp(['myreshape: from ',int2str(K),'x',int2str(L),...
              ' matrix to cell array (method = ',int2str(method),...
              ', sortrows = ',int2str(sortrows),...
              ', largeLimit = ',int2str(largeLimit),').']);
    end
    %
    if (method == 1)  % my own End-Of-Block coding into 5 (or 6) sequences
        % seq 1 : info on W and method
        % seq 2 : DC component
        % seq 3 : All other symbols, except after 0 and 1 (EOB and 0)
        % seq 4-6 : The symbols after 1 (i.e. 0) (no 0 or 1)
        % coded by RLEof1( makePositive(flipud(x(i2:end))) );
        % seq 7 : order (only if sortrows)
        if (Mdc > 0) && (Ndc == 0); Ndc = floor(L/Mdc); end;
        if (Mdc*Ndc == L) % mypred used to predict DC component
            xC = cell(6+sortrows+nofSdc,1);
        else  % DPCM used to predict DC component
            xC = cell(6+sortrows,1);
        end
        xC{1} = [K, L, method, sortrows];
        % next lines should work fine both for sparse and full W
        [wk,wl,w] = find(W);
        nnz = numel(w);        % number of non-zeros
        nzRows = hist(wk,1:K);
        if sortrows
            [temp, order] = sort(nzRows, 'descend');
            [temp,reverseorder] = sort(order);
            xC{7} = order(:);
            W = W(order,:);    % sort rows in W
            nzRows = nzRows(order);  % and non-zeros in each row
            for i=1:nnz; wk(i) = reverseorder(wk(i)); end; % and wk
            % else
            % order = 1:K;
            % reverseorder = 1:K;
        end
        %
        % how to make next section effective is difficult
        EOB = 0;   % the End-Of-Block symbol
        posEOB = zeros(1,L);   % index of last non-zero in each column
        for i=1:nnz
            posEOB(wl(i)) = max(posEOB(wl(i)), wk(i));
        end
        %
        % we may do more advanced prediction of DC, xC{2} may be changed
        if (Mdc*Ndc == L)
            if verbose
                disp('myreshape: use mypred for prediction in DC.');
            end
            %
            dc = full(W(1,:));
            dc(dc >= 0) = dc(dc >= 0) + 1;
            dc(posEOB == 0) = EOB;   % = 0, two symbols for 0, EOB and 1
            xCdc = mypred( reshape(dc, Mdc, Ndc), 'nofS',nofSdc);
            xC{2} = xCdc{1};   % short sequence with side info
            for i=2:numel(xCdc)
                xC{5+sortrows+i} = xCdc{i};
            end
        else   % DPCM for DC
            dc = full(W(1,:));
            dc(dc >= 0) = dc(dc >= 0) + 1;
            dc(posEOB == 0) = EOB;   % = 0, two symbols for 0, EOB and 1
            xC{2} = thisDPCM( dc );
        end
        %
        numelx1  = nnz - nzRows(K);    
        numelx2  = sum(posEOB) - numelx1 - nzRows(K);  
        x1 = ones(numelx1, 1); i1 = 0;  % elements after non-zeros
        x2 = ones(numelx2, 1); i2 = 0;  % elements after zeros
        for i=find(dc)
            if (posEOB(i) == 1)
                i1=i1+1; x1(i1) = EOB;
            else
                last = W(1,i);
                for k=2:posEOB(i)    % to last non-zero in column
                    this = W(k,i);
                    if (last == 0)   % W(k,i) --> x2
                        i2 = i2 + 1;     % if W(k,i)==0, we let x2 be 1
                        if (this > 0)
                            x2(i2) = W(k,i) + 1;
                        elseif (this < 0)    
                            x2(i2) = W(k,i);
                        end
                    else
                        i1 = i1 + 1;     % if W(k,i)==0, we let x1 be 1
                        if (this > 0)
                            x1(i1) = W(k,i) + 1;
                        elseif (this < 0)    
                            x1(i1) = W(k,i);
                        end
                    end
                    last = this;
                end
                if (posEOB(i) < K) % need to write EOB, 
                    i1=i1+1; x1(i1) = EOB; % to x1 since last is non-zero!
                end
            end
        end
        % disp([numelx1, i1, numelx2, i2]);
        % 
        xC{3} = makePositive(x1); % usually best
        [xC{4}, xC{5}, xC{6}] = RLEof1( makePositive(x2) );
        %
    end
    %
    if (method == 2)  % coding of sparse W (non-zeros randomly  placed) 
        % into 4 sequences: info, values, indexes, large-indexes
        % A 'largeLimit' is used, values in seq 3 >= largeLimit are copied
        % to seq 4 and largeLimit is written into these positions in seq 3
        % seq 1 : info on W and method
        % seq 2 : the non-zero values
        % seq 3 : number of zeros before each value and after last
        % seq 4 : large values from sequence 2
        xC = cell(4,1);
        xC{1} = [K, L, method, largeLimit];
        [index, temp, value] = find(W(:));
        index(2:end) = diff(index);   % DPCM of index
        I = find(index >= largeLimit);
        large = index(I);
        index(I) = largeLimit;
        xC{2} = value;
        xC{3} = index-1;
        xC{4} = large - (largeLimit-32);  % code large numbers ok in Huff06
    end
    %
    if (method == 3)  % coding of subbbands (from 2-D wavelet)
        if (Mdc > 0) && (Ndc == 0); Ndc = floor(L/Mdc); end;
        % K, or size(W,1), should be 16, 64, 256, 1024,
        level = floor(log2(K)/2);
        if (4^level ~= K)
            error('myreshape, method 3: K ~= 4^level.');
        end
        b = band2Dwave(level);
        [temp,rowOrder] = sort(b(:),'descend');
        % [temp,invOrder] = sort(rowOrder);
        if m3variant == 1 % the easy obvious way
            xC = myreshape(W(rowOrder,:), 'm',1, 'sDC',[Mdc,Ndc], 'sortRows',0);
            xC{1}(3) = method;
        else  % Some more involved ways, use low (including dc) and high 
            disp('myreshape method 3, these variants are not ready nor good.');
            Ihi = rowOrder( (K/4+1):end );
            Ilo = rowOrder( 1:(K/4) );
            xClo = myreshape(W(Ilo,:), 'm',1, 'sDC',[Mdc,Ndc], 'sortRows',0);
            if m3variant == 2 %  alternative 2:  often GOOD
                x = W(Ihi,:);  [I,J,val] = find(x(:));
                xC = {[I(1);diff(I)]-1, val};
            elseif m3variant == 3 %  alternative 3   (ser ikke like godt ut)
                nnz = sum( W(Ilo,:) ~= 0 );
                x = W(Ihi,(nnz <= 5));  [I1,J,x1] = find(x(:)); I1 = [I1(1); diff(I1)]-1;
                x = W(Ihi,(nnz > 5) & (nnz <= 20));  [I2,J,x2] = find(x(:)); I2 = [I2(1); diff(I2)]-1;
                x = W(Ihi,(nnz > 20));  [I3,J,x3] = find(x(:)); I3 = [I3(1); diff(I3)]-1;
                xC = {I1, I2, I3, x1, x2, x3};
            elseif m3variant == 4 % alternativ 4, heller ikke så godt
                nnz = sum( W(Ihi,:) ~= 0 );
                x = W(Ihi,(nnz > 0) & (nnz <= 4));  [I1,J,x1] = find(x(:)); I1 = [I1(1); diff(I1)]-1;
                x = W(Ihi,(nnz > 4) & (nnz <= 16));  [I2,J,x2] = find(x(:)); I2 = [I2(1); diff(I2)]-1;
                x = W(Ihi,(nnz > 16));  [I3,J,x3] = find(x(:)); I3 = [I3(1); diff(I3)]-1;
                xC = {[nnz(1), diff(nnz)], x1, x2, x3, I1, I2, I3};
            elseif m3variant == 5 % alternativ 5, best but complex
                nnz = sum( W(Ihi,:) ~= 0 );
                nz = (nnz > 0);   % this is a binary 32x32 image (1024 bit)
                x = W(Ihi,nz);  [I,J,x] = find(x(:)); I = [I(1); diff(I)]-1;
                Ilarge = I(I>2);
                J = find(I<=2); 
                Ismall = I(J);
                J = [J(1);diff(J)]-1;
                xChi = {binpred(nz), x, Ilarge, Ismall, J};
                xC = cell( numel(xChi)+numel(xClo),1 );
                for i=1:numel(xClo); xC{i}=xClo{i}; end;
                for j=1:numel(xChi); xC{i+j}=xChi{j}; end;
            end
        end
        % 
    end
    %
    ut = xC;
elseif iscell(inn)   % decoding
    xC = inn;     
    clear inn;
    K = xC{1}(1);      % the first three elements have a special use
    L = xC{1}(2);
    method = xC{1}(3);
    W = zeros(K,L);   % or sparse(K, L);    
    if verbose
        disp(['myreshape: from cell array (',int2str(numel(xC)),...
              ' seqences) to ',int2str(K),'x',int2str(L),...
              ' matrix (method = ',int2str(method),').']);
    end
    if (method == 1)  % my own End-Of-Block coding into 5 (or 6) sequences
        sortrows = xC{1}(4);  % 
        if (numel(xC{2}) ~= L)   % mypred was used to predict DC
            xCdc = cell( (numel(xC)-5-sortrows), 1 );
            xCdc{1} = xC{2}; 
            for i=2:numel(xCdc)
                xCdc{i} = xC{5+sortrows+i};
            end
            dcr = mypred(xCdc);
            W(1,:) = dcr(:);
        else
            W(1,:) = invThisDPCM( xC{2} ); 
        end
        x1 = unmakePositive( xC{3} ); 
        x2 = unmakePositive( invRLEof1( xC{4}, xC{5}, xC{6} ) );
        i1 = 0; i2 = 0;
        for i=1:L
            for k=2:K
                if W(k-1,i) == 0
                    break;  % EOB
                elseif W(k-1,i) == 1
                    i2 = i2+1; W(k,i) = x2(i2);
                else % not (0 or 1)
                    i1 = i1+1; W(k,i) = x1(i1);
                end    
            end
        end
        % 
        W(W > 0) = W(W > 0) - 1;   % remove EOB-symbol (0)
        if sortrows
            W(xC{7}(:),:) = W;
        end
    end
    if (method == 2)  % coding of sparse W (non-zeros randomly  placed) 
        try
            largeLimit = xC{1}(4);
            index = xC{3} + 1;
            index(index == largeLimit) = xC{4} + (largeLimit-32);
            for i=2:numel(index)
                index(i) = index(i-1)+index(i);
            end
            W(index) = xC{2};
        catch ME
            disp('myreshape: ERROR decoding using method 2.');
            disp('Perhaps argument 1 (xC) was not made by myreshape and method 2.');
            disp(ME.message);
        end
    end
    if (method == 3)  % 2D wavelet order of rows in W
        level = floor(log2(K)/2);
        b = band2Dwave(level);
        [temp,rowOrder] = sort(b(:),'descend');
        [temp,invOrder] = sort(rowOrder);
        xC{1}(3) = 1;
        W = myreshape(xC);
        W = W(invOrder,:);
    end
    %
    ut = W;   % return matrix W (in variable xC)
else
    disp(['myreshape: input argument is ',class(inn),', see help.']);
    ut = -1;
end

return   % end of myreshape.m

% ********************************************* S U B F U N C T I O N S *******************

function y = thisDPCM(x)
x = x(:);
y  = [x(1); x(2:end)-x(1:(end-1))];
return;

function x = invThisDPCM(y)
x = zeros(numel(y),1);
x(1) = y(1);
for i=2:numel(y);
    x(i) = x(i-1) + y(i);
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;

%
% RLE is a more complicated function, splitting x into three sequences
%
function [x1,x2,x3] = RLEof1(x)
% Run-Length-Encoding of ones (minimum value in sequence)
% 0 is RLE symbol indicating 2 or more ones
% x1 is RLE encoded sequence for symbols following a 2 or greater
% x2 is RLE encoded sequence for symbols following 0 or 1
% x3 is number of ones for each 0 (RLE symbol)
if (min(x) < 1)
    disp('myreshape.RLEof1: sequence x not as it should be.');
end
N = numel(x);
x1 = zeros(N,1);
x2 = zeros(N,1);
x3 = zeros(N,1);
n1 = 0;
n2 = 0;
n3 = 0;
zc = 0;  % count of ones
for n=1:N
    if x(n) == 1
        zc = zc+1;
    else
        if (zc > 0)  
            if (zc > 1)
                n1 = n1+1; x1(n1) = 0;
                n3 = n3+1; x3(n3) = zc;
            else  %  (zc == 1)
                n1 = n1+1; x1(n1) = 1;
            end
            n2 = n2+1; x2(n2) = x(n);
            zc = 0;  
        else
            n1 = n1+1; x1(n1) = x(n);
        end
    end
end
if (zc > 1)
    n1 = n1+1; x1(n1) = 0;
    n3 = n3+1; x3(n3) = zc;
elseif (zc == 1)
    n1 = n1+1; x1(n1) = 1;
end
x1 = x1(1:n1);
x2 = x2(1:n2)-2;  
x3 = x3(1:n3)-2;
return

function x = invRLEof1(x1,x2,x3)
N = numel(x1)+numel(x2)+sum(x3)+numel(x3);
x = ones(N,1);
n = 0;
n2 = 0;
n3 = 0;
for n1=1:numel(x1)
    if (x1(n1) < 2)
        if (x1(n1) == 1)
            n = n+1;   % x(n) is already 1
        else % == 0
            n3 = n3+1; n = n+x3(n3)+2; 
        end
        if (n < N)
            n = n+1; n2 = n2+1; x(n) = x2(n2)+2;
        end
    else
        n = n+1; x(n) = x1(n1);  % >= 2
    end
end
return;

% a 2D wavelet, as in mylwt2 or lwt2, divides image into block of size NxN
% where N = 2^level. The returned value, b, gives the subband number for
% each of the elementes in the block, b is NxN.
% Here band is numbered from high-high to low-low-...-low
function b = band2Dwave(level)
N = 2^level;
A = [4,2;3,1];   % the first order
K = N/2;
b = kron(ones(K), A);
A = A-1;
step = 2;
while (K > 1)
    K = K/2;
    i = 1:step:N;
    b(i,i) = b(i,i) + kron(ones(K), A);
    step = step*2;
end
return;

% B er logisk array, x er hvor mange riktige prediksjoner det er foran
% hver som er feil. Antall element i x er da lik antall element feil
% prediktert. Enkel prediksjonsregel
% Aritmetisk koding (med ulike context) vil generelt være bedre, i stedet
% for å prediktere verdi, estimeres sannsynlighet for '1'
function x = binpred(B)
x = zeros(numel(B),1);
xv=0; xi=0;
p=0; 
[M,N] = size(B);
i=1; j=1;
if (B(i,j)==p); xv=xv+1; else xi=xi+1; x(xi)=xv; xv=0; end;
for j = 2:N
    p = B(i,j-1);
    if (B(i,j)==p); xv=xv+1; else xi=xi+1; x(xi)=xv; xv=0; end;
end
for i = 2:M
    j = 1;
    p = B(i-1,j);
    if (B(i,j)==p); xv=xv+1; else xi=xi+1; x(xi)=xv; xv=0; end;
    for j = 2:N
        a=B(i-1,j); b=B(i,j-1); 
        if (a==b)
            p=a;
        else
            if (i>2); a=B(i-2,j); end;
            if (j>2); b=B(i,j-2); end;
            if (a==b); p=a; else p=B(i-1,j-1); end;
        end
        if (B(i,j)==p); xv=xv+1; else xi=xi+1; x(xi)=xv; xv=0; end;
    end
end
x = x(1:xi);
return;

% %% TEST code ******************************************* T E S T *****
% % special test: ut = myreshape(25);
% if (nargin == 1) && (numel(inn) == 1)
%     del = inn;
%     A = double(imread('lena.bmp')) - 128;
%     X = myim2col(A, 't','dct', 's',8);
%     W = UniQuant(X, del, del, 2000);
%     xC1 = Mat2Vec(W, 5);
%     xC2 = myreshape(W, 'sortRows',1);   % important to sort!
%     ut = cell(numel(xC1)+numel(xC2)+1,1);
%     disp('Test of myreshape, first sequences from Mat2Vec.');
%     for i=1:numel(xC1)
%         ut{i} = checkSeq(xC1{i}, 1);
%     end
%     disp([blanks(30),'sum all = ',...
%         int2str( ut{1}.b4+ut{2}.b4+ut{3}.b4 )]);
%     disp('And sequences from myreshape (except first).');
%     for i=2:numel(xC2)
%         ut{numel(xC1)+i-1} = checkSeq(xC2{i}, 0);
%         disp([' ',int2str(i),': ',ut{numel(xC1)+i-1}.text]);
%     end
%     disp([blanks(30),'sum 4, 5 and 6 = ',...
%         int2str( ut{6}.b4+ut{7}.b4+ut{8}.b4 )]);
%     disp([blanks(30),'sum 2 to 6     = ',...
%         int2str( ut{4}.b4 + ut{5}.b4 + ut{6}.b4 + ...
%         ut{7}.b4 + ut{8}.b4 + ut{9}.b4 )]);
%     disp('   ---  ');
%     Wr = myreshape(xC2);
%     ut{numel(xC1)+numel(xC2)} = Wr;
%     ut{numel(xC1)+numel(xC2)+1} = xC2{4};
%     disp(['Check inverse of myreshape ',num2str(norm(W-Wr))]);  % should be zero
%     %     disp('Look more at DC from myreshape, xC2{3}.');
%     %     rC = checkSubfun(xC2{3});
%     %     for i=(2:6)
%     %         r = checkSeq(rC{i}, 0);
%     %         disp(['Seq ',int2str(i),' (err=',int2str(rC{i+5}),') ',r.text]);
%     %     end
%     return;
% end
% 
% function r = checkSubfun(x)
% y1 = x(:);
% y2 = makePositive(y1);
% y3 = thisDPCM(y1);
% y4 = makePositive(y3);
% y5 = DDPCM(y1);
% y6 = makePositive(y5);
% y2err = sum(abs(y1 - unmakePositive(y2)));
% y3err = sum(abs(y1 - invThisDPCM(y3)));
% y4err = sum(abs(y3 - unmakePositive(y4)));
% y5err = sum(abs(y1 - invDDPCM(y5)));
% y6err = sum(abs(y5 - unmakePositive(y6)));
% r = {y1,y2,y3,y4,y5,y6, y2err, y3err, y4err, y5err, y6err};
% return;
% 
% function res = checkSeq(x, verbose)
% N = numel(x);               % antall element
% if (N<5)
%     res = struct('text', ['Short sequence N=',int2str(N),], 'N',N);
% else
%     x1 = min(x);
%     x2 = max(x);
%     counts = hist(x, x1:x2);    % antall av hvert symbol
%     [cmax,imax] = max(counts);
%     B = nonzeros(counts);      % antall av hvert brukt symbol
%     a = numel(B);              % antall brukte symbol
%     H = log2(N) - sum(B.*log2(B))/N;  % zero order entropy
%     NH = ceil(N*H);
%     b1 = NH + ceil(a*log2(a));
%     HL = HuffLen(counts);
%     b2 = sum(HL.*counts);
%     HLlen = HuffTabLen(HL);
%     b3 = b2 + HLlen;
%     y = Huff06({x},5);
%     b4 = 8*numel(y);
%     res = struct('text', [int2str(N),' [',int2str(x1),',',int2str(x2),']  ',...
%         int2str(x1+imax-1),':',int2str(cmax),'  ',...
%         sprintf('%5.2f',H),' -> ',int2str(NH),', ',...
%         int2str(b1),'  ',...
%         int2str(b2),'+',int2str(HLlen),'=',int2str(b3),'  ',...
%         int2str(b4)], ...
%         'N',N, 'min',x1, 'max',x2, 'counts',counts, 'maxcount',cmax, ...
%         'maxsymbol',x1+imax-1, 'entropy',H, 'b1',b1, ...
%         'HuffmanLengths',HL, 'b2',b2, 'HtabLen',HLlen', 'b3',b3, ...
%         'y',y, 'b4',b4 );
% end
% if verbose; disp(res.text); end;
% return
% 
% function y = DDPCM(x)
% x = x(:);
% y  = [x(1); x(2)-x(1); x(3:end)-2*x(2:(end-1))+x(1:(end-2))];
% return;
% 
% function x = invDDPCM(y)
% x = zeros(numel(y),1);
% x(1) = y(1);
% x(2) = y(2) + x(1);
% for i=3:numel(y);
%     x(i) = y(i) + 2*x(i-1) - x(i-2);
% end
% return
% 

% % test array
% K = 8; L = 30;  Z = floor(1.25*rand(K,L)); 
% W = Z; I = find(Z);
% W(I) = floor(7*randn(numel(I),1)+2);
% %
% [wk,wl,w] = find(W);
% nnz = numel(w);        % number of non-zeros
% nzRows = hist(wk,1:K);
% [temp, order] = sort(nzRows, 'descend');
% [temp,reverseorder] = sort(order);
% for i=1:nnz; wk(i) = reverseorder(wk(i)); end;
% %
% posEOB = zeros(1,L);   % index of last non-zero in each column
% for i=1:nnz
%     posEOB(wl(i)) = max(posEOB(wl(i)), wk(i));
% end
% %
% dc = full(W(order(1),:));            % assume first row is DC
% dc(dc >= 0) = dc(dc >= 0)+1;  % make room for EOB (=0)
% dc(posEOB==0) = 0;          % some columns be be only zeros
% %
% x1 = ones(nnz - nzRows(order(K)), 1);   % the elements succeeding a non-zeros
% x2 = ones(sum(posEOB)-numel(x1), 1);    % elemets after a zeros
% i1 = 0;
% i2 = 0;
% col = wl(1);
% first = 1;
% last = 0;
% for i=1:nnz
%     if (wl(i) > col) % next colum
%         % ind = first:last is previous column
%         if (first <= last)
%             [temp,ii] = sort(wk(first:last));
%             values = w(first-1+ii);
%             wki = wk(ii(wk(ii)>1));  %
%             % and what to do here must be correct
%         end
%         first = i;
%         col = wl(i);
%     end
%     last = i;
% end
% 
