s05_concat_data.htb 22.5 KB
Newer Older
1
2
3
4
5
6
7
%% SCRIPT DESCRIPTION
% This script combineds the results of the amicas A,B and C. It takes this information and using ICAtest it determines which components were
% and were not replicated across and flags them, along with time sections with high standard deviation. Beta,theta and Alpha power is calculated 
% arcoss time and is marked but not primed for removal.
%
%% From Config          key_strswap     Description
%-----------------------------------------------------------------------------------------------------------------------------
8
%  in_path = 		[in_path]  	    Relative path to input data files assuming cd = work_path
9
10
%  recur_sec = 		[recur_sec]  	Recurrence (sec) for the current artifact detection epoching (e.g. 1)
%  limit_sec = 		[limit_sec]  	Limits (sec) for the current artifact detection epoching (e.g. [0 1])
11
12
13
14
15
16
%  sd_t_meth =      [sd_t_meth]     Method used for flagging epochs (e.g. 'q' (quantiles), or 'na' (default))
%  sd_t_vals  =     [sd_t_vals]     Percentage trim for confidence intervals during epoch standard deviation criteria (e.g. [.3 .7])
%  sd_t_o =         [sd_t_o]        z threshold for flagging epochs during standard deviation criteria (e.g. 6)
%  sd_t_f_meth =    [sd_t_f_meth]   Fixed method used for flagging epochs (should be 'fixed')
%  sd_t_f_vals  =   [sd_t_f_vals]   Percentage trim for confidence intervals during epoch standard deviation criteria (e.g. [.3 .7], leave empty for 'fixed')
%  sd_t_f_o =       [sd_t_f_o]      z threshold for flagging epochs during fixed standard deviation criteria (e.g. .2)
17
18
%  min_gap_ms = 	[min_gap_ms]  	Minimum time (ms) to allow between periods marked for rejection  
%  icr_crit =       [icr_crit]      False discovery rate of component replicatablilty, measure of sensitivity %
19
20
21
%  dip_vol =        [dip_vol]       file path/name of the standard volume file used in dipole fitting (e.g. derivatives/lossless/code/misc/standard_vol.mat).
%  dip_mri =        [dip_mri]       file path/name of the standard mri file used in dipole fitting (e.g. derivatives/lossless/code/misc/standard_mri.mat).
%  dip_elc =        [dip_elc]       file path/name of the standard electrode file used in dipole fitting (derivatives/lossless/code/misc/standard_1020.elc).
22
23
%  dip_transmat =   [dip_transmat]  Transformation matrix to apply during dipole fit (e.g. [0 0 0 0 0 0 1 1 1])
%  cor_transmat =   [cor_transmat]  The transformation matrix to apply between dipole fitting to the standard head and performing ICLabel to the EEGLAB default orientation (e.g.[0,0,0,0,0,1.57,1,1,1])
24
25
%  cor_nosedir =    [cor_nosedir]   Nose direction of montage to apply between dipole fitting to the standard head and performing ICLabel to the EEGLAB default orientation (e.g. +Y)
%  out_path = 		[out_path]  	Relative path to output data files assuming cd = work_path
26
27
28
29
30
31
32

%% LOAD COMPART DATASET
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loads the data file and the first amica model.

logging_log('INFO', 'Loading set file: [batch_dfn,_,-1]_compart_data.set...');
tic;
33
EEG = pop_loadset('filepath','[in_path]','filename','[batch_dfn,_,-1]_compart_data.set');
34
35
36
37
38
39
40
41
42
43
EEG = eeg_checkset( EEG );
logging_log('INFO', 'TIME TO: LOAD DATASET...');
toc

% Set up the ICA channels
EEG.icachansind = marks_label2index(EEG.marks.chan_info,{'manual','rank'},'indices','invert','on');

% Load Amica A
logging_log('INFO', 'Loading Amica ABC...');
tic;
44
EEG.amica(2).models=loadmodout15('[in_path]/[batch_dfn,_,-1]_amicaout_A');
45
46
47
EEG.icaweights=EEG.amica(2).models.W;
EEG.icasphere=EEG.amica(2).models.S(1:EEG.amica(2).models.num_pcs,:);

