function varargout = LIS_expand(Ri, Li, varargin)
% LIS_expand       alpha-expansion (and beta-shrink) on label image using max-flow-algorithm 
%
%  Results are bot returned in scaled form, each term in range 0 to 1, as
%     Etot = Edata + lambda*(E4 + E8); 
%  or as (large) integers for sum in Ri (R scaled and represented as int16)
%     Ctot = Cdata + C4 + C8;
%  where Cdata is sum of the cost in Ri for labeling Li
%  and C4 = n4*bF4 and C8 = n8*bF8
%  and n4 is number of horizontal and vertical edges between labeled pixels
%  and n8 is number of diagonal 'edges' between labeled pixels
%
%   Use of function:
%   ----------------
%   Le = LIS_expand(R, L);       % all options are default
%   [Le, res] = LIS_expand(R, L, 'option',value, ...)
%   [Le, res] = LIS_expand(R, L, Ltrue, 'option',value, ...); % including true labeling
%
%     Le  is resulting labeling after alpha expansion
%     res  is a struct with more results
%     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)
%     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
%     Additional arguments may be given as pairs: 'option',value, ...
%       (or options may be given in one (or more) struct or cell-array)
% 
%       Options used when calculating the cost for a labeling are:
%   '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
% 
%       Options controlling how many times and how alpha-expansion is done:
%   'maxLoop' maximum number of loops to do, each loop may do several 
%             alpha-expansions. If no pixels are swapped in a loop, no more
%             loops are done.
%   'alpha'   the alpha value(s) to use for each alpha-expansion done. It may
%             be a matrix of maxLoop columns. The non-zeros elements in each 
%             column is processed in each loop. 
%             If alpha is not given, (or 'alpha',-1), random values are used.
%   'beta'    0 (default) do not include beta-shrink modification
%             >0 use this particular beta (maxLoop = 1, and numel(alpha) = 1)
%             -1 chose beta at random each time
%             -2 a 'promising' beta is selected each time
%   'brems'   a brake may be included to slow down the 'update speed',
%             default is 0 (no brake), or it may be up to 1 (much brake). 
%             may be given as a row vector, one value for each loop.
%   'bw' or 'borderWidth' restrict alpha extension to be done within this
%             distance from border of alpha region. Default 0 means that 'bw'
%             restriction is not used, all of the image is considered.
%             It may be given as a row vector, one value for each loop.
%
%        Other options
%   'scaleRi' may be given, if not a suitable value is found if needed
%   'FigNo'   number for figure window to display results in, see LIS_plotcol.m
%             default: [0,0,0] (no plot for input or output or expand region)
%   'v' or 'verbose' may be 0/false (default), 1/true or 2 (very verbose).
% 
%   Examples:
%   ---------
%   clear all; tim = 3;  Ltrue = teximage(tim+200);
%   R = LIS_getR(tim);  [I,J,M] = size(R);  
%   sigma = 5; Lg = LIS_smooth(R, sigma);
%   tic; [Le, res] = LIS_expand(R, Lg, Ltrue, 'v',1, 'FigNo',[1,2,0] ); toc;
%   ePar1 = struct('maxLoop',6, 'bw',[2,4,6,10,15,25], 'brems',0, 'alpha',-1);
%   [Le1, res1] = LIS_expand(R, Lg, Ltrue, ePar1, 'v',1, 'FigNo',[1,2,0] ); 
%   ePar2 = struct('maxLoop',1, 'bw',20, 'brems',0, 'alpha',3);
%   [Le2, res2] = LIS_expand(R, Lg, Ltrue, ePar2, 'v',2, 'FigNo',[1,2,3] ); 


%----------------------------------------------------------------------
% Copyright (c) 2012.  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  01.11.2012  Try to make function
% Ver. 1.1  12.11.2012  named it as LIS_expand (in ver. 1.0 it was ex26expand)
% Ver. 1.2  15.11.2012  changed parameters
%----------------------------------------------------------------------

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

mfile = 'LIS_expand';

%% Check if input arguments are given
if (nargin < 2)  % just check number of input arguments 
    t = [mfile,': arguments must be given, see help.'];
    disp(t);
    if nargout >= 1
        varargout{1} = [];
    end
    if nargout >= 2
        varargout{2} = struct('Error',t);
    end
    return
end

