function s = EMP8_argin(varargin)
% EMP8_argin       Arrange input arguments in a struct for use in EMP8_xxx.m 
%                  EMP8: Energy Minimization using Potts model and 8 neighbors
% 
% This function may use:
%   LIS_getR(..)       to get the data term matrix
%   LIS_smooth(..)     to get input labeling L by Gaussian smoothing of R
%   teximage(..)       to get the true labeling, or image in Uim 
%   EMP8_cost_mex(..)  to find cost for input labeling
% The result s may be used as input in EMP8_cost, EMP8_erode, EMP8_swap
%
% Use of function:
% ----------------
%   s = EMP8_argin(R, L);      % all options as default values
%   s = EMP8_argin(R, L, 'option',value, ...);
%   s = EMP8_argin(R, L, Ltrue, 'option',value, ...);
%   s = EMP8_argin(s, 'option',value, ...);    % change/add options
%
%   s   a struct containing fields
%      Ri   MIJx1 (int16) array
%      Li,Lt   IJx1 (uint8) array (input and true labeling)
%      ...   and many more fields
%   If error occurs (input not as expected), s is just the error message
%
%  Input arguments:
%     R is data cost matrix, size IxJxM, class of R should be (double)
%       or (int16) then it is interpreted as Ri (int16) of size MIJx1 (or MxIxJ)
%       or an integer (1-13) for a stored matrix, see LIS_getR.m
%     L is input label matrix, size IxJ, or Li (uint8) of size IJx1 (or IxJ)
%       0 may be used for minimizing data term, or a single number interpreted as
%       sigma to be used as: L = LIS_smooth(R, sigma);
%     Ltrue is true labeling (if available), may be an integer 1-13, 
%       or it may be skipped as it is only used when displaying results
%       see: teximage.m
%  Additional arguments (options), given as pairs:, ..., 'option',value, ... 
%   'lambda'  the parameter lambda to use for Etot (and Ctot), default 2
%   'bF8percent'  how many percent of bF4 the bF8 value should be, default 70
%             this is used if bF4 is not given (bF4 < 0)
%   'bF4'     border factor for 4-neighbors pixels, default -1
%             if not given (< 0) the bF4 value is found from data and lambda
%   'bF8'     if bF4 is given bF8 should also be given
%
%   'swUpq'   1: use only bF4 and bF8 (which may be given by lambda and bF8percent)
%             2: also use Upq which gives a faktor for each pixel to 4 neighbors
%             (the three to the left and the one above, in order: AL, L, BL, A)
%             3: also use Uim and Utab (256 integers)
%   'Upq'    -Is a matrix with all U-factors in the Potts (8-neighbor) model.
%             Since these factors are symmetric only 4 are needed for each
%             pixel, thus size is IxJx4: (i,j,1) is for pixel left and
%             above, (i,j,2) is for pixel left, (i,j,3) is for pixel left
%             and below, (i,j,4) is for pixel above. The normal range for
%             these factors is from 0 to 1, and type is double.
%            -Alternative is to give Upq as integers, (int16) and size
%             4IJx1 (or 4xIxJ). The range should now be 0:PU1, where PU1 is
%             defined in EMP8_common.h (and is 4096)
%            -ex: Uim = teximage(104); [I,J] = size(Uim); Upq = zeros(I,J,4);
%             Upq(:,:,1) = Uim - [Uim(1,:); Uim(2:I,1),Uim(1:(I-1),1:(J-1))];
%             Upq(:,:,2) = Uim - [Uim(:,1),Uim(:,1:(J-1))];
%             Upq(:,:,3) = Uim - [Uim(:,1),[Uim(2:I,1:(J-1));Uim(I,2:J)]];
%             Upq(:,:,4) = Uim - [Uim(1,:);Uim(1:(I-1),:)];
%             Upq = abs(Upq); Upq = 1 - Upq/max(Upq(:));
%             tim = 4;  sigma = 5;  
%             s = EMP8_argin(tim, sigma, tim, 'swUpq',2, 'Upq',Upq, 'v',2);
%   'Uim'    -Is a matrix with an 'image' of size IxJ, each element an
%             integer in range 0 to 255. The Upq value between two pixels
%             is now found by a table lookup where the table index is given
%             by 1+abs(Uim(p)-Uim(q)), range 1 to 256.
%             ex: ..., 'Uim',teximage(104), ...
%   'Utab'   -A table with 256 elements, if type is double range is
%             assumed to be from 0 to 1, if type is int16 range is assumed
%             to be from 0 to PU1 (PU1 is defined as 4096 i EMP8_common.h)
%             Two shortcuts are defined:
%             1: (linear)   Utab = [1,x,y],  ex: x=150;y=0.3;
%               -> Utab = [1-((0:x)*((1-y)/x)),y*ones(1,255-x)];
%             2: (cosine)   Utab = [2,x,y],  ex: x=150;y=0.3;
%               -> Utab = [(1+y)/2+((1-y)/2)*cos((0:x)*(pi/x)),y*ones(1,255-x)];              
%             ex: ..., 'Utab',[2,120,0.2], ...
% 
%   'scaleRi' may be given, if not a suitable value is found if needed
%   'dR',  'minsumRi' and 'maxsumRi' may also be given as input
%   'v' or 'verbose' may be 0/false (default), 1/true or 2 (very verbose).
%
%   'lowerSize' or 'numPix1' components smaller than this (in number of pixels)
%       will always be eroded, larger components will only be eroded if the
%       total cost is decreased. Default is 10.
%   'upperSize' or 'numPix2' components larger than this (in number of pixels)
%       will not be eroded. Default is 2500. If this is 0 no components will be
%       eroded, border region may be eroded though.
%   'shapeFactor' or 'sf'  used for priority of components, range 0 to 1
%       default value 0.6 is ok. Small value means than shape count less, 
%       i.e. only number of pixels counts when shapeFactor is 0. 
%   'borderWidth' og 'bw' gives the width of the border region to be eroded. 
%       Default 0 means that no erosion along borders is done.
%       Alternative -1 use bw = [3,5,8,13,21,34] in EMP8_swap
%   'maxLoop' maximum number of loops to do, each loop may do several 
%             alpha-beta-swaps. If no pixel is swapped in a loop, no more
%             loops are done.
%   'alpha'   the alpha value(s) to use in each swap. 
%   'beta'    the beta value(s) to use in each swap. 
%             Both alpha and beta may be a matrices of maxLoop columns. 
%             The non-zeros elements in each column are combined by each other, 
%             and the pairs (where beta > alpha) are processed in each loop. 
%             -1 may be given to select all combinations in random order
%   'onlyNeighbors' true (default) does only process pair where alpha and 
%             beta are neighbors, i.e. have common border.
%   'useRand' used by ICM to randomize pixel order (or process pixels sequentially)
%             default: true.
%
% Example:
% ---------------------
%   tim = 4;  sigma = 5;  
%   s = EMP8_argin(tim, sigma);   
%   s = EMP8_argin(tim, sigma, 'v',2);  
%   s = EMP8_argin(tim, sigma, tim, 'v',2);   % use Lt as well
%   s = EMP8_argin(tim, sigma, tim, 'swUpq',3, 'Uim',tim, 'Utab',[2,120,0.2], 'v',2);
%   % and plot results in a figure
%   cc = conncomp(reshape(s.Li,s.I,s.J));
%   figure(1); clf; 
%   figTitle = ['Labeling of test image ',int2str(tim),' using sigma = ',num2str(sigma)];
%   LIS_plotcol(reshape(s.Li,s.I,s.J), figTitle, s.Ci.errorRate, ...
%               cc.NumObjects, s.Ci.Etot, s.Ci.Edata, s.Ci.E4+s.Ci.E8);