48
%% PASS EEG DATA THROUGH THE AMICA MATRIX
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loads the B and C amica models. Amica A is used to fill in the ICA information.

% Instead of EEG=eeg_checkset(EEG);
tmpdata = EEG.data(EEG.icachansind,:);
tmpindices = find(~sum(isnan(tmpdata)));
EEG.icaact = zeros(size(EEG.icaweights,1), size(tmpdata,2)); EEG.icaact(:) = NaN;
EEG.icaact(:,tmpindices) = (EEG.icaweights*EEG.icasphere)*tmpdata(:,tmpindices);
clear tmpdata;

% Create winv for topo plots manually
EEG.icawinv = [];
EEG.icawinv = pinv(EEG.icaweights * EEG.icasphere);

% Load Amica B
64
EEG.amica(3).models=loadmodout15('[in_path]/[batch_dfn,_,-1]_amicaout_B');
65
% Load Amica C
66
EEG.amica(4).models=loadmodout15('[in_path]/[batch_dfn,_,-1]_amicaout_C');
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
logging_log('INFO', 'TIME TO: LOAD ICA ABC INFO...');
toc

%% UPDATE REJECTION STRUCTURE...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Creating flag marks to based on the log liklyhodd models found for each of the ICA models

% Flags for amicas A B and C
logging_log('INFO', 'Updating latflaginfo structure...');
tic;

% A
logl_A_flags=zeros(size(EEG.data(1,:)));
logl_A_flag_inds=marks_label2index(EEG.marks.time_info,{'manual'},'indices','invert','on');

% Check to see if the overwrite will be successful. If the overwirte is a different size then some of the old data will be retained. 
if length(EEG.amica(2).models.Lt) ~= length(logl_A_flag_inds)
    EEG.amica(2).models.Lt = EEG.amica(2).models.Lt(1:length(logl_A_flag_inds));
    warning('YOU FORGOT TO DELETE YOUR OLD AMICA FILES AND THE NEW CHANNEL SIZES DO NOT MATCH UP');
end

logl_A_flags(logl_A_flag_inds)=(EEG.amica(2).models.Lt-min(EEG.amica(2).models.Lt))/(max(EEG.amica(2).models.Lt)-min(EEG.amica(2).models.Lt));
EEG.marks=marks_add_label(EEG.marks,'time_info', ...
	{'logl_A',[1,0.45,0],logl_A_flags});
clear logl_A_flags;
% B
logl_B_flags=zeros(size(EEG.data(1,:)));
logl_B_flag_inds=marks_label2index(EEG.marks.time_info,{'manual'},'indices','invert','on');
logl_B_flags(logl_B_flag_inds)=(EEG.amica(3).models.Lt-min(EEG.amica(3).models.Lt))/(max(EEG.amica(3).models.Lt)-min(EEG.amica(3).models.Lt));
EEG.marks=marks_add_label(EEG.marks,'time_info', ...
	{'logl_B',[1,0.45,0],logl_B_flags});
clear logl_B_flags;
% C
logl_C_flags=zeros(size(EEG.data(1,:)));
logl_C_flag_inds=marks_label2index(EEG.marks.time_info,{'manual'},'indices','invert','on');
logl_C_flags(logl_C_flag_inds)=(EEG.amica(4).models.Lt-min(EEG.amica(4).models.Lt))/(max(EEG.amica(4).models.Lt)-min(EEG.amica(4).models.Lt));
EEG.marks=marks_add_label(EEG.marks,'time_info', ...
	{'logl_C',[1,0.45,0],logl_C_flags});
clear logl_C_flags;

logging_log('INFO', 'TIME TO: UPDATE REJECTION STRUCTURE...');
toc

%% IDENTIFY RELIABLE COMPONENTS USING ISCTEST
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Gives the A variable of each of the amicas to the function icatest. Returned are the indecies of the 
% Components that were not replicated across all 3 of the models.

tic;
% Clear a in case of local models are saved already
a = [];
for i=2:length(EEG.amica);
    a(:,:,i-1)=EEG.amica(i).models(1).A;
end