%% defaults, initial values
lambda = 2;
bF4 = -1;
bF8 = -1;
bF8percent = 70;
dR = -1;
minsumRi = -1;
maxsumRi = -1;
figno = [0, 0, 0];
maxLoop = 6;
brems = 0;
beta_in = 0;
bw = [3,5,8,13,21,34];
alpha_in = -1;
verbose = 0;
returnA = false;
returnT = false;
    
%% get the options
if (nargin >= 3)
    nofOptions = nargin-2;  % options start at argument 3 
    optionNumber = 1;
    if isnumeric(varargin{optionNumber})  
        Ltrue = varargin{optionNumber};    % first option was Ltrue
        optionNumber = 2;       % next option
    end
end
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,'M')
        if isnumeric(opVal); M = floor(opVal(1)); end;
    end
    if strcmpi(opName,'I')
        if isnumeric(opVal); I = floor(opVal(1)); end;
    end
    if strcmpi(opName,'J')
        if isnumeric(opVal); J = floor(opVal(1)); end;
    end
    if strcmpi(opName,'scaleRi')
        if isnumeric(opVal); scaleRi = opVal(1); end;
    end
    if strcmpi(opName,'lambda')
        if isnumeric(opVal); lambda = opVal(1); end;
    end
    if strcmpi(opName,'bF4')
        if isnumeric(opVal); bF4 = opVal(1); end;
    end
    if strcmpi(opName,'bF8')
        if isnumeric(opVal); bF8 = opVal(1); end;
    end
    if strcmpi(opName,'bF8percent')
        if isnumeric(opVal); bF8percent = opVal(1); end;
    end
    if strcmpi(opName,'dR')
        if isnumeric(opVal); dR = opVal(1); end;
    end
    if strcmpi(opName,'minsumRi')
        if isnumeric(opVal); minsumRi = opVal(1); end;
    end
    if strcmpi(opName,'maxsumRi')
        if isnumeric(opVal); maxsumRi = opVal(1); end;
    end
    if strcmpi(opName,'alpha')
        if isnumeric(opVal); alpha_in = opVal(1); end;
    end
    if strcmpi(opName,'beta')
        if isnumeric(opVal); beta_in = opVal(1); end;
    end
    if strcmpi(opName,'maxLoop')
        if isnumeric(opVal); maxLoop = opVal(1); end;
    end
    if strcmpi(opName,'bw') || strcmpi(opName,'borderWidth')
        if isnumeric(opVal); bw = reshape(opVal(:),1,numel(opVal)); end;
    end
    if strcmpi(opName,'brems') 
        if isnumeric(opVal); brems = reshape(opVal(:),1,numel(opVal)); end;
    end
    if strcmpi(opName,'FigNo')
        if (isnumeric(opVal) && (numel(opVal) >= 3)); figno = opVal(1:3); end;
        if (isnumeric(opVal) && (numel(opVal) == 2)); figno = [opVal(1),opVal(2),0]; end;
        if (isnumeric(opVal) && (numel(opVal) == 1)); figno = [0,opVal(1),0]; end;
    end
    if strcmpi(opName,'verbose') || strcmpi(opName,'v')
        if (islogical(opVal) && opVal); verbose = 1; end;
        if isnumeric(opVal); verbose = opVal(1); end;
    end
    if strcmpi(opName,'returnA') 
        if islogical(opVal); returnA = opVal(1); end;
        if isnumeric(opVal); returnA = (opVal(1) > 0); end;
    end
    if strcmpi(opName,'returnT') 
        if islogical(opVal); returnT = opVal(1); end;
        if isnumeric(opVal); returnT = (opVal(1) > 0); end;
    end
end

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

%% check input L, the label matrix
if (isnumeric(Li) && (numel(Li) == 1))
    if (Li == 0)
        [~,Li] = min(reshape(Ri, M, I*J));   % integers in Ri may give a little
        % bit different labeling than using R (double) direcly from LIS_getR()
        Li = reshape(Li,I,J);
    else
        Li = Li*ones(I,J);
    end
end
if ~( isa(Li,'uint8') || isa(Li,'double') )
    disp([mfile,' second argument L is not uint8 or double, this may cause problem.']);
