Merge branch 'nguyen-td:master' into BS_dev

README.md

@@ -36,7 +36,7 @@ Next, clone the repository
 ```
 git clone https://github.com/sccn/roiconnect.git
 ```
-That's it! If you want to run the plugin, please start [EEGLAB](https://github.com/sccn/eeglab#to-use-eeglab) first. You may need to add EEGLAB to the [MATLAB path](https://de.mathworks.com/help/matlab/ref/addpath.html).  
+That's it! If you want to run the plugin, please start [EEGLAB](https://github.com/sccn/eeglab#to-use-eeglab) first. You may need to add EEGLAB to the [MATLAB path](https://de.mathworks.com/help/matlab/ref/addpath.html). Some functions may require the additional installation of  [FieldTrip (lite or normal)](https://www.fieldtriptoolbox.org) and [Brainmovie](https://github.com/arnodelorme/brainmovie).
 
 Why is this EEGLAB plugin not in the EEGLAB plugin manager? The plugin is beta. Once it is completed, it will be added to the EEGLAB plugin manager. 
 

libs/nolte/data2bs_univar.m

@@ -21,7 +21,7 @@
 %   respective channel combination. For example, if para.chancomb = [1,2,1],
 %   the output will be cs([1,2,1], f1, f2) = <x(f1)_1*x(f2)_2*conj(x(f1+f2-1)_1)>
 %   (for a comparison, check the default output). In this example, at least two channels 
-%   are required.
+%   are required. Throws an error if there are not enough channels.
 %
 % output: 
 % cs: nchan  by nf by nf tensor for nf frequencies (i.e. nf=maxfreqbins)   

pop_roi_activity.m

@@ -31,8 +31,8 @@
 %                       specified, the signal is zero padded accordingly.
 %                       Default is 0 (means no padding).
 %  'chansel'          - [cell array of string] channel selection. Default is all.
-%  'lowmemory'        - ['on'|'off'] options to run the code with low memory, though, it might take significantly longer to complete.
-%
+%  'lowmemory'        - ['on'|'off'] Option to run the code with low memory, though, it might take significantly longer to complete. When turned on, the estimation of voxel-wise spectral power 
+%                       will require less memory.
 %
 % Other optional inputs:
 %  All ROI_ACTIVITY parameters are accepted as input and passed on.
@@ -207,7 +207,7 @@
         {} { 'style' 'text' 'string' 'Model parameters:' } { 'style' 'edit' 'string' '0.05'  'tag' 'modelparams' } {}  ...
         {} ...
         { 'style' 'text' 'string' 'Other parameters' 'fontweight' 'bold'} ...