% Calculating the replicated component flag
% Replicate criteria changed from 0.05,0.05 to 0.1,0.1 to be less sensitive
ca=isctest(a,[icr_crit],[icr_crit],'mixing');
if isempty(ca)
    mess = ['ISCTest failed to return reasonable values. ISCTest is ' ...
        'meant to find duplicate channels in amica_{A,B,C} if this ' ...
        'failed consider investigating the datafile'];
    logging_log('ERROR', mess);
    error(mess);
end
icr_inds=setdiff(1:min(size(EEG.amica(2).models(1).A)),ca(:,1));

% Creating the replicated component flag
logging_log('INFO', 'Updating chflaginfo structure...');
icr_comp_flags=zeros(min(size(EEG.amica(2).models(1).A)),1);
icr_comp_flags(icr_inds)=1;
EEG.marks=marks_add_label(EEG.marks,'comp_info', ...
	{'ic_rt',[0,0,1],[0,0,1],-1,icr_comp_flags});
logging_log('INFO', 'TIME TO: EDIT FIND REPLICATED COMPONENTS AND UPDATE COMPINFO STRUCT...');
toc

% Window the continuous data
logging_log('INFO', 'Windowing the continous data...');
tic;
EEG=marks_continuous2epochs(EEG,'recurrence',[recur_sec],'limits',[limit_sec]);
logging_log('INFO', 'TIME TO: WINDOW THE CONTINUOUS DATA...');
toc

%% CALCULATE ICAACT SD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Similarly to the IC_SD flag created in compart. Marks sections of time where the components
% have an abnormally high standard deviation.

% Calculating the replicated component flag
logging_log('INFO', 'Calculating the IC sd array for window criteria...');
tic;
epoch_inds=marks_label2index(EEG.marks.time_info,{'manual'},'indexes','invert','on');
159
[EEG,icaact_sd2_t]=chan_variance(EEG,'data_field','icaact', ...
160
161
162
163
164
165
166
		'epoch_inds',epoch_inds, ...
		'varmeasure','absmean', ...
        'detrend','off', ...
        'plot_figs','off');

% Creating the ic_sd component flag
logging_log('INFO', 'Assessing window icaact sd distributions...')
167
[~,flag_sd_t_inds]=marks_array2flags(icaact_sd2_t, ...
168
    'flag_dim','col', ...
169
    'init_method','[sd_t_meth]', ...
170
171
    'init_vals',[sd_t_vals], ...
    'init_crit',[sd_t_o], ...
172
173
174
    'flag_method','[sd_t_f_meth]', ...
    'flag_vals',[ [sd_t_f_vals] ], ...
    'flag_crit',[sd_t_f_o], ...
175
176
177
178
179
180
181
182
183
184
185
186
187
    'plot_figs','off');

% Updating rejection structure
logging_log('INFO', 'Updating latflaginfo structure...');

icsd_epoch_flags=zeros(size(EEG.data(1,:,:)));
icsd_epoch_flags(1,:,epoch_inds(flag_sd_t_inds))=1;
EEG.marks=marks_add_label(EEG.marks,'time_info', ...
	{'ic_sd2',[.4,.6,.6],icsd_epoch_flags});
logging_log('INFO', 'TIME TO: EDIT FIND IC_SD AND UPDATE TIMEINFO REJECTION STRUCT...');
toc
clear icsd_epoch_flags;

188
%% CALCULATE DELTA/THETA ICAACT POWER
189
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
190
% Calculates and flags times of high Theta/Alpha Power
191

192
193
% Calculating the Delta/Theta component flag
logging_log('INFO', 'Calculating the delta/theta IC power array for window criteria...');
194
tic;
195
[EEG,icaact_dt_t]=chan_variance(EEG,'data_field','icaact', ...
196
197
198
		'epoch_inds',epoch_inds, ...
		'detrend','on', ...
        'varmeasure','spect', ...
199
        'spectrange',[0 7], ...
200
201
        'plot_figs','off');

202
% Creating the Delta/Theta component flag
203
logging_log('INFO', 'Assessing window icaact sd distributions...')
204
[~,flag_dt_t_inds]=marks_array2flags(icaact_dt_t, ...
205
    'flag_dim','col', ...
206
    'init_method','[sd_t_meth]', ...
207
208
    'init_vals',[sd_t_vals], ...
    'init_crit',[sd_t_o], ...
209
210
211
    'flag_method','[sd_t_f_meth]', ...
    'flag_vals',[ [sd_t_f_vals] ], ...
    'flag_crit',[sd_t_f_o], ...
212
213
214
215
    'plot_figs','off');