end
if exist('M','var')
    if ~((numel(Ri)/numel(Li)) == M)
        t = [mfile,': sizes of R and L do not match, see help.'];
        disp(t);
        if nargout >= 1
            varargout{1} = struct('Error',t);
        end
        return
    end
end
if size(Li,2) == 1
    if ~(exist('I','var') && exist('J','var'))
        t = [mfile,': label image size IxJ is not availble, see help.'];
        disp(t);
        if nargout >= 1
            varargout{1} = struct('Error',t);
        end
        return
    end
    Li = reshape(Li, I, J);
else
    [I,J] = size(Li);
end
M = floor(numel(Ri)/(I*J));
Li = uint8(Li(:));

%% check input Ltrue, the label matrix
if exist('Ltrue','var')
    if (isnumeric(Ltrue) && (numel(Ltrue) == 1))
        Ltrue = teximage(Ltrue+200);
    end
    Ltrue = uint8(Ltrue(:));
    if (numel(Ltrue) ~= (I*J))
        t = [mfile,': Ltrue is not expected size, ignore it.'];
        disp(t);
        clear Ltrue
    end
end

%% check some important parameters
if (maxLoop < 1); maxLoop = 1; end;
n = min(maxLoop,numel(bw));
bw = [reshape(bw(1:n),1,n), zeros(1,maxLoop-n)];
n = min(maxLoop,numel(brems));
brems = [reshape(brems(1:n),1,n), zeros(1,maxLoop-n)];
if (alpha_in(1) < 0)
    alpha_in = zeros(M,maxLoop);
    for i=1:maxLoop
        alpha_in(:,i) = reshape(randperm(M),M,1);
    end
end

%% Display info for input
if (verbose > 1)  % very verbose
    disp(' ');
    disp([mfile,' started ',datestr(now),'  (lambda=',num2str(lambda),')']);
    disp([' Use Ri (',class(Ri),', ',int2str(size(Ri,1)),'x',int2str(size(Ri,2)),...
          'x',int2str(size(Ri,3)),') and Li (',class(Li),', ',int2str(size(Li,1)),...
          'x',int2str(size(Li,2)),') M=',int2str(M),', I=',int2str(I),', J=',int2str(J),'.']);
    if exist('Ltrue','var')
        disp([' We also have Ltrue (',class(Ltrue),', ',int2str(size(Ltrue,1)),...
            'x',int2str(size(Ltrue,2)),')']);
    end
    if (maxLoop > 1)
        disp([' Do up to ',int2str(maxLoop),' loops.']);
        if nnz(brems)
            disp([' (brems: ',nnz(brems),' non-zeros from ',num2str(min(nonzeros(brems))), ...
                ' to ',num2str(max(brems)),')']);
        end
        if nnz(bw)
            disp([' (bw: ',nnz(bw),' non-zeros from ',num2str(min(nonzeros(bw))), ...
                ' to ',num2str(max(bw)),')']);
        end
    elseif (numel(alpha_in) == 1)
        disp([' Do one alpha-expansion (alpha=',int2str(alpha_in),...
            ', bw=',num2str(bw(1)),', brems=',num2str(brems(1)),').']);
    end
end

%% do the calulations
paramCost = struct('M',M, 'I',I, 'J',J, 'bF4',bF4, 'bF8',bF8,...
    'dR',dR, 'minsumRi',minsumRi, 'maxsumRi',maxsumRi, 'verbose',0, 'FigNo',0, ...
    'bF8percent',bF8percent, 'lambda',lambda, 'useMex',true);
if exist('Ltrue','var') 
    costLtrue = LIS_cost(Ri, Ltrue, Ltrue, paramCost);      % cost of input labeling
    costLtrue.text(10:17) = '   Ltrue';
    paramCost.bF4 = costLtrue.bF4; 
    paramCost.bF8 = costLtrue.bF8; 
    paramCost.lambda = costLtrue.lambda; 
    paramCost.dR = costLtrue.dR; 
    paramCost.minsumRi = costLtrue.minsumRi; 
    paramCost.maxsumRi = costLtrue.maxsumRi;
    paramCost.FigNo = figno(1);
    costLi = LIS_cost(Ri, Li, Ltrue, paramCost);      % cost of input labeling