-        {} { 'style' 'text' 'string' 'Number of dimentions per ROI:'  } { 'style' 'edit' 'string' '3' 'tag' 'pca' } {}  ...
+        {} { 'style' 'text' 'string' 'Number of dimensions per ROI:'  } { 'style' 'edit' 'string' '3' 'tag' 'pca' } {}  ...
         };
 
     [result,usrdat,~,out] = inputgui('geometry', uigeom, 'uilist', uilist, 'helpcom', 'pophelp(''pop_roi_activity'')', ...

pop_roi_connect.m

@@ -17,9 +17,9 @@
 %  'naccu'          - [integer]  Number of accumulation for stats. Default is 0.
 %  'methods'        - [cell] Cell of strings corresponding to methods.
 %                       'CS'    : Cross spectrum
-%                       'aCOH'  : Coherence
-%                       'cCOH'  : (complex-valued) coherency
-%                       'iCOH'  : absolute value of the imaginary part of coherency
+%                       'aCOH'  : Coherence 
+%                       'cCOH'  : (Complex-valued) Coherency
+%                       'iCOH'  : Absolute value of the imaginary part of Coherency
 %                       'GC'    : Granger Causality
 %                       'TRGC'  : Time-reversed Granger Causality
 %                       'wPLI'  : Weighted Phase Lag Index

pop_roi_connectplot.m

@@ -14,12 +14,14 @@
 %  'measure'              - ['psd'|'roipsd'|'trgc'|'crossspecimag'|'crossspecpow'|'mic'|'mim']
 %                           'psd'   : Source power spectrum
 %                           'psdroi': ROI based power spectrum
-%                           'trgc'  : Time-reversed granger causality
-%                           'gc'    : Granger causality
+%                           'TRGC'  : Time-reversed granger causality
+%                           'GC'    : Granger causality
 %                           'crossspecimag': Imaginary part of coherence from cross-spectrum
 %                           'crossspecpow' : Average cross-spectrum power for each ROI
-%                           'mic' : Maximized Imaginary Coherency for each ROI
-%                           'mim' : Multivariate Interaction Measure for each ROI
+%                           'aCOH': Coherence 
+%                           'iCOH': Absolute value of the imaginary part of Coherency
+%                           'MIC' : Maximized Imaginary Coherency for each ROI
+%                           'MIM' : Multivariate Interaction Measure for each ROI
 %                           'pac' : Phase-amplitude coupling for a certain frequency (band) combination based on bicoherence
 %                           'pac_anti': Phase-amplitude coupling for a certain frequency (band) combination based on the antisymmetrized bicoherence
 %  'freqrange'            - [min max] frequency range or [integer] single frequency in Hz. Default is to plot broadband power.
@@ -146,11 +148,33 @@
         splot(end  ).plot3d = plot3dFlag;
     end
 
-    if isfield(EEG.roi, 'COH')
+    if isfield(EEG.roi, 'aCOH')
         splot(end+1).label    = 'ROI to ROI coherence';
         splot(end  ).labelshort = 'Coherence';
-        splot(end  ).acronym  = 'Coh';
-        splot(end  ).unit   = 'Coh';
+        splot(end  ).acronym  = 'aCOH';
+        splot(end  ).unit   = 'aCOH';
+        splot(end  ).cortex = cortexFlag;
+        splot(end  ).matrix = 1;
+        splot(end  ).psd    = -1;
+        splot(end  ).plot3d = plot3dFlag;
+    end
+
+    if isfield(EEG.roi, 'cCOH')
+        splot(end+1).label    = 'ROI to ROI coherency';
+        splot(end  ).labelshort = 'Coherency';
+        splot(end  ).acronym  = 'cCoh';
+        splot(end  ).unit   = 'cCOH';
+        splot(end  ).cortex = cortexFlag;
+        splot(end  ).matrix = 1;
+        splot(end  ).psd    = -1;
+        splot(end  ).plot3d = plot3dFlag;
+        end
+
+    if isfield(EEG.roi, 'iCOH')
+        splot(end+1).label    = 'ROI to ROI absolute value of the imaginary part of coherency';
+        splot(end  ).labelshort = 'Img. part of Coherency';
+        splot(end  ).acronym  = 'iCOH';
+        splot(end  ).unit   = 'iCOH';
         splot(end  ).cortex = cortexFlag;
         splot(end  ).matrix = 1;
         splot(end  ).psd    = -1;
@@ -411,6 +435,17 @@
                 end
                 matrix = squeeze(mean(MI(frq_inds, :, :),1));
                 cortexPlot = mean(matrix, 2);
+
+            case { 'acoh' 'ccoh' 'icoh'}
+                if strcmpi(g.measure, 'aCOH')
+                    matrix = squeeze(mean(S.aCOH(frq_inds, :, :), 1));
+                elseif strcmpi(g.measure, 'cCOH')
+                    error(['Complex values are not supported. To plot the absolute values, compute "aCOH", ' ...
+                        'to plot the imaginary part, compute "iCOH".'])
+                else
+                    matrix = squeeze(mean(S.iCOH(frq_inds, :, :), 1));
+                end
+                cortexPlot = mean(matrix, 2);
 
             case { 'crossspecpow' 'coh' 'crossspecimag' }
                 if strcmpi(g.measure, 'coh')

pop_roi_statsplot.m

@@ -68,7 +68,7 @@
         matrix = squeeze(S.(g.measure)(frq_inds, :, :, :));
     end
 
-    % generate p-values by comparing the true FC (first shuffle) to null distribution
+    % average over one dimension to obtain net FC, then generate p-values by comparing the true FC (first shuffle) to null distribution
     netFC = squeeze(mean(matrix, 2));
     FC_pn = sum(netFC(:, 1) < netFC(:, 2:end), 2)./(size(matrix, 3) - 1);
 

roi_activity.m

@@ -41,7 +41,8 @@
 %  'freqresolution'   - [integer] Desired frequency resolution (in number of frequencies). If
 %                       specified, the signal is zero padded accordingly.
 %                       Default is 0 (means no padding).
-%  'lowmemory'        - ['on'|'off'] options to run the code with low memory, though, it might take significantly longer to complete.
+%  'lowmemory'        - ['on'|'off'] Option to run the code with low memory, though, it might take significantly longer to complete. When turned on, the estimation of voxel-wise spectral power 
+%                       will require less memory.
 %
 % Output:
 %  EEG - EEGLAB dataset with field 'roi' containing connectivity info.

test_pipes/pipeline_connectivity.m

@@ -49,6 +49,6 @@
     EEG = pop_roi_connect(EEG, 'methods', measures(iMeasure));
     t(iMeasure) = toc;
 end
-pop_roi_connectplot(EEG, 'measure', 'iCOH', 'plotcortex', 'on', 'plotmatrix', 'off');  
+pop_roi_connectplot(EEG, 'measure', 'MIM', 'plotcortex', 'on', 'plotmatrix', 'on');  
 pop_roi_connectplot(EEG, 'measure', 'MIM', 'plot3d', 'on',  'plotcortex', 'off', 'plot3dparams', {'thresholdper', 0.05, 'brainmovieopt' {'nodeColorDataRange', [], 'nodeSizeLimits', [0 0.2]}});
 

test_pipes/test_brainplots_loretaatlas.m

@@ -19,7 +19,7 @@
 EEG = pop_leadfield(EEG, 'sourcemodel',fullfile(eeglabp,'plugins','dipfit','LORETA-Talairach-BAs.mat'), ...
     'sourcemodel2mni',[0 -24 -45 0 0 -1.5708 1000 1000 1000] ,'downsample',1);
 
-EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3);
+EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3, 'chansel', EEG.dipfit.chansel);
 EEG = pop_roi_connect(EEG, 'methods', {'MIM'});
 
 %% Plot brain plot with different parameters