% Updating rejection structure
logging_log('INFO', 'Updating latflaginfo structure...');
216
217
icdt_epoch_flags=zeros(size(EEG.data(1,:,:)));
icdt_epoch_flags(1,:,epoch_inds(flag_dt_t_inds))=1;
218
EEG.marks=marks_add_label(EEG.marks,'time_info', ...
219
220
	{'ic_dt',[1,1,0],icdt_epoch_flags});
logging_log('INFO', 'TIME TO: CALCULATE DELTA/THETA ICAACT POWER AND UPDATE MARKS...');
221
toc
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
clear icdt_epoch_flags;

%% CALCULATE ALPHA ICAACT POWER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculates and flags times of high Alpha Power

% Calculating the Alpha component flag
logging_log('INFO', 'Calculating the alpha IC power array for window criteria...');
tic;
[EEG,icaact_a_t]=chan_variance(EEG,'data_field','icaact', ...
		'epoch_inds',epoch_inds, ...
		'detrend','on', ...
        'varmeasure','spect', ...
        'spectrange',[8 13], ...
        'plot_figs','off');

% Creating the Alpha component flag
logging_log('INFO', 'Assessing window icaact sd distributions...')
[~,flag_a_t_inds]=marks_array2flags(icaact_a_t, ...
    'flag_dim','col', ...
    'init_method','[sd_t_meth]', ...
    'init_vals',[sd_t_vals], ...
    'init_crit',[sd_t_o], ...
    'flag_method','[sd_t_f_meth]', ...
    'flag_vals',[ [sd_t_f_vals] ], ...
    'flag_crit',[sd_t_f_o], ...
    'plot_figs','off');

% Updating rejection structure
logging_log('INFO', 'Updating latflaginfo structure...');
ica_epoch_flags=zeros(size(EEG.data(1,:,:)));
ica_epoch_flags(1,:,epoch_inds(flag_a_t_inds))=1;
EEG.marks=marks_add_label(EEG.marks,'time_info', ...
	{'ic_a',[1,1,0],ica_epoch_flags});
logging_log('INFO', 'TIME TO: CALCULATE ALPHA ICAACT POWER AND UPDATE MARKS...');
toc
clear ica_epoch_flags;
259
260
261

%% CALCULATE BETA ICAACT POWER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
262
% Calculates and flags times of high Beta Power. 
263
264
265
266
267
% Note that at the moment this mark is used for removals.

% Calculating the Beta component flag
logging_log('INFO', 'Calculating the beta IC power array for window criteria...');
tic;
268
[EEG,icaact_b_t]=chan_variance(EEG,'data_field','icaact', ...
269
270
271
272
273
274
275
		'epoch_inds',epoch_inds, ...
		'detrend','on', ...
        'varmeasure','spect', ...
        'spectrange',[14 30], ...
        'plot_figs','off');
logging_log('INFO', 'TIME TO: CALCULATE BETA ICAACT POWER...');

276
% Creating the Beta component flag
277
logging_log('INFO', 'Assessing window icaact sd distributions...')
278
[~,flag_b_t_inds]=marks_array2flags(icaact_b_t, ...
279
    'flag_dim','col', ...
280
    'init_method','[sd_t_meth]', ...
281
282
    'init_vals',[sd_t_vals], ...
    'init_crit',[sd_t_o], ...
283
284
285
    'flag_method','[sd_t_f_meth]', ...
    'flag_vals',[ [sd_t_f_vals] ], ...
    'flag_crit',[sd_t_f_o], ...
286
287
288
289
290
291
292
293
294
295
296
297
298
299
    'plot_figs','off');

logging_log('INFO', 'TIME TO: CREATE WINDOW CRITERIA VECTOR...');

% Updating rejection structure
logging_log('INFO', 'Updating latflaginfo structure...');
icb_epoch_flags=zeros(size(EEG.data(1,:,:)));
icb_epoch_flags(1,:,epoch_inds(flag_b_t_inds))=1;
EEG.marks=marks_add_label(EEG.marks,'time_info', ...
	{'ic_b',[0,1,1],icb_epoch_flags});