else
    paramCost.FigNo = figno(1);
    costLi = LIS_cost(Ri, Li, paramCost);      % cost of input labeling
    paramCost.bF4 = costLi.bF4; 
    paramCost.bF8 = costLi.bF8; 
    paramCost.lambda = costLi.lambda;
    paramCost.dR = costLi.dR; 
    paramCost.minsumRi = costLi.minsumRi; 
    paramCost.maxsumRi = costLi.maxsumRi;
end
bF4 = costLi.bF4;
bF8 = costLi.bF8;
lambda = costLi.lambda;
if (verbose >= 1)
    disp([' The used values for bF4 and bF8 and lambda are:  bF4=',num2str(bF4),...
            ', bF8=',num2str(bF8),' and lambda=',num2str(lambda),'.']);
    if exist('Ltrue','var') 
        disp([' ',costLtrue.text]);    
    end
    disp([' ',costLi.text]);    
end

L = reshape(Li,I,J);   % IxJ matrix of uint8
if (bF8 > 0)
    nhood = strel(ones(3,3));
else
    nhood = strel([0,1,0;1,1,1;0,1,0]);
end

% alpha-expansion for all labels, one by one 
for loop = 1:maxLoop
    totalSwapped = 0;
    if (bw(loop) > 0)       
        nh1 =  strel('disk',bw(loop),0);
    end
    for alpha = reshape(alpha_in(:,loop),1,size(alpha_in,1));
        if ((alpha < 1) || (alpha > M))
            continue; 
        end
        num_alpha = nnz(L == alpha);   
        if (num_alpha == (I*J))    % no expansion is possible
            continue; 
        end
        
        if (numel(beta_in) == 1)
            if (beta_in >= 0)
                beta = beta_in; 
                if (beta > M); beta = M; end;
                if (beta == alpha); beta = ceil((M-1)*rand(1)); end;
                if (beta == alpha); beta = M; end;
            elseif (beta_in == -1)
                beta = ceil((M-1)*rand(1)); 
                if (beta == alpha); beta = M; end;
            else
                beta = 0; % -2 should select a 'promising' beta
            end
        else
            beta = 0;  % not ready for these cases yet
        end
     
        if (bw(loop) <= 0) || (bw(loop) > max(I,J))
            mask = true(I,J);   % use all pixels
        else
            if (beta == 0)
                mask = xor(imdilate(L == alpha, nh1), imerode(L == alpha, nhood));
            else
                mask = xor(imdilate(L == alpha, nh1), imerode(L == alpha, nh1));
            end
        end
        num_alpha = nnz(L(mask) == alpha);
        
        if (verbose >= 2)
            disp(['  Loop ',int2str(loop),': process ',int2str(num_alpha),...
                ' alpha pixels and ',int2str(nnz(mask)-num_alpha),...
                ' not-alpha pixels ',int2str(alpha),' (bw=',...
                num2str(bw(loop)),', brems=',num2str(brems(loop)),...
                ', alpha=',int2str(alpha),', beta=',int2str(beta),')']);
        end
        
        [A, T] = alphaExpand_makeGraph(Ri, L, mask, bF4, bF8, alpha, beta, 1+brems(loop) );
        
        [~, labels] = LIS_maxflow(A, T);
        
        pix = find(mask);
        num_pix = nnz(mask);
        if (beta > 0)
            num_2beta = nnz( (L(pix) == alpha) & (labels(1:num_pix) == 0) );
            L(pix( (L(pix) == alpha) & (labels(1:num_pix) == 0) )) = beta; 
        else
            num_2beta = 0;
        end
        L(pix(labels(1:num_pix) == 1)) = alpha;   
        
        num_2alpha = nnz(labels(1:num_pix));
        nofSwapped = 2*num_2beta + num_2alpha - num_alpha;
        if (verbose >= 2)
            disp(['  number of pixels swapped is ',int2str(nofSwapped), ...
                ' (',int2str(num_2alpha-num_alpha+num_2beta), ...
                ' -> ',int2str(alpha),' and ',int2str(num_2beta),...
                ' -> ',int2str(beta),').']);
        end
        totalSwapped = totalSwapped + nofSwapped;
    end
    if (verbose >= 1)
        disp(['  Total number of pixels swapped in loop ',int2str(loop),...
            '    is    ',int2str(totalSwapped)]);
    end
    if (totalSwapped == 0)
        break;
    end
end