%----------------------------------------------------------------------
% Copyright (c) 2013.  Karl Skretting.  All rights reserved.
% University of Stavanger, Signal Processing Group
% Mail:  karl.skretting@uis.no   Homepage:  http://www.ux.uis.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  10.12.2013  Make function
% Ver. 1.1  14.10.2014  May also call function with s from last time
%                       also remove field 'param'
%----------------------------------------------------------------------

% some identities that may be helpful 
% i = rem(pix-1,I)+1;  j = floor((pix-1)/I)+1;  pix = (j-1)*I+i;

%% Check if input arguments are given
if (nargin < 1)  % just check number of input arguments 
    s = [mfilename,': arguments must be given, see help.'];
    disp(s);
    return
end

%% defaults, initial values
s = struct( 'Ri',int16([]), 'Li',uint8([]), 'Lt',uint8([]), ...
            'Upq', int16([]), 'Uim',uint8([]), 'Utab',int16([]), ...
            'M',-1, 'I',-1, 'J',-1, 'P',-1, ...
            'formatCost',' Labeling %8s --> E = %7.5f + %5.2f*(%7.5f + %7.5f) = %7.5f, err. %6.3f%% ', ...
            'PU1',4096, 'swUpq',1, ...
            'lambda',4, 'bF4',-1, 'bF8',-1, 'bF8percent',70, ...
            'scaleRi',-1, 'dR',-1, 'minsumRi',-1, 'maxsumRi',-1, ...
            'numPix1',10, 'numPix2',2500, 'shapeFactor',0.6, ...
            'borderWidth',0, ...   % -1 ---> [3,5,8,13,21,34], ...
            'alpha',-1, 'beta',-1, ...
            'maxLoop',6, 'onlyNeighbors',true, 'useRand',true, ...
            'verbose',0 );
        
