diff --git a/toolbox/forward/bst_headmodeler.m b/toolbox/forward/bst_headmodeler.m index 6af5e4dca..9af1a6298 100644 --- a/toolbox/forward/bst_headmodeler.m +++ b/toolbox/forward/bst_headmodeler.m @@ -22,8 +22,13 @@ % .EEGMethod: Method used to compute the forward model for EEG sensors. % - 'eeg_3sphereberg' : EEG forward modeling with a set of 3 concentric spheres (Scalp, Skull, Brain/CSF) % - 'openmeeg' : OpenMEEG forward model -% .SEEGMethod: 'openmeeg' and 'duneuro' -% .ECOGMethod: 'openmeeg' and 'duneuro' +% .SEEGMethod: Method used to compute the forward model for SEEG sensors. +% - 'homogeneous' : Infinite homogeneous medium +% - 'openmeeg' : OpenMEEG forward model +% - 'duneuro' : DUNEuro forward model +% .ECOGMethod: Method used to compute the forward model for ECOG sensors. +% - 'openmeeg' : OpenMEEG forward model +% - 'duneuro' : DUNEuro forward model % .NIRSMethod: 'import' % % ======= METHODS OPTIONS ============================================= @@ -536,6 +541,13 @@ Gain(~isnan(Gain_dn)) = Gain_dn(~isnan(Gain_dn)); end +%% ===== COMPUTE: SEEG INFINITE HOMOGENEOUS MEDIUM HEADMODEL ===== +if strcmp('homogeneous', OPTIONS.SEEGMethod) + Gain(iSeeg, :) = bst_seeg_homogeneous(OPTIONS.GridLoc, OPTIONS.Channel(iSeeg), sSurfInner, OPTIONS); + strHistory = [strHistory, ' | ', 'Homogeneous medium (SEEG)', ... + ' | ', sprintf('Conductivity: %1.3f S/m', OPTIONS.Conductivity), ... + ' | ', sprintf('Min dist. betweem SEEG contact and dipoles: %1.3f mm', OPTIONS.MinSeegDipoleDist)]; +end %% ===== COMPUTE: BRAINSTORM HEADMODELS ===== Param = []; diff --git a/toolbox/forward/bst_seeg_homogeneous.m b/toolbox/forward/bst_seeg_homogeneous.m new file mode 100644 index 000000000..c56ecf56f --- /dev/null +++ b/toolbox/forward/bst_seeg_homogeneous.m @@ -0,0 +1,113 @@ +function G = bst_seeg_homogeneous(GridLoc, sChannel, sInnerSkull, Options) +% bst_seeg_homogeneous: Calculate the electric potential for infinite homogeneous medium +% +% USAGE: G = bst_seeg_homogeneous(GridLoc, sChannel, sInnerSkull, Options) +% +% INPUT: +% - GridLoc : Dipole locations (in meters) [nDipoles x 3] +% - sChannel : Channel structure [nSensors] +% - sInnerSkull : Inner skull surface structure +% - Options structure +% - Options.Conductivity : Conductivity (S/m) +% - Options.MinSeegDipoleDist : Minimum distance between SEEG and dipoles +% OUTPUTS: +% - G : SEEG forward model gain matrix [nSensors x (3*nDipoles)] +% +% DESCRIPTION: sEEG single layer forward model +% This function computes the voltage potential forward gain matrix for an array of +% sEEG electrodes inside the brain. The conductivity is assumed to be uniform and isotropoic +% inside the medium (that is assumed to be infinite). +% +% For electrodes outside of the brain, the grain is set to 0. +% +% Ref: +% + Grova, C., Aiguabella, M., Zelmann, R., Lina, J.-M., Hall, J.A. and Kobayashi, E. (2016), +% Intracranial EEG potentials estimated from MEG sources: A new approach to correlate MEG and iEEG data in epilepsy. +% Hum. Brain Mapp., 37: 1661-1683. https://doi.org/10.1002/hbm.23127 +% +% + Næss, S., Halnes, G., Hagen, E., Hagler Jr, D. J., Dale, A. M., Einevoll, G. T., & Ness, T. V. (2021). +% Biophysically detailed forward modeling of the neural origin of EEG and MEG signals. NeuroImage, 225, 117467. +% +% dot(n_i, u_ij) +% V(E_j) = -------------------------- +% 4 * pi * sigma0 * (r_ij)^2 +% +% V(E_j) = Electric potential at sensor j +% n_i = Vector, current dipole for source i +% u_ij = Unit vector, oriented from source i to sensor j +% sigma0 = Conductivity of infinite homogeneous medium +% r_ij = Euclidean distance between source i and sensor j +% +% Written as matrix multiplication: V = G * N +% V = Electric potential at contacts [nSensors, nTime] +% G = Gain matrix [nSensors, nDipoles] +% N = Dipole activation currents [nDipoles, nTime] +% + +% @============================================================================= +% This function is part of the Brainstorm software: +% https://neuroimage.usc.edu/brainstorm +% +% Copyright (c) University of Southern California & McGill University +% This software is distributed under the terms of the GNU General Public License +% as published by the Free Software Foundation. Further details on the GPLv3 +% license can be found at http://www.gnu.org/copyleft/gpl.html. +% +% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE +% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY +% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF +% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY +% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. +% +% For more information type "brainstorm license" at command prompt. +% =============================================================================@ +% +% Authors: Edouard Delaire, 2026 + + % Add default options + Options = struct_copy_fields(struct('Conductivity', 0.25, 'MinSeegDipoleDist', 3/1000), Options, 1); + min_distance = Options.MinSeegDipoleDist; + sigma0 = Options.Conductivity; + + NbElectrodes = length(sChannel); + NbVertices = size(GridLoc, 1); + + % Find electrodes that are inside the inner skull + SEEG_Loc = [sChannel.Loc]'; + isSEEGInsideSkull = inpolyhd(SEEG_Loc, sInnerSkull.Vertices, sInnerSkull.Faces); + + % Compute the leadfield + bst_progress('start', 'Computing head model', sprintf('Computing head model for %d contacts...', NbElectrodes), 0, NbElectrodes); + G = zeros(NbElectrodes , 3*NbVertices); + for iContact = 1:NbElectrodes + % Ignore contacts outside of the inner skull + if ~isSEEGInsideSkull(iContact) + continue + end + + % Compute unit vectors from SEEG contact to source points (u_j) + VectorDipolesToSEEG = SEEG_Loc(iContact, :) - GridLoc; + DistanceToDipoles = sqrt(sum(VectorDipolesToSEEG.^2,2)); + VectorDipolesToSEEG = VectorDipolesToSEEG ./ repmat(DistanceToDipoles, 1, 3); + + % Replace short distance (< 'min_distance') with 'min_distance' to prevent instabilities in the model by too close measurements + iShort = find(DistanceToDipoles < min_distance); + if ~isempty(iShort) + % fprintf(' %d vertex had distance to the cortex smaller than %.2f mm to electrodes %s \n', length(iShort), min_distance*1000, sChannel(iContact).Name); + DistanceToDipoles(iShort) = min_distance; + end + + % Compute the leadfield (u_j / (r_j)^2) + scaledVector = VectorDipolesToSEEG ./ repmat(DistanceToDipoles.^2, 1, 3); + + % Organize the matrix as x,y,z + G(iContact, :) = reshape(scaledVector', 1, []); + + bst_progress('inc', 1); + end + + % Add normalization constant + G = G / (4 * pi * sigma0); + + bst_progress('stop'); +end diff --git a/toolbox/forward/panel_headmodel.m b/toolbox/forward/panel_headmodel.m index 21575109f..d70ad294c 100644 --- a/toolbox/forward/panel_headmodel.m +++ b/toolbox/forward/panel_headmodel.m @@ -121,6 +121,7 @@ jCheckMethodSEEG.setSelected(1); % Combobox jComboMethodSEEG = gui_component('ComboBox', jPanelMethod, 'tab hfill', [], [], [], @UpdateComment, []); + jComboMethodSEEG.addItem(BstListItem('homogeneous', '', 'Homogeneous medium', [])); jComboMethodSEEG.addItem(BstListItem('openmeeg', '', 'OpenMEEG BEM', [])); jComboMethodSEEG.addItem(BstListItem('duneuro', '', 'DUNEuro FEM', [])); jComboMethodSEEG.setSelectedIndex(0); @@ -437,6 +438,7 @@ function UpdateComment(varargin) sMethod.Comment = strrep(sMethod.Comment, 'os_meg', 'Overlapping spheres'); sMethod.Comment = strrep(sMethod.Comment, 'meg_sphere', 'Single sphere'); sMethod.Comment = strrep(sMethod.Comment, 'eeg_3sphereberg', '3-shell sphere'); + sMethod.Comment = strrep(sMethod.Comment, 'homogeneous', 'Homogeneous medium'); sMethod.Comment = strrep(sMethod.Comment, 'openmeeg', 'OpenMEEG BEM'); sMethod.Comment = strrep(sMethod.Comment, 'duneuro', 'DUNEuro FEM'); % Grid type @@ -446,6 +448,7 @@ function UpdateComment(varargin) sMethod.Comment = [sMethod.Comment ' (mixed)']; end end + isHomogeneous = any(strcmpi(allMethods, 'homogeneous')); isOpenMEEG = any(strcmpi(allMethods, 'openmeeg')); isDuneuro = any(strcmpi(allMethods, 'duneuro')); % Get protocol description @@ -597,6 +600,38 @@ function UpdateComment(varargin) bst_error('Please import a cortex surface for this subject.', 'Head model', 0); return end + + % ===== INFINITE HOMOGENEOUS MEDIUM ===== + if isHomogeneous + % Interactive interface to set the options + if OPTIONS.Interactive + OPTIONS.Conductivity = []; + OPTIONS.MinSeegDipoleDist = []; + % Get options + optionsStr = {'Brain conductivity (S/m):','Minimum distance between SEEG and dipoles (mm):'}; + res = java_dialog('input', optionsStr, 'SEEG head model options', [], {'0.25','3'}); + if isempty(res) + return + end + resDouble = str2double(res); + % Validate inputs + for iRes = 1 : length(resDouble) + if isnan(resDouble(iRes)) + errMessage = sprintf('Invalid input "%s" for option:\n%s', res{iRes}, optionsStr{iRes}); + return + end + end + OPTIONS.Conductivity = resDouble(1); + OPTIONS.MinSeegDipoleDist = resDouble(2)/1000; + end + % Check that inner skull exists + if ~isempty(sSubject.iInnerSkull) + OPTIONS.InnerSkullFile = file_fullpath(sSubject.Surface(sSubject.iInnerSkull(1)).FileName); + else + errMessage = 'No inner skull surface for this subject.'; + return + end + end % ===== OPENMEEG ===== if isOpenMEEG