if figno(3) > 0
    figure(figno(3)); spy(mask);
    title([strrep(mfile,'_','-'),': Pixels processed in loop ',int2str(loop),...
        ', \alpha=',int2str(alpha)]);
end

paramCost.FigNo = figno(2);
if exist('Ltrue','var') 
    costLe = LIS_cost(Ri, L(:), Ltrue, paramCost);  
else
    costLe = LIS_cost(Ri, L(:), paramCost);  
end
costLe.text(10:17) = ' Lexpand';
if (verbose >= 1)
    disp([' ',costLe.text]);    
end

%% return results
res = struct('costLi',costLi, 'costLe',costLe,'paramCost',paramCost);

if (exist('costTrue','var')); res.costTrue = costTrue; end;         
if (returnA && exist('A','var')); res.A = A; end;
if (returnT && exist('T','var')); res.T = full(T); res.labels = labels; end;

if nargout >= 1
    varargout{1} = L;   % IxJ
end
if nargout >= 2
    varargout{2} = res;
end
if nargout >= 3
    for i=3:nargout
        varargout{i} = [];
    end
end

return

%%  make the graph
function [A,T] = alphaExpand_makeGraph(Ri, L, mask, bF4, bF8, alpha, beta, T1fac)
% alphaExpand_makeGraph   return the graph (A and T) for the max-flow problem
%  
% Ri        MIJx1 (double/int32/int16)
% L           IxJ (double/uint8, 1..M) 
% mask        IxJ (logical)
% bF4,bF8     1x1 (double/int32/int16)
% alpha       Ix1 (double/uint8, 1..M) 
%
% ex:
% clear all; L = [1,1,1,1; 1,2,2,2; 3,3,2,2; 3,3,3,2]; [I,J] = size(L);
% M = 3; Ri = ceil(50*rand(M*I*J,1)); bF4 = 6; bF8 = 5; alpha = 2;
% mask = true(I,J); mask([1,I],[1,J]) = false(2);
% [A,T] = alphaExpand_setAandT(Ri, L, mask, bF4, bF8, alpha);

[I,J] = size(L);
M = numel(Ri)/(I*J);
largeCapacity = 2*ceil(35000 + 4*(bF4+bF8) + 1);

nnzPix = nnz(mask);         % number of pixels to use, each is a node 
pix = zeros(I,J);
pix(mask) = 1:nnzPix;       % renumber the selected pixels (nodes)

% identify edges in L (label image) where both pixels are in pix
Ea = (L(1:(I-1),:) ~= L(2:I,:));   % | edge above, (I-1)xJ
nEa = ~Ea & mask(1:(I-1),:) & mask(2:I,:);
Ea = Ea & mask(1:(I-1),:) & mask(2:I,:);
nnzEa = nnz(Ea);  nnznEa = nnz(nEa);
El = (L(:,1:(J-1)) ~= L(:,2:J));   % - edge left,  Ix(J-1)
nEl = ~El & mask(:,1:(J-1)) & mask(:,2:J);
El = El & mask(:,1:(J-1)) & mask(:,2:J);
nnzEl = nnz(El); nnznEl = nnz(nEl);
num_nodes = nnzPix + nnzEa + nnzEl;
num_edges = nnznEa + nnznEl + 2*(nnzEa + nnzEl);
if (bF8 > 0)
    Eal = (L(1:(I-1),1:(J-1)) ~= L(2:I,2:J)); % \ edge above left, (I-1)x(J-1)
    nEal = ~Eal & mask(1:(I-1),1:(J-1)) & mask(2:I,2:J);
    Eal = Eal & mask(1:(I-1),1:(J-1)) & mask(2:I,2:J);
    nnzEal = nnz(Eal);  nnznEal = nnz(nEal);
    Ebl = (L(1:(I-1),2:J) ~= L(2:I,1:(J-1))); % / edge below left, (I-1)x(J-1)
    nEbl = ~Ebl & mask(1:(I-1),2:J) & mask(2:I,1:(J-1));
    Ebl = Ebl & mask(1:(I-1),2:J) & mask(2:I,1:(J-1));
    nnzEbl = nnz(Ebl);  nnznEbl = nnz(nEbl);
    num_nodes = num_nodes + nnzEal + nnzEbl;
    num_edges = num_edges + nnznEal + nnznEbl + 2*(nnzEal + nnzEbl);
