function varargout = EMP8_ex5(no, K, vecMet, lambda, swUpq, makeFig)
% EMP8_ex5         Complete segmentation experiments, FTCM + EM, 
%                  Outex Contrib test images (same as Randen, Husøy)
%
% Training images, test images and ground truth image are all stored in
% catalog for Outex Contrib test images, in the case below each of the 12 test
% images is in a subcatalog, '000' to '011', under my_matlab_path('ContribSS');
% for Outex test suite see:  http://www.outex.oulu.fi/index.php?page=segmentation
%
% Here all the different steps of FTCM are done
%  1.  Learn dictionaries, using                                 texdictlearn.m 
%  2.  Find class distance for all pixels in test image, using   texclassify.m 
%  3.  Smooth the resulting distance images to get a labeling, 
%  3a.  Gaussian smoothing, (fig 3)                              LIS_smooth.m
%  3b.  Swap moves (fig 4)                                       EMP8_swap.m 
%  3c.  Erode moves (fig 6), and edge smoothing (fig 5 and 7)    EMP8_erode.m 
% 
% use: EMP8_ex5(no, K, vecMet, lambda, swUpq, makeFig);
%      txt = EMP8_ex5(no, K, vecMet, lambda, swUpq, makeFig);
% -------------------------------------------------------------------
%  no       which test image from the Contrib SS test suite to use; 0 to 11
%  K        number of atoms in dictionaries, typically from 50 to 200
%  vecMet   method in makevector.m,  'a' to 'j'  (best is probably 'g')
%           no, K, or vecmet may be row vectors to process several values
%  lambda   the lambda to use in cost function, typically 2 
%  swUpq    1 standard, or 3 use test image
%           lambda or swUpq may be row vectors to process several values
%  makeFig  true: display some figures (1-6, or 7) as well
%  txt      the resulting text
% -------------------------------------------------------------------
% ex:  txt = EMP8_ex5(2, 160, 'g', 3.0, 3, true);
% ex:  res = EMP8_ex5(0:4, 160, 'g', 2:0.5:5, [1,3], false);

%----------------------------------------------------------------------
% Copyright (c) 2014.  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  14.10.2014  Made this experiment based on ex02.m
% Ver. 1.1  16.10.2014  Make it possible to use vectors for input 
%----------------------------------------------------------------------


%%  Check input arguments
if (nargin < 1); no = 0; end;
if (nargin < 2); K = 75; end;
if (nargin < 3); vecMet = 'g'; end;
if (nargin < 4); lambda = 2; end;
if (nargin < 5); swUpq = 1; end;
if (nargin < 6); makeFig = true; end;

imrange = [-1,1];
Mno = [5,5,5,5,5,16,16,10,10,2,2,2];
sigmaTab = [2,3,5,6,8,12,4];

% make it possible to process several cases
tab_no = no;
tab_K = K;
tab_vecMet = vecMet;
dim1 = numel(tab_no)*numel(tab_K)*numel(tab_vecMet);
tab_lambda = lambda;
tab_swUpq = swUpq;
dim2 = numel(tab_lambda)*numel(tab_swUpq);
if dim1*dim2 > 1
    tab_result = cell(dim1,dim2);
end

dim1_index = 0;
for no = tab_no
    for K = tab_K
        for vecMet = tab_vecMet;
            dim1_index = dim1_index+1;
            dim2_index = 0;

%% 1. Dictionary learning,    dictionary size is NxK
N = size(makevector(ones(19),vecMet,10,10,0),1);
n2k = 2500;  % gives 2000*n2k iterations (number of training vectors to process)

learn_opt = struct('imrange',imrange, 'vecMet',vecMet, 'K',K, ...
     'samet','mexOMP', ...   % alternative 'javaormp' or 'ORMP' (slow)
     'saopt',  struct('tnz',3, 'verbose',0), ...    % tnz usualy 3, but may be 4
     'ptc', {{'mtvp','r2avg','traceA','fbB','betamse','d2min','nofr','rcondA','nofo'}}, ...
     'minibatch', repmat([n2k,1],5,1).*[1,20; 2,40; 3,100; 3,200; 2,500], ... 
     'lam0', 0.996, 'lam1', 1.1, ...
     'checkrate', 5000, ...
     'dlim', 2.5, 'outlim', 3, ...
     'saveDs', 3, ...   % give 'dictfile' to store dictionary in
     'verbose', 1 );
    

if isunix(); cat_sep = '/'; else cat_sep = '\'; end;
ou_cat = my_matlab_path('ContribSS');
gtfile = sprintf('%s%s%03i%sground_truth.ras', ou_cat, cat_sep, no, cat_sep);

