%%% Function to find all subjects and respective channels that contribute
%%% to a specified region of interest

% input: all_subjects e.g. [1:22 24:28 30:37 39:49]
% atlasType: 'BA' 
% atlasIdx: [17 18 19],[37],[20],[39 40] %BA codes for Occ, pTC, aTC, LPC
% hemisphere: 'lh','rh','both' %which hemisphere to sample
% bp_chans: 'all' (row+col of grids), 'row' (only row of grids)

function [subjects,numChans,ChanIndices] = findSubs(all_subjects,atlasType,atlasIdx,hemi,bp_chans)
    
    baseDir = 'D:\Elias\iEEG_PFC\ECoG_MemDev\';

    %set dir and fs
    dataDir = [baseDir 'DATA'];
    fs = '\';
    %load subj Information
    load([baseDir 'subject_list.mat']);
    
    %find subjects with channels and collect them
    subjects = [];
    numChans = [];
    ChanIndices = {};
    %for subject = all_subjects
    for snum = 1:length(all_subjects)
        subject = all_subjects(snum);
        load([dataDir fs subNames{subject} fs 'coordinates' fs 'electrodes_bp_atlas_extended.mat'])

        %get idx of electrodes that fall into specified roi
        roi_idx = find(ismember(electrodes.BA,atlasIdx));

        %get idx of electrodes from specified hemisphere
        if strcmp(hemi,'lh')
            hemi_idx = find(strcmp(electrodes.hemisphere,'lh'));
        elseif strcmp(hemi,'rh')
            hemi_idx = find(strcmp(electrodes.hemisphere,'rh'));
        elseif strcmp(hemi,'both')
            hemi_idx = 1:length(electrodes.hemisphere);
        end

        %get idx of electrodes according to bp_chantype (row +/ col of
        %grids
        if strcmp(bp_chans,'all')
            bp_type_idx = 1:length(electrodes.hemisphere);
        elseif strcmp(bp_chans,'row')
            bp_type_idx = find(strcmp(electrodes.chantype,'row'));
        end


        %%% now get idx of all channels that fall into specified ROI
        %%% and that fulfill specified hemisphere and bp_type criteria
        chanIdx = intersect(intersect(bp_type_idx,hemi_idx),roi_idx);
  
        if ~isempty(chanIdx)
            subjects(end+1,1) = subject;
            numChans(end+1,1) = length(chanIdx);
            ChanIndices{end+1,1} = chanIdx;
        end
    end
    %display
    for k = 1:length(subjects)
        subidx = subjects(k);
        disp([subNames{subidx} ' contributes ' num2str(numChans(k)) ' electrodes'])
    end
    disp([num2str(length(subjects)) ' contribute ' num2str(sum(numChans)) ' electrodes'])
end
