% Script generated by Brainstorm (11-Mar-2024)

% Input files
cd('/users/bkavanau/data/bkavanau/TEMP_T999_dots/')
EEGFiles = dir('*.set');
filelist = extractfield(EEGFiles,'name')
ID_List=strrep(filelist,'TB2_','');
ID_List=strrep(ID_List,'TEMP_','');
ID_List = extractBefore(ID_List,4)

ID_List = cellfun(@str2num,ID_List)
[bs, vec] = sort(ID_List(:).');
uvec(vec) = [true, diff(bs) ~= 0];
ID_List_nodup    = ID_List(uvec);

ID_List_cell_nodup = num2cell(ID_List_nodup)'

%%
sFiles = [];

SubjectNames = ID_List_cell_nodup;

  % Loop on subjects
for subj_i = 1:length(SubjectNames)
    subj_ID = cell2mat(SubjectNames(subj_i))
    MRI_Files = sprintf('/users/bkavanau/data/temp/derivatives/fmriprep/sourcedata/freesurfer/sub-%03d', subj_ID)

    % Process: Import anatomy folder. Each participant has 1 MRI (but
    % multiple EEG files).
    sFiles = bst_process('CallProcess', 'process_import_anatomy', sFiles, [], ...
        'subjectname', SubjectNames{1}, ...
        'mrifile',     {MRI_Files, 'FreeSurfer+Thick'}, ...
        'nvertices',   15000, ...
        'nas',         [0, 0, 0], ...
        'lpa',         [0, 0, 0], ...
        'rpa',         [0, 0, 0], ...
        'ac',          [0, 0, 0], ...
        'pc',          [0, 0, 0], ...
        'ih',          [0, 0, 0]);
    
    % Process: MNI normalization
    sFiles = bst_process('CallProcess', 'process_mni_normalize', sFiles, [], ...
        'subjectname', SubjectNames{1}, ...
        'method',      'maff8', ...  % maff8:Affine registration using SPM mutual information algorithm.Estimates a simple 4x4 linear transformation to the MNI space.Included in Brainstorm.
        'uset2',       0)

    
% here is when loop EEGs. 
    cd('/users/bkavanau/data/bkavanau/TEMP_T999_dots/')
    EEGFiles = dir('*.set');
    files = fullfile('*', sprintf('%03d', subj_ID), '*.set')
    files_replace=strrep(files,'/','')
    eeg_dir = dir(files_replace)
    eeg_list = extractfield(eeg_dir,'name')
    eeg_path = extractfield(eeg_dir,'folder')
    eeg_list_path = fullfile(eeg_path,eeg_list)

  for eeg_i = 1:length(eeg_list) %now here hoping to loop through each of the participant's EEG sessions while keeping active their single MRI.

        % Process: Create link to raw file
        sFiles = bst_process('CallProcess', 'process_import_data_raw', sFiles, [], ...
            'subjectname',    SubjectNames{1}, ...
            'datafile',       {eeg_list_path(eeg_i), 'EEG-EEGLAB'}, ...
            'channelreplace', 0, ...
            'channelalign',   0, ...
            'evtmode',        'value');

        % Process: Import MEG/EEG: Existing epochs. The EEG data files are
        % epoched (and already EEGLAB preprocessed) trials from a sternberg
        % working memory task (each trial is ~6 seconds).
        sFiles = bst_process('CallProcess', 'process_import_data_epoch', sFiles, [], ...
            'subjectname',   SubjectNames{1}, ...
            'condition',     '', ...
            'iepochs',       [], ...
            'eventtypes',    '', ...
            'createcond',    0, ...
            'usectfcomp',    1, ...
            'usessp',        1, ...
            'freq',          [], ...
            'baseline',      [], ...
            'blsensortypes', 'MEG, EEG');
        
        % Process: Add EEG positions
        sFiles = bst_process('CallProcess', 'process_channel_addloc', sFiles, [], ...
            'channelfile', {'', ''}, ...
            'usedefault',  'ICBM152: ANT Waveguard 64', ...  % ICBM152: ANT Waveguard 64
            'fixunits',    0, ...
            'vox2ras',     0, ...
            'mrifile',     {'', ''}, ...
            'fiducials',   []);
        
        % Process: Remove head points
        sFiles = bst_process('CallProcess', 'process_headpoints_remove', sFiles, [], ...
            'zlimit', 1000);
        
        % Process: Project electrodes on scalp
        sFiles = bst_process('CallProcess', 'process_channel_project', sFiles, [], ...
            'sensortypes', 'EEG');
        
        % Process: Refine registration
        sFiles = bst_process('CallProcess', 'process_headpoints_refine', sFiles, [], ...
            'tolerance', 0);
        
        % Process: Generate BEM surfaces
        sFiles = bst_process('CallProcess', 'process_generate_bem', sFiles, [], ...
            'subjectname', SubjectNames{1}, ...
            'nscalp',      1922, ...
            'nouter',      1922, ...
            'ninner',      1922, ...
            'thickness',   4, ...
            'method',      'brainstorm');  % Brainstorm
        
        % Process: Compute head model
        sFiles = bst_process('CallProcess', 'process_headmodel', sFiles, [], ...
            'Comment',     '', ...
            'sourcespace', 1, ...  % Cortex surface
            'meg',         3, ...  % Overlapping spheres
            'eeg',         3, ...  % OpenMEEG BEM
            'ecog',        2, ...  % OpenMEEG BEM
            'seeg',        2, ...  % OpenMEEG BEM
            'openmeeg',    struct(...
                 'BemFiles',     {{}}, ...
                 'BemNames',     {{'Scalp', 'Skull', 'Brain'}}, ...
                 'BemCond',      [1, 0.0125, 1], ...
                 'BemSelect',    [1, 1, 1], ...
                 'isAdjoint',    0, ...
                 'isAdaptative', 1, ...
                 'isSplit',      0, ...
                 'SplitLength',  4000), ...
            'channelfile', '');
        
        % Process: Compute covariance (noise or data). I'm not sure if this
        % is correct, particularly my 'datatimewindow', as my data is
        % structured in 6-second epoched files.
        sFiles = bst_process('CallProcess', 'process_noisecov', sFiles, [], ...
            'baseline',       [], ...
            'datatimewindow', [0, 984.972], ...
            'sensortypes',    'MEG, EEG, SEEG, ECOG', ...
            'target',         1, ...  % Noise covariance     (covariance over baseline time window)
            'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
            'identity',       0, ...
            'copycond',       0, ...
            'copysubj',       0, ...
            'copymatch',      0, ...
            'replacefile',    1);  % Replace
        
        % Process: Compute sources [2018]
        sFiles = bst_process('CallProcess', 'process_inverse_2018', sFiles, [], ...
            'output',  2, ...  % Kernel only: one per file
            'inverse', struct(...
                 'Comment',        'MN: EEG', ...
                 'InverseMethod',  'minnorm', ...
                 'InverseMeasure', 'amplitude', ...
                 'SourceOrient',   {{'fixed'}}, ...
                 'Loose',          0.2, ...
                 'UseDepth',       1, ...
                 'WeightExp',      0.5, ...
                 'WeightLimit',    10, ...
                 'NoiseMethod',    'reg', ...
                 'NoiseReg',       0.1, ...
                 'SnrMethod',      'fixed', ...
                 'SnrRms',         1e-06, ...
                 'SnrFixed',       3, ...
                 'ComputeKernel',  1, ...
                 'DataTypes',      {{'EEG'}}));
        
        % Process: Scout time series: A39rd_L A39rd_R A40rd_L A40rd_R
        % A9/46d_L A9/46d_R A9/46v_L A9/46v_R. specific areas of the
        % brainnetome atlas to extract.
        sFiles = bst_process('CallProcess', 'process_extract_scout', sFiles, [], ...
            'timewindow',     [], ...
            'scouts',         {'Brainnetome', {'A39rd_L', 'A39rd_R', 'A40rd_L', 'A40rd_R', 'A9/46d_L', 'A9/46d_R', 'A9/46v_L', 'A9/46v_R'}}, ...
            'flatten',        1, ...
            'scoutfunc',      'mean', ...  % Mean
            'pcaedit',        struct(...
                 'Method',         'pca', ...
                 'Baseline',       [-0.1, 0], ...
                 'DataTimeWindow', [0, 1], ...
                 'RemoveDcOffset', 'file'), ...
            'isflip',         1, ...
            'isnorm',         0, ...
            'concatenate',    1, ...
            'save',           1, ...
            'addrowcomment',  1, ...
            'addfilecomment', []);
        
        % here i'm hoping to export the raw scout time series for each
        % scout (8 total scouts by time). I believe the structure of the
        % .mat file is rows are 8 scouts per epoch and columns are time
        % points.
        export_name = fullfile(eeg_list_path(eeg_i),'_BRAINSTORM')
        export_name_name=strrep(export_name,'.set/','')
        sFiles = bst_process('CallProcess', 'process_export_file', sFiles, [], ...
        'exportmatrix', {export_name_name{1}, 'BST'});

  end
end

