diff --git a/toolbox/forward/bst_duneuro2026.m b/toolbox/forward/bst_duneuro2026.m new file mode 100644 index 0000000000..f823f1b92a --- /dev/null +++ b/toolbox/forward/bst_duneuro2026.m @@ -0,0 +1,617 @@ +function [Gain, errMsg] = bst_duneuro2026(cfg) +% BST_DUNEURO: Call Duneuro to compute a FEM solution for Brainstorm. +% +% USAGE: [Gain, errMsg] = bst_duneuro(cfg) +% NOTE: online viewer of hdf5 files : https://myhdf5.hdfgroup.org/ + +% @============================================================================= +% 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: Takfarinas Medani, 2019-2026 +% Juan Garcia-Prieto, 2019-2021 +% Francois Tadel 2020-2023 + +% Initialize returned values +Gain = []; +%% ===== INSTALL AND GET THE EXECUTABLE ===== +% ***** ToDo: Define how to check the container ******* +% Install/load duneuro container +% Install/load duneuro plugin TODO: +% Note: For the current version of this function, we assume that the appropriate container is installed. +% Three possible runners are possible: +% - via docker +% - via podman +% - via apptainer +runner = cfg.containerRunner; +% Get DUNEuro container executable +% NOTE: temporary place for the containers +% https://github.com/users/MalteHoel/packages/container/package/duneuro_in_docker_testing +bst_plugin('SetProgressLogo', 'duneuro'); + +%% ===== SENSORS ===== +% Select modality for DUNEuro +isEeg = strcmpi(cfg.EEGMethod, 'duneuro2026') && ~isempty(cfg.iEeg); +isMeg = strcmpi(cfg.MEGMethod, 'duneuro2026') && ~isempty(cfg.iMeg); +isEcog = strcmpi(cfg.ECOGMethod, 'duneuro2026') && ~isempty(cfg.iEcog); +isSeeg = strcmpi(cfg.SEEGMethod, 'duneuro2026') && ~isempty(cfg.iSeeg); + +% Get the modality +if ((isEeg || isEcog || isSeeg) && isMeg) + dnModality = 'meeg'; +elseif (isEeg || isEcog || isSeeg) + dnModality = 'eeg'; +elseif isMeg + dnModality = 'meg'; % from DUNEuro side, EEG, sEEG, ECOG uses the same process +else + errMsg = 'No valid modality available.'; + return; +end + +% Get EEG positions +% Combined modalities for EEG/sEEG/EcoG as EEG +if (isEeg || isEcog || isSeeg) + cfg.iEeg = [cfg.iEeg, cfg.iSeeg, cfg.iEcog]; + EegLoc = cat(2, cfg.Channel(cfg.iEeg).Loc); +end + +% Get MEG positions/orientations +if isMeg + MegChannels = []; + for iChan = 1:length(cfg.iMeg) + sChan = cfg.Channel(cfg.iMeg(iChan)); + for iInteg = 1:size(sChan.Loc, 2) + MegChannels = [MegChannels; iChan, sChan.Loc(:,iInteg)', sChan.Orient(:,iInteg)', sChan.Weight]; + end + end + % In the case where the MEG integration points are used + if cfg.UseIntegrationPoint == 0 + % loop over the integration Points + % chan_loc = figure_3d('GetChannelPositions', cfg, cfg.iMeg); % <= this function is not sufficient, we need also the weights. + MegChannelsTemp = []; + for iChan = 1 : length(cfg.iMeg) + group = MegChannels(MegChannels(:,1) == iChan,:); + groupPositive = group(group(:,end)>0,:); + groupNegative = group(group(:,end)<0,:); + if ~isempty(groupPositive) + %equivalentPositionPostive = sum(repmat(abs(groupPositive(:,end)),[1 3]) .* groupPositive(:,2:4)); + equivalentPositionPostive = mean(groupPositive(:,2:4)); + MegChannelsTemp = [MegChannelsTemp; iChan equivalentPositionPostive groupPositive(1,5:7) sum(groupPositive(:,end))]; + end + if ~isempty(groupNegative) + %equivalentPositionNegative = sum(repmat(abs(groupNegative(:,end)),[1 3]) .* groupNegative(:,2:4)); + equivalentPositionNegative = mean(groupNegative(:,2:4)); + MegChannelsTemp = [MegChannelsTemp; iChan equivalentPositionNegative groupNegative(1,5:7) sum(groupNegative(:,end))]; + end + end + MegChannels = MegChannelsTemp; + end +end + +%% ===== HEAD MODEL ===== +% Load FEM mesh +FemMat = load(file_fullpath(cfg.FemFile)); +% Remove unselected tissue from the head model (MEG only), this could be also done for sEEG later +if strcmp(dnModality, 'meg') + % Remove the elements corresponding to the unselected tissues + iRemove = find(~ismember(FemMat.Tissue, find(cfg.FemSelect))); + if ~isempty(iRemove) + FemMat = fem_remove_elem(FemMat, iRemove); + end +elseif strcmp(dnModality,'meeg') && (sum(cfg.FemSelect) ~= length(unique(FemMat.Tissue))) + errMsg = 'Reduced head model cannot be used when computing MEG+EEG simultaneously.'; + return; +end + +% Get temp folder +TmpDir = bst_get('BrainstormTmpDir', 0, 'duneuro'); + +% Display message +bst_progress('text', 'DUNEuro: Writing temporary files...'); +disp(['DUNEURO> Writing temporary files to: ' TmpDir]); + +%% ====== SOURCE SPACE ===== +% Source space type +switch (cfg.HeadModelType) + case 'volume' + % Nothing to do as for now: + % TODO or keep it as it's now.... + case 'surface' + bst_progress('text', 'DUNEuro: Fixing source space...'); + % Read cortex file + sCortex = bst_memory('LoadSurface', cfg.CortexFile); + cfg.GridLoc = sCortex.Vertices; + % Shrink the cortex surface by XX mm + if (cfg.SrcShrink > 0) + % Get spherical coordinates of the surface normals + [azimuth, elevation] = cart2sph(sCortex.VertNormals(:,1), sCortex.VertNormals(:,2), sCortex.VertNormals(:,3)); + % Find components to shrink the surface in the three dimensions + depth = cfg.SrcShrink ./ 1000 .* [cos(elevation) .* cos(azimuth), cos(elevation) .* sin(azimuth), sin(elevation)]; + % Apply to the cortex surface + cfg.GridLoc = cfg.GridLoc - depth; + end + % Force all the dipoles within the GM layer + iGM = find(panel_duneuro('CheckType', FemMat.TissueLabels, 'gray'), 1); + iWM = find(panel_duneuro('CheckType', FemMat.TissueLabels, 'white'), 1); + if cfg.SrcForceInGM && ~isempty(iGM) + % Install/load iso2mesh plugin + [isInstalled, errMsg] = bst_plugin('Install', 'iso2mesh', cfg.Interactive); + if ~isInstalled + % Delete the temporary files + file_delete(TmpDir, 1, 1); + return; + end + % Extract GM vertices and elements + [gmVert, gmElem] = removeisolatednode(FemMat.Vertices, FemMat.Elements(FemMat.Tissue == iGM,:)); + % Compute the centroid of the GM elements + nElem = size(gmElem, 1); + nMesh = size(gmElem, 2); + ElemCenter = zeros(nElem, 3); + for i = 1:3 + ElemCenter(:,i) = sum(reshape(gmVert(gmElem,i), nElem, nMesh)')' / nMesh; + end + + % Extract GM envelope + envFaces = volface(FemMat.Elements(FemMat.Tissue <= iGM,:)); + [envVert, envFaces] = removeisolatednode(FemMat.Vertices, envFaces); + % Find the dipoles outside of the GM envelope + iVertOut = find(~inpolyhedron(envFaces, envVert, cfg.GridLoc)); + if ~isempty(iVertOut) + disp(['DUNEURO> Warning: ' num2str(length(iVertOut)) ' dipole(s) outside of the GM.']); + end + + % If there is a white matter layer: find the dipoles inside the WM and move them outside to the GM + if ~isempty(iWM) + % Extract GM envelope + wmFaces = volface(FemMat.Elements(FemMat.Tissue == iWM,:)); + [wmVert, wmFaces] = removeisolatednode(FemMat.Vertices, wmFaces); + % Find the dipoles outside of the GM envelope + iVertWM = find(inpolyhedron(wmFaces, wmVert, cfg.GridLoc)); + if ~isempty(iVertWM) + disp(['DUNEURO> Warning: ' num2str(length(iVertWM)) ' dipole(s) inside the WM.']); + iVertOut = union(iVertOut, iVertWM); + end + end + + % Move each vertex towards the centroid of the closest GM element + % view_surface_matrix(cfg.GridLoc, sCortex.Faces) + for i = 1:length(iVertOut) + bst_progress('text', sprintf('DUNEuro: Fixing dipole %d/%d...', i, length(iVertOut))); + % Find the closest GM centroid + iTarget = dsearchn(ElemCenter, cfg.GridLoc(iVertOut(i),:)); + targetFaces = volface(gmElem(iTarget,:)); + + % OPTION #1: Replace the vertex position directly with the centroid. + % => Problem: might project multiple vertices on the same centroid... + % cfg.GridLoc(iVertOut ,:) = ElemCenter(iTarget,:); + + % OPTION #2: Gradually move the vertex towards the center of the centroid, until it is located inside the element + % Move the vertex towards the center until it is inside the element + %nFix = 10; + %for iFix = 1:nFix + % tmpVert = (nFix - iFix)/nFix * cfg.GridLoc(iVertOut(i),:) + iFix/nFix * ElemCenter(iTarget,:); + % if inpolyhedron(targetFaces, gmVert, tmpVert) + % distMove = sqrt(sum((cfg.GridLoc(iVertOut(i),:) - tmpVert) .^ 2)) * 1000; + % disp(sprintf('DUNEURO> Dipole #%d moved inside the GM (%1.2fmm)', iVertOut(i), distMove)); + % cfg.GridLoc(iVertOut(i),:) = tmpVert; + % break; + % end + %end + + % OPTION #3: move the vertex towards the centroid of the element, and then place the final dipole + % in the symetric point to the center, as the image of the computed vertex + % x-----o-----x' + % ^ ^ ^____ : x' the image of x, or the final dipole position + % | |_________ : o is the center of the elem, and middle of [x,x'] + % |_____________ : the point inside the element determined by the tmpVert in the following equation + nFix = 10; % divid into 10 segments + iFix = 7; % ratio of the distance tmpVert from the centroid + tmpVert = (nFix - iFix)/nFix * cfg.GridLoc(iVertOut(i),:) + iFix/nFix * ElemCenter(iTarget,:); + distPoint = ElemCenter(iTarget,:) - tmpVert; + newPoint = ElemCenter(iTarget,:) + (distPoint); + distMove = sqrt(sum((cfg.GridLoc(iVertOut(i),:) - newPoint) .^ 2)) * 1000; + % use the image/symeric point if it's inside GM + if inpolyhedron(targetFaces, gmVert, newPoint) + disp(sprintf('DUNEURO> iDipole %d/%d : Dipole #%d moved inside the GM (%1.2fmm) (option 3 :as image)',i,length(iVertOut), iVertOut(i), distMove)); + cfg.GridLoc(iVertOut(i),:) = newPoint; + elseif inpolyhedron(targetFaces, gmVert, tmpVert) % use the original point if it's inside GM + disp(sprintf('DUNEURO> iDipole %d/%d : Warning Dipole #%d moved outside the GM (%1.2fmm) (option 3 :as image)', i,length(iVertOut),iVertOut(i), distMove)); + % use the original distance unstead of the image + distMove = sqrt(sum((cfg.GridLoc(iVertOut(i),:) - tmpVert) .^ 2)) * 1000; + disp(sprintf('DUNEURO> iDipole %d/%d : Correction 1: Dipole #%d moved inside the GM (%1.2fmm) (option 3: not image)', i,length(iVertOut),iVertOut(i), distMove)); + cfg.GridLoc(iVertOut(i),:) = tmpVert; + else % Use the option 2 defined by Francois + disp(sprintf('DUNEURO> iDipole %d/%d : Warning Dipole #%d moved outside the GM (%1.2fmm) (option 3 :as image)', i,length(iVertOut),iVertOut(i), distMove)); + nFix = 20; % with 10 it's not working for some extrem case, then I upgrade it to 20 + for iFix = 1: nFix + tmpVert = (nFix - iFix)/nFix * cfg.GridLoc(iVertOut(i),:) + iFix/nFix * ElemCenter(iTarget,:); + if inpolyhedron(targetFaces, gmVert, tmpVert) + distMove = sqrt(sum((cfg.GridLoc(iVertOut(i),:) - tmpVert) .^ 2)) * 1000; + disp(sprintf('DUNEURO> iDipole %d/%d : Correction 2: Dipole #%d moved inside the GM (%1.2fmm) (option2)', i,length(iVertOut),iVertOut(i), distMove)); + cfg.GridLoc(iVertOut(i),:) = tmpVert; + break; + end + end + end + end + % view_surface_matrix(cfg.GridLoc, sCortex.Faces) + end + case 'mixed' + % TODO : not used ? +end + +%% ===== SOURCE MODEL ===== +bst_progress('text', 'DUNEuro: Writing temporary files...'); +bstdn_write_source_space(TmpDir, cfg.GridLoc); + +%% ===== SENSOR MODEL ===== +% Write the EEG electrode file +if isEeg || isEcog || isSeeg + bstdn_write_pem_electrodes(TmpDir, EegLoc', 'measurement'); +end +% Write the MEG sensors data +if isMeg + % Write coil data + % coil_to_channel_transform = eye(length(MegChannels)); + coil_to_channel_transform = MegChannels(:,8:end); + dnbst_write_magnetometers(TmpDir, MegChannels(:,2:4), MegChannels(:,5:7), coil_to_channel_transform); + build_coil_to_channel_transform_matrix(MegChannels); +end + +%% ===== CONDUCTIVITY MODEL ===== +% Isotropic without tensor +if ~cfg.UseTensor + % Create tensor per tissue format from the isoconductivity + tensors = zeros(length(cfg.FemCond), 3, 3); + for iTissue = 1 : length(FemMat.TissueLabels) + tensors(iTissue,1,1) = cfg.FemCond(iTissue); + tensors(iTissue,2,2) = cfg.FemCond(iTissue); + tensors(iTissue,3,3) = cfg.FemCond(iTissue); + end +% NOTE: online viewer of hdf5 files : https://myhdf5.hdfgroup.org/ +else % With tensor (isotropic or anisotropic) + % Transformation matrix and tensor mapping on each direction + tensors = zeros(length(FemMat.Elements), 3,3); + for iElem = 1 : length(FemMat.Elements) + temp0 = reshape(FemMat.Tensors(iElem,:),3,[]); + T1 = temp0(:,1:3); % get the 3 eigen vectors + l = diag(temp0(:,4)); % get the eigen value as 3x3 + temp = T1 * l * T1'; % reconstruct the tensors + tensors(iElem,:,:) = temp; + end + % write the tensors + % TODO: Need to double check if it is correct with Malte. + % I use all the size of the elem == size of the tensors + % it seems that we can squeeze it down to set only one value for the + % iso tissues ==> smaller I/O file +end +dnbst_write_volume_conductor(TmpDir, FemMat, tensors); + +%% ===== TRANSFER MATRIX CONFIGURATION ===== + transfer_matrix_config = []; +transfer_matrix_config.name = 'compute_transfer_matrix'; +if strcmp(dnModality, 'eeg') + transfer_matrix_config.do_eeg = 'True'; + transfer_matrix_config.do_meg = 'False'; +end +if strcmp(dnModality, 'meg') + transfer_matrix_config.do_meg = 'True'; + transfer_matrix_config.do_eeg = 'False'; +end +if strcmp(dnModality, 'meeg') + transfer_matrix_config.do_meg = 'True'; + transfer_matrix_config.do_eeg = 'True'; +end +transfer_matrix_config.residual_reduction = cfg.residual_reduction; +transfer_matrix_config.nr_threads = num2str(cfg.NbOfThread); + +%% ===== LEADFIELD MATRIX CONFIGURATION ===== +leadfield_config = []; +leadfield_config.name = 'compute_leadfield'; +if strcmp(dnModality, 'eeg') + leadfield_config.do_meg = 'False'; + leadfield_config.do_eeg = 'True'; +end +if strcmp(dnModality, 'meg') + leadfield_config.do_meg = 'True'; + leadfield_config.do_eeg = 'False'; +end +if strcmp(dnModality, 'meeg') + leadfield_config.do_meg = 'True'; + leadfield_config.do_eeg = 'True'; +end +% set final hard code value +leadfield_config.eeg_scaling = cfg.eeg_scaling; % check with Malte if those value are optimised +leadfield_config.meg_scaling = cfg.meg_scaling; % Malte to check and get final value for MKSA system. +leadfield_config.sourcemodel = cfg.SrcModel2026; % [select from the interface: 'multipolar_venant', 'local_subtraction', 'partial_integration'] +leadfield_config.nr_threads = num2str(cfg.NbOfThread); % same as above +% ToDo from the interface and add it to the duneuro def options +if length(cfg.GridLoc) >= 100000 + % use the chunk dipoles / make it as inputs from the def optons + UI + leadfield_config.do_sanity_check = 'True'; + leadfield_config.chunked_computation = 'True'; + leadfield_config.chunk_size = '100000'; +end + +%% ===== RUN DUNEURO ====== +bst_progress('text', 'DUNEuro: Computing leadfield...'); +% disp(['DUNEURO> System call: ' callStr]); +tic; +% Call DUNEuro and Compute transfer matrix +[status, errMsg ] = bst_run_duneuro_task(TmpDir, transfer_matrix_config, runner); +transferMatrix_time = toc; + +tic; +% Call DUNEuro and Compute the Leadfield +[status, errMsg ] = bst_run_duneuro_task(TmpDir, leadfield_config, runner); +leadField_time = toc; + +if (status ~= 0) + errMsg = 'Error during the DUNEuro computation, see logs in the command window.'; + return; +end +disp(['DUNEURO> FEM computation completed in: ' num2str(toc) 's']); + +%% ===== READ LEADFIELD ====== +bst_progress('text', 'DUNEuro: Reading leadfield...'); +% For checking: read the used source space +% source_positions = transpose(h5read(fullfile(TmpDir, 'duneuro_io.hdf5'), '/measurement/source_space/positions')); +% EEG +if (isEeg || isEcog || isSeeg) + GainEeg = transpose(h5read(fullfile(TmpDir, 'duneuro_io.hdf5'), '/measurement/source_space/leadfield/EEG')); +end +if (isMeg) + GainMeg = transpose(h5read(fullfile(TmpDir, 'duneuro_io.hdf5'), '/measurement/source_space/leadfield/MEG')); +end + +% Fill the unused channels with NaN +Gain = NaN * zeros(length(cfg.Channel), 3 * length(cfg.GridLoc)); +if (isEeg || isEcog || isSeeg) + Gain(cfg.iEeg,:) = GainEeg; + % ToDo: add ref electrode / identify the reference from the channel + % file and post process the final Leadfield +end + +if isMeg % The MEG is not fully tested as we need to convert from integration points to final position + Gain(cfg.iMeg,:) = GainMeg; +end + +% Remove logo +bst_plugin('SetProgressLogo', []); +end + +%% ================================================================================= +% === SUPPORT FUNCTIONS ========================================================= +% ================================================================================= +function [iOk, errMsg] = dnbst_write_volume_conductor(duneuro_io_dir, FemMat, tensors) + iOk = 0; errMsg = ''; + nodes = FemMat.Vertices; + elements = FemMat.Elements - 1; + labels = FemMat.Tissue - 1; + io_file_path = fullfile(duneuro_io_dir, 'duneuro_io.hdf5'); + + % get information about the input data + nr_nodes = size(nodes, 1); + nr_elements = size(elements, 1); + nr_labels = size(labels, 1); + nr_unique_tensors = size(tensors,1); + + dim = 3; + nr_nodes_per_tetrahedron = 4; + if nr_elements ~= nr_labels + errMsg = 'Number of elements does not match number of labels'; return; + end + if size(nodes, 2) ~= dim + errMsg = 'Number of columns of nodes array must be 3'; return; + end + if size(elements, 2) ~= nr_nodes_per_tetrahedron + errMsg = 'The DUNEuro container interface currently only supports tetrahedral meshes (or the number of column of the elements array is wrong)'; return; + end + + if (size(tensors, 2) ~= dim) || (size(tensors, 3) ~= dim) + errMsg = 'The shape of the tensors array must be (K, 3, 3)'; return; + end + % create output file + % DUNEuro is a C++ toolbox, it utilizes row-major ordering. Matlab, on the other hand, + % uses colum-major ordering. For a consistent array shape, we thus need to transpose + % the arrays (resp. apply a permutation for the tensor array) + h5create(io_file_path, "/volume_conductor/nodes", [dim nr_nodes ], Datatype="double"); + h5create(io_file_path, "/volume_conductor/elements", [nr_nodes_per_tetrahedron nr_elements], Datatype="int32"); + h5create(io_file_path, "/volume_conductor/labels", nr_elements, Datatype="int32"); + h5create(io_file_path, "/volume_conductor/tensors", [dim dim nr_unique_tensors], Datatype="double"); + + h5write(io_file_path, "/volume_conductor/nodes", nodes'); + h5write(io_file_path, "/volume_conductor/elements", elements'); + h5write(io_file_path, "/volume_conductor/labels", labels'); + h5write(io_file_path, "/volume_conductor/tensors", permute(tensors, [3 2 1])); + + h5writeatt(io_file_path, "/volume_conductor", 'type', 'fitted'); + h5writeatt(io_file_path, "/volume_conductor", 'element_type', 'tetrahedron'); + + % if all goes well, return 1 + iOk = 1; +end + +function [iOk, errMsg] = dnbst_write_magnetometers(duneuro_io_dir, coil_positions, coil_orientations, coil_to_channel_transform) + iOk = 0; errMsg = ''; + + io_file_path = fullfile(duneuro_io_dir, 'duneuro_io.hdf5'); + + nr_magnetometers = size(coil_positions, 1); + nr_channels = size(coil_to_channel_transform, 1); + + dim = 3; + + if size(coil_positions, 2) ~= dim + errMsg = 'Number of columns of coil position array must be 3'; return; + end + + if size(coil_orientations, 2) ~= dim + errMsg = 'Number of columns of coil orientation array must be 3'; return; + end + + if size(coil_orientations, 1) ~= nr_magnetometers + error('Position and orientation arrays must have matching number of rows'); return; + end + + % if size(coil_to_channel_transform, 2) ~= nr_magnetometers + % errMsg = 'Number of columns of transformation matrix must match number of coils'; return; + % end + + h5create(io_file_path, "/measurement/sensors/magnetometers/positions", [dim nr_magnetometers], Datatype="double"); + h5create(io_file_path, "/measurement/sensors/magnetometers/orientations", [dim nr_magnetometers], Datatype="double"); + % h5create(io_file_path, "/measurement/sensors/magnetometers/coil_to_channel_transform", [nr_magnetometers, nr_channels], Datatype="double"); + + h5write(io_file_path, "/measurement/sensors/magnetometers/positions", coil_positions'); + h5write(io_file_path, "/measurement/sensors/magnetometers/orientations", coil_orientations'); + % h5write(io_file_path, "/measurement/sensors/magnetometers/coil_to_channel_transform", coil_to_channel_transform); + + % if all goes well, return 1 + iOk = 1; +end + +function MEGChannelTransform = build_coil_to_channel_transform_matrix(MegChannels) + nb_magnetic_field_values = size(MegChannels, 1); + + % Extract channel indices and weights + chan_idx = MegChannels(:,1); + w = MegChannels(:,end); + + % Find unique channels + channels = unique(chan_idx); + nb_chan = length(channels); + + % Preallocate channel transform + MEGChannelTransform = zeros(nb_chan, nb_magnetic_field_values); + + % Build transform channel by channel + for k = 1:nb_chan + % Get indices of projected magnetic field values that contribute + % to the current channel + rows_k = find(chan_idx == channels(k)); + nr_fields_in_channel = numel(rows_k); + + for l=1:nr_fields_in_channel + MEGChannelTransform(k, rows_k(l)) = w(rows_k(l)); + end + end + % not sure about this + h5create(io_file_path, "/measurement/sensors/magnetometers/coil_to_channel_transform", [nb_chan, nb_magnetic_field_values], Datatype="double"); + h5write(io_file_path, "/measurement/sensors/magnetometers/coil_to_channel_transform", MEGChannelTransform); +end + +function [iOk, errMsg] = bstdn_write_pem_electrodes(duneuro_io_dir, electrode_positions, electrode_type_flag) + iOk = 0; errMsg = ''; + + io_file_path = fullfile(duneuro_io_dir, 'duneuro_io.hdf5'); + + nr_electrodes = size(electrode_positions, 1); + dim = 3; + + if size(electrode_positions, 2) ~= dim + errMsg = 'The electrode array must have 3 columns'; return; + end + + if strcmp(electrode_type_flag, 'measurement') + h5create(io_file_path, "/measurement/sensors/electrodes/positions", [dim nr_electrodes], Datatype="double"); + h5write(io_file_path, "/measurement/sensors/electrodes/positions", electrode_positions'); + h5writeatt(io_file_path, "/measurement/sensors/electrodes", 'mode', 'PEM'); + elseif strcmp(electrode_type_flag, 'stimulation') + h5create(io_file_path, "/stimulation/TDCS/electrodes/positions", [dim nr_electrodes], Datatype="double"); + h5write(io_file_path, "/stimulation/TDCS/electrodes/positions", electrode_positions'); + h5writeatt(io_file_path, "/stimulation/TDCS/electrodes", 'mode', 'PEM'); + else + errMsg = 'Electrode type flag needs to be either "measurement" or "stimulation"'; return; + end + iOk = 1; +end + +function [iOk, errMsg] = bstdn_write_source_space(duneuro_io_dir, dipole_positions) + iOk = 0; errMsg = ''; + + io_file_path = fullfile(duneuro_io_dir, 'duneuro_io.hdf5'); + + nr_dipoles = size(dipole_positions, 1); + % no_dipoles = size(dipole_orientations, 1); + dim = 3; + + if size(dipole_positions, 2) ~= dim + errMsg = 'The dipole position array must have 3 columns'; return; + end + + h5create(io_file_path, "/measurement/source_space/positions", [dim nr_dipoles], Datatype="double"); + h5write(io_file_path, "/measurement/source_space/positions", dipole_positions'); + + iOk = 1; +end + + +function [status, errMsg ] = bst_run_duneuro_task(duneuro_io_dir, config, runner) + status = 0; errMsg = ''; + task_name = config.name; + + config_file_path = fullfile(duneuro_io_dir, 'config.ini'); + io_file_path = fullfile(duneuro_io_dir, 'duneuro_io.hdf5'); + + % first write config file + file_handle = fopen(config_file_path, 'wt'); + + fprintf(file_handle, ['[task_info]\ntask_list=' task_name '\n\n[' task_name '_config]\n']); + + config_keys = fieldnames(config); + for i = 1:length(config_keys) + current_key = config_keys{i}; + current_value = config.(current_key); + + % write everything except the task name + if ~strcmp(current_key, 'name') + fprintf(file_handle, [current_key '=' current_value '\n']); + end + end + + fclose(file_handle); + + % now execute system call to start the container + if strcmp(runner, 'docker') + + runner_system_call = ['docker run -t --rm -v ' duneuro_io_dir ':/duneuro/external_mount ghcr.io/maltehoel/duneuro_in_docker_testing:wip']; + elseif strcmp(runner, 'podman') + % if we run podman in rootless mode, we need to make the IO directory writable + % for the container user + duneuro_container_uid = '50000'; % this uid is explicitely set in the Dockerfile + change_ownership_in_container_command = ['podman unshare chown ' duneuro_container_uid ':' duneuro_container_uid ' -R ' duneuro_io_dir]; + container_command = ['podman run -t --rm -v ' duneuro_io_dir ':/duneuro/external_mount ghcr.io/maltehoel/duneuro_in_docker_testing:wip']; + change_ownership_back_command = ['podman unshare chown 0:0 -R ' duneuro_io_dir]; + runner_system_call = [change_ownership_in_container_command ' && ' container_command ' && ' change_ownership_back_command]; + elseif strcmp(runner, 'apptainer') + % runner_system_call = ['apptainer run --bind ' duneuro_io_dir ':/duneuro/external_mount docker://ghcr.io/maltehoel/duneuro_in_docker_testing:wip']; + %runner_system_call = ['apptainer run --bind ' duneuro_io_dir ':/duneuro/external_mount duneuro_testing.sif']; + + %system('bash -lc "module load apptainer && which apptainer"') + runner_system_call = ... + ['bash -lc "module load apptainer && apptainer run --bind ' ... + duneuro_io_dir ':/duneuro/external_mount duneuro_testing.sif"']; + +% status = system(runner_system_call); + + else + errMsg = 'unknown runner'; + return; + end + + status = system(runner_system_call); +end \ No newline at end of file diff --git a/toolbox/forward/bst_headmodeler.m b/toolbox/forward/bst_headmodeler.m index 004556c98d..cd54bde308 100644 --- a/toolbox/forward/bst_headmodeler.m +++ b/toolbox/forward/bst_headmodeler.m @@ -504,12 +504,17 @@ %% ===== COMPUTE: DUNEURO ===== -if ismember('duneuro', {OPTIONS.MEGMethod, OPTIONS.EEGMethod, OPTIONS.ECOGMethod, OPTIONS.SEEGMethod}) +if ismember('duneuro', {OPTIONS.MEGMethod, OPTIONS.EEGMethod, OPTIONS.ECOGMethod, OPTIONS.SEEGMethod}) ... + || ismember('duneuro2026', {OPTIONS.MEGMethod, OPTIONS.EEGMethod, OPTIONS.ECOGMethod, OPTIONS.SEEGMethod}) % Start progress bar bst_progress('start', 'Head modeler', 'Starting Duneuro...'); bst_progress('setimage', 'plugins/duneuro_logo.png'); % Run duneuro FEM computation - [Gain_dn, errMessage] = bst_duneuro(OPTIONS); + if ismember('duneuro', {OPTIONS.MEGMethod, OPTIONS.EEGMethod, OPTIONS.ECOGMethod, OPTIONS.SEEGMethod}) + [Gain_dn, errMessage] = bst_duneuro(OPTIONS); + elseif ismember('duneuro2026', {OPTIONS.MEGMethod, OPTIONS.EEGMethod, OPTIONS.ECOGMethod, OPTIONS.SEEGMethod}) + [Gain_dn, errMessage] = bst_duneuro2026(OPTIONS); + end % Comment in history field dn_types = {'MEG', 'EEG', 'ECOG', 'SEEG'}; Lia = ismember({OPTIONS.MEGMethod, OPTIONS.EEGMethod, OPTIONS.ECOGMethod, OPTIONS.SEEGMethod}, 'duneuro'); @@ -539,8 +544,8 @@ %% ===== COMPUTE: BRAINSTORM HEADMODELS ===== Param = []; -if (~isempty(OPTIONS.MEGMethod) && ~ismember(OPTIONS.MEGMethod, {'openmeeg', 'duneuro'})) || ... - (~isempty(OPTIONS.EEGMethod) && ~ismember(OPTIONS.EEGMethod, {'openmeeg', 'duneuro'})) +if (~isempty(OPTIONS.MEGMethod) && ~ismember(OPTIONS.MEGMethod, {'openmeeg', 'duneuro', 'duneuro2026'})) || ... + (~isempty(OPTIONS.EEGMethod) && ~ismember(OPTIONS.EEGMethod, {'openmeeg', 'duneuro', 'duneuro2026'})) % ===== DEFINE SPHERES FOR EACH SENSOR ===== Param = repmat(struct('Center', [], 'Radii', []), 1, length(OPTIONS.Channel)); diff --git a/toolbox/forward/duneuro_defaults.m b/toolbox/forward/duneuro_defaults.m index 7359d4239d..c2f3596460 100644 --- a/toolbox/forward/duneuro_defaults.m +++ b/toolbox/forward/duneuro_defaults.m @@ -75,11 +75,22 @@ cfgDef.BstMegTransferFile = 'meg_transfer.dat'; cfgDef.BstEegLfFile = 'eeg_lf.dat'; cfgDef.BstMegLfFile = 'meg_lf.dat'; - % [MEG computation Options] cfgDef.UseIntegrationPoint = 1; cfgDef.EnableCacheMemory = 0; cfgDef.MegPerBlockOfSensor = 0; % ToDo +% DN2026 +cfgDef.containerRunner = 'docker'; % ToDo +[cfgDef.NbOfThread, cfgDef.NbOfThreadMax] = GetThreads(); +cfgDef.SrcModel2026 = 'multipolar_venant'; % 'multipolar_venant', 'local_subtraction', 'partial_integration' +cfgDef.residual_reduction = '1e-16'; % for the transfer matrix +%From Malte: +% 1) a derivation of the EEG and MEG scaling factors to get the output in volt and tesla. The short version is the following: +% If you measure distances in meter, conductivites in Siemens/meter, and dipole moments in Amperemeter, then the following setting +% results in the duneuro EEG leadfied being given in volt and the duneuro MEG leadfield being given in tesla. +cfgDef.eeg_scaling = '1e0'; % +cfgDef.meg_scaling = '1e-7'; % + % Use default values if not set if (nargin == 0) || isempty(cfg) cfg = cfgDef; @@ -87,9 +98,6 @@ end % Add missing values cfg = struct_copy_fields(cfg, cfgDef, 0); - - - % % The reste is not needed... we keep it just in case % % subpart [analytic_solution] % cfg.minifile.solution.analytic_solution.radii = [1 2 3 4 ]; @@ -134,7 +142,14 @@ % fprintf(fid, 'reference = %s\n',cfg.minifile.wrapper.outputtreecompare.reference); % fprintf(fid, 'type = %s\n',cfg.minifile.wrapper.outputtreecompare.type); % fprintf(fid, 'absolute = %s\n',cfg.minifile.wrapper.outputtreecompare.absolute); +end - - - +%% ===== THREAD DETECTION FUNCTION ===== +function [nThreadOptimal, nPhysical] = GetThreads() + try + nPhysical = feature('numcores'); + nThreadOptimal = max(1, round(nPhysical/2)); + catch + nThreadOptimal = 1; % Safe fallback + end +end \ No newline at end of file diff --git a/toolbox/forward/panel_duneuro2026.m b/toolbox/forward/panel_duneuro2026.m new file mode 100644 index 0000000000..2c18f6bc77 --- /dev/null +++ b/toolbox/forward/panel_duneuro2026.m @@ -0,0 +1,429 @@ +function varargout = panel_duneuro2026(varargin) +% PANEL_DUNEURO2026: DUNEuro2026 options +% +% USAGE: bstPanel = panel_duneuro2026('CreatePanel', OPTIONS) : Call from the interactive interface +% bstPanel = panel_duneuro2026('CreatePanel', sProcess, sFiles) : Call from the process editor +% s = panel_duneuro2026('GetPanelContents') + +% @============================================================================= +% 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: Francois Tadel, 2020 +% Takfarinas Medani, 2026 : adapted from panel_duneuro + +eval(macro_method); +end + + +%% ===== CREATE PANEL ===== +function [bstPanelNew, panelName] = CreatePanel(sProcess, sFiles) %#ok + panelName = 'DuneuroOptions2026'; + % Java initializations + import java.awt.*; + import javax.swing.*; + + % GUI CALL: panel_duneuro('CreatePanel', OPTIONS) + if (nargin == 1) + OPTIONS = sProcess; + % Check if there is only MEG, for simplified model by default + isMegOnly = ~ismember('duneuro2026', {OPTIONS.EEGMethod, OPTIONS.ECOGMethod, OPTIONS.SEEGMethod}); + isMeg = isequal(OPTIONS.MEGMethod, 'duneuro2026'); + % PROCESS CALL: panel_duneuro('CreatePanel', sProcess, sFiles) + else + OPTIONS = sProcess.options.duneuro.Value; + % List of sensors + Modalities = {'MEG MAG', 'MEG GRAD', 'MEG', 'EEG', 'ECOG', 'SEEG'}; + if ~isempty(sFiles(1).ChannelTypes) + Modalities = intersect(sFiles(1).ChannelTypes, Modalities); + if any(ismember({'MEG MAG','MEG GRAD'}, Modalities)) + Modalities = setdiff(Modalities, 'MEG'); + end + end + isMeg = any(ismember({'MEG', 'MEG MAG', 'MEG GRAD'}, Modalities)); + isMegOnly = all(ismember(Modalities, {'MEG', 'MEG MAG', 'MEG GRAD'})); + % Get FEM files + sSubject = bst_get('Subject', sFiles(1).SubjectFile); + if isempty(sSubject.iFEM) + error('No available FEM mesh file for this subject.'); + end + OPTIONS.FemFile = sSubject.Surface(sSubject.iFEM).FileName; + end + % Default options + % defOPTIONS = bst_get('DuneuroOptions'); + defOPTIONS = duneuro_defaults(); + OPTIONS = struct_copy_fields(OPTIONS, defOPTIONS, 0); + + % ==== GET MESH INFO ==== + % Load tissue labels + FemMat = load(file_fullpath(OPTIONS.FemFile), 'TissueLabels'); + % Get default conductivities + OPTIONS.FemNames = FemMat.TissueLabels; + OPTIONS.FemCond = GetDefaultCondutivity(OPTIONS.FemNames); + % EEG: Select all layers; MEG: Select only the innermost layer + if isMegOnly + OPTIONS.FemSelect = zeros(size(OPTIONS.FemCond)); + OPTIONS.FemSelect(1) = 1; + else + OPTIONS.FemSelect = ones(size(OPTIONS.FemCond)); + end + % Get size of Tensors matrix + Tensors = whos('-file', file_fullpath(OPTIONS.FemFile), 'Tensors'); + OPTIONS.UseTensor = (~isempty(Tensors) && all(Tensors.size > 0)); + + % ==== FRAME STRUCTURE ==== + % Create main panel: split in top (interface) / left(standard options) / right (details) + jPanelNew = gui_component('panel'); + jPanelNew.setBorder(BorderFactory.createEmptyBorder(12,12,12,12)); + % Create top panel + jPanelTop = gui_river([1,1], [0,6,6,6]); + jPanelNew.add(jPanelTop, BorderLayout.NORTH); + % Create left panel + jPanelLeft = java_create('javax.swing.JPanel'); + jPanelLeft.setLayout(GridBagLayout()); + jPanelNew.add(jPanelLeft, BorderLayout.WEST); + % Create right panel + jPanelRight = java_create('javax.swing.JPanel'); + jPanelRight.setLayout(GridBagLayout()); + jPanelRight.setBorder(BorderFactory.createEmptyBorder(0,12,0,0)); + jPanelNew.add(jPanelRight, BorderLayout.EAST); + % Default constrains + c = GridBagConstraints(); + c.fill = GridBagConstraints.HORIZONTAL; + c.weightx = 1; + c.weighty = 0; + + % ===== PANEL LEFT: FEM LAYERS ===== + jPanelLayers = gui_river([2,2], [0,6,6,6], 'FEM layers & conductivities'); + nLayers = length(OPTIONS.FemNames); + jCheckLayer = javaArray('javax.swing.JCheckBox', nLayers); + jTextCond = javaArray('javax.swing.JComponent', nLayers); + % Loop on each layer + for i = 1:nLayers + % Add layer + jCheckLayer(i) = gui_component('checkbox', jPanelLayers, 'br', OPTIONS.FemNames{i}, [], [], @(h,ev)UpdatePanel(), []); + if ~OPTIONS.UseTensor + jTextCond(i) = gui_component('texttime', jPanelLayers, 'tab', num2str(OPTIONS.FemCond(i), '%g'), [], [], [], []); + else + jTextCond(i) = gui_component('label', jPanelLayers, 'tab', '(using tensors)', [], [], [], []); + end + % Default selection of layers + jCheckLayer(i).setSelected(OPTIONS.FemSelect(i)); + end + c.gridy = 1; + jPanelLeft.add(jPanelLayers, c); + + % ==== PANEL RIGHT: FEM SOURCE MODEL ==== + jPanelSrcModel = gui_river([1,1], [0,6,6,6], 'FEM source model'); + jGroupSrcModel = ButtonGroup(); + jRadioSrcModelVen = gui_component('radio', jPanelSrcModel, 'br', 'Multipolar Venant (optimal, faster)', jGroupSrcModel, '', @(h,ev)UpdatePanel(1), []); + jRadioSrcModelSub = gui_component('radio', jPanelSrcModel, 'br', 'Local Subtraction (recomended, slower)', jGroupSrcModel, '', @(h,ev)UpdatePanel(1), []); + jRadioSrcModelPar = gui_component('radio', jPanelSrcModel, 'br', 'Partial Integration (basic, fastest)', jGroupSrcModel, '', @(h,ev)UpdatePanel(1), []); + switch lower(OPTIONS.SrcModel) + case 'multipolar_venant', jRadioSrcModelVen.setSelected(1); + case 'local_subtraction', jRadioSrcModelSub.setSelected(1); + case 'partial_integration', jRadioSrcModelPar.setSelected(1); + end + c.gridy = 3; + jPanelLeft.add(jPanelSrcModel, c); + +% ==== PANEL RIGHT: INPUT OPTIONS ==== +jPanelInput = gui_river([1,1], [0,6,6,6], 'Source space'); + % Shrink source space + gui_component('label', jPanelInput, '', 'Shrink source space: ', [], '', [], []); + jTextSrcShrink = gui_component('texttime', jPanelInput, '', '', [], '', [], []); + gui_validate_text(jTextSrcShrink, [], [], {0,100,100}, '', 0, OPTIONS.SrcShrink, []); + gui_component('label', jPanelInput, '', ' mm'); + % Force source space inside grey matter + iGM = find(CheckType(OPTIONS.FemNames, 'gray'), 1); + if ~isempty(iGM) + jCheckSrcForceInGM = gui_component('checkbox', jPanelInput, 'br', ['Force source space inside layer "' OPTIONS.FemNames{iGM} '"'], [], '', [], []); + else + jCheckSrcForceInGM = []; + end +c.gridy = 3; +jPanelRight.add(jPanelInput, c); + +% ==== PANEL RIGHT: SOLVER OPTIONS ==== +jPanelOptSolver = gui_river([3,3], [0,6,6,6], 'Solver options'); + % Number of thread + gui_component('label', jPanelOptSolver, [], 'Number of Thread : ', [], '', [], []); + jTextNbOfThread = gui_component('texttime', jPanelOptSolver, 'tab', '', [], '', [], []); + gui_validate_text(jTextNbOfThread, [], [], 0:OPTIONS.NbOfThreadMax, '', 0, OPTIONS.NbOfThread, []); + gui_component('label', jPanelOptSolver, '[]', ['[1, ' num2str(OPTIONS.NbOfThreadMax) ']']); +c.gridy = 2; +jPanelRight.add(jPanelOptSolver, c); + + % ==== PANEL RIGHT: MEG COMPUTATIONS OPTIONS ==== + jPanelMegComputationOption = gui_river([1,1], [0,6,6,6], 'MEG computation options'); + if isMeg + % Use integration Points, recommended for high mesh density + jCheckUseIntegrationPoint = gui_component('checkbox', jPanelMegComputationOption, 'br', 'Use MEG integration points', [], '', [], []); + % Enable MEG cache memory for high mesh density if users do not + % high memory, or want to use the integration points + jCheckEnableCacheMemory = gui_component('checkbox', jPanelMegComputationOption, 'br', 'Enable cache memory', [], '', [], []); + % Enable the MEG Computation per block of sensors + ... jCheckMegPerBlockOfSensor = gui_component('checkbox', jPanelMegComputationOption, 'br', 'Compute per block of sensors [Todo]', [], '', [], []); + % Set jCheckUseIntegrationPoint to 1 as default option + if (OPTIONS.UseIntegrationPoint) + jCheckUseIntegrationPoint.setSelected(1); + end + c.gridy = 4; + jPanelRight.add(jPanelMegComputationOption, c); + else + jCheckUseIntegrationPoint = []; + jCheckEnableCacheMemory = []; + jCheckMegPerBlockOfSensor = []; + end + + % ==== PANEL RIGHT: OUTPUT OPTIONS ==== + jPanelOutput = gui_river([1,1], [0,6,6,6], 'Output options'); + % Save transfer matrix + jCheckSaveTransfer = gui_component('checkbox', jPanelOutput, '', 'Save transfer matrix', [], '', [], []); + if (OPTIONS.BstSaveTransfer) + jCheckSaveTransfer.setSelected(1); + end + c.gridy = 5; + jPanelRight.add(jPanelOutput, c); + + % ===== GLUE ===== + % Add a glue at the bottom of the left panel (for appropriate scaling with the right) + c.gridy = 6; + c.weighty = 1; + jPanelLeft.add(Box.createVerticalGlue(), c); + % Add a glue at the bottom of the right panel (for appropriate scaling with the left) + c.gridy = 5; + c.weighty = 1; + jPanelRight.add(Box.createVerticalGlue(), c); + + % ===== VALIDATION BUTTONS ===== + jPanelValid = gui_river([1,1], [12,6,6,6]); + % Expert/normal mode + jButtonExpert = gui_component('button', jPanelValid, [], 'Show details', [], [], @SwitchExpertMode_Callback, []); + gui_component('label', jPanelValid, 'hfill', ' '); + % Ok/Cancel + gui_component('Button', jPanelValid, 'right', 'Cancel', [], [], @ButtonCancel_Callback, []); + gui_component('Button', jPanelValid, [], 'OK', [], [], @ButtonOk_Callback, []); + c.gridy = 4; + c.weighty = 0; + jPanelLeft.add(jPanelValid, c); + + % ===== HELP BUTTON ===== + jPanelHelp = gui_river([1,1], [12,6,6,6]); + jButtonHelp = gui_component('Button', jPanelHelp, 'right', 'Help', [], [], @(h,ev)web('https://neuroimage.usc.edu/brainstorm/Tutorials/Duneuro', '-browser'), []); + jButtonHelp.setForeground(Color(.7, 0, 0)); + c.gridy = 6; + c.weighty = 0; + jPanelRight.add(jPanelHelp, c); + + % ===== PANEL CREATION ===== + % Return a mutex to wait for panel close + bst_mutex('create', panelName); + % Create the BstPanel object that is returned by the function + ctrl = struct('jCheckLayer', jCheckLayer, ... + 'jTextCond', jTextCond, ... + 'jRadioSrcModelVen', jRadioSrcModelVen, ... + 'jRadioSrcModelSub', jRadioSrcModelSub, ... + 'jRadioSrcModelPar', jRadioSrcModelPar, ... + 'jTextSrcShrink', jTextSrcShrink, ... + 'jTextNbOfThread', jTextNbOfThread, ... + 'jCheckSrcForceInGM', jCheckSrcForceInGM, ... + 'jCheckSaveTransfer', jCheckSaveTransfer, ... + 'jCheckUseIntegrationPoint', jCheckUseIntegrationPoint,... + 'jCheckEnableCacheMemory', jCheckEnableCacheMemory,... + ...'jCheckMegPerBlockOfSensor', jCheckMegPerBlockOfSensor,... + 'UseTensor', OPTIONS.UseTensor); + ctrl.FemNames = OPTIONS.FemNames; + % Create the BstPanel object that is returned by the function + bstPanelNew = BstPanel(panelName, jPanelNew, ctrl); + % Update comments + UpdatePanel(1); + +%% ================================================================================= +% === LOCAL CALLBACKS ============================================================ +% ================================================================================= + %% ===== BUTTON: CANCEL ===== + function ButtonCancel_Callback(varargin) + % Close panel + gui_hide(panelName); + end + + %% ===== BUTTON: OK ===== + function ButtonOk_Callback(varargin) + % Release mutex and keep the panel opened + bst_mutex('release', panelName); + end + + %% ===== SWITCH EXPERT MODE ===== + function SwitchExpertMode_Callback(varargin) + % Toggle expert mode + ExpertMode = bst_get('ExpertMode'); + bst_set('ExpertMode', ~ExpertMode); + % Update comment + UpdatePanel(1); + end + + %% ===== UPDATE PANEL ====== + % USAGE: UpdatePanel(isForced = 0) + function UpdatePanel(isForced) + % Default values + if (nargin < 1) || isempty(isForced) + isForced = 0; + end + % Expert mode / Normal mode + if isForced + ExpertMode = bst_get('ExpertMode'); + % Show/hide panels + jPanelRight.setVisible(ExpertMode); + jPanelSrcModel.setVisible(ExpertMode); + jPanelInput.setVisible(ExpertMode); + jPanelOutput.setVisible(ExpertMode); + jPanelOptSolver.setVisible(ExpertMode); + jPanelMegComputationOption.setVisible(ExpertMode); + % Update expert button + if ExpertMode + jButtonExpert.setText('Hide details'); + else + jButtonExpert.setText('Show details'); + end + % Get old panel + [bstPanelOld, iPanel] = bst_get('Panel', 'DuneuroOptions2026'); + container = get(bstPanelOld, 'container'); + % Re-pack frame + if ~isempty(container) + jFrame = container.handle{1}; + if ~isempty(jFrame) + jFrame.pack(); + end + end + end + % FEM Layers + for j = 1:nLayers + jTextCond(j).setEnabled(jCheckLayer(j).isSelected()); + end + end +end + + +%% ================================================================================= +% === EXTERNAL CALLBACKS ========================================================= +% ================================================================================= +%% ===== GET PANEL CONTENTS ===== +function s = GetPanelContents() %#ok + % Get panel controls handles + ctrl = bst_get('PanelControls', 'DuneuroOptions2026'); + if isempty(ctrl) + s = []; + return; + end + % Get default duneuro options + s = duneuro_defaults(); + + % FEM layers + for i = 1:length(ctrl.jCheckLayer) + s.FemSelect(i) = ctrl.jCheckLayer(i).isSelected(); + if ~ctrl.UseTensor + s.FemCond(i) = str2double(char(ctrl.jTextCond(i).getText())); + end + end + s.UseTensor = ctrl.UseTensor; + % Source model + if ctrl.jRadioSrcModelVen.isSelected() + s.SrcModel = 'multipolar_venant'; + elseif ctrl.jRadioSrcModelSub.isSelected() + s.SrcModel = 'local_subtraction'; + elseif ctrl.jRadioSrcModelPar.isSelected() + s.SrcModel = 'partial_integration'; + end + % Input options + s.SrcShrink = str2double(ctrl.jTextSrcShrink.getText()); + s.NbOfThread = str2double(ctrl.jTextNbOfThread.getText()); + if ~isempty(ctrl.jCheckSrcForceInGM) + s.SrcForceInGM = ctrl.jCheckSrcForceInGM.isSelected(); + else + s.SrcForceInGM = 0; + end + % Output options + s.BstSaveTransfer = ctrl.jCheckSaveTransfer.isSelected(); + % MEG Computation options + if ~isempty(ctrl.jCheckUseIntegrationPoint) + s.UseIntegrationPoint = ctrl.jCheckUseIntegrationPoint.isSelected(); + else + s.UseIntegrationPoint = 1; + end + + if ~isempty(ctrl.jCheckEnableCacheMemory) + s.EnableCacheMemory = ctrl.jCheckEnableCacheMemory.isSelected(); + else + s.EnableCacheMemory = 0; + end +end + + +%% ===== GET DEFAULT CONDUCTIVITIES ===== +function FemCond = GetDefaultCondutivity(FemNames, Reference) + % Default reference + if (nargin < 2) || isempty(Reference) + Reference = 'simbio'; + end + % Default conductivity values + switch (Reference) + case 'simbio' + conductivity = [0.14, 0.33, 1.79, 0.025, 0.008, 0.43]; + case 'simnibs' + conductivity = [0.126, 0.275, 1.654, 0.010, 0.465]; % SimNIBS paper & soft + end + % By default: conductivity of the grey matter + FemCond = conductivity(2) * ones(1, length(FemNames)); + % Detect the conductivity layer by name + for i = 1:length(FemNames) + if CheckType(FemNames{i}, 'white') + FemCond(i) = conductivity(1); + elseif CheckType(FemNames{i}, 'gray') + FemCond(i) = conductivity(2); + elseif CheckType(FemNames{i}, 'csf') + FemCond(i) = conductivity(3); + elseif CheckType(FemNames{i}, 'skull') + FemCond(i) = conductivity(5); + elseif CheckType(FemNames{i}, 'scalp') + FemCond(i) = conductivity(6); + end + end +end + + + % Check the type of a layer based on its name +function isType = CheckType(strName, strType) + if iscell(strName) + isType = cellfun(@(c)CheckType(c, strType), strName); + else + strName = lower(strName); + switch strType + case 'white' + isType = ~isempty(strfind(strName, 'white')) || ~isempty(strfind(strName, 'wm')); + case 'gray' + isType = ~isempty(strfind(strName, 'brain')) || ~isempty(strfind(strName, 'grey')) || ~isempty(strfind(strName, 'gray')) || ~isempty(strfind(strName, 'gm')) || ~isempty(strfind(strName, 'cortex')); + case 'csf' + isType = ~isempty(strfind(strName, 'csf')) || ~isempty(strfind(strName, 'inner')); + case 'skull' + isType = ~isempty(strfind(strName, 'spong')) || ... % 'Skull spongia' + ~isempty(strfind(strName, 'bone')) || ~isempty(strfind(strName, 'skull')) || ~isempty(strfind(strName, 'outer')); % 'Skull compacta' + case 'scalp' + isType = ~isempty(strfind(strName, 'skin')) || ~isempty(strfind(strName, 'scalp')) || ~isempty(strfind(strName, 'head')); + end + end +end \ No newline at end of file diff --git a/toolbox/forward/panel_headmodel.m b/toolbox/forward/panel_headmodel.m index 21575109f5..750d03e9a0 100644 --- a/toolbox/forward/panel_headmodel.m +++ b/toolbox/forward/panel_headmodel.m @@ -80,6 +80,7 @@ jComboMethodMEG.addItem(BstListItem('os_meg', '', 'Overlapping spheres', [])); jComboMethodMEG.addItem(BstListItem('openmeeg', '', 'OpenMEEG BEM', [])); jComboMethodMEG.addItem(BstListItem('duneuro', '', 'DUNEuro FEM', [])); + jComboMethodMEG.addItem(BstListItem('duneuro2026', '', 'DUNEuro FEM 2026', [])); jComboMethodMEG.setSelectedIndex(1); else jCheckMethodMEG = []; @@ -95,6 +96,7 @@ jComboMethodEEG.addItem(BstListItem('eeg_3sphereberg', '', '3-shell sphere', [])); jComboMethodEEG.addItem(BstListItem('openmeeg', '', 'OpenMEEG BEM', [])); jComboMethodEEG.addItem(BstListItem('duneuro', '', 'DUNEuro FEM', [])); + jComboMethodEEG.addItem(BstListItem('duneuro2026', '', 'DUNEuro FEM 2026', [])); jComboMethodEEG.setSelectedIndex(1); else jCheckMethodEEG = []; @@ -109,6 +111,7 @@ jComboMethodECOG = gui_component('ComboBox', jPanelMethod, 'tab hfill', [], [], [], @UpdateComment, []); jComboMethodECOG.addItem(BstListItem('openmeeg', '', 'OpenMEEG BEM', [])); jComboMethodECOG.addItem(BstListItem('duneuro', '', 'DUNEuro FEM', [])); + jComboMethodECOG.addItem(BstListItem('duneuro2026', '', 'DUNEuro FEM 2026', [])); jComboMethodECOG.setSelectedIndex(0); else jCheckMethodECOG = []; @@ -123,6 +126,7 @@ jComboMethodSEEG = gui_component('ComboBox', jPanelMethod, 'tab hfill', [], [], [], @UpdateComment, []); jComboMethodSEEG.addItem(BstListItem('openmeeg', '', 'OpenMEEG BEM', [])); jComboMethodSEEG.addItem(BstListItem('duneuro', '', 'DUNEuro FEM', [])); + jComboMethodSEEG.addItem(BstListItem('duneuro2026', '', 'DUNEuro FEM 2026', [])); jComboMethodSEEG.setSelectedIndex(0); else jCheckMethodSEEG = []; @@ -439,6 +443,7 @@ function UpdateComment(varargin) sMethod.Comment = strrep(sMethod.Comment, 'eeg_3sphereberg', '3-shell sphere'); sMethod.Comment = strrep(sMethod.Comment, 'openmeeg', 'OpenMEEG BEM'); sMethod.Comment = strrep(sMethod.Comment, 'duneuro', 'DUNEuro FEM'); + sMethod.Comment = strrep(sMethod.Comment, 'duneuro2026', 'DUNEuro FEM 2026'); % Grid type if strcmpi(sMethod.HeadModelType, 'volume') sMethod.Comment = [sMethod.Comment ' (volume)']; @@ -448,6 +453,7 @@ function UpdateComment(varargin) end isOpenMEEG = any(strcmpi(allMethods, 'openmeeg')); isDuneuro = any(strcmpi(allMethods, 'duneuro')); + isDuneuro2026 = any(strcmpi(allMethods, 'duneuro2026')); % Get protocol description ProtocolInfo = bst_get('ProtocolInfo'); @@ -689,7 +695,7 @@ function UpdateComment(varargin) OPTIONS.FemFile = sSubject.Surface(sSubject.iFEM(1)).FileName; % Interactive interface to set the OpenMEEG options if OPTIONS.Interactive - DuneuroOptions = gui_show_dialog('DUNEuro options', @panel_duneuro, 1, [], OPTIONS); + DuneuroOptions = gui_show_dialog('DUNEuro Options', @panel_duneuro, 1, [], OPTIONS); if isempty(DuneuroOptions) bst_progress('stop'); return; @@ -699,6 +705,25 @@ function UpdateComment(varargin) end end + % ===== DUNEURO ===== + if isDuneuro2026 + % Get default FEM head model + if isempty(sSubject.iFEM) + errMessage = 'No FEM head model available for this subject.'; + return; + end + OPTIONS.FemFile = sSubject.Surface(sSubject.iFEM(1)).FileName; + % Interactive interface to set the OpenMEEG options + if OPTIONS.Interactive + DuneuroOptions = gui_show_dialog('DUNEuro options', @panel_duneuro2026, 1, [], OPTIONS); + if isempty(DuneuroOptions) + bst_progress('stop'); + return; + end + % Copy the selected options to the OPTIONS structure + OPTIONS = struct_copy_fields(OPTIONS, DuneuroOptions, 1); + end + end % ===== COMPUTE HEADMODEL ===== % Start process [OPTIONS, errMessage] = bst_headmodeler(OPTIONS); diff --git a/toolbox/process/functions/process_fem_tensors.m b/toolbox/process/functions/process_fem_tensors.m index d0274fc01c..d97be24840 100644 --- a/toolbox/process/functions/process_fem_tensors.m +++ b/toolbox/process/functions/process_fem_tensors.m @@ -496,6 +496,4 @@ function ClearTensors(FemFile) term3(:,3) = axis(:,3).*temp; v_rot = term1 + term2 + term3; -end - - +end \ No newline at end of file