% learn a dictionary for each training image (ras file)
M = Mno(no+1);         % number of training images (candidate textures)
Dcell = cell(M,1);
for m = 1:M
    % 
    % 1.) check that the training image exist
    fn = sprintf('%s%s%03i%strain_%02i.ras', ou_cat, cat_sep, no, cat_sep, m);
    if ~exist(fn,'file')
        disp([mfilename,': can not find ',fn]);
        continue;
    end
    % figure(m); clf; imagesc(imread(fn)); colormap(gray); title(['Original image ',int2str(m)]);
    % A = teximage(fn,imrange); 
    % the learning function just need the file name of the training image
    %
    % 2.) learn a dictionary for fn and store it in dictfile
    dictfile = sprintf('%s%s%03i%sD%02i%s%03i.mat', ou_cat, cat_sep, no, cat_sep, m, vecMet, K);
    Dcell{m} = dictfile;
    if exist(dictfile,'file')
        disp(' ');
        disp([mfilename,': Dictionary ',dictfile,' exists. ']);
        disp('  Delete it before you learn it again. Now it is used as it was.');
        Ds = load(dictfile);
    else
        Ds = texdictlearn('imfile',fn, learn_opt, 'dictfile',dictfile);
    end
    if isfield(Ds,'text')
        disp('Field ''text'' in the mat-file is:'); 
        disp(Ds.text);
        disp(' ');
    end
    %  texdictfig(Ds, 1);    % display figure of learning data
end

%% 2. texture classification of each pixel based using the dictionaries
disp(' ');
testim = sprintf('%s%s%03i%sproblem.ras', ou_cat, cat_sep, no, cat_sep);
% next two lines to store R in a quite large file
% resfile = sprintf('%s%s%03i%sR%s%03i.mat', ou_cat, cat_sep, no, cat_sep, vecMet, K);
% classify_opt = struct('resfile',resfile, 'gt',gtfile, 'verbose',1);
% or next line does not store R (it should take less than one minute to generate it)
classify_opt = struct('gt',gtfile, 'verbose',1);
R = texclassify(testim, Dcell, classify_opt);
disp(' ');

if makeFig
    figure(1); clf; 
    imagesc(teximage(testim,imrange)); colormap(gray);
    title(['Test image: ',testim],'interpreter','none');
end

for lambda = tab_lambda
    for swUpq = tab_swUpq
        dim2_index = dim2_index+1;

