function [output] = DemoModels
% Run: DemoModels

%% GOAL: One script to offer the following options...
%        For any image in the original dataset:
%         (1) Display original image
%         (2) Display observer fixations  [new dialog box-- which fixations to show]
%         (3) Show context oracle map
%         (4) Show computational SALIENCY map
%         (5) Show computational CONTEXT map
%         (6) Show computational TARGET FEATURES map
%         (7) Show COMBINED model maps  [new dialog box-- which guidance sources to combine]
%         (8) Display ROC curve (& Behavioral Summary)

%% SET UP:
warning off
close all

%Folders:
HOME = pwd;
init(HOME)

% Dialog Box Input Variables
prompt={'Image (random or enter name- eg. U0888_0000.jpg)'; 'Calculate Human agreement? (y or n)'; 'Fixations to DEFINE region (1:5)';'Fixations to TEST region (1:5)';...
    'Select 2 thresholds to compare (eg. ''5  25'' indicates 5% and 25%)'; ...
    'Show: Saliency Map (y or n)';  'Show: Target features Map (y or n)'; 'Show: Context Map (y or n)'; 'Show: Context Oracle Map (y or n)'; ...
    'Show 2-source model: CONTEXT x SALIENCY';'Show 2-source model: CONTEXT x TARGET';'Show 2-source model: TARGET x SALIENCY';...
    'Show Triple model: '; 'Overlay fixations inside/outide each region?';'Calculate ROC for selected models? (y or n)';};
def={'random';'n'; '1 2 3 4 5'; '1 2 3 4 5'; '5  25'; 'y'; 'y'; 'y'; 'n';'n';'n';'n';'y';'y';'y';};
Boxtitle='Input Variables';
lineNo=1; userinput=inputdlg(prompt,Boxtitle,lineNo,def);        %Input dialog box
imgName = (userinput{1});                           %Subject initials
flag = (char(userinput(2)));  [doUs.DispObsFixations] = convYN(flag);  %Display window with observer's fixtions
LookWhichFixs = str2num(userinput{3});
TestWithFixs = str2num(userinput{4});
thresh = str2num(userinput{5}); thresh = thresh/100;                        %Region size = Fraction of image area
flag = (char(userinput(6)));  [doUs.ShowSaliencyMap] = convYN(flag);        %Display One-source map: Saliency model
flag = (char(userinput(7)));  [doUs.ShowTargetFeaturesMap] = convYN(flag);  %Display One-source map: Target features model
flag = (char(userinput(8)));  [doUs.ShowContextMap] = convYN(flag);         %Display One-source map: Context model
flag = (char(userinput(9)));  [doUs.ShowContextOracleMap] = convYN(flag);   %Display One-source map: Context oracle
flag = (char(userinput(10))); [doUs.dbl_Context_Sal] = convYN(flag);        %Display Two-source map: CONTEXT x SALIENCY
flag = (char(userinput(11))); [doUs.dbl_Context_Tgt] = convYN(flag);        %Display Two-source map: CONTEXT x TARGET
flag = (char(userinput(12))); [doUs.dbl_Tgt_Sal] = convYN(flag);            %Display Two-source map: TARGET x SALIENCY
flag = (char(userinput(13))); [doUs.triple] = convYN(flag);                 %Display Two-source map: TARGET x SALIENCY
flag = (char(userinput(14))); [doUs.overlayFixations] = convYN(flag);       %Overlay fixations inside/outside each region?
flag = (char(userinput(15))); [doUs.calcROC] = convYN(flag); 

%Load Image:
if strmatch(imgName,'random');
    stimList_ALL = dir(fullfile(pwd,'Dataset_STIMULI','*.jpg'));
    stimList_ALL = {stimList_ALL.name}';
    namesIdx = stimList_ALL(randperm(length(stimList_ALL)));
    imgName = namesIdx{1};
end
imgPath = which(imgName);
if isempty(imgPath);
    disp(['ERROR:  ' imgName '  was not found ']); return;
