%%%
%%%    This file is part of the evaluation framework for the SABS dataset:
%%%    http://www.vis.uni-stuttgart.de/index.php?id=sabs
%%%
%%%    DISCLAIMER OF WARRANTY
%%%    This source code is provided "as is" and without warranties.
%%%    If you experience any difficulties, please drop be an email:
%%%    benjamin.hoeferlin@vis.uni-stuttgart.de
%%%
function EvaluateForegroundMasks()

    %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Users have to define variables in this block
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Number of evaluations to run (e.g., for each parameter tested)
    numberOfEvaluations = 1;
    
    % Choose evaluation scenario
    % 1) Basic
    % 2) Bootstrap(short)  
    % 3) Uninteresting
    % 4) Sudden Scence Change
    % 5) Gradual Scene Change
    % 6) Shadow
    % 7) Camouflage
    % 8) No Camouflage
    % 9) Noise 
    % 10) MPEG4
    % 11) Bootstrap (long)
    loadScenario = 1;
    testString = {'Basic'; 'Bootstrap(kurz)'; 'Uninteresting'; 'Sudden Scene Change'; 'Gradual Scene Change'; 'Shadow'; 'Camouflage'; 'No Camouflage'; 'Noise'; 'Mpeg4'; 'Bootstrap(lang)'};
    disp('====================');
    disp(['Evaluating Scenario: ' testString{loadScenario, 1}]);
    disp('====================');
    
    
    % Path, folder and filename construction
    % A typical foreground mask filename would e.g., look like:
    % ./EvalMasks/Out03/fg_0201.png
    
    fgPath = './SubtractionMasks'; % root path to subfolders containing the forground masks for evaluation
    
    fgFolderPrefix = 'Out';   % prefix of all folders containing the fg masks - if fgFolderHasIndex is set to true, this is the folder name
    fgFolderHasIndex = true;  % switch to false if only a single folder should be tested
    fgFolderStartIndex = 1;
    fgFolderNumLeadingZeros = 0;
    
    fgFilePrefix = 'fg_';      % prefix of all forground mask filenames
    fgFileStartIndex = 1;
    fgFileNumLeadingZeros = 0;
    fgFileImageType = '.png';
     
    
    
    %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Further variables. NO user interaction is required in this section
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % path to ground-truth data
    gtPath = './GT/';
    gtPrefix = 'GT';

    
    
    % set variables to load gt data
    numRows = 600;
    numColumns = 800;
    
    
    % set evaluation parameters
    if( (loadScenario == 1) || (loadScenario == 2) || (loadScenario == 4) || ...
        (loadScenario == 5) || (loadScenario == 6) || (loadScenario == 9) || ...
        (loadScenario==10) || (loadScenario == 11))
        evaluationParameters = [numRows,numColumns,1,numRows,1,numColumns];
    elseif( (loadScenario==3) ) % dynamic background
        evaluationParameters = [numRows,numColumns,200,560,100,380];
    elseif( (loadScenario==7) || (loadScenario==8) ) % camouflage / no camouflage
        evaluationParameters = [numRows,numColumns,150,500,400,600];
    end
    % set frame numbers for evaluation
    if((loadScenario == 1)||(loadScenario == 2)||(loadScenario == 3)|| ...
       (loadScenario == 4)||(loadScenario == 6)||(loadScenario == 9)|| ...
       (loadScenario == 10))
        firstFrameNumber = 801;
        lastFrameNumber = 1400;
        numberOfFrames = 600;
    elseif((loadScenario == 5))
        firstFrameNumber = 1;
        lastFrameNumber = 1400;
        numberOfFrames = 1400;
    elseif((loadScenario == 7)||(loadScenario == 8))
        firstFrameNumber = 1;
        lastFrameNumber = 600;
        numberOfFrames = 600;
    elseif((loadScenario == 11))
        firstFrameNumber = 1;
        lastFrameNumber = 1400;
        numberOfFrames = 1400;
    end
        
    fgFileOffset = fgFileStartIndex - firstFrameNumber;
    
    metrics = zeros(numberOfFrames,8);
    meanMetrics = zeros(numberOfEvaluations,4);
    meanBasicMetrics = zeros(numberOfEvaluations,4);
    rocValues = zeros(numberOfEvaluations,2);
    prValues = zeros(numberOfEvaluations,2);
    
    
    %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % COMPUTATIONS - NO user interaction required beyond this point
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    for currentEvaluation = 1 : numberOfEvaluations
        
        for f = firstFrameNumber : lastFrameNumber
            
            % create index for ground-truth data
            gtIndex = sprintf('%04u', f);
                       
            % create index for forground masks
            if(fgFileNumLeadingZeros == 0)
                fgFileIndex = int2str(f + fgFileOffset);
            else    
                fgFileIndex = sprintf(['%0' int2str(fgFileNumLeadingZeros) 'u'], ...
                    (f + fgFileOffset));
            end
            
            % current evaluation folder index
            evaluationIndex = (fgFolderStartIndex + (currentEvaluation-1));
            
            % create ground truth image name and read ground truth image
            gtName = [gtPath gtPrefix gtIndex '.png'];
            gtImage = imread(gtName);
            
            % create filename of foreground mask and read the image
            if(fgFolderHasIndex)
                if(fgFolderNumLeadingZeros == 0)
                    fgFilename = fullfile(fgPath, [fgFolderPrefix int2str(evaluationIndex)], [fgFilePrefix fgFileIndex fgFileImageType]);
                else
                    fgFolderName = sprintf(['%s%0' int2str(fgFolderNumLeadingZeros) 'u'], ...
                        fgFolderPrefix, evaluationIndex );
                    fgFilename = fullfile(fgPath, fgFolderName, [fgFilePrefix fgFileIndex fgFileImageType]);
                end
            else
                fgFilename = fullfile(fgPath, fgFolderPrefix, [fgFilePrefix fgFileIndex fgFileImageType]);
            end
            fgImage = imread(fgFilename);
            
            % display processing information
            if( mod(f-firstFrameNumber, 100) == 0)
                disp(fgFilename);
            end
            
            [tp,tn,fp,fn] = ComputeBasicMetrics(fgImage,gtImage,evaluationParameters);
            
            metrics(f,1) = tp;
            metrics(f,2) = tn;
            metrics(f,3) = fp;
            metrics(f,4) = fn;
            
        end
                
        plotMetrics(metrics,firstFrameNumber,lastFrameNumber,currentEvaluation);
        
        % compute mean metrics
        [meanTP,meanTN,meanFP,meanFN,meanTpr,meanFpr,meanPrecission,meanFscore] = meanMetricValues(metrics,firstFrameNumber);
        meanBasicMetrics(currentEvaluation,1) = meanTP;
        meanBasicMetrics(currentEvaluation,2) = meanTN;
        meanBasicMetrics(currentEvaluation,3) = meanFP;
        meanBasicMetrics(currentEvaluation,4) = meanFN;
        meanMetrics(currentEvaluation,1) = meanTpr;
        meanMetrics(currentEvaluation,2) = meanFpr;
        meanMetrics(currentEvaluation,3) = meanPrecission;
        meanMetrics(currentEvaluation,4) = meanFscore;
        prValues(currentEvaluation,1) = meanPrecission;
        prValues(currentEvaluation,2) = meanTpr;
        rocValues(currentEvaluation,1) = meanTpr;
        rocValues(currentEvaluation,2) = meanFpr;
        
    end
    
    % plot the residual plots
    disp('    TPR       FPR       Prec      F-Score');
    disp(meanMetrics);
    plotFailureMetrics(meanBasicMetrics);
    plotRoc(rocValues);
    plotPR(prValues);