logging_log('INFO', 'TIME TO: CALCULATE BETA ICAACT POWER AND UPDATE MARKS...');
toc
clear icb_epoch_flags;

300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
%% CALCULATE LOW GAMMA ICAACT POWER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculates and flags times of high Low-Gamma Power. 
% Note that at the moment this mark is used for removals.

% Calculating the Beta component flag
logging_log('INFO', 'Calculating the low-gamma IC power array for window criteria...');
tic;
[EEG,icaact_lg_t]=chan_variance(EEG,'data_field','icaact', ...
		'epoch_inds',epoch_inds, ...
		'detrend','on', ...
        'varmeasure','spect', ...
        'spectrange',[31 60], ...
        'plot_figs','off');
logging_log('INFO', 'TIME TO: CALCULATE LOW-GAMMA ICAACT POWER...');

% Creating the LOW-GAMMA component flag
logging_log('INFO', 'Assessing window icaact sd distributions...')
[~,flag_lg_t_inds]=marks_array2flags(icaact_lg_t, ...
    'flag_dim','col', ...
    'init_method','[sd_t_meth]', ...
    'init_vals',[sd_t_vals], ...
    'init_crit',[sd_t_o], ...
    'flag_method','[sd_t_f_meth]', ...
    'flag_vals',[ [sd_t_f_vals] ], ...
    'flag_crit',[sd_t_f_o], ...
    'plot_figs','off');

logging_log('INFO', 'TIME TO: CREATE WINDOW CRITERIA VECTOR...');

% Updating rejection structure
logging_log('INFO', 'Updating latflaginfo structure...');
iclg_epoch_flags=zeros(size(EEG.data(1,:,:)));
iclg_epoch_flags(1,:,epoch_inds(flag_lg_t_inds))=1;
EEG.marks=marks_add_label(EEG.marks,'time_info', ...
	{'ic_lg',[0,1,1],iclg_epoch_flags});
logging_log('INFO', 'TIME TO: CALCULATE LOW-GAMMA ICAACT POWER AND UPDATE MARKS...');
toc
clear iclg_epoch_flags;


%% CALCULATE HIGH GAMMA ICAACT POWER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculates and flags times of high High-Gamma Power. 
% Note that at the moment this mark is used for removals.

% Calculating the Beta component flag
logging_log('INFO', 'Calculating the high-gamma IC power array for window criteria...');
tic;
[EEG,icaact_hg_t]=chan_variance(EEG,'data_field','icaact', ...
		'epoch_inds',epoch_inds, ...
		'detrend','on', ...
        'varmeasure','spect', ...
        'spectrange',[61 100], ...
        'plot_figs','off');
logging_log('INFO', 'TIME TO: CALCULATE HIGH-GAMMA ICAACT POWER...');

% Creating the HIGH-GAMMA component flag
logging_log('INFO', 'Assessing window icaact sd distributions...')
[~,flag_hg_t_inds]=marks_array2flags(icaact_hg_t, ...
    'flag_dim','col', ...
    'init_method','[sd_t_meth]', ...
    'init_vals',[sd_t_vals], ...
    'init_crit',[sd_t_o], ...
    'flag_method','[sd_t_f_meth]', ...
    'flag_vals',[ [sd_t_f_vals] ], ...
    'flag_crit',[sd_t_f_o], ...
    'plot_figs','off');

logging_log('INFO', 'TIME TO: CREATE WINDOW CRITERIA VECTOR...');

% Updating rejection structure
logging_log('INFO', 'Updating latflaginfo structure...');
ichg_epoch_flags=zeros(size(EEG.data(1,:,:)));
ichg_epoch_flags(1,:,epoch_inds(flag_hg_t_inds))=1;
EEG.marks=marks_add_label(EEG.marks,'time_info', ...
	{'ic_hg',[0,1,1],ichg_epoch_flags});
logging_log('INFO', 'TIME TO: CALCULATE HIGH-GAMMA ICAACT POWER AND UPDATE MARKS...');
toc
clear ichg_epoch_flags;
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406

