diff --git a/gui/mainGUI.m b/gui/mainGUI.m index c566f86..99fe3eb 100755 --- a/gui/mainGUI.m +++ b/gui/mainGUI.m @@ -929,6 +929,7 @@ function createbutton_Callback(hObject, eventdata, handles) handles.params.ChannelReductionParams.tobeExcludedChans = ... str2num(get(handles.excludeedit, 'String')); + handles.params.ChannelReductionParams.readdExcludedChans = true; % hard coded for now! else handles.params.ChannelReductionParams = struct([]); end diff --git a/preprocessing/RecommendedParameters.m b/preprocessing/RecommendedParameters.m index 030cd6d..e548d50 100755 --- a/preprocessing/RecommendedParameters.m +++ b/preprocessing/RecommendedParameters.m @@ -77,7 +77,8 @@ DetrendingParams = struct(); - ChannelReductionParams = struct('tobeExcludedChans', []); + ChannelReductionParams = struct('tobeExcludedChans', [], ... + 'readdExcludedChans', false); EEGSystem = struct('name', 'Others',... diff --git a/preprocessing/performETguidedICA.m b/preprocessing/performETguidedICA.m index 173e866..5614a0b 100644 --- a/preprocessing/performETguidedICA.m +++ b/preprocessing/performETguidedICA.m @@ -1,4 +1,4 @@ -function [EEG] = performETguidedICA(EEG, params) +function [EEG, ET] = performETguidedICA(EEG, params) % For details see the underlying publication: Dimigen, 2020, NeuroImage @@ -116,7 +116,12 @@ % % % EEG = pop_detecteyemovements(EEG,[LX LY],[RX RY],THRESH,MINDUR,DEG_PER_PIXEL,SMOOTH,0,25,2,PLOTFIG,WRITESAC,WRITEFIX); -%% Create optimized data for ICA training (OPTICAT, Dimigen, 2018) +%% keep ET data for reattaching it at the end of the preprocessing stream +nbchan = length(ETparams.sync.importColumns); +ET = EEG; +ET.data = ET.data(end-nbchan+1:end, :); % keep only ET data +ET.nbchan = nbchan; +ET.chanlocs = ET.chanlocs(end-nbchan+1:end); OW_PROPORTION = 1.0; % overweighting proportion SACCADE_WINDOW = [str2double(params.from_edit) str2double(params.to_edit)]; % time window to overweight (-20 to 10 ms is default) diff --git a/preprocessing/performICLabel.m b/preprocessing/performICLabel.m index 7c5d39d..2a617e2 100755 --- a/preprocessing/performICLabel.m +++ b/preprocessing/performICLabel.m @@ -1,4 +1,4 @@ -function EEG = performICLabel(EEG, varargin) +function [EEG, ET] = performICLabel(EEG, varargin) % performICLabel perform Independent Component Analysis (ICA) on the high % passsed EEG and classifies bad components using ICLabel. % This function applies a high pass filter before the ICA. But the output @@ -114,8 +114,9 @@ EEG_orig.automagic.iclabel.ETguidedICA.performed = 'no'; try if ETguidedICA - [~, EEG] = evalc('performETguidedICA(EEG, addETdataParams)'); - EEG.data = EEG.data(1:size(EEG_orig.data, 1), :); % remove ET data + [~, EEG, ET] = evalc('performETguidedICA(EEG, addETdataParams)'); + % remove ET data + EEG.data = EEG.data(1:size(EEG_orig.data, 1), :); EEG.nbchan = EEG_orig.nbchan; EEG.chanlocs = EEG_orig.chanlocs; EEG_orig.automagic.iclabel.ETguidedICA.performed = 'yes'; @@ -125,6 +126,10 @@ fprintf('\n ET guided ICA skipped. Continue with the standard ICA... \n') end +if ~exist('ET','var') + ET = []; +end + %% Run ICA [~, EEG, ~] = evalc('pop_runica(EEG, ''icatype'',''runica'',''chanind'',EEG.icachansind)'); diff --git a/preprocessing/preprocess.m b/preprocessing/preprocess.m index c933bc0..ca4f92a 100755 --- a/preprocessing/preprocess.m +++ b/preprocessing/preprocess.m @@ -129,7 +129,6 @@ ICLabelParams.addETdataParams = p.Results.addETdataParams; end - if isempty(Settings) Settings = Recs.Settings; end @@ -144,12 +143,20 @@ data_orig = data; % for plot with original data -% Set system dependent parameters and eeparate EEG from EOG -[EEG, EOG, EEGSystem, MARAParams] = ... +% Set system dependent parameters and separate EEG from EOG +[EEG, EOG, EXCLUDED, EEGSystem, MARAParams] = ... systemDependentParse(data, EEGSystem, ChannelReductionParams, ... MARAParams, ORIGINAL_FILE); EEGRef = EEG; +% free memory if EXCLUDED channel data is not supposed to be reattached +% later +if ~isempty(ChannelReductionParams) + if ~ChannelReductionParams.readdExcludedChans + clear EXCLUDED + end +end + % Trim data EEG.automagic.TrimData.performed = 'no'; if isfield(TrimDataParams, 'changeCheck') @@ -305,7 +312,7 @@ EEG.automagic.iclabel.performed = 'no'; if ( ~isempty(ICLabelParams) ) try - EEG = performICLabel(EEG, ICLabelParams); + [EEG, ET] = performICLabel(EEG, ICLabelParams); EEG.icachansind = find(~EEG.automagic.preprocessing.removedMask); catch ME message = ['ICLabel is not done on this subject, continue with the next steps: ' ... @@ -413,7 +420,7 @@ clear chan_nb re_chan; end -% Write back output +% Write back output (ref channel) if ~isempty(EEGSystem.refChan) removedChans(removedChans >= EEGSystem.refChan.idx)=removedChans(removedChans >= EEGSystem.refChan.idx)+1; EEG.automagic.autoBadChans = removedChans; @@ -421,6 +428,33 @@ else EEG.automagic.autoBadChans = removedChans; end + +% Put back excluded channels (only if wanted) +if ~isempty(EEG.automagic.channelReduction.params) + if (EEG.automagic.channelReduction.params.readdExcludedChans && ... + ~isempty(EEG.automagic.channelReduction.excludedChannels)) + excludedChans = sort(EEG.automagic.channelReduction.excludedChannels); % should be ordered (for step below) + for i_ch_excl = 1:length(excludedChans) + excluded_chan = excludedChans(i_ch_excl); + EEG.data = [EEG.data(1:excluded_chan-1,:); ... + EXCLUDED.data(i_ch_excl,:);... + EEG.data(excluded_chan:end,:)]; + EXCLUDED.chanlocs(i_ch_excl).maraLabel = []; + EEG.chanlocs = [EEG.chanlocs(1:excluded_chan-1), ... + EXCLUDED.chanlocs(i_ch_excl), ... + EEG.chanlocs(excluded_chan:end)]; + end + EEG.nbchan = size(EEG.data,1); + clear excluded_chan i_ch_excl + + % Write back output (excluded channels) + for excluded_chan = excludedChans + removedChans(removedChans >= excluded_chan)=removedChans(removedChans >= excluded_chan)+1; + end + end +end +EEG.automagic.autoBadChans = removedChans; + EEG.automagic.params = params; EEG.automagic = rmfield(EEG.automagic, 'preprocessing'); @@ -429,6 +463,26 @@ allSteps.EEGFinal = EEG; end +% attach Eyetracker data (if wanted) +addETdata = true; +EEG.automagic.channelReduction.params.addETdata = addETdata; % hard coded for now +%addETevents = true; % might be seperately set, if preferred +if exist('ET', 'var') + if ~isempty(ET) && EEG.automagic.channelReduction.params.addETdata + EEG.data = [EEG.data; ET.data]; % add to the very end, such that it doesnt affect any index + f = fieldnames(ET.chanlocs); + toKeep = fieldnames(EEG.chanlocs); + toRemove = f(~ismember(f,toKeep)); + EEG.chanlocs = [EEG.chanlocs, rmfield(ET.chanlocs, toRemove)]; + EEG.nbchan = size(EEG.data,1); +% end +% if ~isempty(ET) && addETevents + EEG.event = ET.event; + EEG.urevent = ET.urevent; + %EEG.eventdescription = ET.eventdescription; + end +end + %% Creating the final figure to save %%% Find the colormap selected @@ -491,6 +545,8 @@ else final_idx = 1:size(EEG.data, 1); end +% only keep original EEG channels (no ET, no excluded channels) +final_idx = final_idx(ismember(final_idx, EEG.automagic.channelReduction.usedEEGChannels)); %eeg figure subplot(13,1,2:3) @@ -726,6 +782,9 @@ else final_idx = 1:size(EEG.data, 1); end +% only keep original EEG channels (no ET, no excluded channels) +final_idx = final_idx(ismember(final_idx, EEG.automagic.channelReduction.usedEEGChannels)); + fig2 = figure('visible', 'off'); ax = gca; % outerpos = ax.OuterPosition; diff --git a/preprocessing/systemDependentParse.m b/preprocessing/systemDependentParse.m index 31ca2f0..123a172 100755 --- a/preprocessing/systemDependentParse.m +++ b/preprocessing/systemDependentParse.m @@ -1,4 +1,4 @@ -function [EEG, EOG, EEGSystem, MARAParams] = ... +function [EEG, EOG, EXCLUDED, EEGSystem, MARAParams] = ... systemDependentParse(EEG, EEGSystem, ... ChannelReductionParams, MARAParams, ORIGINAL_FILE) %#ok % systemDependentParse parse EEG input depending on the EEG system. @@ -449,9 +449,13 @@ clear all_chans; end +% Store excluded channel data for later retrieval +[~, EXCLUDED] = evalc('pop_select( EEG , ''channel'', tobeExcludedChans)'); + % Seperate EEG channels from EOG channels [~, EOG] = evalc('pop_select( EEG , ''channel'', eog_channels)'); [~, EEG] = evalc('pop_select( EEG , ''channel'', eeg_channels)'); + % Map original channel lists to new ones after the above separation if ~isempty(EEGSystem.refChan) [~, idx] = ismember(EEGSystem.refChan.idx, eeg_channels); @@ -463,6 +467,7 @@ if chanSize(2) == 1 EEG.chanlocs = EEG.chanlocs'; EOG.chanlocs = EOG.chanlocs'; + EXCLUDED.chanlocs = EXCLUDED.chanlocs'; end clear chanSize; diff --git a/src/Block.m b/src/Block.m index 764aa95..60f1a05 100755 --- a/src/Block.m +++ b/src/Block.m @@ -528,8 +528,17 @@ function excludeChannelsFromRBC(self, exclude_chans) if(isempty(EEG)) return; end + % index EEG channels for quality asssessment + eeg_chans = 1:EEG.nbchan; + if isfield(EEG.automagic, 'channelReduction') + if EEG.automagic.channelReduction.params.readdExcludedChans % in this case take original eeg indices. ET is also excluded by default + eeg_chans = EEG.automagic.channelReduction.usedEEGChannels; + elseif EEG.automagic.channelReduction.params.addETdata % in this case, original eeg indices dont work, but we have to excluce the ET channels at the end + eeg_chans = 1:EEG.nbchan-6; % hard coded for now + end + end rating_badchans = []; - qScore = calcQuality(EEG, rating_badchans, ... + qScore = calcQuality(EEG, eeg_chans, rating_badchans, ... self.project.qualityThresholds); qScoreIdx.OHA = arrayfun(@(x) ceil(length(x.OHA)/2), qScore); qScoreIdx.THV = arrayfun(@(x) ceil(length(x.THV)/2), qScore); @@ -610,29 +619,43 @@ function interpolate(self) InterpolationParams = self.params.InterpolationParams; end - %save old ica data which gets corrupted in eeg_interp method: - orig_icasphere=EEG.icasphere; - orig_icachansind=EEG.icachansind; - orig_icaweights=EEG.icaweights; - orig_icawinv=EEG.icawinv; + % to avoid excluded misc and ET channels going into the interpolation, + % we also let the misc channels be interpolated (temporarily). + % draw info from automagic params + EEGtoInterp = EEG; + misc_chans = []; + if isfield(automagic, 'channelReduction') + usedEEGChannels = automagic.channelReduction.usedEEGChannels; + %misc_chans = [1:EEG.nbchan]; + misc_chans = find(~cellfun('isempty', { EEGtoInterp.chanlocs.theta })); % only those with loc info, as others are ignored by interpolation function anyway.. + misc_chans = misc_chans(~ismember(misc_chans, usedEEGChannels)); + end + + % do interpolation if size(EEG.data,1)==length(interpolate_chans) disp('All channels are bad. Skipping interpolation...'); filenamE = strsplit(EEG.comments,filesep); filenamE = filenamE{end}; automagic.error_msg = ['Interpolation skipped because all channels are bad: ',filenamE]; else - EEG = eeg_interp(EEG ,interpolate_chans , InterpolationParams.method); + EEGtoInterp = eeg_interp(EEGtoInterp ,[interpolate_chans misc_chans], InterpolationParams.method); end - %put the original icadata back into the structure - EEG.icasphere=orig_icasphere; - EEG.icachansind=orig_icachansind; - EEG.icaweights= orig_icaweights; - EEG.icawinv=orig_icawinv; + % now add only the interpolated EEG channels back to EEG struct + EEG.data(interpolate_chans,:) = EEGtoInterp.data(interpolate_chans,:); rating_badchans = unique([self.finalBadChans interpolate_chans]); rating_badchans = setdiff(rating_badchans, ... self.project.manuallyExcludedRBCChans); - qScore = calcQuality(EEG, rating_badchans, ... + % index EEG channels for quality asssessment + eeg_chans = 1:EEG.nbchan; + if isfield(automagic, 'channelReduction') + if automagic.channelReduction.params.readdExcludedChans % in this case take original eeg indices. ET is also excluded by default + eeg_chans = automagic.channelReduction.usedEEGChannels; + elseif automagic.channelReduction.params.addETdata % in this case, original eeg indices dont work, but we have to excluce the ET channels at the end + eeg_chans = [1:EEG.nbchan-6]; % hard coded for now + end + end + qScore = calcQuality(EEG, eeg_chans, rating_badchans, ... self.project.qualityThresholds); qScoreIdx.OHA = arrayfun(@(x) ceil(length(x.OHA)/2), qScore); qScoreIdx.THV = arrayfun(@(x) ceil(length(x.THV)/2), qScore); diff --git a/src/calcQuality.m b/src/calcQuality.m index 9a7c46a..c6e24bf 100755 --- a/src/calcQuality.m +++ b/src/calcQuality.m @@ -1,4 +1,4 @@ -function Q = calcQuality(EEG, bad_chans, varargin) +function Q = calcQuality(EEG, eeg_chans, bad_chans, varargin) % Calculates quality measures of a dataset based on the following metrics: % % - The ratio of data points that exceed the absolute value a certain @@ -61,8 +61,13 @@ disp('No bad channel information...') end %% Data preparation -% Data +% only rate data that was actually prepocessed (without excluded and +% readded channels, or added ET channels) + +% remove readded EXCLUDED and ET channels from data (have to remove it from EEG, as calcCHV uses it) +EEG.data = EEG.data(eeg_chans,:); % i think its not necessary to shift around bad_chans indices, as they are not used specifically?! X = EEG.data; + % Get dimensions of data t = size(X,2); c = size(X,1);