end




%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% helper functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%
% function [tp,tn,fp,fn] = ComputeBasicMetrics(subImage,gtImage)
% this computes the basic metrics tp,tn,fp,fn by comparing the subtraction
% image of a video frame against the corresponding ground truth data
% input : 
%   subImage             - subtraction image
%   gtImage              - ground truth data
%   evaluationParameters - array containing numRows,numColumns,sr,er,sc,ec
% output :
%   basic metrics tp,tn,fp,fn
function [tp,tn,fp,fn] = ComputeBasicMetrics(subImage,gtImage,evaluationParameters)

    % get evaluationParameters
    sr = evaluationParameters(1,3);
    er = evaluationParameters(1,4);
    sc = evaluationParameters(1,5);
    ec = evaluationParameters(1,6);
    
    % additional parameters
    tp = 0;
    tn = 0;
    fp = 0;
    fn = 0;
    
    % convert subImage to binary image
    subImage = im2bw(subImage);
    
    % define structure element for erode image
    se = [1,1,1;1,1,1;1,1,1];
    
    % compute erode image
    binaryGT = im2bw(gtImage,0.001);
    erodeImage = imerode(binaryGT,se);
    

    % loop through frame and compare with ground truth
    for r = sr : er
        for c = sc : ec            
            
            % increase the basic metric for the current pixel if binaryGT
            % == erodeImage (no contour pixel)
            if(binaryGT(r,c) == erodeImage(r,c))
                if( (subImage(r,c) == 1) && (binaryGT(r,c) ~= 0) )  % ( (gtImage(r,c,1) ~= 0 ) || (gtImage(r,c,2) ~= 0 ) || (gtImage(r,c,3) ~= 0 ) )
                    tp = tp + 1;
                elseif( (subImage(r,c) == 1) && (binaryGT(r,c) == 0) ) %( (gtImage(r,c,1) == 0 ) && (gtImage(r,c,2) == 0 ) && (gtImage(r,c,3) == 0 ) ) 
                    fp = fp + 1;
                elseif( (subImage(r,c) == 0) && (binaryGT(r,c) == 0) ) % ( (gtImage(r,c,1) == 0 ) && (gtImage(r,c,2) == 0 ) && (gtImage(r,c,3) == 0 ) ) 
                    tn = tn + 1;
                elseif( (subImage(r,c) == 0) && (binaryGT(r,c) ~= 0) ) % ( (gtImage(r,c,1) ~= 0 ) || (gtImage(r,c,2) ~= 0 ) || (gtImage(r,c,3) ~= 0 ) ) 
                    fn = fn + 1;
                end

            end
            
        end
    end
    

