function varargout = process_mfip_ExcitabilityIndex( varargin )
% PROCESS compute excitability measure based on phase synchrony
% Giovanni Pellegrino, MD, PhD, Giulio Ferrazzi, PhD, Venice, 2021
eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>

sProcess = db_template('ProcessDesc');
% Description the process
sProcess.Comment     = 'Excitability Index';
sProcess.Category    = 'Custom';
sProcess.SubGroup    = 'mfip';
sProcess.Index       = 6703;
sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Connectivity';
% Definition of the input accepted by this process
sProcess.InputTypes  = {'results'};
sProcess.OutputTypes = {'results'};
sProcess.nInputs     = 1;
sProcess.nMinFiles   = 1;

% === TIME WINDOW
sProcess.options.timewindow.Comment = 'Time window:';
sProcess.options.timewindow.Type    = 'timewindow';
sProcess.options.timewindow.Value   = [];
sProcess.options.timewindow.Group   = 'input';
% === Scouts
sProcess.options.scouts.Comment = '';
sProcess.options.scouts.Type    = 'scout_confirm';
sProcess.options.scouts.Value   = {};
sProcess.options.scouts.InputTypes = {'results'};
% === FREQ BANDS
sProcess.options.freqbands.Comment = 'Frequency bands for the Hilbert transform:';
sProcess.options.freqbands.Type    = 'groupbands';
sProcess.options.freqbands.Value   = bst_get('DefaultFreqBands');
% === AVERAGE SIGNALS BEFORE COMPUTING MEASURE
sProcess.options.AveFirst.Comment    = 'Average signals before computing measure (default yes!)';
sProcess.options.AveFirst.Type       = 'checkbox';
sProcess.options.AveFirst.Value      = 1;
% == Size of moving cyrcle
sProcess.options.circlesize.Comment = 'Size of the moving circle in vertices (only if not working with scout):';
sProcess.options.circlesize.Type    = 'value';
sProcess.options.circlesize.Value   = {0,'neighbours vertices', 0};

end


%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
Comment = sProcess.Comment;
end

function OutputFiles = Run(sProcess, sInput)
OutputFiles = {};

%WOrk on cortex and compute  Geodesic Distances if needed
sSource = load(file_fullpath(sInput(1).FileName)); % load input file
Cortex = load(file_fullpath(sSource.SurfaceFile)); %load cortex related to input file

if isempty(sProcess.options.scouts.Value)
    
    disp('Computing geodesic distances between vertices ...')
    tGeodesic=tic;
    
    Atlas_side = Cortex.Atlas(strcmp({Cortex.Atlas.Name},'Structures')); % Load data identifing right and left hemispheres
    VerticesLeft = Atlas_side.Scouts(strcmp({Atlas_side.Scouts.Label},'Cortex L')); %Identify vertices of left hemisphere
    VerticesRight = Atlas_side.Scouts(strcmp({Atlas_side.Scouts.Label},'Cortex R')); %Identify vertices of right hemisphere
    if ~isempty(VerticesLeft) % This if is because cortex may not comprise both hemispheres
        VerticesLeft = VerticesLeft.Vertices;
    end
    if ~isempty(VerticesRight)
        VerticesRight = VerticesRight.Vertices;
    end
    % Geodesic Distance between vertices % belonging to different hemispheres will be inf.
    VertConnTriu = triu(Cortex.VertConn);
    [x,y] = find(VertConnTriu);
    D = sparse(size(VertConnTriu,1),size(VertConnTriu,2));
    ind = sub2ind(size(VertConnTriu),x,y);
    D(ind) = sqrt((Cortex.Vertices(x,1) - Cortex.Vertices(y,1)).^2 + ...
        (Cortex.Vertices(x,2) - Cortex.Vertices(y,2)).^2 + ...
        (Cortex.Vertices(x,3) - Cortex.Vertices(y,3)).^2);
    D = D + D';
    [Geodesic_Distance] = graphallshortestpaths(D)*1000;
    tEndGeodesic = toc(tGeodesic);
    disp(['It took ' num2str(tEndGeodesic) ' seconds to compute geodesic distances between vertices']);
    clear x y VertConnTriu tGeodesic tEndGeodesic D ind
end
% end geodesic

% Load data (vertices x time or vertices x time x trials)
if sProcess.options.AveFirst.Value == 1
    %[Stat, Messages, iOutFiles, AllEvents] = bst_avg_files(FilesListA, FilesListB, Function, isVariance, isWeighted, isMatchRows, isZeroBad, isPercent)
    DataInputAve = bst_avg_files({sInput.FileName}, [], 'mean', 0, 0, 0, 1, 0);
    DataInput = DataInputAve.mean;
else
    % Loop over input files (fix this that takes too much memory)
    % DataInput has the following dimensions: vertices x time x trials
    for iFile = 1:length(sInput)
        % Load reference signal
        sInputA = bst_process('LoadInputFile', sInput(iFile).FileName);
        if iFile == 1
            DataInput = sInputA.Data;
        else
            DataInput (:,:,1+iFile) = sInputA.Data;
        end
        
    end
end

% Take only the time of interest
sInputA = bst_process('LoadInputFile', sInput(1).FileName);

% if ~isempty(sProcess.options.timewindow.Value{1})
%     iTime = panel_time('GetTimeIndices', sInputA.Time, cell2mat(sProcess.options.timewindow.Value(1)));
%     DataInput = DataInput(:,iTime,:);
% else
    iTime = 1:size(DataInput,2);