else
    img = imread([imgName]);
    load(fullfile('Dataset_EyeData\eyeData',imgName(1:end-4)))
    eval(['stimFile = ' imgName(1:end-4) ';'])
end;
x=strfind(imgName,'_'); imgName2 = [imgName(1:x) '_' imgName(x+1:end)];
imgName

%Initialize stuff:
stimSize = size(img);  stimSize = stimSize(1:2);  %stimulus dimensions
Tsal_px = thresh * (stimSize(1)*stimSize(2));
whiteMat = zeros(stimSize);                   %blank matrix
LookWhichFixs = LookWhichFixs + 1; 

% GAMMA (for multisource models)
exponentSAL = 0.025;                % modulates: SALIENCY MAP
exponentINRIA = 0.025;              % modulates: OBJECT MAP (INRIA)
gamEx.ContextXSaliency = 0.675;
gamEx.ContextXObject = 0.995;       % Better performance for values closer to 1
gamEx.SaliencyXObject = 0.005;      % Better performance for values closer to 0
gamEx.Triple1 = 0.1;
gamEx.Triple2 = 0.85;

%Figure stuff:
scnsize = get(0,'ScreenSize');
fullscn=[5 40  scnsize(3)-10 scnsize(4)-120];

%% ==================================== %%
%% (2) OBSERVER FIXATIONS
%Eliminate trials with lost eyetracking (or wrong response):
criteria = {'badNaN';'badRT';'accuracy';};
setBadRT = 5000; %cut for badRT trials (ms)
badCriteria = {criteria{1} criteria{3}};
doUs.inclTfixs = 1;

% -A- CHECK TRIAL VALIDITY (no lost eyetracking, Correct Response, Etc)...
Nsubs = size(stimFile.subdata,2); goodTrials = [];
for s = 1:Nsubs
    badNaN = double(stimFile.subdata(s).badNaN);
    badRT = double(stimFile.subdata(s).RT > setBadRT);
    accuracy = isempty(strmatch(stimFile.subdata(s).accuracy,'correct'));
    trial_Is_Ok = 1;
    for c=1:length(badCriteria); if eval(badCriteria{c}); trial_Is_Ok = 0; end;  end
    if doUs.inclTfixs
        hasFixs = find(stimFile.subdata(s).fixType~=0);  %(anything beside INITIAL fixation)
    else hasFixs = find(stimFile.subdata(s).fixType==1);  %(only SEARCH Fixs)
    end
    durCutoff = find(stimFile.subdata(s).fixDur >= 120);
    hasFixs = intersect(hasFixs,durCutoff);
    if isempty(hasFixs); trial_Is_Ok = 0; end
    if trial_Is_Ok; goodTrials = cat(2,goodTrials,s); end
end
Nsubs = length(goodTrials);