else
    nnzEal = 0; nnznEal = 0; nnzEbl = 0; nnznEbl = 0;
end
% num_nodes (except s and t), num_edges (except to s and t)
% num_edges = num_edges + I*J + num_nodes;  % add edges to s and t
% num_nodes = num_nodes + 2;                % add s and t
% disp([mfile,': nume_nodes = ',int2str(num_nodes),', num_edges = ',int2str(num_edges)]);
        
% divide the edges into 8 groups here,
n5 = nnzPix + (1:nnzEa)';   % number for new nodes by split of | edges
n6 = nnzPix + nnzEa + (1:nnzEl)'; % and or split of - edges
if (bF8 > 0)
    n7 = nnzPix + nnzEa + nnzEl +(1:nnzEal)'; %  \ edges
    n8 = nnzPix + nnzEa + nnzEl + nnzEal +(1:nnzEbl)'; %  \ edges
    
    i = [ pix([nEa; false(1,J)]); ...   %  | edges
        pix([nEl, false(I,1)]); ...   %  - edges
        pix([nEal,false(I-1,1); false(1,J)]); ...  % \ edges
        pix([false(I-1,1),nEbl; false(1,J)]); ...  % / edges
        pix([Ea; false(1,J)]); n5; ...    % split of | edges
        pix([El, false(I,1)]); n6; ...    % split of -edges
        pix([Eal,false(I-1,1); false(1,J)]); n7; ... % split of \ edges
        pix([false(I-1,1),Ebl; false(1,J)]); n8 ]; % split of / edges
    
    j = [ pix([false(1,J); nEa]); ...   %  | edges
        pix([false(I,1), nEl]); ...   %  - edges
        pix([false(1,J); false(I-1,1),nEal]); ...  % \ edges
        pix([false(1,J); nEbl,false(I-1,1)]); ...  % / edges
        n5; pix([false(1,J); Ea]); ...    % split of | edges
        n6; pix([false(I,1), El]); ...    % split of - edges
        n7; pix([false(1,J); false(I-1,1),Eal]); ... % split of \ edges
        n8; pix([false(1,J); Ebl,false(I-1,1)]) ]; % split of / edges
    
    s = [bF4*ones(nnznEa+nnznEl,1); bF8*ones(nnznEal+nnznEbl,1); ...
        bF4*ones(2*(nnzEa+nnzEl),1); bF8*ones(2*(nnzEal+nnzEbl),1)];
    
else
    i = [ pix([nEa; false(1,J)]); ...   %  | edges
        pix([nEl, false(I,1)]); ...   %  - edges
        pix([Ea; false(1,J)]); n5; ...    % split of | edges
        pix([El, false(I,1)]); n6 ];      % split of -edges
    
    j = [ pix([false(1,J); nEa]); ...   %  | edges
        pix([false(I,1), nEl]); ...   %  - edges
        n5; pix([false(1,J); Ea]); ...    % split of | edges
        n6; pix([false(I,1), El]) ];      % split of - edges
    
    s = bF4*ones(num_edges,1);
end
        
% the matrix A should be symmetric (undirected capacity)
A = sparse([i;j], [j;i], [s;s], num_nodes, num_nodes, 2*num_edges);

% T should contain the additional edges from the nodes to s and t
i = find(mask);
j = (L(i) == alpha);
% T1 are edges towards source (alpha), contains cost for assigning a node/label
% as alpha, 
T1 = [T1fac*double(Ri((i-1)*M+alpha));   % edges for pixels nodes to S (source)
      zeros(nnzEa + nnzEl + nnzEal + nnzEbl,1)]; % , ... % and extra nodes to S
% T2 are edges towards sink (not alpha), cost for keeping the value, or assign
% it to beta (if it was alpha, here indicated by j)
if (beta == 0)
    T2 = [double(Ri(M*(i-1)+double(L(i)))) + ... % edges for pixels nodes to T (sink)
          largeCapacity*double(j); ...  % add for those labeled alpha 
          bF4*ones(nnzEa+nnzEl,1); ...  % extra nodes on | and - edges to T
          bF8*ones(nnzEal+nnzEbl,1)];  % extra nodes on \ and / edges to T