% end

Output = zeros(size(DataInput(:,:,1)));
% End load data

% Sampling frequency and round
sfreq = 1 ./ (sInputA.Time(2) - sInputA.Time(1));
sfreq = round(sfreq * 1e6) * 1e-6;

% loop over frequency bands
for we = 1 : size(sProcess.options.freqbands.Value,1)
    
    % Here band by band, the measure is computed
    
    % if an atlas was selected, then goes by ROIs
    if ~isempty(sProcess.options.scouts.Value)
        
        scoutName = sProcess.options.scouts.Value{1,1}; %identify Atlas Name
        All_ROIs = Cortex.Atlas(strcmp({Cortex.Atlas.Name},scoutName)); %find all ROIs of that Atlas
        
        for wsd = 1:length (All_ROIs.Scouts) %Loop over parcels
            
            disp(['Computing measure for -->    ' num2str(wsd) '/' num2str(length (All_ROIs.Scouts)) '  -->  ' All_ROIs.Scouts(wsd).Label])
            
            VerticesROI = All_ROIs.Scouts(strcmp({All_ROIs.Scouts.Label},All_ROIs.Scouts(wsd).Label)).Vertices; %Find vertices belonging to ROI
            SignalROI = DataInput(VerticesROI,:,:);
            %[Nodes Time Trial] = size(SignalROI);
            E = [];
            for t = 1 : size(SignalROI,3)
                E = [E; my_phasesync(SignalROI(:,:,t),sfreq,str2num(sProcess.options.freqbands.Value{we,2}))];
            end
            E = mean(E,1); % mean trial
            % R = mean(E,2); % Uncomment this if you wanna average over
            % time. Remember also to change output structure (no time, only
            % one map)
            
            Output(VerticesROI,:) = repmat(E,[size(SignalROI,1) 1]); %Project same value to all vertices of ROI
            
        end
    end
    
    
    % if no atlas is selected start from a vertex and move forward
    if isempty(sProcess.options.scouts.Value)
        
        
        Size_moving_circle = sProcess.options.circlesize.Value{1,1};
        Atlas_side = Cortex.Atlas(strcmp({Cortex.Atlas.Name},'Structures')); % Load data identifing right and left hemispheres
        VerticesLeft = Atlas_side.Scouts(strcmp({Atlas_side.Scouts.Label},'Cortex L')); %Identify vertices of left hemisphere
        VerticesRight = Atlas_side.Scouts(strcmp({Atlas_side.Scouts.Label},'Cortex R')); %Identify vertices of right hemisphere
        
        %Left hemisphere
        for iLeft=1:length(VerticesLeft.Vertices)
            
            %disp(['Doing L Hemisphere --->   ' num2str(iLeft) '/' num2str(length(VerticesLeft.Vertices))])
            disp(['Doing L Hemisphere'])
            
            
            %find row with iLeft Vertex, sort and take the smallest vertices
            [xs, index] = sort(Geodesic_Distance(VerticesLeft.Vertices(iLeft),:));
            index = index(1:Size_moving_circle);
            SignalROI = DataInput(index,:,:);
            
            [Nodes Time Trial] = size(SignalROI);
            E = [];
            for t = 1 : Trial
                E = [E; my_phasesync(SignalROI(:,:,t),sfreq,str2num(sProcess.options.freqbands.Value{we,2}))];
            end
            E = mean(E,1); % mean trial
            R = mean(E,2);
            
            Output(index(1),:) = E;
            
        end
        
        %Right hemisphere
        for iRight=1:length(VerticesRight.Vertices)
            
            %disp(['Doing R Hemisphere --->   ' num2str(iRight) '/' num2str(length(VerticesRight.Vertices))])
            disp(['Doing R Hemisphere'])
            
            %find row with iRight Vertex, sort and take the smallest vertices
            [xs, index] = sort(Geodesic_Distance(VerticesRight.Vertices(iRight),:));
            index = index(1:Size_moving_circle);
            SignalROI = DataInput(index,:,:);
            
            [Nodes Time Trial] = size(SignalROI);
            E = [];
            for t = 1 : Trial
                E = [E; my_phasesync(SignalROI(:,:,t),sfreq,str2num(sProcess.options.freqbands.Value{we,2}))];
            end
            E = mean(E,1); % mean trial
            R = mean(E,2);
            
            Output(index(1),:) = E;
            
        end
    end
    
    Output(isnan(Output)) = 0;
    
    sResult = sSource;
    sResult.ImageGridAmp=Output;
    sResult.Time = sInputA.Time(iTime);
    sResult.ImagingKernel = [];
    sResult.Whitener = [];
    sResult.History = [];
    sResult.OriginalOptions = sProcess.options;
    sResult.Comment = [sInput(1).SubjectName '_' sInput(1).Condition '_ExcitabilityIndex | ' sProcess.options.freqbands.Value{we,1}];
    db_add(sInput(1).iStudy, sResult, 0);
    db_reload_studies(sInput(1).iStudy);
    panel_protocols('UpdateTree');
    
end
end

function E = my_phasesync(X,fs,iband)

Wn=iband/(fs/2);
[b1,b2]=butter(2,Wn);

[verticies time] = size(X);

XF = zeros(size(X));
H = zeros(size(X));
for v = 1 : verticies
    x = X(v,:);
    x = filtfilt(b1,b2,x);
    XF(v,:) = x;
    H(v,:) = angle(hilbert(x));
end
E = abs(sum(exp(1i.*H),1))./verticies;
end