%% get the arguments and options
nofNumArg = 0;
if isnumeric(varargin{nofNumArg+1})  
    nofNumArg = nofNumArg+1;
    s.Ri = varargin{nofNumArg};    % first option was Ri
end
if isnumeric(varargin{nofNumArg+1})  
    nofNumArg = nofNumArg+1;
    s.Li = varargin{nofNumArg};    % next argument was Li
end
if isnumeric(varargin{nofNumArg+1})  
    nofNumArg = nofNumArg+1;
    s.Lt = varargin{nofNumArg};    % and third possible Lt
end
%
nofOptions = nargin;  
optionNumber = nofNumArg+1; % options start at argument nofNumArg+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,'R') || strcmpi(opName,'Ri')
        if isnumeric(opVal); s.Ri = opVal; end;
    end
    if strcmpi(opName,'L') || strcmpi(opName,'Li')
        if isnumeric(opVal); s.Li = opVal; end;
    end
    if strcmpi(opName,'Ltrue') || strcmpi(opName,'Lt')
        if isnumeric(opVal); s.Lt = opVal; end;
    end
    if strcmpi(opName,'M')
        if isnumeric(opVal); s.M = floor(opVal(1)); end;
    end
    if strcmpi(opName,'I')
        if isnumeric(opVal); s.I = floor(opVal(1)); end;
    end
    if strcmpi(opName,'J')
        if isnumeric(opVal); s.J = floor(opVal(1)); end;
    end
    if strcmpi(opName,'P')
        if isnumeric(opVal); s.P = floor(opVal(1)); end;
    end
    if strcmpi(opName,'bF4')
        if isnumeric(opVal); s.bF4 = opVal(1); end;
    end
    if strcmpi(opName,'bF8')
        if isnumeric(opVal); s.bF8 = opVal(1); end;
    end
    if strcmpi(opName,'bF8percent')
        if isnumeric(opVal); s.bF8percent = opVal(1); end;
    end
    if strcmpi(opName,'swUpq')
        if isnumeric(opVal); s.swUpq = floor(opVal(1)); end;
        if ((s.swUpq < 1) || (s.swUpq > 3)); s.swUpq = 1; end;
    end
    if strcmpi(opName,'lambda')
        if isnumeric(opVal); s.lambda = opVal(1); end;
    end
    if strcmpi(opName,'scaleRi')
        if isnumeric(opVal); s.scaleRi = opVal(1); end;
    end
    if strcmpi(opName,'dR')
        if isnumeric(opVal); s.dR = opVal(1); end;
    end
    if strcmpi(opName,'minsumRi')
        if isnumeric(opVal); s.minsumRi = opVal(1); end;
    end
    if strcmpi(opName,'maxsumRi')
        if isnumeric(opVal); s.maxsumRi = opVal(1); end;
    end
    if strcmpi(opName,'Upq')
        if isnumeric(opVal); s.Upq = opVal; end;
    end
    if strcmpi(opName,'Uim')
        if isnumeric(opVal); s.Uim = opVal; end;
    end
    if strcmpi(opName,'Utab')
        if isnumeric(opVal); s.Utab = opVal; end;
    end
    if strcmpi(opName,'lowerSize') || strcmpi(opName,'numPix1')
        if isnumeric(opVal); s.numPix1 = opVal(1); end;
    end
    if strcmpi(opName,'upperSize') || strcmpi(opName,'numPix2')
        if isnumeric(opVal); s.numPix2 = opVal(1); end;
    end
    if strcmpi(opName,'shapeFactor') || strcmpi(opName,'sf')
        if isnumeric(opVal); s.shapeFactor = opVal(1); end;
    end
    if strcmpi(opName,'bw') || strcmpi(opName,'borderWidth')
        if isnumeric(opVal); s.borderWidth = opVal; end;
    end
    if strcmpi(opName,'alpha')
        if isnumeric(opVal); s.alpha = opVal; end;
    end
    if strcmpi(opName,'beta')
        if isnumeric(opVal); s.beta = opVal; end;
    end
    if strcmpi(opName,'maxLoop')
        if isnumeric(opVal); s.maxLoop = opVal(1); end;
    end
    if strcmpi(opName,'onlyNeighbors') || strcmpi(opName,'onlyN')
        if isnumeric(opVal); s.onlyNeighbors = (opVal(1) ~= 0); end;
        if islogical(opVal); s.onlyNeighbors = opVal(1); end;
    end
    if strcmpi(opName,'useRand') 
        if isnumeric(opVal); s.useRand = (opVal(1) ~= 0); end;
        if islogical(opVal); s.useRand = opVal(1); end;
    end
    if strcmpi(opName,'verbose') || strcmpi(opName,'v')
        if (islogical(opVal) && opVal); s.verbose = 1; end;
        if isnumeric(opVal); s.verbose = opVal(1); end;
    end