end

%%%
% this function returns the mean metric values computed for a set of
% subtraction images. the mean of TPR,FPR,Precision,F-Score are
% returned (in this order, mean metrics is a (1x4) matrix)
% input :
%   metric     - matrix containing basic metrics of the frames
%              - startFrame, index of the first frame, number of frames for
%                computing mean values is defined by this value
% output :
%              - (1x4) matrix containing the mean metric values as
%                discussed above
function [meanTP,meanTN,meanFP,meanFN,meanTPR,meanFPR,meanPrecission,meanFScore] = meanMetricValues(metrics, firstFrame)

    lastFrame = size(metrics,1);
    numberOfFrames = lastFrame - firstFrame + 1;
    meanTP = 0; 
    meanTN = 0;
    meanFP = 0;
    meanFN = 0;
    
    for f = firstFrame : lastFrame
       meanTP = meanTP + metrics(f,1);
       meanTN = meanTN + metrics(f,2); 
       meanFP = meanFP + metrics(f,3); 
       meanFN = meanFN + metrics(f,4); 
    end
    
    % compute the mean metric values and store them in meanMetrics
    meanTP = meanTP/numberOfFrames;
    meanTN = meanTN/numberOfFrames;
    meanFP = meanFP/numberOfFrames;
    meanFN = meanFN/numberOfFrames;
    meanTPR = meanTP / (meanTP + meanFN);
    meanFPR = 1 - (meanTN)/(meanTN + meanFP);
    meanPrecission = meanTP / (meanTP + meanFP);
    meanFScore = (2*meanTPR*meanPrecission)/(meanTPR + meanPrecission);