% -B- TALLY A LIST OF OBSERVER'S FIXATIONS
recordOfFixs = [];
for s=goodTrials
    if doUs.inclTfixs
        hasFixs = find(stimFile.subdata(s).fixType~=0);  %(anything beside INITIAL fixation)
    else hasFixs = find(stimFile.subdata(s).fixType==1);  %(only SEARCH Fixs)
    end
    whichFixs = LookWhichFixs(find(ismember(LookWhichFixs,hasFixs)));
    nfixs = length(whichFixs);
    fixXY = cat(2,cat(2,repmat(s,nfixs,1),whichFixs'),round(stimFile.subdata(s).fixXY(whichFixs,:))); %[s Nfix X Y]
    if nfixs>0
        badX = find(fixXY(:,3)>800); fixXY(badX,3) = 800;
        badY= find(fixXY(:,4)>600);   fixXY(badY,4) = 600;
        recordOfFixs = cat(1,recordOfFixs,fixXY);
    end
    clear fixXY
end
%Correct for initial fixation index (first SEARCH fixation should be labeled "1")
recordOfFixs(:,2) = recordOfFixs(:,2) - 1;
Nfixs = size(recordOfFixs,1)

%Make a matrix (same size as stimulus) with fixation locations indicated:
fixationMat = zeros(stimSize);
Nfixs = size(recordOfFixs,1);
for f=1:Nfixs
    i = round(recordOfFixs(f,4)); j = round(recordOfFixs(f,3));
    fixationMat(i,j) = fixationMat(i,j)+1;
end

output.imgName = stimFile.imgName;
output.NcontribSs = Nsubs;
output.DefineWithFixs = LookWhichFixs-1;
output.TestWithFixs = TestWithFixs;
output.recordOfFixs = [{'OsNo' 'FixNo' 'X' 'Y'}; num2cell(recordOfFixs)];

% -C- (Optional) DISPLAY IMAGE & FIXATIONS
if doUs.DispObsFixations
    if Nfixs ~=0
        figure('Position',fullscn);
        subplot(1,2,1); imagesc(img); axis off equal; hold on; title([ 'ORIGINAL STIMULUS: ' imgName2]);
        subplot(1,2,2); imagesc(img); axis off equal; hold on; title(['Fixations used to define region']);
        subplot(1,2,2); plot(recordOfFixs(:,3),recordOfFixs(:,4),'o','MarkerSize',12,'LineWidth',2,'MarkerEdgeColor','y','MarkerFaceColor','r');   %REGION
        drawnow
    end

    % -D- (Optional) CALCULATE INTER-OBSERVER AGREEMENT
    disp('Calculating human agreement')
    if Nfixs ~=0
        %% Determine agreement among subsets of observers (using all-except-one method described in paper) ==== %%xs
        if Nsubs > 1; % Image must have more than valid observer (otherwise NULL values):
            N = 0;
            for NthSub = goodTrials; %Exclude one subject at a time (subject "N")
                N = N+1;

                %%%% Split up "recordOfFixs" into 'everyone except Nth' & 'Nth'
                includeObs = setdiff(goodTrials,NthSub);
                DefineRegionFixations = recordOfFixs(ismember(recordOfFixs(:,1),includeObs),:);
                TestRegionFixations = ...
                    recordOfFixs(ismember(recordOfFixs(:,1),NthSub) & ismember(recordOfFixs(:,2),TestWithFixs),:);
                Nfixs_define = size(DefineRegionFixations,1);
                Nfixs_test = size(TestRegionFixations,1);

                %%%% Create imgFixRegions matrix for MAKE fixation predictions (sum ea. fixation, blur)
                imgFixRegions = whiteMat;
                for f=1:Nfixs_define
                    i = round(DefineRegionFixations(f,4)); j = round(DefineRegionFixations(f,3));
                    imgFixRegions(i,j) = imgFixRegions(i,j)+1;
                end
                imgFixRegions = gaussian(imgFixRegions, 8);
                %%%% Create agreeTextFix matrix for TESTING fixation predictions
                agreeTestFix = whiteMat;
                if Nfixs_test > 0
                    for f=1:Nfixs_test
                        i = round(TestRegionFixations(f,4)); j = round(TestRegionFixations(f,3));
                        agreeTestFix(i,j) = 1;
                    end
                end

                %% Get ROC curves for within-image human agreement
                [aROC, falseAlarmRate, detectionRate] = areaROC(imgFixRegions, agreeTestFix);
                output.HOagree(N).areaROC = aROC;
                output.HOagree(N).faRate = falseAlarmRate;
                output.HOagree(N).detRate = detectionRate;

                clear *Matrix *RegionFixations
            end % of goodTrials
            output.meanHOagree = [mean([output.HOagree.areaROC])] %'Mean'
        else
            output.HOagree = [];
            output.meanHOagree = [NaN NaN Nsubs] %'Mean' 'Stdev' 'Nsubs'
        end
    else
        output.HOagree = [];
        output.meanHOagree = [NaN NaN Nsubs] %'Mean' 'Stdev' 'Nsubs'
    end
end

%% ==================================== %%
%% (3) CONTEXT ORACLE MAP
if doUs.ShowContextOracleMap
    %% Load data (mean context oracle map)
    imgPath = which([imgName(1:end-4) '_contextOracleMap.mat']);
    if isempty(imgPath);
        disp(['ERROR:  ' imgName '  was not found in the folder /Dataset_ContextOracle']); return;
    else
        load(imgPath);
    end;

    %% DISPLAY:
    map = contextOracleMap;
    figure('Position',fullscn,'Name',['Context Oracle Map for: ' imgName]);
    subplot(2,2,1:2);
    imagesc(img); axis off equal; hold on; title([ 'Original Stimulus: ' imgName2]);
    for t = 1:2
        [imgFixRegions_thresh{t}] = findPredRegion(map,Tsal_px(t),0);
        subplot(2,2,t+2)
        imagesc(imgFixRegions_thresh{t}.*double(img(:,:,1)));
        colormap(gray); axis off equal; hold on; title(['Empirically Defined Context Region: Thresh ' num2str(thresh(t)*100) '%']);
    end
    output.contextOracleMap = contextOracleMap;
end

%% ==================================== %%
%% (4) SALIENCY MAP
if doUs.ShowSaliencyMap
    %% Load data (mean context oracle map)
    imgPath = which([imgName(1:end-4) '_SaliencyMap.mat']);
    if isempty(imgPath); %Calculate saliency map for this image:
        %% Resizing the image speeds things up:
        if stimSize(1) >= 600; salMapDownsize = 2; else;  salMapDownsize = 1; end
        img_resized = imresize(img,1/salMapDownsize);
        [mIMA,nIMA,p]=size(img_resized);
        disp(['Image: ' imgName ]); disp(['Running Saliency Code on image: ' imgName ' ...'])

        %% Set parameters for salient features:
        Nsc = 4; % # of scales
        Nor = 4; % # of orientations
        Rinhibition = 40; % Radium of inhibition (used when identifying top N most salient points)

        [vL, vLind]=features(img, Nsc, Nor);
        l=0;
        X = zeros(mIMA*nIMA,p*Nsc*Nor);
        for n=1:size(img,3)
            for b=2:size(vLind,1)-1
                x =  (pyrBand(vL(:,n),vLind(:,:,n),b));
                x = imresize(x, [mIMA,nIMA], 'nearest');
                x = (x - mean(x(:)))/std(x(:));
                l=l+1;
                X(:,l)=x(:);
            end
        end

        [mu, cY, k] = fitPowerExponential (X,2);
        saliencyMap = powerExponential(X, mu, cY, k);
        saliencyMap = 1./saliencyMap;
        saliencyMap = reshape(saliencyMap, [mIMA nIMA]);

        salMap_Blur = real(gaussian(saliencyMap.^exponentSAL, 8));
        salMap_Blur = imresize(salMap_Blur,salMapDownsize); salMap = salMap_Blur;

        output.salMap = salMap_Blur;
        save(fullfile(HOME,'saliencyMaps',[imgName(1:end-4) '_SaliencyMap.mat']),'salMap') %Save saliency map in /saliencyMaps
    else
        x = load(imgPath);
        output.salMap = x.salMap;
    end;

    %% DISPLAY:
    map = output.salMap;
    figure('Position',fullscn,'Name',['SALIENCY MAP for: ' imgName]);
    subplot(2,2,1:2)
    imagesc(img); axis off equal; hold on; title([ 'Original Stimulus: ' imgName2]);
    for t = 1:2
        [imgFixRegions_thresh{t}] = findPredRegion(map,Tsal_px(t),0);
        subplot(2,2,t+2)
        imagesc(imgFixRegions_thresh{t}.*double(img(:,:,1)));
        colormap(gray); axis off equal; hold on; title(['SALIENCY Defined Region: Thresh ' num2str(thresh(t)*100) '%']);
    end
    if doUs.overlayFixations
        for t = 1:2
            [fixIn fixOut] = checkFixRegions(output,map,Tsal_px(t));
            subplot(2,2,t+2)
            plot(fixIn(:,1),fixIn(:,2),'wo','MarkerSize',10,'LineWidth',2,'MarkerFaceColor','g');
            plot(fixOut(:,1),fixOut(:,2),'wo','MarkerSize',10,'LineWidth',2,'MarkerFaceColor','r');
            text(0,size(img,1)+50,[num2str(size(fixIn,1)) ' out of ' num2str(size(fixIn,1)+size(fixOut,1)) ' fixations IN region'])
        end
    end
end

%% ==================================== %%
%% (5) Computational CONTEXT MAP
if doUs.ShowContextMap
    paramFile =  'personwalking_Nb8_NgF4';
    ContextMat = foa(uint8(img),fullfile('context','contextModels',paramFile));
    output.contextMap = ContextMat;

    %% DISPLAY:
    map = output.contextMap;
    figure('Position',fullscn,'Name',['Computational CONTEXT MAP for: ' imgName]);
    subplot(2,2,1:2)
    imagesc(img); axis off equal; hold on; title([ 'Original Stimulus: ' imgName2]);
    for t = 1:2
        [imgFixRegions_thresh{t}] = findPredRegion(map,Tsal_px(t),0);
        subplot(2,2,t+2)
        imagesc(imgFixRegions_thresh{t}.*double(img(:,:,1)));
        colormap(gray); axis off equal; hold on; title(['Context Defined Region: Thresh ' num2str(thresh(t)*100) '%']);
    end
    if doUs.overlayFixations
        for t = 1:2
            [fixIn fixOut] = checkFixRegions(output,map,Tsal_px(t));
            subplot(2,2,t+2)
            plot(fixIn(:,1),fixIn(:,2),'wo','MarkerSize',10,'LineWidth',2,'MarkerFaceColor','g');
            plot(fixOut(:,1),fixOut(:,2),'wo','MarkerSize',10,'LineWidth',2,'MarkerFaceColor','r');
            text(0,size(img,1)+50,[num2str(size(fixIn,1)) ' out of ' num2str(size(fixIn,1)+size(fixOut,1)) ' fixations IN region'])
        end
    end
end

%% ==================================== %%
%% (6) TARGET FEATURES MAP
if doUs.ShowTargetFeaturesMap
    %% Load data (mean context oracle map)
    imgPath = which([imgName(1:end-4) '_INRIAobjdet.mat']);
    if isempty(imgPath);
        disp(['ERROR:  ' imgName '  was not found in the folder /targetfeatureMaps']); return;
    else
        x = load(imgPath);
    end;
    % The +10 is an arbitrary correction for the INRIA maps, which
    % tend to have negative numbers (but don't go below -10)
    ObjMat_gaus = real(gaussian((x.output+10).^exponentINRIA, 8));
    output.objectMap = ObjMat_gaus;

    %% DISPLAY:
    map = output.objectMap;
    figure('Position',fullscn,'Name',['TARGET FEATURES MAP for: ' imgName]);
    subplot(2,2,1:2)
    imagesc(img); axis off equal; hold on; title([ 'Original Stimulus: ' imgName2]);
    for t = 1:2
        [imgFixRegions_thresh{t}] = findPredRegion(map,Tsal_px(t),0);
        subplot(2,2,t+2)
        imagesc(imgFixRegions_thresh{t}.*double(img(:,:,1)));
        colormap(gray); axis off equal; hold on; title(['Target Features Region: Thresh ' num2str(thresh(t)*100) '%']);
    end
    if doUs.overlayFixations
        for t = 1:2
            [fixIn fixOut] = checkFixRegions(output,map,Tsal_px(t));
            subplot(2,2,t+2)
            plot(fixIn(:,1),fixIn(:,2),'wo','MarkerSize',10,'LineWidth',2,'MarkerFaceColor','g');
            plot(fixOut(:,1),fixOut(:,2),'wo','MarkerSize',10,'LineWidth',2,'MarkerFaceColor','r');
            text(0,size(img,1)+50,[num2str(size(fixIn,1)) ' out of ' num2str(size(fixIn,1)+size(fixOut,1)) ' fixations IN region'])
        end
    end
end

%% ==================================== %%
%% (7) COMBINED SOURCES OF GUIDANCE MAPS

% (i) 2-source model: Context x Saliency
if doUs.dbl_Context_Sal
    %Check for each component in 2-source model:
    outputnames = fieldnames(output);
    if isempty(strmatch('salMap',outputnames)); blah = showSaliencyMap(imgName,[],0); output.salMap = blah.salMap; end
    if isempty(strmatch('contextMap',outputnames)); blah = showContextMap(imgName,[],0); output.contextMap = blah.contextMap;  end

    output.dbl_Context_Saliency_Map = output.contextMap.^(1-gamEx.ContextXSaliency).*output.salMap.^gamEx.ContextXSaliency;

    map = output.dbl_Context_Saliency_Map; MapName = 'TWO-Sources of Guidance: CONTEXT x SALIENCY Models';
    figure('Position',fullscn,'Name',MapName);
    subplot(2,2,1:2)
    imagesc(img); axis off equal; hold on; title([ 'Original Stimulus: ' imgName2]);
    for t = 1:2
        [imgFixRegions_thresh{t}] = findPredRegion(map,Tsal_px(t),0);
        subplot(2,2,t+2)
        imagesc(imgFixRegions_thresh{t}.*double(img(:,:,1)));
        colormap(gray); axis off equal; hold on; title([MapName ': Thresh ' num2str(thresh(t)*100) '%']);
    end
    if doUs.overlayFixations
        for t = 1:2
            [fixIn fixOut] = checkFixRegions(output,map,Tsal_px(t));
            subplot(2,2,t+2)
            plot(fixIn(:,1),fixIn(:,2),'wo','MarkerSize',10,'LineWidth',2,'MarkerFaceColor','g');
            plot(fixOut(:,1),fixOut(:,2),'wo','MarkerSize',10,'LineWidth',2,'MarkerFaceColor','r');
            text(0,size(img,1)+50,[num2str(size(fixIn,1)) ' out of ' num2str(size(fixIn,1)+size(fixOut,1)) ' fixations IN region'])
        end
    end
end

% (ii) 2-source model: Context x Target
if doUs.dbl_Context_Tgt
    %Check for each component in 2-source model:
    outputnames = fieldnames(output);
    if isempty(strmatch('contextMap',outputnames)); blah = showContextMap(imgName,[],0); output.contextMap = blah.contextMap;  end
    if isempty(strmatch('objectMap',outputnames)); blah = showTargetFeaturesMap(imgName,[],0); output.objectMap = blah.objectMap;  end

    output.dbl_Context_Target_Map = output.contextMap.^(1-gamEx.ContextXObject).*output.objectMap.^gamEx.ContextXObject;
    
    map = output.dbl_Context_Target_Map; MapName = 'TWO-Sources of Guidance: CONTEXT x TARGET Models';
    figure('Position',fullscn,'Name',MapName);
    subplot(2,2,1:2)
    imagesc(img); axis off equal; hold on; title([ 'Original Stimulus: ' imgName2]);
    for t = 1:2
        [imgFixRegions_thresh{t}] = findPredRegion(map,Tsal_px(t),0);
        subplot(2,2,t+2)
        imagesc(imgFixRegions_thresh{t}.*double(img(:,:,1)));
        colormap(gray); axis off equal; hold on; title([MapName ': Thresh ' num2str(thresh(t)*100) '%']);
    end
    if doUs.overlayFixations
        for t = 1:2
            [fixIn fixOut] = checkFixRegions(output,map,Tsal_px(t));
            subplot(2,2,t+2)
            plot(fixIn(:,1),fixIn(:,2),'wo','MarkerSize',10,'LineWidth',2,'MarkerFaceColor','g');
            plot(fixOut(:,1),fixOut(:,2),'wo','MarkerSize',10,'LineWidth',2,'MarkerFaceColor','r');
            text(0,size(img,1)+50,[num2str(size(fixIn,1)) ' out of ' num2str(size(fixIn,1)+size(fixOut,1)) ' fixations IN region'])
        end
    end
end

% (iii) 2-source model: Target x Saliency
if doUs.dbl_Tgt_Sal
    %Check for each component in 2-source model:
    outputnames = fieldnames(output);
    if isempty(strmatch('objectMap',outputnames)); blah = showTargetFeaturesMap(imgName,[],0); output.objectMap = blah.objectMap;  end
    if isempty(strmatch('salMap',outputnames)); blah = showSaliencyMap(imgName,[],0); output.salMap = blah.salMap;  end

    output.dbl_Target_Saliency_Map = output.objectMap.^(1-gamEx.SaliencyXObject).*output.salMap.^gamEx.SaliencyXObject;
    
    map = output.dbl_Target_Saliency_Map; MapName = 'TWO-Sources of Guidance: TARGET x SALIENCY Models';
    figure('Position',fullscn,'Name',MapName);
    subplot(2,2,1:2)
    imagesc(img); axis off equal; hold on; title([ 'Original Stimulus: ' imgName2]);
    for t = 1:2
        [imgFixRegions_thresh{t}] = findPredRegion(map,Tsal_px(t),0);
        subplot(2,2,t+2)
        imagesc(imgFixRegions_thresh{t}.*double(img(:,:,1)));
        colormap(gray); axis off equal; hold on; title([MapName ': Thresh ' num2str(thresh(t)*100) '%']);
    end
    if doUs.overlayFixations
        for t = 1:2
            [fixIn fixOut] = checkFixRegions(output,map,Tsal_px(t));
            subplot(2,2,t+2)
            plot(fixIn(:,1),fixIn(:,2),'wo','MarkerSize',10,'LineWidth',2,'MarkerFaceColor','g');
            plot(fixOut(:,1),fixOut(:,2),'wo','MarkerSize',10,'LineWidth',2,'MarkerFaceColor','r');
            text(0,size(img,1)+50,[num2str(size(fixIn,1)) ' out of ' num2str(size(fixIn,1)+size(fixOut,1)) ' fixations IN region'])
        end
    end
end

% (vi) 3-source model: Context x Saliency x Target
if doUs.triple
    %Check for each component in 2-source model:
    outputnames = fieldnames(output);
    if isempty(strmatch('contextMap',outputnames)); blah = showContextMap(imgName,[],0); output.contextMap = blah.contextMap;  end
    if isempty(strmatch('objectMap',outputnames)); blah = showTargetFeaturesMap(imgName,[],0); output.objectMap = blah.objectMap;  end
    if isempty(strmatch('salMap',outputnames)); blah = showSaliencyMap(imgName,[],0); output.salMap = blah.salMap;  end

    output.tripleSourceMap = (output.contextMap.^(1-gamEx.Triple1-gamEx.Triple2)).*(output.salMap.^gamEx.Triple1).*(output.objectMap.^gamEx.Triple2);
    
    map = output.tripleSourceMap; MapName = 'THREE-Sources of Guiance: CONTEXT x TARGET x SALIENCY Models';
    figure('Position',fullscn,'Name',MapName);
    subplot(2,2,1:2)
    imagesc(img); axis off equal; hold on; title([ 'Original Stimulus: ' imgName2]);
    for t = 1:2
        [imgFixRegions_thresh{t}] = findPredRegion(map,Tsal_px(t),0);
        subplot(2,2,t+2)
        imagesc(imgFixRegions_thresh{t}.*double(img(:,:,1)));
        colormap(gray); axis off equal; hold on; title([MapName ': Thresh ' num2str(thresh(t)*100) '%']);
    end
    if doUs.overlayFixations
        for t = 1:2
            [fixIn fixOut] = checkFixRegions(output,map,Tsal_px(t));
            subplot(2,2,t+2)
            plot(fixIn(:,1),fixIn(:,2),'wo','MarkerSize',10,'LineWidth',2,'MarkerFaceColor','g');
            plot(fixOut(:,1),fixOut(:,2),'wo','MarkerSize',10,'LineWidth',2,'MarkerFaceColor','r');
            text(0,size(img,1)+50,[num2str(size(fixIn,1)) ' out of ' num2str(size(fixIn,1)+size(fixOut,1)) ' fixations IN region'])
        end
    end

    %     output.tripleSource_Oracle_Map = (output.contextOracleMap.^(1-gamEx.Triple1-gamEx.Triple2)).*(output.salMap.^gamEx.Triple1).*(output.objectMap.^gamEx.Triple2);
end

%% ==================================== %%
%% (8) SHOW ROC CURVES
%% ==================================== %%
if doUs.calcROC
    %Determine which maps have been requested:
    fldnames = fieldnames(output)';
    whichMapsIdx = strfind(fldnames,'Map');
    for i=1:length(whichMapsIdx); whichMaps(i) = ~isempty(whichMapsIdx{i}); end
    whichMaps = fldnames(find(whichMaps));

    %Calculate ROC for each requested map:
    for i=1:length(whichMaps)
        map = eval(['output.' whichMaps{i}]);
        [aROC, falseAlarmRate, detectionRate] = areaROC(map, fixationMat);
        %Save results in output data struct:
        eval(['output.ROC_' whichMaps{i} '.areaROC = aROC;']);
        eval(['output.ROC_' whichMaps{i} '.faRate = falseAlarmRate;']);
        eval(['output.ROC_' whichMaps{i} '.detRate = detectionRate;']);
    end
end

%% ==================================== %%
%% HELPER FUNCTIONS %%
%% ==================================== %%
function init(HOME)
%% Set up folders in path
addpath(genpath(fullfile(HOME, 'Tools' )))
addpath(fullfile(HOME,'LabelMeToolbox'))
addpath(fullfile(HOME, 'Dataset_STIMULI'))
addpath(fullfile(HOME, 'Dataset_ContextOracle'))
addpath(fullfile(HOME, 'Dataset_EyeData','eyeData'))
addpath(fullfile(HOME,'targetfeatureMaps','INRIAmaps2000ped_64'));
addpath(fullfile(HOME,'saliencyMaps'));

function [fixIn fixOut] = checkFixRegions(output,map,tsal_px)
%Check whether observer fixations <output.recordOfFixs> fall within
% the prediction map <map> at a particular threshold <tsal_px>
i = [output.recordOfFixs{2:end,4}]; j = [output.recordOfFixs{2:end,3}];
Nfixs = length(i);
imgReg = findPredRegion(map,tsal_px,0);
for f=1:Nfixs; 
    fixIdx(f)= imgReg(i(f),j(f)) > 0; ;  
end;
fixIn = find(fixIdx); 
fixOut = setdiff([1:Nfixs],fixIn);
if isempty(fixIn); fixIn = zeros(1,2); fixIn(1,:) = [];
else;              fixIn = [output.recordOfFixs{fixIn+1,3}; output.recordOfFixs{fixIn+1,4};]';
end
if isempty(fixOut); fixOut = zeros(1,2); fixOut(1,:) = [];
else;               fixOut = [output.recordOfFixs{fixOut+1,3}; output.recordOfFixs{fixOut+1,4};]';
end

function [b] = convYN(a); if strmatch(a,'y'); b = 1; else b = 0; end