end

%% check input s.Ri, the data matrix
if (isnumeric(s.Ri) && (numel(s.Ri) == 1))
    s.Ri = LIS_getR(s.Ri);    % actually R (double) is returned
end
if ~( isa(s.Ri,'int16') || isa(s.Ri,'double') )
    disp([mfilename,' first argument R (Ri) is not int16 or double, this may cause problem.']);
end
if isa(s.Ri,'int16') 
    if size(s.Ri,2) > 1
        [s.M,s.I,s.J] = size(s.Ri);    % if not a column vector size should be MxIxJ here!
    end
    s.Ri = s.Ri(:);
else  % process as if s.Ri is double
    [s.I,s.J,s.M] = size(s.Ri);
    if (isnumeric(s.Li) && (numel(s.Li) == 1)) % s.Li is sigma, Ri is R
        s.Li = uint8(reshape(LIS_smooth(s.Ri, s.Li), s.I*s.J, 1));
    end
    if (s.scaleRi < 0)   % if >0 it was given as input parameter
        s.scaleRi = 32000/max(abs(s.Ri(:)));
    end
    s.Ri = int16(s.scaleRi*reshape(permute(s.Ri,[3,1,2]),s.M*s.I*s.J,1));   % and make it int16
end
% s.Ri should now be the (int16) column vector with M*I*J elements

%% check input L, the label matrix
if ~( isa(s.Li,'uint8') || isa(s.Li,'double') )
    disp([mfilename,' second argument L (s.Li) is not uint8 or double, this may cause problem.']);
end
if (numel(s.Li) == 0)  % []    use Lmin
    s.Li = uint8(reshape( LIS_smooth(double(permute(reshape(s.Ri,s.M,s.I,s.J),[2,3,1])), 0), s.I*s.J, 1));
elseif (numel(s.Li) == 1)  % a single number is assumed to be sigma
    s.Li = uint8(reshape( LIS_smooth(double(permute(reshape(s.Ri,s.M,s.I,s.J),[2,3,1])), s.Li), s.I*s.J, 1));