%% ********* LABELING from R and smoothing ***************
txt = [mfilename,': Outex Contrib SS no = ',sprintf('%03i',no), ...
    ', vecMet = ''',vecMet,''', N = ',int2str(N),', K = ',int2str(K), ...
    ', lambda = ',sprintf('%5.2f',lambda),', swUpq = ',int2str(swUpq)];

% find (and display) cost for true labeling
Ltrue = teximage(gtfile);
Ltrue = uint8(Ltrue(:));
if swUpq == 3
    disp('Use swUpq = 3');
    s = EMP8_argin(R, 0, Ltrue, 'lambda',lambda, 'swUpq',3, 'Uim',teximage(testim), 'Utab',[2,160,0.2], ...
        'numPix1',20, 'numPix2',5000, 'verbose',2);
else
    disp('Use swUpq = 1');
    s = EMP8_argin(R, 0, Ltrue, 'lambda',lambda, 'swUpq',1, 'numPix1',20, 'numPix2',5000, 'verbose',2);
end
if dim1*dim2 > 1
    s.verbose = 0;
end

% s includes cost for Ltrue, s.Ct, and for Lmin (L=0 --> s.Li = Lmin), s.Ci.
t1 = sprintf(s.formatCost, 'Ltrue', s.Ct.Edata, s.Ct.lambda, s.Ct.E4, s.Ct.E8, s.Ct.Etot, 0);
disp(t1);
txt = char(txt, t1);
if makeFig
    cc = conncomp(reshape(s.Li,s.I,s.J));
    figTitle = ['True labeling for Contrib SS test number ',sprintf('%03i',no)];
    figure(2); clf;
    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);
end

% minimum labeling          
t1 = sprintf(s.formatCost, 'Lmin', s.Ci.Edata, s.Ci.lambda, s.Ci.E4, s.Ci.E8, s.Ci.Etot, 100*s.Ci.errorRate);
disp(t1);
txt = char(txt, t1);

% Gaussian smoothing with sigmaTab
for sigma = sigmaTab;
    Lg = LIS_smooth(double(permute(reshape(s.Ri,s.M,s.I,s.J),[2,3,1])), sigma);
    s.Li = uint8(reshape( Lg, s.I*s.J, 1 )); 
    s.Ci = EMP8_cost(s); 
    t1 = sprintf(s.formatCost, sprintf('Lg%3.1f',sigma), ...
        s.Ci.Edata, s.Ci.lambda, s.Ci.E4, s.Ci.E8, s.Ci.Etot, 100*s.Ci.errorRate);
    disp(t1); 
    txt = char(txt,t1);
end
if makeFig  
    cc = conncomp(reshape(s.Li,s.I,s.J));
    figTitle = sprintf('Lg%3.1f  labeling for Contrib SS test number %03i',sigma,no);
    figure(3); clf;
    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);
end

% swap (from s.Li = Lg, with sigma as last entry in sigmaTab)
s.borderWidth = -1;    % use a good increasing sequence
s.maxLoop = 10;
[Ls,cs] = EMP8_swap(s);
t1 = sprintf([s.formatCost,', %7.1f [ms]'], 'Lswap', cs.Edata, cs.lambda, ...
    cs.E4, cs.E8, cs.Etot, 100*cs.errorRate, 1000*cs.time);
disp(t1); 
txt = char(txt,t1);
if makeFig  
    cc = conncomp(reshape(Ls,s.I,s.J));
    figTitle = sprintf('Lswap labeling for Contrib SS test number %03i',no);
    figure(4); clf;
    LIS_plotcol(reshape(Ls,s.I,s.J), figTitle, cs.errorRate, ...
            cc.NumObjects, cs.Etot, cs.Edata, cs.E4+cs.E8);
end

% erode and border smoothing (from s.Li = Lg, with sigma as last entry in sigmaTab)
s.borderWidth = 5;
[s.Li, ce] = EMP8_erode(s);
temp = ce.time;
nofComp = ce.erodeNofComp;
[Le, ce] = EMP8_icm(s);
ce.time = ce.time+temp;
%
t1 = sprintf([s.formatCost,', %7.1f [ms]'], 'Lerode+', ce.Edata, ce.lambda, ...
    ce.E4, ce.E8, ce.Etot, 100*ce.errorRate, 1000*ce.time);
disp(t1); 
txt = char(txt,t1);
if makeFig  
    figTitle = sprintf('Lerode+ labeling for Contrib SS test number %03i',no);
    figure(5); clf;
    LIS_plotcol(reshape(Le,s.I,s.J), figTitle, ce.errorRate, ...
            nofComp, ce.Etot, ce.Edata, ce.E4+ce.E8);
end

% simple erode (from s.Li = Lg, with sigma as last entry in sigmaTab)
s.borderWidth = 0;
[Le, ce] = EMP8_erode(s);
t1 = sprintf([s.formatCost,', %7.1f [ms]'], 'Lerode', ce.Edata, ce.lambda, ...
    ce.E4, ce.E8, ce.Etot, 100*ce.errorRate, 1000*ce.time);
disp(t1); 
txt = char(txt,t1);
if makeFig  
    figTitle = sprintf('Lerode labeling for Contrib SS test number %03i',no);
    figure(6); clf;
    LIS_plotcol(reshape(Le,s.I,s.J), figTitle, ce.errorRate, ...
            ce.erodeNofComp, ce.Etot, ce.Edata, ce.E4+ce.E8);
end

% and border smoothing by swap
s.Li = uint8(reshape( Le, s.I*s.J, 1 ));   % start from erode results!
s.borderWidth = [5,10,20];
s.maxLoop = 3;
[Ls,cs] = EMP8_swap(s);
t1 = sprintf([s.formatCost,', %7.1f [ms]'], 'Le+swap', cs.Edata, cs.lambda, ...
    cs.E4, cs.E8, cs.Etot, 100*cs.errorRate, 1000*(ce.time+cs.time));
disp(t1); 
txt = char(txt,t1);
if makeFig  
    cc = conncomp(reshape(Ls,s.I,s.J));
    figTitle = sprintf('Lerode+swap labeling for Contrib SS test number %03i',no);
    figure(7); clf;
    LIS_plotcol(reshape(Ls,s.I,s.J), figTitle, cs.errorRate, ...
            cc.NumObjects, cs.Etot, cs.Edata, cs.E4+cs.E8);
end

disp(' ');
disp(txt);
disp(' ');

if dim1*dim2 > 1
    tab_result{dim1_index, dim2_index} = txt;
end

    end
end
        end
    end
end

if (dim1*dim2 > 1) && (size(tab_result,1) == 12)
    % display summary for images
    j = 1;  % or whatever in dim2
    res = zeros(12,1);
    pos = 80:85;   % error rate
    for line = 11:14;
        disp(tab_result{1,j}(1,:));  % first line
        for i = 1:12
            disp(tab_result{i,j}(line,:));
            res(i) = eval(tab_result{i,j}(line,pos));
        end
        t1 = blanks(numel(tab_result{i,j}(line,:)));
        t1(pos) = sprintf('%6.3f',mean(res));
        disp(t1);
        disp(' ');
    end
end

if nargout >= 1
    if dim1*dim2 > 1
        varargout{1} = tab_result;
    else
        varargout{1} = txt;
    end
end

return

