%% HistoEmotions is function that plots number of fixations in each roi for chosen emotion as well as average duration of fix and statistics
%% Input data
% path: path to fixation data
% emotion: choose the emotion
% mini: start of interval in video in ms
% maxi: end of interval in video in ms
% efed: experiment (EFE1/EFE2; EFED1/EFED2)
% limitdwn: minimum duration of fixation in ms
% limitup: maximum duration of fixation in ms
%% Output data
% statistics: gives percentage of fixation in each roi from the total
% VidInd: Gives index of video/image of the emotion
% NbFx: Number of fixation in each video/image

function [statistics, VidInd, NbFx]=HistoEmotions (path, emotion, mini, maxi, efed, limitdwn, limitup)
%% preset 
if mini==0
    mini=1;
else 
    mini=mini;
end
switch emotion
    case 'Neutral'
        emo= 'N';
    case 'Fear'
        emo= 'F';
    case 'Disgust'
        emo= 'D';
    case 'Surprise'
        emo= 'S';
    case 'Happiness'
        emo= 'G';
end

%% load data
load(path)

%% choosing the particular videos for the chosen emotion from all subjects
FIXperemo=zeros(1,13);   % initiate vector to save fixations
DURperemo=zeros(1,13);   % initiate vector to save durations
KO=zeros(17,13,400);     % 
for i=1:length(fieldnames(EventFix.(efed))) % for all subjects 
    subj=fieldnames(EventFix.(efed));
    subj=subj(i,1);
    subj=subj{1,1};      % picks i subject
    subjDUR=zeros(1,13); % empty vector to save durations in subject
    subjFIX=zeros(1,13); % empty vector to save fixations in subject
    k=1;
    for j=1:length(EventFix.(efed).(subj))
        nature=EventFix.(efed).(subj)(1,j).Nature; 
        if strcmp (nature, emo) % to check for the right video
            if isempty(EventFix.(efed).(subj)(1,j).Roi) %if doesnt contain ROI, then move on
                fix=zeros(1,13);
                duration=zeros(1,13);
            else
                mat=EventFix.(efed).(subj)(1, j).FixInVid; % load the matrix of fixations                 
                NbFix(i,k)=EventFix.(efed).(subj)(1,j).NbFix; % saves total number of fixation for all subjects and all videos
             %% cut to interval and optimalization (takes the whole duration of fixation that belongs in the interval)
                l=1;
                m=1;
                if sum(mat(:,mini))==1 
                    while sum(mat(:,mini-l))==1
                        l=l+1;
                    end
                    minim=mini-l;
                else 
                    minim=mini;
                end
                if sum(mat(:,maxi))==1 
                    while sum(mat(:,maxi+m))==1
                        m=m+1;
                    end
                    maxim=maxi+m;
                else 
                    maxim=maxi;
                end
                mat(:,[1:minim])=0;
                mat(:,[maxim:end])=0; % now we have ones only where the fixation appear in given interval    
                %% Number of fixations in the interval
                fix=zeros(1,13);
                duration=zeros(1,13);
                B.NumObjects=[];
                for pp=1:9     % number of rows (just ROI!!!)
                    B = bwconncomp(mat(pp, [mini:maxi])); % connects components in each row (1 next to each other)
                    pom=0;
                    dur=0;
                    for o=1:B.NumObjects % Check length of every fixation if it fits the condition.
                        if (length(B.PixelIdxList{1,o}))>=limitup || (length(B.PixelIdxList{1,o}))<=limitdwn
                            pom=pom+1;
                         else
                            pom=pom+0;
                            dur=length(B.PixelIdxList{1,o})+dur; % sums all durations in the row that belong to the interval
                        end
                    end
                    fix(pp)=B.NumObjects-pom; % gives number of fixations in each row (sum gives total number of fix)
                    duration(pp)=dur; % save the sum of durations into the belonging roi
                end
                NbFx(i,k)=sum(fix); % matrix for number of fixation, each row = subject. Without out of ROI
                VidInd(i,k)=j;      % gives index of particular videos
                k=k+1;
            end
        else
            fix=zeros(1,13);  
            duration=zeros(1,13);
        end
        subjFIX=subjFIX+fix;
        subjDUR=subjDUR+duration;
    end
    FIXperemo=subjFIX+FIXperemo; % final vector of fixations in each ROI
    DURperemo=subjDUR+DURperemo; % sum of durations in each ROI