else
    if ((s.M > 0) && (numel(s.Ri) ~= (numel(s.Li)*s.M))) 
        s = [mfilename,': sizes of R (s.Ri) and L (s.Li) do not match, see help.'];
        disp(s);
        return
    end
    if (size(s.Li,2) == 1)   % Li is one dimensional, need I or J
        if ((s.I < 0) && (s.J > 0)); s.I = floor(numel(s.Li)/s.J); end;
        if ((s.J < 0) && (s.I > 0)); s.J = floor(numel(s.Li)/s.I); end;
        if ((s.I < 0) && (s.J < 0)) 
            s = [mfilename,': label image size IxJ is not availble, see help.'];
            disp(s);
            return
        end
        if (numel(s.Li) ~= (s.I*s.J))
            s = [mfilename,': label image has not IxJ elements, see help.'];
            disp(s);
            return
        end
    else  % Li is two-dimensional
        [I,J] = size(s.Li);
        if (s.I < 0); s.I = I; end;
        if (s.J < 0); s.J = J; end;
    end
    if (s.M < 0); s.M = floor(numel(s.Ri)/(s.I*s.J)); end;
    if (numel(s.Ri) ~= (s.I*s.J*s.M))
        s = [mfilename,': Ri has not IxJ*M elements, see help.'];
        disp(s);
        return
    end
    s.Li = uint8(s.Li(:));
end


%% check input Ltrue, the label matrix
if numel(s.Lt > 0)
    if (isnumeric(s.Lt) && (numel(s.Lt) == 1))
        s.Lt = teximage(rem(s.Lt,100)+200);
    end
    s.Lt = uint8(s.Lt(:));
    if (numel(s.Lt) ~= (s.I*s.J))
        t = [mfilename,': Ltrue is not expected size, ignore it.'];
        disp(t);
        s.Lt = uint8([]);
    end
end

%% check mex-file
if ~exist('EMP8_erode_mex','file')
    t = [mfilename,': EMP8_erode_mex is not availble, may be trouble.'];
    disp(t);
end
if ~exist('EMP8_cost_mex','file')
    t = [mfilename,': EMP8_cost_mex is not availble, may be trouble.'];
    disp(t);
end

% and check the borderWidth option
if max(s.borderWidth) > (min(s.I,s.J)/4)    
    s.borderWidth(s.borderWidth > min(s.I,s.J)/4) = floor(min(s.I,s.J)/4);
end

%% check swUpq and related arguments
if (s.swUpq == 1)
    % this is ok
elseif (s.swUpq == 2) 
    % check Upq
    if ((numel(s.Upq) == 0) || ~isnumeric(s.Upq))
        s = [mfilename,': Option Upq was not given as expected, (and swUpq == 2).'];
        disp(s);
        return
    end
    if ~( isa(s.Upq,'int16') || isa(s.Upq,'double') )
        disp([mfilename,' Upq is not int16 or double, this may cause problem.']);
    end
    if isa(s.Upq,'int16')  % expected size is 4IJx1 (or 4xIxJ)
        if size(sUpq,2) > 1
            [m,i,j] = size(s.Upq);   
            if ((m ~= 4) || (i ~= s.I) || (j ~= s.J))
                disp([mfilename,' Upq is not expected size, this may cause problem.']);
            end
        end
        s.Upq = s.Upq(:);     % 4IJx1
    else  % we assume that Upq is double, and size is IxJx4
        [i,j,m] = size(s.Upq);
        if ((m ~= 4) || (i ~= s.I) || (j ~= s.J))
            disp([mfilename,' Upq is not expected size, this may cause problem.']);
        end
        s.Upq = int16(s.PU1*reshape(permute(s.Upq,[3,1,2]),4*s.I*s.J,1));
    end