end


%%%
% this function plots the computed metrics in a matlab figure
% input : 
%   metrics -              array containing the basic metrics of the subtraction phase
%   evaluationParameters - contains numRows,numColumns (sr,er,sc,ec)
% output :
%   -
function plotMetrics(metrics,firstFrame,lastFrame,currentEvaluation)
    figure(1);
    % get number of evaluated frames
    numberOfFrames = size(metrics,1);
    
    % create arrays for composed metrics (recall, precision, fpr and
    % f-score)
    recallValues = zeros(1,numberOfFrames);
    precisionValues = zeros(1,numberOfFrames);
    fprValues = zeros(1,numberOfFrames);
    fscoreValues = zeros(1,numberOfFrames);

    % fill these arrays
    for f = 1 : numberOfFrames
        if(f < firstFrame)
            recallValues(1,f) = 0;
            precisionValues(1,f) = 0;
            fprValues(1,f) = 0;
            fscoreValues(1,f) = 0;
        else
            recallValues(1,f) = metrics(f,1)/(metrics(f,1)+metrics(f,4));
            precisionValues(1,f) = metrics(f,1)/(metrics(f,1)+metrics(f,3));
            fprValues(1,f) = 1 - ((metrics(f,2))/(metrics(f,2)+metrics(f,3)));
            fscoreValues(1,f) = (2*recallValues(1,f)*precisionValues(1,f))/(recallValues(1,f)+precisionValues(1,f));
        end
    end
    
    % set color of plot according to currentEvaluation. after 10th plot,
    % all further plots will have the same color. (too many plots will be confusing)
    switch(currentEvaluation)
        case 1
            color = [0,0,0];
        case 2
            color = [1,0,0];
        case 3
            color = [0,1,0];
        case 4
            color = [0,0,1];
        case 5
            color = [0,1,1];
        case 6
            color = [1,0,1];
        case 7
            color = [0.5,0,0];
        case 8
            color = [0,0.5,0];
        case 9
            color = [0,0,0.5];
        case 10
            color = [0.5,0.5,0];
        otherwise
            color = [0.5,0,0.5];
    end

    % plot recall
    hold on
    frameIndexes = 1:1:lastFrame;
    subplot(2,2,1);
    plot(frameIndexes,recallValues(1,:),'Color',color);
    title('Recall - Plot')
    xlabel('Frame')
    ylabel('Recall')
    axis([0,lastFrame,0,1.0])
    hold off
    
    % plot precision
    hold on
    subplot(2,2,2);
    plot(frameIndexes,precisionValues(1,:),'Color',color);
    title('Precision - Plot')
    xlabel('Frame')
    ylabel('Precision')
    axis([0,lastFrame,0,1.0])
    hold off
    
    % plot false positive rate
    hold on
    subplot(2,2,3);
    plot(frameIndexes,fprValues(1,:),'Color',color);
    title('FPR - Plot')
    xlabel('Frame')
    ylabel('FPR')
    axis([0,lastFrame,0,1.0])
    hold off
    
    % plot f-score
    hold on
    subplot(2,2,4);
    plot(frameIndexes,fscoreValues(1,:),'Color',color);
    title('F-Score - Plot')
    xlabel('Frame')
    ylabel('F-Score')
    axis([0,lastFrame,0,1.0])
    hold off
    
end