end
%% 3d matrix where i=subject, j=ROI, k= fixation, in cell its duration
kap(:,1,:)=cat(3,KO(:,1,:),KO(:,2,:));
kap(:,2,:)=cat(3,KO(:,3,:),KO(:,7,:));
kap(:,3:5,:)=cat(3,KO(:,4:6,:),zeros(17,3,400));
kap(:,6:7,:)=cat(3,KO(:,8:9,:),zeros(17,2,400));
kap(kap==0)=NaN;
kap=sort(kap,3);

%% regrouping new: Eyes, Brows, Mouth, Chin, Nose, Corr, Front
emoDUR=zeros(1,6);
emoFIX=zeros(1,6);
emoDUR(1,1)=DURperemo(1,1)+DURperemo(1,2); %Eyes
emoDUR(1,2)=DURperemo(1,3)+DURperemo(1,7); %Brows
emoDUR(1,3:5)=DURperemo(1,4:6);
emoDUR(1,6:7)=DURperemo(1,8:9);
emoFIX(1,1)=FIXperemo(1,1)+FIXperemo(1,2);
emoFIX(1,2)=FIXperemo(1,3)+FIXperemo(1,7);
emoFIX(1,3:5)=FIXperemo(1,4:6);
emoFIX(1,6:7)=FIXperemo(1,8:9);
emoFIX=emoFIX';
emoDUR=emoDUR';
emoDUR=round(emoDUR./emoFIX); % average duration in each ROI

%% histogram
x=1:1:7;
% bar properties
offset = (x(2)-x(1))/5;
width = (x(2)-x(1))/2.5;
colors = {'r','g'};
% plot
[ax, hBar1, hBar2]=plotyy(x-offset,emoDUR,x+offset,emoFIX, @(x,y) bar(x,y,width,colors{1}), @(x,y) bar(x,y,width,colors{2}));
Labels = {'Eyes','Brows','Mouth','Chin','Nose','Corr','Front'};
set(ax, 'XTick', 1:7, 'XTickLabel', Labels, 'fontsize', 14);
title(sprintf('Emotion: %s, %d-%d ms, experiment: %s', emotion, mini, maxi, efed),'FontSize', 18,'fontweight','bold')
xlabel('ROIs', 'fontsize', 16)
ylabel(ax(1),'Average duration of fixation [ms]', 'color', 'r','fontsize', 16)
ylabel(ax(2),'Number of fixations', 'color', 'g','fontsize', 16)
legend('Average duration of 1 fix in ROI','Number of fixations in ROI')

%% normalization of scales to show text the right way
yLimits = get(ax(1),'YLim');
yLimits2 = get(ax(2),'YLim');
a=0;
A=0;
B=yLimits2(2);
b=yLimits(2);
for i=1:length(emoFIX)
apt(i)=a+(emoFIX(i)-A)*(b-a)/(B-A);
end
text(x-offset,emoDUR,num2str(emoDUR),... 
'HorizontalAlignment','center',... 
'VerticalAlignment','bottom')
text(x+offset,apt,num2str(emoFIX),... 
'HorizontalAlignment','center',... 
'VerticalAlignment','bottom')

%% percentage of fixations in ROI
for i=1:length(emoFIX)
    perc(i)=round(emoFIX(i)/sum(emoFIX)*100);
end
statistics={sprintf('Number of fixations: %d in ROI', sum(emoFIX)), sprintf('Eyes: %d %%', perc(1)), sprintf('Brows: %d %%', perc(2)), sprintf('Mouth: %d %%', perc(3)),sprintf('Chin: %d %%', perc(4)), sprintf('Nose: %d %%', perc(5)), sprintf('Corr: %d %%', perc(6)), sprintf('Front: %d %%', perc(7))} ;
