Making bad channel detection portable

@@ -172,12 +172,12 @@ uicontrol(Fmenu,'Style','Text','String',SPMc,...
 %-Data
 %-Data
 %-----------------------------------------------------------------------
 %-----------------------------------------------------------------------
 uicontrol(Fmenu,'Style','PopUp',...
 uicontrol(Fmenu,'Style','PopUp',...
-    'String','Import|  Convert|  Electrodes positions|  Bipolar montage|  Bad channels|  Events|  Set time origin',...
+    'String','Import|  Convert|  Electrodes positions|  Bipolar montage|  Set bad channels|  Detect bad channels|  Events|  Set time origin',...
 	'Position',[035 379 100 030].*WS,...
 	'Position',[035 379 100 030].*WS,...
     'FontSize', FS(9), ...
     'FontSize', FS(9), ...
 	'ToolTipString', 'Convert data in SPM EEG format',...
 	'ToolTipString', 'Convert data in SPM EEG format',...
     'CallBack','spm(''PopUpCB'',gcbo)',...
     'CallBack','spm(''PopUpCB'',gcbo)',...
-	'UserData',{'ImaGIN_spm_eeg_converteeg2mat;','ImaGIN_Electrode;','ImaGIN_BipolarMontage;','ImaGIN_BadChannelSet;','ImaGIN_Events;','ImaGIN_TimeZero;'},...
+	'UserData',{'ImaGIN_spm_eeg_converteeg2mat;','ImaGIN_Electrode;','ImaGIN_BipolarMontage;','ImaGIN_BadChannelSet;','ImaGIN_BadChannel;','ImaGIN_Events;','ImaGIN_TimeZero;'},...
     'Tag', 'EEG');
     'Tag', 'EEG');
 
 
 uicontrol(Fmenu,'Style','PopUp',...
 uicontrol(Fmenu,'Style','PopUp',...

 function D = ImaGIN_BadChannel(S)
 function D = ImaGIN_BadChannel(S)
-% Extract badchannels indices
+% Detect bad channels indices with a trained classifier.
+%
 % -=============================================================================
 % -=============================================================================
 % This function is part of the ImaGIN software: 
 % This function is part of the ImaGIN software: 
 % https://f-tract.eu/
 % https://f-tract.eu/
@@ -16,62 +17,96 @@ function D = ImaGIN_BadChannel(S)
 %
 %
 % Authors: Viateur Tuyisenge & Olivier David
 % Authors: Viateur Tuyisenge & Olivier David
 
 
-try
-    FileIn = S.dataset; % name of cropped iEEG file .dat/.mat
-    %badDir = S.DirFileOut; % directory where iEEG file .dat/.mat with badchannels indices will be stored
-    trainDir = S.trainBase;  % directory where training set is stored
-    try
-        D = spm_eeg_load(FileIn); % Load the cropped meeg object
-    catch
-        FileIn = spm_select(1, '\.mat$', 'Select data file');
-        D=spm_eeg_load(FileIn);
+% Name of input iEEG file .dat/.mat (cropped)
+if (nargin >= 1) && isfield(S, 'dataset') && ~isempty(S.dataset)
+    FileIn = S.dataset; 
+else
+    FileIn = spm_select(1, '\.mat$', 'Select EEG mat file');
+    if isempty(FileIn)
+        return;
     end
     end
+end
+
+% Directory where iEEG file .dat/.mat with badchannels indices will be stored
+if (nargin >= 1) && isfield(S, 'FileOut') && ~isempty(S.FileOut)
+    [badDir, FileOut] = fileparts(S.FileOut);
+    isNewFile = 1;
+else
+    [badDir, FileOut] = fileparts(FileIn);
+    isNewFile = 0;
+end
+
+% Directory where training set is stored (if not defined, use the toolbox folder)
+if (nargin >= 1) && isfield(S, 'trainBase') && ~isempty(S.trainBase)
+    trainDir = S.trainBase;
+else
+    trainDir = fileparts(mfilename('fullpath'));
+end
+
+
+% Detect bad channels and save summary images
+try
+    % Load the cropped meeg object
+    D = spm_eeg_load(FileIn);
     
     
-    [badDir,FileOut,~] = fileparts(S.FileOut);
-    bPrefix = strcat(badDir,'/',FileOut);
-    
+    % Compute signal features (apply artifact correction only for stims)
     clear S2;
     clear S2;
     S2.FileName = FileIn;
     S2.FileName = FileIn;
-
+    if ~isempty(D.events) && ~isempty(D.events.type) && ismember('Stim', {D.events.type})
+        S2.InterpolationFilter = 1;
+    else
+        S2.InterpolationFilter = 0;
+    end
+    T = ImaGIN_FeatureSEEG(S2);
     
     
-%    T = ImaGIN_FeatureSEEG(S2); % function returns a table T of features
+    % If the trained classifier is not available: compute it
+    trainedFile = fullfile(trainDir, 'ImaGIN_trainedClassifier.mat');
+    if ~exist(trainedFile, 'file')
+        % Get train base file
+        trainBaseFile = fullfile(trainDir, 'ImaGIN_trainBaseFeatures.mat');
+        if ~exist(trainBaseFile, 'file')
+            error(['File ' trainBaseFile ' was not found.']);
+        end
+        % Load the train base
+        trainBase = load(trainBaseFile);
+        % Train the classifier
+        trainedClassifier = ImaGIN_trainClassifier(trainBase.predictors, trainBase.response);
+    % Otherwise, load the trained classifier
+    else
+        trainedClassifier = load(trainedFile);
+    end
     
     
-%     Tbase = readtable(strcat(trainDir, '/trainBaseFeatures.csv')); % load training base
-%     trainingData = Tbase(:,2:9);
-%     [trainedClassifier, ~] = ImaGIN_trainClassifier(trainingData); % train the model
-    %Uncomment the following line to use the trained model
-  
-%     load(strcat(trainDir, '/ImaGIN_trainedClassifier.mat'))  %  load the trained classifier
-%     channelClass = trainedClassifier.predictFcn(T(:,2:8)); % predict new dataset
-%     bd = strcmp(channelClass','Bad);
-
-    bd = strcmp('channelClass','Bad');
+    % Predict new dataset
+    channelClass = trainedClassifier.predictFcn(T(:,2:8)); 
+    % Get list of detected bad channels
+    bIdx = find(strcmp(channelClass, 'Bad'));
     
     
-    bIdx = find(bd);
-    badFile = fopen(strcat(bPrefix,'_bIdx.txt'),'w'); % Save badchannel indices in .txt file
-    fprintf(badFile,'%d\n',bIdx(:));
+    % Save bad channel indices in .txt file
+    badFile = fopen(fullfile(badDir, [FileOut, '_bIdx.txt']), 'w'); 
+    fprintf(badFile, '%d\n', bIdx(:));
     fclose(badFile);
     fclose(badFile);
     
     
-%     Tnew = [T channelClass];
-%     Tnew.Properties.VariableNames{'Var9'} = 'Note';
-%     csvfilename = strcat(bPrefix,'.csv'); % Save feature table & badchannels indices
-%     writetable(Tnew,csvfilename,'Delimiter',',');
-    
+    % Add badchannel index in meeg object
     if ~isempty(bIdx)
     if ~isempty(bIdx)
-        D = badchannels(D,bIdx,1); %Add badchannel index in meeg object
+        D = badchannels(D,bIdx,1); 
+    end
+    % Create new file in output (otherwise, simply update the input file)
+    if isNewFile
+        Dbad = clone(D, fullfile(badDir, FileOut), [D.nchannels D.nsamples D.ntrials]); % save meeg with badchannel indices in Badchannel directory
+        Dbad(:,:,:) = D(:,:,:);
     end
     end
-    Dbad = clone(D, bPrefix, [D.nchannels D.nsamples D.ntrials]); % save meeg with badchannel indices in Badchannel directory
-    Dbad(:,:,:) = D(:,:,:);
+    % Save modified list of bad channels
     save(Dbad);
     save(Dbad);
-    %Channel plots and ScreenShots
-    figDir = strcat(badDir, '/ScreenShot');
+    
+    % Channel plots and ScreenShots
+    figDir = fullfile(badDir, 'ScreenShot');
     if ~exist(figDir, 'dir')
     if ~exist(figDir, 'dir')
         mkdir(figDir);
         mkdir(figDir);
     end
     end
     
     
     close all;
     close all;
     
     
-    Size = 8;  %Number of channels per screenshot
+    Size = 8;  % Number of channels per screenshot
     n_c  = size(D,1);
     n_c  = size(D,1);
     if n_c >= Size
     if n_c >= Size
         tmp = floor(n_c/Size);
         tmp = floor(n_c/Size);
@@ -148,8 +183,9 @@ try
         close
         close
     end
     end
     set_final_status('OK')
     set_final_status('OK')
+    
 catch ME
 catch ME
-    fprintf('%s \n',ME.message)
+    fprintf('%s \n', ME.message)
     set_final_status('NOK')
     set_final_status('NOK')
 end 
 end 
     
     

@@ -12,24 +12,29 @@ function T =  ImaGIN_FeatureSEEG(S)
 %
 %
 % Copyright (c) 2000-2018 Inserm U1216
 % Copyright (c) 2000-2018 Inserm U1216
 % =============================================================================-
 % =============================================================================-
-% 
+%
 % Authors: Viateur Tuyisenge & Olivier David
 % Authors: Viateur Tuyisenge & Olivier David
 
 
 sFile = S.FileName;
 sFile = S.FileName;
-[pth, fName, ~] = fileparts(sFile);
-clear S
-S.Fname = sFile;
-S.EventType = 'Stim';
-%S.StartInterpolation= -0.015;
-%S.EndInterpolation  = 0.015;
-S.StartInterpolation = -0.008;
-S.EndInterpolation   = 0.008;
-
-D = ImaGIN_InterpolationFilter(S);
-clear S
-
-iFile = [pth '/i' fName];
-nFile = [pth '/ni' fName];
+[pth, fName] = fileparts(sFile);
+
+% Apply interpolation filter to remove artifact (unless disabled)
+if ~isfield(S, 'InterpolationFilter') || isempty(S.InterpolationFilter) || S.InterpolationFilter
+    clear S
+    S.Fname = sFile;
+    S.EventType = 'Stim';
+    %S.StartInterpolation= -0.015;
+    %S.EndInterpolation  = 0.015;
+    S.StartInterpolation = -0.008;
+    S.EndInterpolation   = 0.008;
+    D = ImaGIN_InterpolationFilter(S);
+    strInterp = 'i';
+else
+    strInterp = '';
+end
+
+iFile = [pth '/' strInterp fName];
+nFile = [pth '/n' strInterp fName];
 P.Fname =  iFile;
 P.Fname =  iFile;
 P.Freq  = 50;
 P.Freq  = 50;
 ImaGIN_NotchFilter(P) % notch filter 50Hz
 ImaGIN_NotchFilter(P) % notch filter 50Hz
@@ -38,19 +43,21 @@ P.Fname =  nFile;
 P.Freq  = 100;
 P.Freq  = 100;
 ImaGIN_NotchFilter(P) % notch filter 100Hz
 ImaGIN_NotchFilter(P) % notch filter 100Hz
 clear P.Freq;
 clear P.Freq;
+if ~isempty(strInterp)
+    delete([iFile '.*']);
+end
 delete([nFile '.*']);
 delete([nFile '.*']);
-delete([iFile '.*']);
-nnFile  =  [pth '/nni' fName];
+nnFile  =  [pth '/nn' strInterp fName];
 P.Fname =  nnFile;
 P.Fname =  nnFile;
 P.Freq  = 150;
 P.Freq  = 150;
 ImaGIN_NotchFilter(P) % notch filter 150Hz
 ImaGIN_NotchFilter(P) % notch filter 150Hz
 clear P.Freq;
 clear P.Freq;
 delete([nnFile,'.*']);
 delete([nnFile,'.*']);
 
 
-nnnFile  =  [pth '/nnni' fName];
+nnnFile  =  [pth '/nnn' strInterp fName];
 P.LFname = nnnFile;
 P.LFname = nnnFile;
 ImaGIN_LowPassFilter(P) % lowpass filter 0.2Hz
 ImaGIN_LowPassFilter(P) % lowpass filter 0.2Hz
-lpf_nFile = [pth '/lpf_nnni' fName];
+lpf_nFile = [pth '/lpf_nnn' strInterp fName];
 delete([nnnFile,'.*'])
 delete([nnnFile,'.*'])
 
 
 D = spm_eeg_load(lpf_nFile);
 D = spm_eeg_load(lpf_nFile);