test_pipes/test_flex_epochlength.m

@@ -19,7 +19,7 @@
 EEG = pop_leadfield(EEG, 'sourcemodel',fullfile(eeglabp,'functions','supportfiles','head_modelColin27_5003_Standard-10-5-Cap339.mat'), ...
     'sourcemodel2mni',[0 -24 -45 0 0 -1.5708 1000 1000 1000] ,'downsample',1);
 
-EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',1);
+EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3, 'chansel', EEG.dipfit.chansel);
 
 measures = { 'MIM' 'wPLI' };
 for iMeasure = 1:length(measures)

test_pipes/test_fooof.m

@@ -20,7 +20,7 @@
     'sourcemodel2mni',[0 -24 -45 0 0 -1.5708 1000 1000 1000] ,'downsample',1);
 
 % compute fooof
-EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3, 'fooof', 'on', 'fooof_frange', [1 30]);
+EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3, 'fooof', 'on', 'fooof_frange', [1 30], 'chansel', EEG.dipfit.chansel);
 
 % access values
 results = EEG.roi.fooof_results;

test_pipes/test_largebarplot_dkatlas.m

@@ -19,7 +19,7 @@
 EEG = pop_leadfield(EEG, 'sourcemodel',fullfile(eeglabp,'functions','supportfiles','head_modelColin27_5003_Standard-10-5-Cap339.mat'), ...
     'sourcemodel2mni',[0 -24 -45 0 0 -1.5708 1000 1000 1000] ,'downsample',1);
 
-EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3);
+EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3, 'chansel', EEG.dipfit.chansel);
 tic
 EEG = pop_roi_connect(EEG, 'methods', {'MIM'});
 toc

test_pipes/test_matrixplots_brainregions.m

@@ -19,7 +19,7 @@
 EEG = pop_leadfield(EEG, 'sourcemodel',fullfile(eeglabp,'plugins','dipfit','LORETA-Talairach-BAs.mat'), ...
     'sourcemodel2mni',[0 -24 -45 0 0 -1.5708 1000 1000 1000] ,'downsample',1);
 
-EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','Brain-Regions','nPCA',3);
+EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','Brain-Regions','nPCA',3, 'chansel', EEG.dipfit.chansel);
 EEG = pop_roi_connect(EEG, 'methods', {'MIM'});
 
 %% Plot matrix with different parameters

test_pipes/test_matrixplots_dkatlas.m