else
    T2 = [double(Ri(M*(i-1)+double(L(i)))).*double(~j) + ... % edges for not alpha pixels nodes to T (sink)
          T1fac*double(Ri(M*(i-1)+beta)).*double(j); ...  % and for those labeled alpha 
          bF4*ones(nnzEa+nnzEl,1); ...  % extra nodes on | and - edges to T
          bF8*ones(nnzEal+nnzEbl,1)];  % extra nodes on \ and / edges to T
end
  
if (nnzPix < (I*J))  % should adjust edge capacities at mask borders
    % disp(['nnzPix = ',int2str(nnzPix)]);
    if (bF8 > 0)
        nhood = strel(ones(3,3));
    else
        nhood = strel([0,1,0;1,1,1;0,1,0]);
    end
    mask = xor(mask, imdilate(mask,nhood));  % the pixels just outside of mask 
    % original mask is in pix, i.e.   orig_mask = (pix > 0);
    Lm = L(mask);
    f = find(pix);
    pe = [zeros(1,J+2); zeros(I,1), pix, zeros(I,1); zeros(1,J+2)]; 
    % bF4 = 2*bF4; bF8 = 2*bF8;   % keep edges more as they were!

    me = [false(I,1),mask,false(I,1); false(2,J+2)];   % above
    pn = pe(me); Ln = zeros(size(Lm)); Ln(pn > 0) = L(f(nonzeros(pn)));
    i = pn(((Ln == Lm) & (Ln ~= alpha)));   % those where neighbor has same label
    if (numel(i) > 0); T1(i) = T1(i) + bF4; end;

    me = [false(1,J+2); mask,false(I,2); false(1,J+2)];   % left
    pn = pe(me); Ln = zeros(size(Lm)); Ln(pn > 0) = L(f(nonzeros(pn)));
    i = pn(((Ln == Lm) & (Ln ~= alpha)));  % those where neighbor has same label
    if (numel(i) > 0); T1(i) = T1(i) + bF4; end;
    
    me = [false(2,J+2); false(I,1),mask,false(I,1)];   % below
    pn = pe(me); Ln = zeros(size(Lm)); Ln(pn > 0) = L(f(nonzeros(pn)));
    i = pn(((Ln == Lm) & (Ln ~= alpha)));  % those where neighbor has same label
    if (numel(i) > 0); T1(i) = T1(i) + bF4; end;
    
    me = [false(1,J+2); false(I,2),mask; false(1,J+2)];   % right
    pn = pe(me); Ln = zeros(size(Lm)); Ln(pn > 0) = L(f(nonzeros(pn)));
    i = pn(((Ln == Lm) & (Ln ~= alpha)));  % those where neighbor has same label
    if (numel(i) > 0); T1(i) = T1(i) + bF4; end;

    if (bF8 > 0)
        me = [mask,false(I,2); false(2,J+2)];   % above and left
        pn = pe(me); Ln = zeros(size(Lm)); Ln(pn > 0) = L(f(nonzeros(pn)));
        i = pn(((Ln == Lm) & (Ln ~= alpha)));   % those where neighbor has same label
        if (numel(i) > 0); T1(i) = T1(i) + bF8; end;
        
        me = [false(2,J+2); mask,false(I,2)];   % below and left
        pn = pe(me); Ln = zeros(size(Lm)); Ln(pn > 0) = L(f(nonzeros(pn)));
        i = pn(((Ln == Lm) & (Ln ~= alpha)));   % those where neighbor has same label
        if (numel(i) > 0); T1(i) = T1(i) + bF8; end;
        
        me = [false(I,2),mask; false(2,J+2)];   % above and right
        pn = pe(me); Ln = zeros(size(Lm)); Ln(pn > 0) = L(f(nonzeros(pn)));
        i = pn(((Ln == Lm) & (Ln ~= alpha)));   % those where neighbor has same label
        if (numel(i) > 0); T1(i) = T1(i) + bF8; end;
        
        me = [false(2,J+2); false(I,2),mask];   % below and right
        pn = pe(me); Ln = zeros(size(Lm)); Ln(pn > 0) = L(f(nonzeros(pn)));
        i = pn(((Ln == Lm) & (Ln ~= alpha)));   % those where neighbor has same label
        if (numel(i) > 0); T1(i) = T1(i) + bF8; end;
        %
    end
        
end

% cancel direct flow from S to pixel/node to T
T = sparse([max(T1-T2,0),max(T2-T1,0)]);

return