function R2 = texclassify(testim, Dcell, varargin)
% texclassify      Texture Classification of a test image using dictionaries  
%                  This is the classification part of the FTCM method, see
% Texture Classification (Segmentation) Tools for Matlab.
% http://www.ux.uis.no/~karlsk/tct/index.html
%
% use: R2 = texclassify(testim, Dcell, varargin);
% ---------------------------------------------------------------------
%  R2      : 2-norm of representation error for each test vector, size is IxJxM
%  testim  : the (grayscale) test image of size IxJ, may be given as
%            - an image (scaled to be in the same range as training images)
%            - a file name, which is read using teximage and option imscale
%  Dcell   : A cell array with M elements, each element generates a layer in R2.  
%            Each cell element should contain a dictionary and may be either:
%            - the dictionary as an NxK numeric array 
%            - a dictionary struct Ds (for example as returned by texdictlearn)
%            - the name of a mat file where a dictionary structure is stored,
%              i.e.  Ds = load(Dcell{m});
%  Additional options may be given as pairs: 'option',value, ...
%  A needed option not given will have the value given in the dictionary (Ds)
%   imrange   may be [-1,1] (default), [0,1], [-128,127], [0,255]
%   vecMet    the neighborhood for each pixel as im makevector.m 'a'-'j'
%   samet     the sparse approximation method to use.
%   saopt     the additional options to use in sparseapprox.m as a struct. 
%   gt        ground truth may be given as a file name or an IxJ array
%   resfile   name of mat file where results are stored. If this file exist R2
%             will be read from this file and  not generated again
%   verbose   verbose level given as 0 (only warnings) 1 (normal) 2 (much)
% ---------------------------------------------------------------------

%----------------------------------------------------------------------
% 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  20.02.2013  Made function based on an older variant (tcimage.m)  
%----------------------------------------------------------------------

%% check input
mfile = 'texclassify'; 

if (nargin < 2)
    error([mfile,': wrong number of arguments, see help.']);
end
if ~(iscell(Dcell))
    error([mfile,': second argument Dcell should be a cell array, see help.']);
end

% check dictionaries first, set Ds as first found
M = numel(Dcell);
for m=1:M
    if isstruct(Dcell{m})
        if ~exist('Ds','var'); Ds = Dcell{m}; end;
    end
    if ischar(Dcell{m})
        fn = Dcell{m};   
        % should be a mat-file, check mat extension
        if ((numel(fn) < 5) || ~strcmpi(fn((end-3):end),'.mat'))
            fn = [Dcell{m},'.mat'];
            Dcell{m} = fn;
        end
        if exist(fn, 'file')
            if ~exist('Ds','var'); Ds = load(fn); end;
        else
            disp([mfile,': can not load dictionary ',int2str(m),' from ',fn]);
        end
    end
end
imrange = [-1,1];
verbose = 0;
if exist('Ds','var')
    if isfield(Ds,'letter'); vecMet = Ds.letter; end;   % obsolete name for vecMet
    if isfield(Ds,'vecMet'); vecMet = Ds.vecMet; end;
    if isfield(Ds,'samet'); samet = Ds.samet; end;
    if isfield(Ds,'saopt'); saopt = Ds.saopt; end;
    if isfield(Ds,'sety'); sety = Ds.sety; end;   % some special samet may use this
end    

%% get options
nofOptions = nargin - 2;
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,'imrange') )
        if isnumeric(opVal) && (numel(imrange)==2)
            imrange = opVal;
        else
            error([mfile,': not numeric for option ',opName]);
        end
    end 
    if strcmpi(opName,'vecMet') || strcmpi(opName,'letter')   % 
        if ischar(opVal)
            vecMet = lower(opVal(1));
            if (sum(vecMet == 'abcdefghij') == 0); vecMet = 'd'; end;
        else
            error([mfile,': not char for option ',opName]);
        end
    end 
    %
    if ( strcmpi(opName,'samet') )
        if ischar(opVal)
            samet = opVal;
        else
            error([mfile,': not character value (string) for option ',opName]);
        end
    end 
    if ( strcmpi(opName,'saopt') )
        if isstruct(opVal)
            saopt = opVal;
        else
            error([mfile,': not struct for option ',opName]);
        end
    end 
    if ( strcmpi(opName,'sety') )
        if isnumeric(opVal) 
            sety = opVal;
        else
            error([mfile,': not numeric array in option ',opName]);
        end
    end 
    %
    if strcmpi(opName,'gt') || strcmpi(opName,'groundTruth') || ...
            strcmpi(opName,'gtfile') || strcmpi(opName,'ground_truth')
        if ischar(opVal)
            B = teximage(opVal);  % ground truth image
        elseif isnumeric(opVal)
            B = opVal;  % ground truth image
        else
            error([mfile,': not char or numeric value for option ',opName]);
        end
    end
    if strcmpi(opName,'resfile')
        if ischar(opVal)
            resfile = opVal;
        else
            error([mfile,': not char for option ',opName]);
        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

