eegmp_preproc/code/a3_ica/a3_ica.m

function a3_ica(id, rootpath)

if ismac % run if function is not pre-compiled
    id = '5'; % test for example subject
    currentFile = mfilename('fullpath');
    [pathstr,~,~] = fileparts(currentFile);
    cd(fullfile(pathstr,'..', '..'))
    rootpath = pwd;
end

% inputs
pn.eeg_BIDS = fullfile(rootpath, 'data', 'inputs', 'rawdata', 'eeg_BIDS');
pn.channel_locations = fullfile(rootpath, 'data', 'inputs', 'rawdata', 'channel_locations');
pn.events = fullfile(rootpath, 'data', 'inputs', 'rawdata', 'events');
pn.tools = fullfile(rootpath, 'tools');
% outputs
pn.eeg_ft = fullfile(rootpath, 'data', 'outputs', 'eeg');
pn.history = fullfile(rootpath, 'data', 'outputs', 'history');

if ismac % run if function is not pre-compiled
     addpath(fullfile(pn.tools, 'fieldtrip')); ft_defaults;
     addpath(fullfile(pn.tools, 'helpers'));
end

%% define IDs for visual screening

% N = 33
IDs = tdfread(fullfile(pn.eeg_BIDS, 'participants.tsv'));
IDs = cellstr(IDs.participant_id);

id = str2num(id);

display(['processing ' num2str(IDs{id})]);

if ~exist(fullfile(pn.eeg_ft, [IDs{id}, '_task-xxxx_eeg_ica.mat']),'file')

    %%  load raw data & exclude parts containing artifacts

    % load config
    config = [];
    load(fullfile(pn.history, [IDs{id}, '_task-xxxx_config.mat']),'config');

    config_tmp.(['trl']) = config.trl_ica1;
    config_tmp.(['visual_inspection']) = config.visual_inspection;

    % define segment(s) to be read by fieldtrip
    cfg.trl = config.trl_ica1;

    % load data
    load(fullfile(pn.eeg_ft, [IDs{id}, '_task-xxxx_eeg_raw.mat']),'data_eeg');

    % remove M2 (collinear with M1 due to avg. mastoid REF: mean(M1,M2)=0)
    cfg_preproc               = [];
    cfg_preproc.channel       = {'all', '-M2'};
    data_eeg = ft_preprocessing(cfg_preproc, data_eeg);
    
    % adjust final length
    cfg.trl(end,2) = size(data_eeg.trial{1},2);

    % "segment" data
    data = cm_segmentation_of_continuous_data_fieldtrip_20150825(data_eeg,cfg);    

    % clear cfg structure
    clear cfg

    %% segmentation for ICA

    % define settings
    cfg.length = 2;
    cfg.n      = 5000;      % keep all possible trials
    cfg.type   = 'rnd';     % select trials randomly if > 1500 trials available
    cfg.seed   = 20170915 + str2num(IDs{id});

    % arbitrary segmentation - segments a 2 sec
    % NOTE original segments will be overwritten
    data = cm_arbitrary_segmentation_fieldtrip_20151002(data,cfg);

    %%  ICA

    % date
    dt = date;

    % ica config
    cfg.method           = 'runica';
    cfg.channel          = {'all'}; % additional channels should be excluded already...
    cfg.trials           = 'all';
    cfg.numcomponent     = 'all';
    cfg.demean           = 'no';
    cfg.runica.extended  = 1;
    cfg.runica.logfile   = fullfile(pn.history, ['log_' IDs{id} '_ica1_' dt '.txt']);

    % run ICA
    icadat = ft_componentanalysis(cfg,data);

    %%  automatic ICA labeling

    [iclabels] = cm_automatic_IC_detection_20170919(data,icadat);

    %%  save data for ICA labeling

    config.trl = config_tmp.trl;

    config.visual_inspection = config_tmp.visual_inspection;

    % - include ICA solution in data
    data.topo 	   = icadat.topo;
    data.unmixing  = icadat.unmixing;
    data.topolabel = icadat.topolabel;
    data.cfg       = icadat.cfg;

    % - include ICA solution in config
    config.ica1.date      = dt;
    config.ica1.topo 	  = icadat.topo;
    config.ica1.unmixing  = icadat.unmixing;
    config.ica1.topolabel = icadat.topolabel;
    config.ica1.cfg       = icadat.cfg;
    config.ica1.iclabels.auto = iclabels;

    % fieldtrip format electrode information
    chanlocs_besa = readtable(fullfile(pn.channel_locations, 'chanlocs_besa.txt'));
    elec.pnt = [chanlocs_besa.Var2, chanlocs_besa.Var3, chanlocs_besa.Var4];
    elec.label = chanlocs_besa.Var1;
    data.elec = elec;
    
    % EEGLAB format electrode information    
    chanlocs_ced = readtable(fullfile(pn.channel_locations, 'chanlocs_ced.txt'));
    chanlocs_ced = table2struct(chanlocs_ced);
    chanlocs_ced = rmfield(chanlocs_ced, {'Number', 'type', 'Var12'});    
    [~, loc] = ismember(data.label, {chanlocs_ced.labels});
    data.chanlocs(1,loc) = chanlocs_ced(loc,1);

    % - include channel information in config
    config.elec     = data.elec;
    config.chanlocs = data.chanlocs;

    % keep ICA labels
    data.iclabels = iclabels;

    % save data
    save(fullfile(pn.eeg_ft, [IDs{id}, '_task-xxxx_eeg_ica.mat']),'data');
    
    % save config
    config.preproc_version = '20170915';
    save(fullfile(pn.history, [IDs{id}, '_task-xxxx_config.mat']),'config');

 end % file available