function plotPR(prValues)

    
    figure(4);

    % get number of evaluations (length of roc values)
    numberOfEvaluations = size(prValues,1);
    
    % init value array
    precissionValues = zeros(1,numberOfEvaluations);
    recallValues = zeros(1,numberOfEvaluations);

    % write rocValues in inverse order -> plot is monotonically increasing
    for e = 1 : numberOfEvaluations
       precissionValues(1,e) = prValues(numberOfEvaluations+1-e,1);
       recallValues(1,e) = prValues(numberOfEvaluations+1-e,2);
    end

    % plot these values. (0,0) and (1,1) are fixed points to get better
    % curve plot curve as a line
    hold on
    plot(recallValues,precissionValues);
    title('PR - Plot')
    xlabel('Recall')
    ylabel('Precission')
    axis([0,1.0,0,1.0])
    hold off
    % emphasize computed values
    hold on
    plot(recallValues,precissionValues,'x');
    axis([0,1.0,0,1.0])
    hold off

end

%%%
% function plotMeanMetrics(metrics)
% this function plots the mean metric values as a bar plot for each
% evaluated folder of subtraction images
function plotFailureMetrics(meanBasicMetrics)
    
    figure(2);
    values = zeros(size(meanBasicMetrics,1),2);
    maxValue = 0;
    for f = 1 : size(meanBasicMetrics,1)
        values(f,1) = meanBasicMetrics(f,3);
        values(f,2) = meanBasicMetrics(f,4);
        maxValue = max(maxValue,(values(f,1)+values(f,2)));
    end
    axis([0  maxValue 0 size(meanBasicMetrics,1)+1]);
    title('False Classification')
    xlabel('Number of false classified pixels')
    ylabel('Evaluation number')
    hold on
    barh(values,'stacked');

end

%%%
% function [tpr,fpr] = computeRocValue(metrics)
% this function coumputes evaluates the result of a bgs method and maps it
% into roc-space. points in this space have the form (tpr,fpr) where tpr
% is the (mean) true positive rate and (mean) true negative rate for an
% evaluated video.
% input : 
%   metrics - array containing the basic metrics for each frame of a bgs
%             methos
% output : 
%   tpr : - mean true positive rate
%   fpr : - mean false positive rate
function [tpr,fpr] = computeRocValue(metrics)

    % compute the mean basic metrics of the sequence needed for tpr and fpr
    meanTp = sum(metrics(:,1));
    meanTn = sum(metrics(:,2));
    meanFp = sum(metrics(:,3));
    meanFn = sum(metrics(:,4));
    
    % compute tpr & fpr
    tpr = meanTp / (meanTp + meanFn);
    fpr = 1 - ((meanTn)/(meanFp+meanTn));
    
end

%%%
% function plotRoc(rocValues)
% this method plots the roc curve with respect to the evaluated tpr and fpr
% metric values for each sequence evaluated
% input:
%   rocValues - a numberOfEvaluationx2 matrix containing (mean) tpr and fpr metric
%               values for each sequence evaluated.
% output :
%   -       
function plotRoc(rocValues)

    figure(3);

    % get number of evaluations (length of roc values)
    numberOfEvaluations = size(rocValues,1);
    
    % init value array
    tprValues = zeros(1,numberOfEvaluations);
    fprValues = zeros(1,numberOfEvaluations);

    % write rocValues in inverse order -> plot is monotonically increasing
    for e = 1 : numberOfEvaluations
       tprValues(1,e) = rocValues(numberOfEvaluations+1-e,1);
       fprValues(1,e) = rocValues(numberOfEvaluations+1-e,2);
    end

    % plot these values. (0,0) and (1,1) are fixed points to get better
    % curve
    % plot curve as a line
    hold on
    plot([0,fprValues,1],[0,tprValues,1]);
    title('ROC - Plot')
    xlabel('False Positive Rate')
    ylabel('True Positive Rate')
    axis([0,1.0,0,1.0])
    hold off
    % emphasize computed values
    hold on
    plot(fprValues,tprValues,'x');
    axis([0,1.0,0,1.0])
    hold off
    
end