%% check if resfile exist
if exist('resfile','var')
    if ((numel(resfile) < 5) || ~strcmpi(resfile((end-3):end),'.mat'))
        resfile = [resfile,'.mat'];
    end
    if exist(resfile,'file')
        load(resfile);
        if exist('R2','var')
            disp([mfile,': load ''R2'' from ',resfile,' and return ''R2''.']);
            if exist('tekst','var'); 
                disp(tekst);  %#ok<NODEF>
            end;
            return;
        else
            disp([mfile,': the results stored in ',resfile,' does not contain ''R2''.']);
            disp(['  make ''R2'' again and store it in ',resfile,'.']);
        end
    end
end
R2 = -1;

%% check for needed arguments
if ~exist('vecMet','var')
    disp([mfile,': option vecMet (as used in makevector) is not given.']);
    return;
end
if ~exist('samet','var')
    disp([mfile,': option samet (sparse approximation method) is not given.']);
    return;
end
if ~exist('saopt','var')
    disp([mfile,': option saopt (sparse approximation options) is not given.']);
    return;
end

%% load test image 
if isnumeric(testim)
    A = testim;
    testim = 'supplied image';
elseif ischar(testim)
    A = teximage(testim, imrange);  
end
[I,J] = size(A);
tekst = char([mfile,': started ',datestr(now)], ...
            ['  Size of ',testim,' is ',int2str(I),'x',int2str(J)],...
            ['  and range is [',num2str(min(A(:))),', ',num2str(max(A(:))),']'] );
if (verbose > 0);  disp(tekst); end;
% expand imagefile A by 5 on all sides
A = [ A(5:(-1):1, :); A; A(I:(-1):(I-4), :) ];
A = [ A(:, 5:(-1):1), A, A(:, J:(-1):(J-4)) ];


%% use each dictionary on the test image
if verbose
    disp([mfile,': processing ',testim,' using ']);
end
R2 = zeros(I,J,M);

for m=1:M
    if isstruct(Dcell{m})
        Ds = Dcell{m}; 
        dictfile = '*dictionary given in struct*';
        if isfield(Ds,'dictfile'); dictfile = Ds.dictfile; end;
    elseif ischar(Dcell{m})
        dictfile = Dcell{m};   
        if exist(dictfile, 'file')
            Ds = load(dictfile); 
        else
            disp([mfile,': can not load dictionary ',int2str(m),' from ',dictfile]);
            continue;
        end
    elseif isnumeric(Dcell{m})
        D = Dcell{m};
        dictfile = '*dictionary given in array*';
    else
        disp([mfile,': dictionary ',int2str(m),' is not given in input, see help.']);
        continue;
    end
    
    if ~isnumeric(Dcell{m})  % so Ds was read here
        if isfield(Ds,'D')
            D = Ds.D;
        else
            disp([mfile,': dictionary ',int2str(m),' has no D field.']);
            continue;
        end
        if isfield(Ds,'letter'); vecMet = Ds.letter; end;   % obsolete name for vecMet
        if isfield(Ds,'vecMet'); vecMet = Ds.vecMet; end;
        if isfield(Ds,'samet'); samet = Ds.samet; end;
        if isfield(Ds,'saopt'); saopt = Ds.saopt; end;
        if isfield(Ds,'sety'); sety = Ds.sety; end;   % some special samet may use this
    end
    
    if verbose
        fprintf('  dictionary %2i, %50s row %3i ',m,dictfile,0);
    end
    
    %
    for i=1:I   % process one row
        if verbose && (rem(i,10) == 0)
            fprintf('\b\b\b\b%3i ',i);
        end
        X = makevector(A, vecMet, 5+i, 5+(1:J), 0);
        % first check for some special sparse approximation methods not (yet)
        % implemented in sparseapprox.m
        if strcmpi(samet, 'cafixed')
            W = cafixed(X, D, sety, 0);
        elseif strcmpi(samet, 'convexsets')
            W = convexsets(X, D, sety, 0);
        elseif strcmpi(samet, 'greedy')
            W = greedy(X, D, saopt.tnz, 0);
        else
            W = sparseapprox(X, D, samet, saopt);
        end
        % R1(i,:,m) = sum(abs( X - ts.D*W ));  % 1-norm
        R2(i,:,m) = sqrt(sum( ( X - D*W ).^2 ));  % 2-norm
    end
    if verbose
        fprintf('\b\b\b\b%3i \n',i);
    end
    %
end
    
%% check number of correct classified pixels
if exist('B','var')
    if (size(B,1) ~= I) || (size(B,2) ~= J) || (min(B(:)) < 1) || (max(B(:)) > M)
        disp([mfile,': ground truth image is not as expected, we ignore it.']);
    else  % B is probably ok
        [~, Y] = min(R2,[],3);
        cc = sum(sum( Y == B ));  % number of Correct Classified pixels
        percenterrors = 100*(1-cc/(I*J));
        temp = sprintf('  Direct classification of %s gives %5.2f percent errors.',testim,percenterrors);
        tekst = char(tekst,temp);  %#ok<NASGU>
        if verbose; disp(temp); end;
    end
end
%
if exist('resfile','var') 
    save(resfile,'R2', 'testim', 'percenterrors', 'mfile', 'tekst', 'vecMet', 'samet', 'saopt');
end

return