@@ -19,7 +19,7 @@
 EEG = pop_leadfield(EEG, 'sourcemodel',fullfile(eeglabp,'functions','supportfiles','head_modelColin27_5003_Standard-10-5-Cap339.mat'), ...
     'sourcemodel2mni',[0 -24 -45 0 0 -1.5708 1000 1000 1000] ,'downsample',1);
 
-EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3);
+EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3, 'chansel', EEG.dipfit.chansel);
 EEG = pop_roi_connect(EEG, 'methods', {'MIM'});
 
 %% Plot matrix with different parameters

test_pipes/test_matrixplots_loretaatlas.m

@@ -19,7 +19,7 @@
 EEG = pop_leadfield(EEG, 'sourcemodel',fullfile(eeglabp,'plugins','dipfit','LORETA-Talairach-BAs.mat'), ...
     'sourcemodel2mni',[0 -24 -45 0 0 -1.5708 1000 1000 1000] ,'downsample',1);
 
-EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3);
+EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3, 'chansel', EEG.dipfit.chansel);
 EEG = pop_roi_connect(EEG, 'methods', {'MIM'});
 
 %% Plot matrix with different parameters

test_pipes/test_mim_stats.m

@@ -22,7 +22,15 @@
 EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3, 'chansel', EEG.dipfit.chansel);
 
 %% Create null distribution
-EEG1_MIM = pop_roi_connect(EEG, 'methods', {'CS' 'cCOH', 'wPLI' 'MIM'}, 'freqresolution', 200, 'roi_selection', {1, 3, 7}, 'conn_stats', 'on', 'nshuf', 1000, 'poolsize', 12); % takes very long!
-EEG2_MIM = pop_roi_connect(EEG, 'methods', {'CS' 'cCOH', 'wPLI' 'MIM'}, 'conn_stats', 'on', 'nshuf', 30); % takes very long!
+EEG1_MIM = pop_roi_connect(EEG, 'methods', {'CS' 'cCOH', 'wPLI' 'MIM'}, 'freqresolution', 200, 'roi_selection', {1, 3, 7}, 'conn_stats', 'on', 'nshuf', 1000, 'poolsize', 12); 
+EEG2_MIM = pop_roi_connect(EEG, 'methods', {'CS' 'cCOH', 'wPLI' 'MIM'}, 'conn_stats', 'on', 'nshuf', 1000); % takes very long!
 % load('test_pipes/MIM_shuf.mat')
-pop_roi_statsplot(EEG1_MIM, 'measure', 'MIM', 'freqrange', [8 13]);
+pop_roi_statsplot(EEG1_MIM, 'measure', 'MIM', 'freqrange', [8 13]);
+
+%% Unit test
+% Test if the first shuffle of the surrogate connectivity matrix in EEG1 is in fact the true matrix
+EEG1 = pop_roi_connect(EEG, 'methods', {'MIM'}, 'conn_stats', 'on', 'nshuf', 3);
+EEG2 = pop_roi_connect(EEG, 'methods', {'MIM'}, 'conn_stats', 'off');
+if ~isequal(squeeze(EEG1.roi.MIM(:, :, :, 1)), EEG2.roi.MIM) 
+    error 'The first shuffle in the surrogate connectivity array is not the true matrix.'
+end

test_pipes/test_roi_selection.m

@@ -21,7 +21,7 @@
 % EEG = pop_leadfield(EEG, 'sourcemodel',fullfile(eeglabp,'plugins','dipfit','LORETA-Talairach-BAs.mat'), ...
 %     'sourcemodel2mni',[0 -24 -45 0 0 -1.5708 1000 1000 1000] ,'downsample',1);
 
-EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3);
+EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3, 'chansel', EEG.dipfit.chansel);
 
 %% Test PAC for a selection of ROIs
 low = 10;

test_pipes/test_zeropadding.m

@@ -23,7 +23,7 @@
 
 %% Test zero padding in roi_activity and plot results
 % sanity check with default frequency resolution (freqresolution = 0)
-EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3);
+EEG = pop_roi_activity(EEG, 'leadfield',EEG.dipfit.sourcemodel,'model','LCMV','modelparams',{0.05},'atlas','LORETA-Talairach-BAs','nPCA',3, 'chansel', EEG.dipfit.chansel);
 
 % subplot(2,1,1)
 % for j = 1:size(EEG.roi.source_roi_power,2)