%% MARKS UPDATE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Merges the epochs back to continuous data. Flags marks gaps and combines time 
% the new time marks into the manual marks.

% Concat epochs back into continuous data
logging_log('INFO', 'Concatenating windowed data...');
tic;
EEG=marks_epochs2continuous(EEG);
EEG=eeg_checkset(EEG,'eventconsistency');
logging_log('INFO', 'TIME TO: CONCATENATE THE WINDOWED DATA INTO CONTINUOUS DATA...');
toc

% Mark flag gaps between time marks
tic;
EEG=pop_marks_flag_gap(EEG,{'manual','mark_gap','ic_sd2'}, ...
                            [min_gap_ms],'mark_gap',[.8,.8,.8],'offsets',[0 0],'ref_point','both');
logging_log('INFO', 'TIME TO: MARK FLAG GAPS THAT ARE < 2 SECONDS...');
toc

% Combine the removal marks into manual
tic;
EEG=pop_marks_merge_labels(EEG,'time_info',{'manual','ic_sd2','mark_gap'},'target_label','manual');
logging_log('INFO', 'TIME TO: COMBINE MARKS STRUCTURE INTO MANUAL FLAGS...');
toc

407
408
409
410
411
% begin prep for dipole fit
logging_log('INFO','Checking for JVM...');
javachk('jvm');

% TODO(brad) this isn't added for some reason
412
addpath([pwd '/derivatives/lossless/code/dependencies/eeglab_asr_amica/plugins/Fieldtrip-lite141209/src'])
413
414
415
416
417
418
419
420

%% DIPFIT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Does the dipole fitting procedures. This takes a substantial amount of
% time.
% TODO: Make sure this is all needed, seems to be repeating things
logging_log('INFO','FITTING DIPOLES...');
tic;
421
EEG = pop_dipfit_settings( EEG, 'hdmfile','[dip_vol]', ... 
422
				'coordformat','MNI', ...
423
424
425
				'mrifile','[dip_mri]', ... 
				'chanfile','[dip_elc]', ...
				'coord_transform',[dip_transmat], ...
426
427
				'chansel',[1:EEG.nbchan] );

428
429
430
431
% OLD SINGLE THREADED SETUP
% EEG = pop_multifit(EEG, [1:EEG.nbchan] ,'threshold',100,'plotopt',{'normlen' 'on'});

% NEW MULTI-THREADED SETUP
432
EEG = pop_multifit_bucanl_parallel(EEG, [1:EEG.nbchan], 8,'threshold',100,'plotopt',{'normlen' 'on'});
433
434
435
436
437
438
439

EEG=eeg_checkset(EEG);
logging_log('INFO','TIME TO: FIT DIPOLES...');
toc