elseif (s.swUpq == 3) 
    % check Uim and Utab
    if ((numel(s.Uim) == 0) || ~isnumeric(s.Uim))
        s = [mfilename,': Option Uim was not given as expected, (and swUpq == 3).'];
        disp(s);
        return
    end
    if (isnumeric(s.Uim) && (numel(s.Uim) == 1))
        s.Uim = teximage(rem(s.Uim,100)+100);        
    end
    if (numel(s.Uim) ~= (s.I*s.J))
        s = [mfilename,': Option Uim is not same size as image.'];
        disp(s);
        return
    end
    s.Uim = uint8(s.Uim(:));
    if (size(s.Uim,1) ~= (s.I*s.J))
        s.Uim = s.Uim(:);
    end
    %
    if ((numel(s.Utab) == 0) || ~isnumeric(s.Utab))
        s = [mfilename,': Option Utab was not given as expected, (and swUpq == 3).'];
        disp(s);
        return
    end
    if (numel(s.Utab) == 3)
        x = s.Utab(2); y = s.Utab(3);
        if (s.Utab(1) == 1)
            s.Utab = [1-((0:x)*((1-y)/x)),y*ones(1,255-x)];
        else    % if (Utab(1) == 2)
            s.Utab = [(1+y)/2+((1-y)/2)*cos((0:x)*(pi/x)),y*ones(1,255-x)];
        end
    end
    if (numel(s.Utab) >= 256)
        s.Utab = s.Utab((1:256)');
    else
        s.Utab = [s.Utab((1:end)'); s.Utab(end)*ones(256-numel(s.Utab),1)];
    end
    if ~isa(s.Utab,'int16')
        s.Utab = int16(s.PU1*s.Utab(:));
    end
else
    s = [mfilename,': illegal value of swUpq (it is ',num2str(s.swUpq),'), see help.'];
    disp(s);
    return
end

%% Display info for input
if (s.verbose > 1)  % very verbose
    disp(' ');
    disp([mfilename,' started ',datestr(now),'.']); 
    disp([' M=',int2str(s.M),', I=',int2str(s.I),', J=',int2str(s.J),', P=',int2str(s.P)]);
    %   (lambda=',num2str(lambda),', swUpq=',int2str(swUpq),')']);
    disp([' Now Ri (',class(s.Ri),') has ',int2str(size(s.Ri,1)),'=',int2str(numel(s.Ri)),' elements.']);
    disp([' and Li (',class(s.Li),') has ',int2str(size(s.Li,1)),'=',int2str(numel(s.Li)),' elements.']);
    if numel(s.Lt) > 0
        disp([' and Lt (',class(s.Lt),') has ',int2str(size(s.Lt,1)),'=',int2str(numel(s.Lt)),' elements.']);
    end
    if numel(s.Upq) > 0
        disp([' and Upq (',class(s.Upq),') has ',int2str(size(s.Upq,1)),'=',int2str(numel(s.Upq)),' elements.']);
    end
    if numel(s.Uim) > 0
        disp([' and Uim (',class(s.Uim),') has ',int2str(size(s.Uim,1)),'=',int2str(numel(s.Uim)),' elements.']);
    end
    if numel(s.Utab) > 0
        disp([' and Utab (',class(s.Utab),') has ',int2str(size(s.Utab,1)),'=',int2str(numel(s.Utab)),' elements.']);
    end
    disp([' bF4 = ',int2str(s.bF4),', bF8 = ',int2str(s.bF8),', lambda = ',num2str(s.lambda), ...
          ', bF8percent = ',int2str(s.bF8percent)  ]);
    disp([' numPix1 = ',int2str(s.numPix1),', numPix2 = ',int2str(s.numPix2), ...
          ', shapeFactor = ',num2str(s.shapeFactor), ', borderWidth = ',num2str(s.borderWidth)]);
    disp([' alpha = ',int2str(s.alpha),', beta = ',int2str(s.beta), ...
          ', maxLoop = ',num2str(s.maxLoop), ...
          ', onlyNeighbors = ',int2str(s.onlyNeighbors),', useRand = ',num2str(s.useRand)]);
end

%% do some calulations to find more parameters
if exist('EMP8_cost_mex','file') && ((s.bF4 < 0) || (s.dR < 0))
    param = [s.M,s.I,s.J,-1, s.bF4,s.bF8,s.swUpq,0, s.lambda, s.dR, s.minsumRi, s.maxsumRi, ...
             s.numPix1, s.numPix2, s.shapeFactor, -1];   % 16 parameters for mex-file
    if (s.bF4 < 0) 
        param(6) = s.bF8percent; 
    elseif (s.bF8 < 0)  % and (bF4 >= 0)
        param(6) = floor(s.bF4*(s.bF8percent/100)+0.5);
    end
    tic;
    if (s.swUpq == 1)
        Ci = EMP8_cost_mex(s.Ri, s.Li, param);
    elseif (s.swUpq == 2)
        Ci = EMP8_cost_mex(s.Ri, s.Li, s.Upq, param);
    elseif (s.swUpq == 3)
        Ci = EMP8_cost_mex(s.Ri, s.Li, s.Uim, s.Utab, param);
    end
    Ci(18) = toc;
    if (s.bF4 < 0)
        s.bF4 = Ci(5);
        s.bF8 = Ci(6);
    end
    s.dR = Ci(14); 
    s.minsumRi = Ci(15); 
    s.maxsumRi = Ci(16);
    if (numel(s.Lt) == (s.I*s.J))
        er = nnz(s.Li ~= s.Lt)/(s.I*s.J);
        Ci(20) = er;
    else
        er = -0.01;
    end
    s.Ci = costTilStruct(Ci);
    if (numel(s.Lt) == (s.I*s.J))
        if (s.swUpq == 1)
            Ct = EMP8_cost_mex(s.Ri, s.Lt, param);
        elseif (s.swUpq == 2)
            Ct = EMP8_cost_mex(s.Ri, s.Lt, s.Upq, param);
        elseif (s.swUpq == 3)
            Ct = EMP8_cost_mex(s.Ri, s.Lt, s.Uim, s.Utab, param);
        end
        s.Ct = costTilStruct(Ct);
    end
    if (s.verbose > 1)  % very verbose
        disp([mfilename,': Find cost for input labeling (mainly to find bF4  and/or dR)']);
        disp([' bF4=',int2str(s.bF4),', bF8=',int2str(s.bF8),', lambda=',num2str(s.lambda), ...
              ', bF8percent=',int2str(s.bF8percent)  ]);
        disp([' dR=',int2str(s.dR),', minsumRi=',int2str(s.minsumRi),...
              ', maxsumRi=',int2str(s.maxsumRi)  ]);
        s.ti = sprintf(s.formatCost, 'Li', Ci([10,13,11,12,9]), 100*er);
        if (er < 0) % it is nicer to print nothing rather the negative error rate
            s.ti = [s.ti(1:72),'.'];
        end
        disp(s.ti);  
        if exist('Ct','var')
            s.tt = sprintf(s.formatCost, 'Ltrue', Ct([10,13,11,12,9]), 0);
            disp(s.tt);  
        end
    end
end

return


function s = costTilStruct(c)
if (numel(c)<20)
    c = [c(:);-ones(20-numel(c),1)];
end
s = struct('Ctot',c(1),'Cdata',c(2),'C4',c(3),'C8',c(4),...
           'bF4',c(5),'bF8',c(6),'n4',c(7),'n8',c(8),...
           'Etot',c(9),'Edata',c(10),'E4',c(11),'E8',c(12),...
           'lambda',c(13),'dR',c(14),'minsumRi',c(15),'maxsumRi',c(16),...
           'updated',c(17),'time',c(18),'unused',c(19),'errorRate',c(20));
return;
       

% % s.param is not used any more, neither as array (p) or struct (s)
% function s = paramTilStruct(p)
% if (numel(p)<16)
%     p = [p(:);-ones(16-numel(p),1)];
% end
% s = struct('M',p(1),'I',p(2),'J',p(3),'P',p(4),...
%            'bF4',p(5),'bF8',p(6),'swUpq',p(7),'verbose',p(8),...
%            'alpha',p(9),'dR',p(10),'minsumRi',p(11),'maxsumRi',p(12),...
%            'lowerSize',p(13),'upperSize',p(14),'shapeFactor',p(15),'unused',p(16));
% return;