%% MARK ICs FOR REJECTION BASED ON RV
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
440
% Creates marks based on dipfit results. Dipfit marks are added
441
% for rejection while ICMARC are just for reference.
442
% IC RESIDUAL VARIANCE FLAG
443
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
444
445
446
447
% Residuals greater than 15% are flagged - can be done here as it only references
% dipoles and their values.
%logging_log('INFO','UPDATE MARKS STRUCTURE WITH IC_RV...');
%tic;
448
% Dipole criteria changed from 0.15 to 0.2 to be less sensitive
449
450
451
452
453
454
%EEG.reject.gcompreject(find([EEG.dipfit.model.rv]>.15))=1;
%EEG.reject.gcompreject(find(isnan([EEG.dipfit.model.rv])))=1;
%EEG.marks=marks_add_label(EEG.marks,'comp_info', ...
%	{'ic_rv',[1,0,0],[1,0,0],-1,EEG.reject.gcompreject'});
%logging_log('INFO','TIME TO: UPDATE MARKS STRUCTURE...');
%toc
455
456

% COMBINE MARKS STRUCTURE INTO MANUAL FLAGS...
457
458
459
460
%tic;
%EEG=pop_marks_merge_labels(EEG,'comp_info',{'manual','ic_rv'},'target_label','manual');
%logging_log('INFO','TIME TO: COMBINE MARKS STRUCTURE INTO MANUAL FLAGS...');
%toc
461

462
463
464
465
466
467
468
469
% Adjust montage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This will pretty much always be dataset specific - so be careful
EEG = warp_locs( EEG, '[dip_elc]', ...
    'transform',[cor_transmat], ...
    'manual','off');
EEG.chaninfo.nosedir='[cor_nosedir]';
EEG=eeg_checkset(EEG);
470

471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
%START TEMPORARY MOVE TO QC SCRIP
% %% Use the ICLabel plugin 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % ICLabel is a database style plugin that classifies ICA components into the
% % following seven categories
% 
% % (Brain)(Muscle)(Eye)(Heart)(Line Noise)(Channel Noise)(Other)
% % (  1  )(  2   )( 3 )(  4  )(    5     )(      6      )(  7  )
% 
% % Colour structure and groupings are as follows:
% % ----------
% % Brain         -> [0, 0.4, 0]
% % ----------
% % Muscle        -> [0.4, 0, 0]
% % Eye           -> [0.8, 0, 0]
% % Heart         -> [0.8, 0.4, 0]
% % ----------
% % Line Noise    -> [0.2, 0.4, 0.8]
% % Channel noise -> [0.4, 0, 0.8]
% % ----------
% % Other         -> [0.8, 0.8, 0]
% % ----------
% 
% % Create lookup table of the above information.
% pairs = {'brain',[0, 0.4, 0];
% 'muscle',[0.4,0,0];
% 'eye',[0.8,0,0];
% 'heart',[0.8,0.4,0];
% 'line_noise',[0.2,0.4,0.8];
% 'chan_noise',[0.4,0,0.8];
% 'other',[0.8,0.8,0]}
% lookupTable = cell2struct(pairs, {'name', 'colour'}, 2);
% 
% % Below code is how to call proper label lookup
% % [prob, classind] = max(EEG.etc.ic_classification.ICLabel.classifications(1, :));
% % Below line should give you number of components
% % length(EEG.etc.ic_classification.ICLabel.classifications);
% 
% % Run ICLabel
% EEG = iclabel(EEG);
% 
% % Extract labels from the structure.
% rawLabels = zeros([1,length(EEG.etc.ic_classification.ICLabel.classifications)]);
% for index = 1:length(EEG.etc.ic_classification.ICLabel.classifications)
%     [p,i] = max(EEG.etc.ic_classification.ICLabel.classifications(index, :));
%     rawLabels(index) = i;
% end
% 
% %ICLabel Exporting to marks
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% logging_log('INFO','UPDATE MARKS STRUCTURE WITH ICLABEL BRAIN COMPS...');
% tic;
% for compStage = 1:7
%      ICLabel_pass= zeros([1,length(rawLabels)]);
%      ICLabel_pass(find(rawLabels==compStage)) = 1;
%      ICLabel_pass = ICLabel_pass';
%      EEG.marks=marks_add_label(EEG.marks,'comp_info', ...
%           {lookupTable(compStage).name,lookupTable(compStage).colour, ...
%           lookupTable(compStage).colour,1,ICLabel_pass});
% end
% logging_log('INFO','TIME TO: UPDATE MARKS STRUCTURE...');
% toc
%END TEMPORARY MOVE TO QC SCRIP
534
535


536
%% SAVE LOSSLESS DERIVATIVE *_eeg_ll.set FILE (same name as the input file from the 
537
% BIDS root directory [but located in the derivatives/lossless directory]).
538
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
539
logging_log('INFO','Saving set file: [out_path]/[batch_dfn,_,-1]_eeg_ll.set...');
540
tic;
541
EEG = pop_saveset( EEG, 'filename','[out_path]/[batch_dfn,_,-1]_eeg_ll.set');
542
logging_log('INFO','TIME TO: SAVE DIP DATA SET...');
543
544
toc

545
546
% Save diagnostic arrays
try OCTAVE_VERSION;
547
    save('-mat7-binary','[out_path]/[batch_dfn,_,-1]_icaact_sd2_t.mat','icaact_sd2_t');
548
catch
549
    save('[out_path]/[batch_dfn,_,-1]_icaact_sd2_t.mat','icaact_sd2_t');
550
551
end

552
553
554
logging_log('INFO', 'Scheduler: [scheduler]');
print_chan_sample(EEG);