From 28a4b24d5fa3b6496e6e67a42b5a7edaeb697ef5 Mon Sep 17 00:00:00 2001 From: tmedani Date: Wed, 4 Feb 2026 14:48:07 -0800 Subject: [PATCH 01/11] Add DUNEuro 2026 FEM backend and UI Introduce support for the DUNEuro 2026 workflow: add new implementation bst_duneuro_2026.m (exports FEM to HDF5, writes sensors/source space, runs containerized DUNEuro tasks, reads leadfields) and a new GUI panel panel_duneuro2026.m. Update duneuro_defaults.m with SrcModel2026 default. Integrate the new backend into bst_headmodeler.m so headmodeler can invoke bst_duneuro_2026 when OPTIONS.*Method == 'duneuro_2026' and include it in the method checks. --- toolbox/forward/bst_duneuro_2026.m | 958 ++++++++++++++++++++++++++++ toolbox/forward/bst_headmodeler.m | 13 +- toolbox/forward/duneuro_defaults.m | 1 + toolbox/forward/panel_duneuro2026.m | 533 ++++++++++++++++ toolbox/forward/panel_headmodel.m | 25 + 5 files changed, 1526 insertions(+), 4 deletions(-) create mode 100644 toolbox/forward/bst_duneuro_2026.m create mode 100644 toolbox/forward/panel_duneuro2026.m diff --git a/toolbox/forward/bst_duneuro_2026.m b/toolbox/forward/bst_duneuro_2026.m new file mode 100644 index 0000000000..3ec7782926 --- /dev/null +++ b/toolbox/forward/bst_duneuro_2026.m @@ -0,0 +1,958 @@ +function [Gain, errMsg] = bst_duneuro_2026(cfg) +% BST_DUNEURO: Call Duneuro to compute a FEM solution for Brainstorm. +% +% USAGE: [Gain, errMsg] = bst_duneuro(cfg) + +% @============================================================================= +% 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, Juan Garcia-Prieto, 2019-2026 +% Francois Tadel 2020-2023 + +% Initialize returned values +Gain = []; +% ***** ToDo: Define how to check the container ******* +% Install/load duneuro plugin +[isInstalled, errMsg, PlugDesc] = bst_plugin('Install', 'duneuro', cfg.Interactive); +if ~isInstalled + return; +end +bst_plugin('SetProgressLogo', 'duneuro'); + +% Get DUNEuro executable +DuneuroExe = bst_fullfile(PlugDesc.Path, PlugDesc.SubFolder, 'bin', ['bst_duneuro_meeg_', bst_get('OsType')]); +if ispc + DuneuroExe = [DuneuroExe, '.exe']; +else + DuneuroExe = [DuneuroExe, '.app']; +end + + + +%% ===== SENSORS ===== +% Select modality for DUNEuro +isEeg = strcmpi(cfg.EEGMethod, 'duneuro_2026') && ~isempty(cfg.iEeg); +isMeg = strcmpi(cfg.MEGMethod, 'duneuro_2026') && ~isempty(cfg.iMeg); +isEcog = strcmpi(cfg.ECOGMethod, 'duneuro_2026') && ~isempty(cfg.iEcog); +isSeeg = strcmpi(cfg.SEEGMethod, 'duneuro_2026') && ~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(iInteg)]; + 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)); +% Get mesh type +switch size(FemMat.Elements,2) + case 4, ElementType = 'tetrahedron'; + case 8, ElementType = 'hexahedron'; +end +% 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 +% Hexa mesh: detect whether the geometry was adapted +if strcmpi(ElementType, 'hexahedron') + GeometryAdapted = []; + % Detect using the options of the Brainstorm process that created the file + if isfield(FemMat, 'History') && ~isempty(FemMat.History) && ~isempty(strfind([FemMat.History{:,3}], 'NodeShift')) + strOptions = [FemMat.History{:,3}]; + iTag = strfind(strOptions, 'NodeShift'); + val = sscanf(strOptions(iTag:end), 'NodeShift=%f'); + if ~isempty(val) + GeometryAdapted = (val > 0); + end + end + % Otherwise, try to guess based on the geometry + if isempty(GeometryAdapted) + % Compute the distance between the first two nodes of each element + dist = sqrt(sum([FemMat.Vertices(FemMat.Elements(:,1),1) - FemMat.Vertices(FemMat.Elements(:,2),1), ... + FemMat.Vertices(FemMat.Elements(:,2),2) - FemMat.Vertices(FemMat.Elements(:,2),2), ... + FemMat.Vertices(FemMat.Elements(:,2),3) - FemMat.Vertices(FemMat.Elements(:,2),3)] .^ 2, 2)); + % If the distance is not constant: then the geomtry is adapted + GeometryAdapted = (max(abs(dist - dist(1))) > 1e-9); + end + % Copy value in DUNEuro options + if ~isempty(GeometryAdapted) + cfg.GeometryAdapted = GeometryAdapted; + disp(['DUNEURO> Detected parameter: GeometryAdapted=', bool2str(GeometryAdapted)]); + end +end +% % Isotropic +% if (cfg.Isotropic) +% % Isotropic without tensor +% if ~cfg.UseTensor +% if strcmp(ElementType,'hexahedron') +% MeshFile = 'head_model.dgf'; +% else +% MeshFile = 'head_model.msh'; +% end +% % Isotropic with tensor +% else +% if strcmp(ElementType,'hexahedron') +% errMsg = 'Using the tensor model with hexahedral mesh is not supported for now.'; +% return; +% end +% MeshFile = 'head_model.geo'; +% end +% % Anisotropic (with tensor) +% else +% if strcmp(ElementType,'hexahedron') +% errMsg = 'Using the anisotropy model with hexahedral mesh is not supported for now.'; +% return; +% end +% MeshFile = 'head_model.geo'; +% cfg.UseTensor = true; +% 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' + % 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) + +% %%%% ============================================= +% % Now similar process for the WM with an extra and unexpedted step ... +% if ~isempty(iWM) +% disp('Checking the the WM ...'); +% %Extract WM surface +% wm_tetra = FemMat.Elements(FemMat.Tissue <= iWM,:); +% wm_face = volface(wm_tetra); +% [nwmf, ewmf] = removeisolatednode(FemMat.Vertices,wm_face); +% % check if any dipoles is inside the WM +% wMfv.vertices = nwmf; +% wMfv.faces = ewmf; +% tic +% wMin = inpolyhedron(wMfv, sCortex.Vertices); +% wMindex_in = find(wMin); +% disp(['There are ' num2str(sum(wMin)) ' dipoles inside the WM']); +% disp('Moving these dipoles to the GM tissues ...'); +% twm = toc; +% if ~isempty(wMindex_in) +% % 1- move the dipole from inside the WM to the GM surface +% % ==> this is for testing, when we use the centroide +% % directely, some dipole remains within the WM ... +% GMcentroide = 0; % just for testing, to use directely the GM centroides +% if GMcentroide == 1 +% k = dsearchn(elem_centroide,sCortex.Vertices(wMindex_in,:)); +% NewVertices(wMindex_in ,:) = ElemCenter(k,:); +% wMoutFinal = inpolyhedron(wMfv, NewVertices); +% disp(['Now, there are ' num2str(sum(wMoutFinal)) ' dipoles inside the WM']); +% else +% k = dsearchn(gMfv.vertices,sCortex.Vertices(wMindex_in,:)); +% NewVertices(wMindex_in ,:) = gMfv.vertices(k,:); +% wMoutFinal = inpolyhedron(wMfv, NewVertices); +% disp(['There are ' num2str(sum(wMoutFinal)) ' dipoles inside the WM']); +% disp(['These dipoles are moved to the GM outer surface nodes']); +% disp(['Moving these dipoles to the nearest centroide of the GM element...']); +% % Now, we move these same dipole from the surface to the +% % nearest centroide ==> to be sure it's inside the GM +% % ==> fill partially the FEM condition +% % move again to the centoide +% k = dsearchn(ElemCenter,NewVertices(wMindex_in,:)); +% NewVertices(wMindex_in ,:) = ElemCenter(k,:); +% % just to check +% wMoutFinal = inpolyhedron(wMfv, NewVertices); +% disp(['Now, there are ' num2str(sum(wMoutFinal)) ' dipoles inside the WM']); +% disp(['All the dipoles are moved to nearest centroid of the GM']); +% end +% end +% end + + end + case 'mixed' + % TODO : not used ? +end + + +%% ===== SOURCE MODEL ===== +bst_progress('text', 'DUNEuro: Writing temporary files...'); +% Write the source/dipole file +% DipoleFile = fullfile(TmpDir, 'dipole_model.txt'); +% Unconstrained orientation for each dipole. ie the output file have 3*N dipoles (3 directions for each) +% [x1,y1,z1, 1, 0, 0,; +% x1,y1,z1, 0, 1, 0; +% x1,y1,z1, 0, 0, 1; +% x2,y2,z2, 1, 0, 0; +% .... ]; +% dipoles = [kron(cfg.GridLoc, ones(3,1)), kron(ones(size(cfg.GridLoc,1), 1), eye(3))]; +% fid = fopen(DipoleFile, 'wt+'); +% fprintf(fid, '%d %d %d %d %d %d \n', dipoles'); +% fclose(fid); +[iOk, errMsg] = 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 + % ********* TODO *********** + % Get the correct weigh from the original channel file + coil_to_channel_transform = eye(length(MegChannels)); + coil_to_channel_transform = coil_to_channel_transform.* MegChannels(:,end); + dnbst_write_magnetometers(TmpDir, MegChannels(:,2:4), MegChannels(:,5:7), coil_to_channel_transform); +end + +%% ===== CONDUCTIVITY MODEL ===== +% % Isotropic without tensor +% if ~cfg.UseTensor +% CondFile = fullfile(TmpDir, 'conductivity_model.con'); +% fid = fopen(CondFile, 'w'); +% fprintf(fid, '%d\t', cfg.FemCond); +% fclose(fid); +% % With tensor (isotropic or anisotropic) +% else +% CondFile = fullfile(TmpDir, 'conductivity_model.knw'); +% % Transformation matrix and tensor mapping on each direction +% CondTensor = zeros(length(FemMat.Elements),6) ; +% for ind =1 : length(FemMat.Elements) +% temp0 = reshape(FemMat.Tensors(ind,:),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 +% CondTensor(ind,:) = [temp(1) temp(5) temp(9) temp(4) temp(8) temp(7)]; % this is the right order +% end +% % write the tensors +% out_fem_knw(FemMat, CondTensor, CondFile); +% end +% % Write mesh model +% MeshFile = fullfile(TmpDir, MeshFile); +% out_fem(FemMat, MeshFile); +%% DN2026 +% export volume conductor model +% write_volume_conductor(TmpDir, nodes, elements, labels, tensors); +% write hdf5 files : +% online viewer of hdf5 files : https://myhdf5.hdfgroup.org/ +dnbst_write_volume_conductor(TmpDir, FemMat, cfg.FemCond); + + +% run DUNEuro container +% Source Space configuration +sourcespace_config = []; +sourcespace_config.name = 'create_sourcespace'; +% sourcespace_config.sourcecompartment = '1'; +% sourcespace_config.gridsize = '2.0'; +% sourcespace_config.enforceDistanceCondition = 'True'; +% sourcespace_config.distanceThreshold = '0.5'; +sourcespace_config.enforceVenantCondition = 'True'; +% sourcespace_config.visualizeSourceSpace = 'True'; + + +% 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 = '1e-16'; +transfer_matrix_config.nr_threads = '-1'; + +% Leadfield configuration +leadfield_config = []; +leadfield_config.name = 'compute_leadfield'; +if strcmp(dnModality, 'eeg') + leadfield_config.do_meg = 'False'; + leadfield_config.do_eeg = 'True'; + leadfield_config.eeg_scaling = '1e0'; + leadfield_config.meg_scaling = '1e5'; +end +if strcmp(dnModality, 'meg') + leadfield_config.do_meg = 'True'; + leadfield_config.do_eeg = 'False'; + leadfield_config.eeg_scaling = '1e0'; + leadfield_config.meg_scaling = '1e5'; +end +if strcmp(dnModality, 'meeg') + leadfield_config.do_meg = 'True'; + leadfield_config.do_eeg = 'True'; + leadfield_config.eeg_scaling = '1e0'; + leadfield_config.meg_scaling = '1e5'; +end +leadfield_config.sourcemodel = cfg.SrcModel2026; % [select from the interface: 'multipolar_venant', 'local_subtraction', 'partial_integration'] +leadfield_config.nr_threads = '-1'; + +%% ===== RUN DUNEURO ====== +% Assemble command line +% callStr = ['"' DuneuroExe '"' ' ' '"' IniFile '"']; +bst_progress('text', 'DUNEuro: Computing leadfield...'); +runner = 'docker'; +% disp(['DUNEURO> System call: ' callStr]); +tic; +% Call DUNEuro +if 0 + [status, errMsg ] = bst_run_duneuro_task(TmpDir, sourcespace_config, runner); % wheew is the outputs? +end +% Compute transfer matrix +[status, errMsg ] = bst_run_duneuro_task(TmpDir, transfer_matrix_config, runner); +% Compute the Leadfield +[status, errMsg ] = bst_run_duneuro_task(TmpDir, leadfield_config, runner); + +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...'); +% % EEG +% if (isEeg || isEcog || isSeeg) +% GainEeg = in_duneuro_bin(fullfile(TmpDir, cfg.BstEegLfFile))'; +% end +% read data +% 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 isMeg + Gain(cfg.iMeg,:) = GainMeg; + % scaling the MEG Gain matrix + Gain(cfg.iMeg,:) = GainMeg/1000; +end +if (isEeg || isEcog || isSeeg) + Gain(cfg.iEeg,:) = GainEeg; +end + +% Delete the temporary files +% file_delete(TmpDir, 1, 1); + + +%% ===== SAVE TRANSFER MATRIX ====== +disp('DUNEURO> TODO: Save transferOut.dat to database.') + +% Remove logo +bst_plugin('SetProgressLogo', []); + + +% %% ===== WRITE MINI FILE ===== +% % Open the mini file +% IniFile = fullfile(TmpDir, 'duneuro_minifile.mini'); +% fid = fopen(IniFile, 'wt+'); +% % General setting +% fprintf(fid, '__name = %s\n\n', IniFile); +% if strcmp(cfg.SolverType, 'cg') +% fprintf(fid, 'type = %s\n', cfg.FemType); +% end +% fprintf(fid, 'element_type = %s\n', ElementType); +% fprintf(fid, 'solver_type = %s\n', cfg.SolverType); +% fprintf(fid, 'geometry_adapted = %s\n', bool2str(cfg.GeometryAdapted)); +% fprintf(fid, 'tolerance = %d\n', cfg.Tolerance); +% % [electrodes] +% if isEcog || isSeeg +% % Instead of selecting the electrode on the outer surface, +% % uses the nearest FEM node as the electrode location +% cfg.ElecType = 'closest_subentity_center'; +% end +% if strcmp(dnModality, 'eeg') || strcmp(dnModality, 'meeg') +% fprintf(fid, '[electrodes]\n'); +% fprintf(fid, 'filename = %s\n', fullfile(TmpDir, ElecFile)); +% fprintf(fid, 'type = %s\n', cfg.ElecType); +% fprintf(fid, 'codims = %s\n', '3'); +% end +% % [meg] +% if strcmp(dnModality, 'meg') || strcmp(dnModality, 'meeg') +% fprintf(fid, '[meg]\n'); +% fprintf(fid, 'intorderadd = %d\n', cfg.MegIntorderadd); +% fprintf(fid, 'type = %s\n', cfg.MegType); +% fprintf(fid, 'cache.enable = %s\n',bool2str(cfg.EnableCacheMemory) ); +% % [coils] +% fprintf(fid, '[coils]\n'); +% fprintf(fid, 'filename = %s\n', CoilFile); +% % [projections] +% fprintf(fid, '[projections]\n'); +% fprintf(fid, 'filename = %s\n', ProjFile); +% end +% % [dipoles] +% fprintf(fid, '[dipoles]\n'); +% fprintf(fid, 'filename = %s\n', DipoleFile); +% % [volume_conductor.grid] +% fprintf(fid, '[volume_conductor.grid]\n'); +% fprintf(fid, 'filename = %s\n', MeshFile); +% % [volume_conductor.tensors] +% fprintf(fid, '[volume_conductor.tensors]\n'); +% fprintf(fid, 'filename = %s\n', CondFile); +% % [solver] +% fprintf(fid, '[solver]\n'); +% fprintf(fid, 'solver_type = %s\n', cfg.SolvSolverType); +% fprintf(fid, 'preconditioner_type = %s\n', cfg.SolvPrecond); +% if strcmp(cfg.SolverType, 'cg') +% fprintf(fid, 'cg_smoother_type = %s\n', cfg.SolvSmootherType); +% end +% fprintf(fid, 'intorderadd = %d\n', cfg.SolvIntorderadd); +% % Discontinuous Galerkin +% if strcmp(cfg.SolverType, 'dg') +% fprintf(fid, 'dg_smoother_type = %s\n', cfg.DgSmootherType); +% fprintf(fid, 'scheme = %s\n', cfg.DgScheme); +% fprintf(fid, 'penalty = %d\n', cfg.DgPenalty); +% fprintf(fid, 'edge_norm_type = %s\n', cfg.DgEdgeNormType); +% fprintf(fid, 'weights = %s\n', bool2str(cfg.DgWeights)); +% fprintf(fid, 'reduction = %s\n', bool2str(cfg.DgReduction)); +% end +% % [solution] +% fprintf(fid, '[solution]\n'); +% fprintf(fid, 'post_process = %s\n', bool2str(cfg.SolPostProcess)); % true/false +% fprintf(fid, 'subtract_mean = %s\n', bool2str(cfg.SolSubstractMean)); % boolean +% % [solution.solver] +% fprintf(fid, '[solution.solver]\n'); +% fprintf(fid, 'reduction = %d\n', cfg.SolSolverReduction); +% % [solution.source_model] +% fprintf(fid, '[solution.source_model]\n'); +% fprintf(fid, 'type = %s\n', cfg.SrcModel); +% fprintf(fid, 'intorderadd = %d\n', cfg.SrcIntorderadd); +% fprintf(fid, 'intorderadd_lb = %d\n', cfg.SrcIntorderadd_lb); +% fprintf(fid, 'numberOfMoments = %d\n', cfg.SrcNbMoments); +% fprintf(fid, 'referenceLength = %d\n', cfg.SrcRefLen); +% fprintf(fid, 'weightingExponent = %d\n', cfg.SrcWeightExp); +% fprintf(fid, 'relaxationFactor = %e\n', 10^(-cfg.SrcRelaxFactor)); +% fprintf(fid, 'mixedMoments = %s\n', bool2str(cfg.SrcMixedMoments)); +% fprintf(fid, 'restrict = %s\n', bool2str(cfg.SrcRestrict)); +% fprintf(fid, 'initialization = %s\n', cfg.SrcInit); +% % [brainstorm] +% fprintf(fid, '[brainstorm]\n'); +% fprintf(fid, 'modality = %s\n', dnModality); +% fprintf(fid, 'output_folder = %s\n', [TmpDir, filesep]); +% fprintf(fid, 'save_eeg_transfer_file = %s\n', bool2str(cfg.BstSaveTransfer)); +% fprintf(fid, 'save_meg_transfer_file = %s\n', bool2str(cfg.BstSaveTransfer)); +% fprintf(fid, 'save_meeg_transfer_file = %s\n', bool2str(cfg.BstSaveTransfer)); +% fprintf(fid, 'eeg_transfer_filename = %s\n', cfg.BstEegTransferFile); +% fprintf(fid, 'meg_transfer_filename = %s\n', cfg.BstMegTransferFile); +% fprintf(fid, 'eeg_leadfield_filename = %s\n', cfg.BstEegLfFile); +% fprintf(fid, 'meg_leadfield_filename = %s\n', cfg.BstMegLfFile); +% % Close file +% fclose(fid); + + +% %% ===== RUN DUNEURO ====== +% % Assemble command line +% callStr = ['"' DuneuroExe '"' ' ' '"' IniFile '"']; +% bst_progress('text', 'DUNEuro: Computing leadfield...'); +% disp(['DUNEURO> System call: ' callStr]); +% tic; +% % Call DUNEuro +% status = system(callStr) +% 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...'); +% % EEG +% if (isEeg || isEcog || isSeeg) +% GainEeg = in_duneuro_bin(fullfile(TmpDir, cfg.BstEegLfFile))'; +% end +% +% %MEG +% To check with Malte: is the MEG leadfield is the finla or just the +% % seconday? +% if isMeg +% GainMeg = in_duneuro_bin(fullfile(TmpDir, cfg.BstMegLfFile))'; +% % === POST-PROCESS MEG LEADFIELD === +% % Compute the total magnetic field +% dipoles_pos_orie = [kron(cfg.GridLoc,ones(3,1)), kron(ones(length(cfg.GridLoc),1), eye(3))]; +% +% % a- Compute the MEG Primary Magnetic B-field analytically (formula of Sarvas) +% dip_pos = dipoles_pos_orie(:,1:3); +% dip_mom = dipoles_pos_orie(:,4:6); +% Bp = zeros(size(MegChannels,1), size(dip_pos,1)); +% for i = 1:size(CoilsLoc,1) +% for j = 1 : size(dip_pos,1) +% R = CoilsLoc(i,:); +% R_0 = dip_pos(j,:); +% A = R - R_0; +% a = norm(A); +% aa = A./(a^3); +% BpHelp = cross(dip_mom(j,:),aa); +% Bp(i,j) = BpHelp * CoilsOrient(i, :)'; % projection of the primary B-field along the coil orientations +% end +% end +% +% % b- The total magnetic field B = Bp + Bs; +% % full B-field +% Bs = GainMeg; +% mu = 4*pi*1e-4; % check the value of the units maybe it needs to be mu = 4*pi*1e-7 +% Bfull = (mu/(4*pi)) * (Bp - Bs); +% +% % c- Apply the weight : +% [channelIndex] = unique(MegChannels(:,1)); +% nbChannel = length(channelIndex); +% weighted_B = zeros(nbChannel,size(Bfull,2)); +% for iCh = 1 : nbChannel +% communChannel = find(iCh==MegChannels(:,1)); +% BcommunChannel = Bfull(communChannel(:),:); +% WcommunChannel = MegChannels(communChannel(:), 8: end); +% weighted_B(iCh,:) = sum (BcommunChannel.*WcommunChannel,1); +% end +% GainMeg = weighted_B; +% end + +% % Fill the unused channels with NaN +% Gain = NaN * zeros(length(cfg.Channel), 3 * length(cfg.GridLoc)); +% if isMeg +% Gain(cfg.iMeg,:) = GainMeg; +% % scaling the MEG Gain matrix +% Gain(cfg.iMeg,:) = GainMeg/1000; +% end +% if (isEeg || isEcog || isSeeg) +% Gain(cfg.iEeg,:) = GainEeg; +% end +% +% % Delete the temporary files +% file_delete(TmpDir, 1, 1); +% +% +% %% ===== SAVE TRANSFER MATRIX ====== +% disp('DUNEURO> TODO: Save transferOut.dat to database.') +% +% % Remove logo +% bst_plugin('SetProgressLogo', []); + +end + + + + +%% ================================================================================= +% === SUPPORT FUNCTIONS ========================================================= +% ================================================================================= + +%% ===== BOOL => STR ===== +function str = bool2str(bool) + if bool + str = 'true'; + else + str = 'false'; + end +end + + +function [iOk, errMsg] = dnbst_write_volume_conductor(duneuro_io_dir, FemMat, isoCond) + iOk = 0; errMsg = ''; + + nodes = FemMat.Vertices; + elements = FemMat.Elements - 1; + labels = FemMat.Tissue - 1; + + io_file_path = fullfile(duneuro_io_dir, 'duneuro_io.hdf5'); + + % % we assume that this function create the IO file. We abort the execution + % % if the file already exists + % if isfile(io_file_path) + % errMsg = ['There already exists a file at the location ', io_file_path]; + % end + + % get information about the input data + nr_nodes = size(nodes, 1); + nr_elements = size(elements, 1); + nr_labels = size(labels, 1); + nr_unique_tensors = length(isoCond); + + 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 + % Create tensor format from the isoconductivity + tensors = zeros(length(isoCond), 3, 3); + for iTissue = 1 : length(FemMat.TissueLabels) + tensors(iTissue,1,1) = isoCond(iTissue); + tensors(iTissue,2,2) = isoCond(iTissue); + tensors(iTissue,3,3) = isoCond(iTissue); + 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'); + + % % we assume that the output file already exists + % if ~isfile(io_file_path) + % error('The DUNEuro IO file does not exist, please call write_volume_conductor first'); + % end + + 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 [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'); + + % % we assume that the output file already exists + % if ~isfile(io_file_path) + % error('The DUNEuro IO file does not exist, please call write_volume_conductor first'); + % end + + 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 + % if size(electrode_positions, 2) ~= dim + % errMsg = 'The dipole orientation array must have 3 columns'; + % end + + % if strcmp(electrode_type_flag, 'measurement') + h5create(io_file_path, "/measurement/source_space/positions", [dim nr_dipoles], Datatype="double"); + h5write(io_file_path, "/measurement/source_space/positions", dipole_positions'); + + % else + % errMsg = 'Dipole flag needs to be either "measurement"'; return; + % end + 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']; +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..927f13f7c6 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('duneuro_2026', {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('duneuro_2026', {OPTIONS.MEGMethod, OPTIONS.EEGMethod, OPTIONS.ECOGMethod, OPTIONS.SEEGMethod}) + [Gain_dn, errMessage] = bst_duneuro_2026(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..7b4051624f 100644 --- a/toolbox/forward/duneuro_defaults.m +++ b/toolbox/forward/duneuro_defaults.m @@ -60,6 +60,7 @@ cfgDef.SolSolverReduction = 1e-10; % [solution.source_model] cfgDef.SrcModel = 'venant'; % partial_integration, venant, subtraction +cfgDef.SrcModel2026 = 'multipolar_venant'; % 'multipolar_venant', 'local_subtraction', 'partial_integration' cfgDef.SrcIntorderadd = 0; cfgDef.SrcIntorderadd_lb = 2; cfgDef.SrcNbMoments = 3; diff --git a/toolbox/forward/panel_duneuro2026.m b/toolbox/forward/panel_duneuro2026.m new file mode 100644 index 0000000000..384f716819 --- /dev/null +++ b/toolbox/forward/panel_duneuro2026.m @@ -0,0 +1,533 @@ +function varargout = panel_duneuro2026(varargin) +% PANEL_DUNEURO: DUNEuro options +% +% USAGE: bstPanel = panel_duneuro('CreatePanel', OPTIONS) : Call from the interactive interface +% bstPanel = panel_duneuro('CreatePanel', sProcess, sFiles) : Call from the process editor +% s = panel_duneuro('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 + +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('duneuro_2026', {OPTIONS.EEGMethod, OPTIONS.ECOGMethod, OPTIONS.SEEGMethod}); + isMeg = isequal(OPTIONS.MEGMethod, 'duneuro_2026'); + % 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 LEFT: FEM METHOD TYPE ==== +% jPanelType = gui_river([1,1], [0,6,6,6], 'FEM method type'); +% jGroupFemType = ButtonGroup(); +% jRadioFemTypeFit = gui_component('radio', jPanelType, 'br', 'Fitted', jGroupFemType, '', [], []); +% jRadioFemTypeUnfit = gui_component('radio', jPanelType, 'br', 'Unfitted', jGroupFemType, '', [], []); +% switch lower(OPTIONS.FemType) +% case 'fitted', jRadioFemTypeFit.setSelected(1); +% case 'unfitted', jRadioFemTypeUnfit.setSelected(1); +% end +% c.gridy = 2; +% jPanelLeft.add(jPanelType, c); + + % ==== PANEL LEFT: FEM SOLVER TYPE ==== +% jPanelSolverType = gui_river([1,1], [0,6,6,6], 'FEM solver type'); +% jGroupSolverType = ButtonGroup(); +% jRadioSolverTypeCg = gui_component('radio', jPanelSolverType, 'br', 'CG: Continuous Galerkin', jGroupSolverType, '', [], []); +% jRadioSolverTypeDg = gui_component('radio', jPanelSolverType, 'br', 'DG: Discontinuous Galerkin', jGroupSolverType, '', [], []); +% switch lower(OPTIONS.SolverType) +% case 'cg', jRadioSolverTypeCg.setSelected(1); +% case 'dg', jRadioSolverTypeDg.setSelected(1); +% end +% c.gridy = 2; +% jPanelLeft.add(jPanelSolverType, 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: VENANT OPTIONS ==== +% jPanelOptVen = gui_river([3,3], [0,6,6,6], 'Venant options'); +% % Number of moments +% gui_component('label', jPanelOptVen, [], 'Number of moments (1-5): ', [], '', [], []); +% jTextNbMoments = gui_component('texttime', jPanelOptVen, 'tab', '', [], '', [], []); +% gui_validate_text(jTextNbMoments, [], [], 1:5, '', 0, OPTIONS.SrcNbMoments, []); +% % Reference length +% gui_component('label', jPanelOptVen, 'br', 'Reference length (1-100): ', [], '', [], []); +% jTextRefLen = gui_component('texttime', jPanelOptVen, 'tab', '', [], '', [], []); +% gui_validate_text(jTextRefLen, [], [], 1:100, '', 0, OPTIONS.SrcRefLen, []); +% % Weighting exponent +% gui_component('label', jPanelOptVen, 'br', 'Weighting exponent (1-3): ', [], '', [], []); +% jTextWeightExp = gui_component('texttime', jPanelOptVen, 'tab', '', [], '', [], []); +% gui_validate_text(jTextWeightExp, [], [], 1:3, '', 0, OPTIONS.SrcWeightExp, []); +% % Relaxation Factor +% gui_component('label', jPanelOptVen, 'br', 'Relaxation factor exponent (3-9): ', [], '', [], []); +% jTextRelaxFactor = gui_component('texttime', jPanelOptVen, 'tab', sprintf('%e', OPTIONS.SrcRelaxFactor), [], '', [], []); +% gui_validate_text(jTextRelaxFactor, [], [], 3:9, '', 0, OPTIONS.SrcRelaxFactor, []); +% % Mixed moments +% jCheckMixedMoments = gui_component('checkbox', jPanelOptVen, 'br', 'Mixed moments', [], '', [], []); +% if (OPTIONS.SrcMixedMoments == 1) +% jCheckMixedMoments.setSelected(1); +% end +% % Restrict +% jCheckRestrict = gui_component('checkbox', jPanelOptVen, 'br', 'Restrict', [], '', [], []); +% if (OPTIONS.SrcRestrict == 1) +% jCheckRestrict.setSelected(1); +% end +% c.gridy = 1; +% jPanelRight.add(jPanelOptVen, c); + + % ==== PANEL RIGHT: SUBTRACTION OPTIONS ==== +% jPanelOptSub = gui_river([3,3], [0,6,6,6], 'Subtraction options'); +% % Number of moments +% gui_component('label', jPanelOptSub, [], 'intorderadd (1-5): ', [], '', [], []); +% jTextIntorderadd = gui_component('texttime', jPanelOptSub, 'tab', '', [], '', [], []); +% gui_validate_text(jTextIntorderadd, [], [], 0:5, '', 0, OPTIONS.SrcIntorderadd, []); +% % Number of moments +% gui_component('label', jPanelOptSub, 'br', 'intorderadd_lb (1-5): ', [], '', [], []); +% jTextIntorderadd_lb = gui_component('texttime', jPanelOptSub, 'tab', '', [], '', [], []); +% gui_validate_text(jTextIntorderadd_lb, [], [], 0:5, '', 0, OPTIONS.SrcIntorderadd_lb, []); +% c.gridy = 2; +% jPanelRight.add(jPanelOptSub, 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: 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, ... + ... 'jRadioFemTypeFit', jRadioFemTypeFit, ... + ... 'jRadioFemTypeUnfit', jRadioFemTypeUnfit, ... + ... 'jRadioSolverTypeCg', jRadioSolverTypeCg, ... + ... 'jRadioSolverTypeDg', jRadioSolverTypeDg, ... + 'jRadioSrcModelVen', jRadioSrcModelVen, ... + 'jRadioSrcModelSub', jRadioSrcModelSub, ... + 'jRadioSrcModelPar', jRadioSrcModelPar, ... + ... 'jTextNbMoments', jTextNbMoments, ... + ... 'jTextRefLen', jTextRefLen, ... + ... 'jTextWeightExp', jTextWeightExp, ... + ... 'jTextRelaxFactor', jTextRelaxFactor, ... + ... 'jCheckMixedMoments', jCheckMixedMoments, ... + ... 'jCheckRestrict', jCheckRestrict, ... + ... 'jTextIntorderadd', jTextIntorderadd, ... + ... 'jTextIntorderadd_lb', jTextIntorderadd_lb, ... + 'jTextSrcShrink', jTextSrcShrink, ... + '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); +% jPanelType.setVisible(ExpertMode); +% jPanelSolverType.setVisible(ExpertMode); + jPanelSrcModel.setVisible(ExpertMode); +% jPanelOptVen.setVisible(ExpertMode && jRadioSrcModelVen.isSelected()); +% jPanelOptSub.setVisible(ExpertMode && jRadioSrcModelSub.isSelected()); + jPanelInput.setVisible(ExpertMode); + jPanelOutput.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; +% % FEM method type +% if ctrl.jRadioFemTypeFit.isSelected() +% s.FemType = 'fitted'; +% elseif ctrl.jRadioFemTypeUnfit.isSelected() +% s.FemType = 'unfitted'; +% end + % FEM solver type + % if ctrl.jRadioSolverTypeCg.isSelected() + % s.SolverType = 'cg'; + % elseif ctrl.jRadioSolverTypeDg.isSelected() + % s.SolverType = 'dg'; + % end + % 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 + % % Venant options + % if strcmpi(s.SrcModel, 'venant') + % s.SrcNbMoments = str2double(ctrl.jTextNbMoments.getText()); + % s.SrcRefLen = str2double(ctrl.jTextRefLen.getText()); + % s.SrcWeightExp = str2double(ctrl.jTextWeightExp.getText()); + % s.SrcRelaxFactor = str2double(ctrl.jTextRelaxFactor.getText()); + % s.SrcMixedMoments = ctrl.jCheckMixedMoments.isSelected(); + % s.SrcRestrict = ctrl.jCheckRestrict.isSelected(); + % % Subtraction options + % elseif strcmpi(s.SrcModel, 'subtraction') + % s.SrcIntorderadd = str2double(ctrl.jTextIntorderadd.getText()); + % s.SrcIntorderadd_lb = str2double(ctrl.jTextIntorderadd_lb.getText()); + % end + % Input options + s.SrcShrink = str2double(ctrl.jTextSrcShrink.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 + +% if ~isempty(ctrl.jCheckMegPerBlockOfSensor) +% s.MegPerBlockOfSensor = ctrl.jCheckMegPerBlockOfSensor.isSelected(); +% else +% s.MegPerBlockOfSensor = 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 + + +%% ===== DETECTION FUNCTION ===== +% 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 diff --git a/toolbox/forward/panel_headmodel.m b/toolbox/forward/panel_headmodel.m index 21575109f5..11d3f450d7 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('duneuro_2026', '', '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('duneuro_2026', '', 'DUNEuro FEM2026', [])); 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('duneuro_2026', '', 'DUNEuro FEM2026', [])); 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('duneuro_2026', '', 'DUNEuro FEM2026', [])); 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, 'duneuro_2026', 'DUNEuro FEM2026'); % 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, 'duneuro_2026')); % Get protocol description ProtocolInfo = bst_get('ProtocolInfo'); @@ -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); From 92952ae268bdaf8a92d6ec53376e08a5eda5a5e5 Mon Sep 17 00:00:00 2001 From: tmedani Date: Sat, 7 Feb 2026 18:07:51 -0800 Subject: [PATCH 02/11] Refactor DUNEuro container workflow and I/O Rework bst_duneuro_2026 to use the DUNEuro container workflow and simplify I/O. Removed the old and unused sections. --- toolbox/forward/bst_duneuro_2026.m | 547 ++++------------------------- 1 file changed, 78 insertions(+), 469 deletions(-) diff --git a/toolbox/forward/bst_duneuro_2026.m b/toolbox/forward/bst_duneuro_2026.m index 3ec7782926..0cb7f27942 100644 --- a/toolbox/forward/bst_duneuro_2026.m +++ b/toolbox/forward/bst_duneuro_2026.m @@ -21,29 +21,25 @@ % For more information type "brainstorm license" at command prompt. % =============================================================================@ % -% Authors: Takfarinas Medani, Juan Garcia-Prieto, 2019-2026 +% 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 plugin -[isInstalled, errMsg, PlugDesc] = bst_plugin('Install', 'duneuro', cfg.Interactive); -if ~isInstalled - return; -end +% 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 = 'docker'; +% Get DUNEuro container executable bst_plugin('SetProgressLogo', 'duneuro'); -% Get DUNEuro executable -DuneuroExe = bst_fullfile(PlugDesc.Path, PlugDesc.SubFolder, 'bin', ['bst_duneuro_meeg_', bst_get('OsType')]); -if ispc - DuneuroExe = [DuneuroExe, '.exe']; -else - DuneuroExe = [DuneuroExe, '.app']; -end - - - %% ===== SENSORS ===== % Select modality for DUNEuro isEeg = strcmpi(cfg.EEGMethod, 'duneuro_2026') && ~isempty(cfg.iEeg); @@ -109,7 +105,6 @@ % Get mesh type switch size(FemMat.Elements,2) case 4, ElementType = 'tetrahedron'; - case 8, ElementType = 'hexahedron'; end % Remove unselected tissue from the head model (MEG only), this could be also done for sEEG later if strcmp(dnModality, 'meg') @@ -122,62 +117,10 @@ errMsg = 'Reduced head model cannot be used when computing MEG+EEG simultaneously.'; return; end -% Hexa mesh: detect whether the geometry was adapted -if strcmpi(ElementType, 'hexahedron') - GeometryAdapted = []; - % Detect using the options of the Brainstorm process that created the file - if isfield(FemMat, 'History') && ~isempty(FemMat.History) && ~isempty(strfind([FemMat.History{:,3}], 'NodeShift')) - strOptions = [FemMat.History{:,3}]; - iTag = strfind(strOptions, 'NodeShift'); - val = sscanf(strOptions(iTag:end), 'NodeShift=%f'); - if ~isempty(val) - GeometryAdapted = (val > 0); - end - end - % Otherwise, try to guess based on the geometry - if isempty(GeometryAdapted) - % Compute the distance between the first two nodes of each element - dist = sqrt(sum([FemMat.Vertices(FemMat.Elements(:,1),1) - FemMat.Vertices(FemMat.Elements(:,2),1), ... - FemMat.Vertices(FemMat.Elements(:,2),2) - FemMat.Vertices(FemMat.Elements(:,2),2), ... - FemMat.Vertices(FemMat.Elements(:,2),3) - FemMat.Vertices(FemMat.Elements(:,2),3)] .^ 2, 2)); - % If the distance is not constant: then the geomtry is adapted - GeometryAdapted = (max(abs(dist - dist(1))) > 1e-9); - end - % Copy value in DUNEuro options - if ~isempty(GeometryAdapted) - cfg.GeometryAdapted = GeometryAdapted; - disp(['DUNEURO> Detected parameter: GeometryAdapted=', bool2str(GeometryAdapted)]); - end -end -% % Isotropic -% if (cfg.Isotropic) -% % Isotropic without tensor -% if ~cfg.UseTensor -% if strcmp(ElementType,'hexahedron') -% MeshFile = 'head_model.dgf'; -% else -% MeshFile = 'head_model.msh'; -% end -% % Isotropic with tensor -% else -% if strcmp(ElementType,'hexahedron') -% errMsg = 'Using the tensor model with hexahedral mesh is not supported for now.'; -% return; -% end -% MeshFile = 'head_model.geo'; -% end -% % Anisotropic (with tensor) -% else -% if strcmp(ElementType,'hexahedron') -% errMsg = 'Using the anisotropy model with hexahedral mesh is not supported for now.'; -% return; -% end -% MeshFile = 'head_model.geo'; -% cfg.UseTensor = true; -% 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]); @@ -306,75 +249,13 @@ end end % view_surface_matrix(cfg.GridLoc, sCortex.Faces) - -% %%%% ============================================= -% % Now similar process for the WM with an extra and unexpedted step ... -% if ~isempty(iWM) -% disp('Checking the the WM ...'); -% %Extract WM surface -% wm_tetra = FemMat.Elements(FemMat.Tissue <= iWM,:); -% wm_face = volface(wm_tetra); -% [nwmf, ewmf] = removeisolatednode(FemMat.Vertices,wm_face); -% % check if any dipoles is inside the WM -% wMfv.vertices = nwmf; -% wMfv.faces = ewmf; -% tic -% wMin = inpolyhedron(wMfv, sCortex.Vertices); -% wMindex_in = find(wMin); -% disp(['There are ' num2str(sum(wMin)) ' dipoles inside the WM']); -% disp('Moving these dipoles to the GM tissues ...'); -% twm = toc; -% if ~isempty(wMindex_in) -% % 1- move the dipole from inside the WM to the GM surface -% % ==> this is for testing, when we use the centroide -% % directely, some dipole remains within the WM ... -% GMcentroide = 0; % just for testing, to use directely the GM centroides -% if GMcentroide == 1 -% k = dsearchn(elem_centroide,sCortex.Vertices(wMindex_in,:)); -% NewVertices(wMindex_in ,:) = ElemCenter(k,:); -% wMoutFinal = inpolyhedron(wMfv, NewVertices); -% disp(['Now, there are ' num2str(sum(wMoutFinal)) ' dipoles inside the WM']); -% else -% k = dsearchn(gMfv.vertices,sCortex.Vertices(wMindex_in,:)); -% NewVertices(wMindex_in ,:) = gMfv.vertices(k,:); -% wMoutFinal = inpolyhedron(wMfv, NewVertices); -% disp(['There are ' num2str(sum(wMoutFinal)) ' dipoles inside the WM']); -% disp(['These dipoles are moved to the GM outer surface nodes']); -% disp(['Moving these dipoles to the nearest centroide of the GM element...']); -% % Now, we move these same dipole from the surface to the -% % nearest centroide ==> to be sure it's inside the GM -% % ==> fill partially the FEM condition -% % move again to the centoide -% k = dsearchn(ElemCenter,NewVertices(wMindex_in,:)); -% NewVertices(wMindex_in ,:) = ElemCenter(k,:); -% % just to check -% wMoutFinal = inpolyhedron(wMfv, NewVertices); -% disp(['Now, there are ' num2str(sum(wMoutFinal)) ' dipoles inside the WM']); -% disp(['All the dipoles are moved to nearest centroid of the GM']); -% end -% end -% end - end case 'mixed' % TODO : not used ? end - %% ===== SOURCE MODEL ===== bst_progress('text', 'DUNEuro: Writing temporary files...'); -% Write the source/dipole file -% DipoleFile = fullfile(TmpDir, 'dipole_model.txt'); -% Unconstrained orientation for each dipole. ie the output file have 3*N dipoles (3 directions for each) -% [x1,y1,z1, 1, 0, 0,; -% x1,y1,z1, 0, 1, 0; -% x1,y1,z1, 0, 0, 1; -% x2,y2,z2, 1, 0, 0; -% .... ]; -% dipoles = [kron(cfg.GridLoc, ones(3,1)), kron(ones(size(cfg.GridLoc,1), 1), eye(3))]; -% fid = fopen(DipoleFile, 'wt+'); -% fprintf(fid, '%d %d %d %d %d %d \n', dipoles'); -% fclose(fid); [iOk, errMsg] = bstdn_write_source_space(TmpDir, cfg.GridLoc); %% ===== SENSOR MODEL ===== @@ -393,52 +274,11 @@ end %% ===== CONDUCTIVITY MODEL ===== -% % Isotropic without tensor -% if ~cfg.UseTensor -% CondFile = fullfile(TmpDir, 'conductivity_model.con'); -% fid = fopen(CondFile, 'w'); -% fprintf(fid, '%d\t', cfg.FemCond); -% fclose(fid); -% % With tensor (isotropic or anisotropic) -% else -% CondFile = fullfile(TmpDir, 'conductivity_model.knw'); -% % Transformation matrix and tensor mapping on each direction -% CondTensor = zeros(length(FemMat.Elements),6) ; -% for ind =1 : length(FemMat.Elements) -% temp0 = reshape(FemMat.Tensors(ind,:),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 -% CondTensor(ind,:) = [temp(1) temp(5) temp(9) temp(4) temp(8) temp(7)]; % this is the right order -% end -% % write the tensors -% out_fem_knw(FemMat, CondTensor, CondFile); -% end -% % Write mesh model -% MeshFile = fullfile(TmpDir, MeshFile); -% out_fem(FemMat, MeshFile); -%% DN2026 -% export volume conductor model -% write_volume_conductor(TmpDir, nodes, elements, labels, tensors); -% write hdf5 files : -% online viewer of hdf5 files : https://myhdf5.hdfgroup.org/ dnbst_write_volume_conductor(TmpDir, FemMat, cfg.FemCond); +% NOTE: online viewer of hdf5 files : https://myhdf5.hdfgroup.org/ - -% run DUNEuro container -% Source Space configuration -sourcespace_config = []; -sourcespace_config.name = 'create_sourcespace'; -% sourcespace_config.sourcecompartment = '1'; -% sourcespace_config.gridsize = '2.0'; -% sourcespace_config.enforceDistanceCondition = 'True'; -% sourcespace_config.distanceThreshold = '0.5'; -sourcespace_config.enforceVenantCondition = 'True'; -% sourcespace_config.visualizeSourceSpace = 'True'; - - -% Transfer Matrix configuration -transfer_matrix_config = []; +%% ===== TRANSFER MATRIX CONFIGURATION ===== + transfer_matrix_config = []; transfer_matrix_config.name = 'compute_transfer_matrix'; if strcmp(dnModality, 'eeg') transfer_matrix_config.do_eeg = 'True'; @@ -453,47 +293,43 @@ transfer_matrix_config.do_eeg = 'True'; end transfer_matrix_config.residual_reduction = '1e-16'; -transfer_matrix_config.nr_threads = '-1'; +transfer_matrix_config.nr_threads = '-1'; % can be used as user parameters +% Check with Malte is there is an optimised number without overwhelming the +% user computer. -% Leadfield configuration +%% ===== LEADFIELD MATRIX CONFIGURATION ===== leadfield_config = []; leadfield_config.name = 'compute_leadfield'; if strcmp(dnModality, 'eeg') leadfield_config.do_meg = 'False'; leadfield_config.do_eeg = 'True'; - leadfield_config.eeg_scaling = '1e0'; - leadfield_config.meg_scaling = '1e5'; end if strcmp(dnModality, 'meg') leadfield_config.do_meg = 'True'; leadfield_config.do_eeg = 'False'; - leadfield_config.eeg_scaling = '1e0'; - leadfield_config.meg_scaling = '1e5'; end if strcmp(dnModality, 'meeg') leadfield_config.do_meg = 'True'; leadfield_config.do_eeg = 'True'; - leadfield_config.eeg_scaling = '1e0'; - leadfield_config.meg_scaling = '1e5'; end +leadfield_config.eeg_scaling = '1e0'; % check with Malte if those value are optimised +leadfield_config.meg_scaling = '1e5'; + leadfield_config.sourcemodel = cfg.SrcModel2026; % [select from the interface: 'multipolar_venant', 'local_subtraction', 'partial_integration'] -leadfield_config.nr_threads = '-1'; +leadfield_config.nr_threads = '-1'; % same as above %% ===== RUN DUNEURO ====== -% Assemble command line -% callStr = ['"' DuneuroExe '"' ' ' '"' IniFile '"']; bst_progress('text', 'DUNEuro: Computing leadfield...'); -runner = 'docker'; % disp(['DUNEURO> System call: ' callStr]); tic; -% Call DUNEuro -if 0 - [status, errMsg ] = bst_run_duneuro_task(TmpDir, sourcespace_config, runner); % wheew is the outputs? -end -% Compute transfer matrix +% Call DUNEuro and Compute transfer matrix [status, errMsg ] = bst_run_duneuro_task(TmpDir, transfer_matrix_config, runner); -% Compute the Leadfield +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.'; @@ -503,11 +339,7 @@ %% ===== READ LEADFIELD ====== bst_progress('text', 'DUNEuro: Reading leadfield...'); -% % EEG -% if (isEeg || isEcog || isSeeg) -% GainEeg = in_duneuro_bin(fullfile(TmpDir, cfg.BstEegLfFile))'; -% end -% read data +% 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) @@ -526,238 +358,24 @@ end 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 -% Delete the temporary files -% file_delete(TmpDir, 1, 1); - - -%% ===== SAVE TRANSFER MATRIX ====== -disp('DUNEURO> TODO: Save transferOut.dat to database.') - % Remove logo bst_plugin('SetProgressLogo', []); - - -% %% ===== WRITE MINI FILE ===== -% % Open the mini file -% IniFile = fullfile(TmpDir, 'duneuro_minifile.mini'); -% fid = fopen(IniFile, 'wt+'); -% % General setting -% fprintf(fid, '__name = %s\n\n', IniFile); -% if strcmp(cfg.SolverType, 'cg') -% fprintf(fid, 'type = %s\n', cfg.FemType); -% end -% fprintf(fid, 'element_type = %s\n', ElementType); -% fprintf(fid, 'solver_type = %s\n', cfg.SolverType); -% fprintf(fid, 'geometry_adapted = %s\n', bool2str(cfg.GeometryAdapted)); -% fprintf(fid, 'tolerance = %d\n', cfg.Tolerance); -% % [electrodes] -% if isEcog || isSeeg -% % Instead of selecting the electrode on the outer surface, -% % uses the nearest FEM node as the electrode location -% cfg.ElecType = 'closest_subentity_center'; -% end -% if strcmp(dnModality, 'eeg') || strcmp(dnModality, 'meeg') -% fprintf(fid, '[electrodes]\n'); -% fprintf(fid, 'filename = %s\n', fullfile(TmpDir, ElecFile)); -% fprintf(fid, 'type = %s\n', cfg.ElecType); -% fprintf(fid, 'codims = %s\n', '3'); -% end -% % [meg] -% if strcmp(dnModality, 'meg') || strcmp(dnModality, 'meeg') -% fprintf(fid, '[meg]\n'); -% fprintf(fid, 'intorderadd = %d\n', cfg.MegIntorderadd); -% fprintf(fid, 'type = %s\n', cfg.MegType); -% fprintf(fid, 'cache.enable = %s\n',bool2str(cfg.EnableCacheMemory) ); -% % [coils] -% fprintf(fid, '[coils]\n'); -% fprintf(fid, 'filename = %s\n', CoilFile); -% % [projections] -% fprintf(fid, '[projections]\n'); -% fprintf(fid, 'filename = %s\n', ProjFile); -% end -% % [dipoles] -% fprintf(fid, '[dipoles]\n'); -% fprintf(fid, 'filename = %s\n', DipoleFile); -% % [volume_conductor.grid] -% fprintf(fid, '[volume_conductor.grid]\n'); -% fprintf(fid, 'filename = %s\n', MeshFile); -% % [volume_conductor.tensors] -% fprintf(fid, '[volume_conductor.tensors]\n'); -% fprintf(fid, 'filename = %s\n', CondFile); -% % [solver] -% fprintf(fid, '[solver]\n'); -% fprintf(fid, 'solver_type = %s\n', cfg.SolvSolverType); -% fprintf(fid, 'preconditioner_type = %s\n', cfg.SolvPrecond); -% if strcmp(cfg.SolverType, 'cg') -% fprintf(fid, 'cg_smoother_type = %s\n', cfg.SolvSmootherType); -% end -% fprintf(fid, 'intorderadd = %d\n', cfg.SolvIntorderadd); -% % Discontinuous Galerkin -% if strcmp(cfg.SolverType, 'dg') -% fprintf(fid, 'dg_smoother_type = %s\n', cfg.DgSmootherType); -% fprintf(fid, 'scheme = %s\n', cfg.DgScheme); -% fprintf(fid, 'penalty = %d\n', cfg.DgPenalty); -% fprintf(fid, 'edge_norm_type = %s\n', cfg.DgEdgeNormType); -% fprintf(fid, 'weights = %s\n', bool2str(cfg.DgWeights)); -% fprintf(fid, 'reduction = %s\n', bool2str(cfg.DgReduction)); -% end -% % [solution] -% fprintf(fid, '[solution]\n'); -% fprintf(fid, 'post_process = %s\n', bool2str(cfg.SolPostProcess)); % true/false -% fprintf(fid, 'subtract_mean = %s\n', bool2str(cfg.SolSubstractMean)); % boolean -% % [solution.solver] -% fprintf(fid, '[solution.solver]\n'); -% fprintf(fid, 'reduction = %d\n', cfg.SolSolverReduction); -% % [solution.source_model] -% fprintf(fid, '[solution.source_model]\n'); -% fprintf(fid, 'type = %s\n', cfg.SrcModel); -% fprintf(fid, 'intorderadd = %d\n', cfg.SrcIntorderadd); -% fprintf(fid, 'intorderadd_lb = %d\n', cfg.SrcIntorderadd_lb); -% fprintf(fid, 'numberOfMoments = %d\n', cfg.SrcNbMoments); -% fprintf(fid, 'referenceLength = %d\n', cfg.SrcRefLen); -% fprintf(fid, 'weightingExponent = %d\n', cfg.SrcWeightExp); -% fprintf(fid, 'relaxationFactor = %e\n', 10^(-cfg.SrcRelaxFactor)); -% fprintf(fid, 'mixedMoments = %s\n', bool2str(cfg.SrcMixedMoments)); -% fprintf(fid, 'restrict = %s\n', bool2str(cfg.SrcRestrict)); -% fprintf(fid, 'initialization = %s\n', cfg.SrcInit); -% % [brainstorm] -% fprintf(fid, '[brainstorm]\n'); -% fprintf(fid, 'modality = %s\n', dnModality); -% fprintf(fid, 'output_folder = %s\n', [TmpDir, filesep]); -% fprintf(fid, 'save_eeg_transfer_file = %s\n', bool2str(cfg.BstSaveTransfer)); -% fprintf(fid, 'save_meg_transfer_file = %s\n', bool2str(cfg.BstSaveTransfer)); -% fprintf(fid, 'save_meeg_transfer_file = %s\n', bool2str(cfg.BstSaveTransfer)); -% fprintf(fid, 'eeg_transfer_filename = %s\n', cfg.BstEegTransferFile); -% fprintf(fid, 'meg_transfer_filename = %s\n', cfg.BstMegTransferFile); -% fprintf(fid, 'eeg_leadfield_filename = %s\n', cfg.BstEegLfFile); -% fprintf(fid, 'meg_leadfield_filename = %s\n', cfg.BstMegLfFile); -% % Close file -% fclose(fid); - - -% %% ===== RUN DUNEURO ====== -% % Assemble command line -% callStr = ['"' DuneuroExe '"' ' ' '"' IniFile '"']; -% bst_progress('text', 'DUNEuro: Computing leadfield...'); -% disp(['DUNEURO> System call: ' callStr]); -% tic; -% % Call DUNEuro -% status = system(callStr) -% 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...'); -% % EEG -% if (isEeg || isEcog || isSeeg) -% GainEeg = in_duneuro_bin(fullfile(TmpDir, cfg.BstEegLfFile))'; -% end -% -% %MEG -% To check with Malte: is the MEG leadfield is the finla or just the -% % seconday? -% if isMeg -% GainMeg = in_duneuro_bin(fullfile(TmpDir, cfg.BstMegLfFile))'; -% % === POST-PROCESS MEG LEADFIELD === -% % Compute the total magnetic field -% dipoles_pos_orie = [kron(cfg.GridLoc,ones(3,1)), kron(ones(length(cfg.GridLoc),1), eye(3))]; -% -% % a- Compute the MEG Primary Magnetic B-field analytically (formula of Sarvas) -% dip_pos = dipoles_pos_orie(:,1:3); -% dip_mom = dipoles_pos_orie(:,4:6); -% Bp = zeros(size(MegChannels,1), size(dip_pos,1)); -% for i = 1:size(CoilsLoc,1) -% for j = 1 : size(dip_pos,1) -% R = CoilsLoc(i,:); -% R_0 = dip_pos(j,:); -% A = R - R_0; -% a = norm(A); -% aa = A./(a^3); -% BpHelp = cross(dip_mom(j,:),aa); -% Bp(i,j) = BpHelp * CoilsOrient(i, :)'; % projection of the primary B-field along the coil orientations -% end -% end -% -% % b- The total magnetic field B = Bp + Bs; -% % full B-field -% Bs = GainMeg; -% mu = 4*pi*1e-4; % check the value of the units maybe it needs to be mu = 4*pi*1e-7 -% Bfull = (mu/(4*pi)) * (Bp - Bs); -% -% % c- Apply the weight : -% [channelIndex] = unique(MegChannels(:,1)); -% nbChannel = length(channelIndex); -% weighted_B = zeros(nbChannel,size(Bfull,2)); -% for iCh = 1 : nbChannel -% communChannel = find(iCh==MegChannels(:,1)); -% BcommunChannel = Bfull(communChannel(:),:); -% WcommunChannel = MegChannels(communChannel(:), 8: end); -% weighted_B(iCh,:) = sum (BcommunChannel.*WcommunChannel,1); -% end -% GainMeg = weighted_B; -% end - -% % Fill the unused channels with NaN -% Gain = NaN * zeros(length(cfg.Channel), 3 * length(cfg.GridLoc)); -% if isMeg -% Gain(cfg.iMeg,:) = GainMeg; -% % scaling the MEG Gain matrix -% Gain(cfg.iMeg,:) = GainMeg/1000; -% end -% if (isEeg || isEcog || isSeeg) -% Gain(cfg.iEeg,:) = GainEeg; -% end -% -% % Delete the temporary files -% file_delete(TmpDir, 1, 1); -% -% -% %% ===== SAVE TRANSFER MATRIX ====== -% disp('DUNEURO> TODO: Save transferOut.dat to database.') -% -% % Remove logo -% bst_plugin('SetProgressLogo', []); - end - - - %% ================================================================================= % === SUPPORT FUNCTIONS ========================================================= % ================================================================================= - -%% ===== BOOL => STR ===== -function str = bool2str(bool) - if bool - str = 'true'; - else - str = 'false'; - end -end - - function [iOk, errMsg] = dnbst_write_volume_conductor(duneuro_io_dir, FemMat, isoCond) iOk = 0; errMsg = ''; - nodes = FemMat.Vertices; elements = FemMat.Elements - 1; labels = FemMat.Tissue - 1; - io_file_path = fullfile(duneuro_io_dir, 'duneuro_io.hdf5'); - - % % we assume that this function create the IO file. We abort the execution - % % if the file already exists - % if isfile(io_file_path) - % errMsg = ['There already exists a file at the location ', io_file_path]; - % end - + % get information about the input data nr_nodes = size(nodes, 1); nr_elements = size(elements, 1); @@ -766,15 +384,12 @@ 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 - + 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 @@ -814,11 +429,6 @@ io_file_path = fullfile(duneuro_io_dir, 'duneuro_io.hdf5'); - % % we assume that the output file already exists - % if ~isfile(io_file_path) - % error('The DUNEuro IO file does not exist, please call write_volume_conductor first'); - % end - nr_magnetometers = size(coil_positions, 1); nr_channels = size(coil_to_channel_transform, 1); @@ -910,49 +520,48 @@ 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']; -else - errMsg = 'unknown runner'; - return; -end - -status = system(runner_system_call) + 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']; + else + errMsg = 'unknown runner'; + return; + end + + status = system(runner_system_call); end \ No newline at end of file From 8edc89498ea408fb60d38182d6d0d20109cfc3a4 Mon Sep 17 00:00:00 2001 From: tmedani Date: Sat, 7 Feb 2026 18:38:33 -0800 Subject: [PATCH 03/11] Refactor DUNEuro2026 panel and forward code Rename and clean up DUNEuro2026 integration: - Switch identifiers from 'duneuro_2026' to 'duneuro2026' in UI, - Update panel name/comments, - Streamline bst_duneuro_2026 - Minor UI string tweaks and other small whitespace/comment cleanups. --- toolbox/forward/bst_duneuro_2026.m | 37 ++------ toolbox/forward/panel_duneuro2026.m | 132 ++-------------------------- toolbox/forward/panel_headmodel.m | 14 +-- 3 files changed, 24 insertions(+), 159 deletions(-) diff --git a/toolbox/forward/bst_duneuro_2026.m b/toolbox/forward/bst_duneuro_2026.m index 0cb7f27942..e0eb10bb19 100644 --- a/toolbox/forward/bst_duneuro_2026.m +++ b/toolbox/forward/bst_duneuro_2026.m @@ -102,10 +102,6 @@ %% ===== HEAD MODEL ===== % Load FEM mesh FemMat = load(file_fullpath(cfg.FemFile)); -% Get mesh type -switch size(FemMat.Elements,2) - case 4, ElementType = 'tetrahedron'; -end % 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 @@ -129,6 +125,7 @@ % 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...'); @@ -256,7 +253,7 @@ %% ===== SOURCE MODEL ===== bst_progress('text', 'DUNEuro: Writing temporary files...'); -[iOk, errMsg] = bstdn_write_source_space(TmpDir, cfg.GridLoc); +bstdn_write_source_space(TmpDir, cfg.GridLoc); %% ===== SENSOR MODEL ===== % Write the EEG electrode file @@ -266,8 +263,6 @@ % Write the MEG sensors data if isMeg % Write coil data - % ********* TODO *********** - % Get the correct weigh from the original channel file coil_to_channel_transform = eye(length(MegChannels)); coil_to_channel_transform = coil_to_channel_transform.* MegChannels(:,end); dnbst_write_magnetometers(TmpDir, MegChannels(:,2:4), MegChannels(:,5:7), coil_to_channel_transform); @@ -317,7 +312,6 @@ leadfield_config.sourcemodel = cfg.SrcModel2026; % [select from the interface: 'multipolar_venant', 'local_subtraction', 'partial_integration'] leadfield_config.nr_threads = '-1'; % same as above - %% ===== RUN DUNEURO ====== bst_progress('text', 'DUNEuro: Computing leadfield...'); % disp(['DUNEURO> System call: ' callStr]); @@ -351,17 +345,16 @@ % Fill the unused channels with NaN Gain = NaN * zeros(length(cfg.Channel), 3 * length(cfg.GridLoc)); -if isMeg - Gain(cfg.iMeg,:) = GainMeg; - % scaling the MEG Gain matrix - Gain(cfg.iMeg,:) = GainMeg/1000; -end 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 @@ -466,12 +459,7 @@ iOk = 0; errMsg = ''; io_file_path = fullfile(duneuro_io_dir, 'duneuro_io.hdf5'); - - % % we assume that the output file already exists - % if ~isfile(io_file_path) - % error('The DUNEuro IO file does not exist, please call write_volume_conductor first'); - % end - + nr_electrodes = size(electrode_positions, 1); dim = 3; @@ -505,18 +493,11 @@ if size(dipole_positions, 2) ~= dim errMsg = 'The dipole position array must have 3 columns'; return; end - % if size(electrode_positions, 2) ~= dim - % errMsg = 'The dipole orientation array must have 3 columns'; - % end - % if strcmp(electrode_type_flag, 'measurement') h5create(io_file_path, "/measurement/source_space/positions", [dim nr_dipoles], Datatype="double"); h5write(io_file_path, "/measurement/source_space/positions", dipole_positions'); - - % else - % errMsg = 'Dipole flag needs to be either "measurement"'; return; - % end - iOk = 1; + + iOk = 1; end diff --git a/toolbox/forward/panel_duneuro2026.m b/toolbox/forward/panel_duneuro2026.m index 384f716819..fc3ffed8c2 100644 --- a/toolbox/forward/panel_duneuro2026.m +++ b/toolbox/forward/panel_duneuro2026.m @@ -1,9 +1,9 @@ function varargout = panel_duneuro2026(varargin) -% PANEL_DUNEURO: DUNEuro options +% PANEL_DUNEURO2026: DUNEuro2026 options % -% USAGE: bstPanel = panel_duneuro('CreatePanel', OPTIONS) : Call from the interactive interface -% bstPanel = panel_duneuro('CreatePanel', sProcess, sFiles) : Call from the process editor -% s = panel_duneuro('GetPanelContents') +% 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: @@ -24,6 +24,7 @@ % =============================================================================@ % % Authors: Francois Tadel, 2020 +% Takfarinas Medani, 2026 : adapted from panel_duneuro eval(macro_method); end @@ -40,8 +41,8 @@ if (nargin == 1) OPTIONS = sProcess; % Check if there is only MEG, for simplified model by default - isMegOnly = ~ismember('duneuro_2026', {OPTIONS.EEGMethod, OPTIONS.ECOGMethod, OPTIONS.SEEGMethod}); - isMeg = isequal(OPTIONS.MEGMethod, 'duneuro_2026'); + 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; @@ -125,30 +126,6 @@ end c.gridy = 1; jPanelLeft.add(jPanelLayers, c); - -% % ==== PANEL LEFT: FEM METHOD TYPE ==== -% jPanelType = gui_river([1,1], [0,6,6,6], 'FEM method type'); -% jGroupFemType = ButtonGroup(); -% jRadioFemTypeFit = gui_component('radio', jPanelType, 'br', 'Fitted', jGroupFemType, '', [], []); -% jRadioFemTypeUnfit = gui_component('radio', jPanelType, 'br', 'Unfitted', jGroupFemType, '', [], []); -% switch lower(OPTIONS.FemType) -% case 'fitted', jRadioFemTypeFit.setSelected(1); -% case 'unfitted', jRadioFemTypeUnfit.setSelected(1); -% end -% c.gridy = 2; -% jPanelLeft.add(jPanelType, c); - - % ==== PANEL LEFT: FEM SOLVER TYPE ==== -% jPanelSolverType = gui_river([1,1], [0,6,6,6], 'FEM solver type'); -% jGroupSolverType = ButtonGroup(); -% jRadioSolverTypeCg = gui_component('radio', jPanelSolverType, 'br', 'CG: Continuous Galerkin', jGroupSolverType, '', [], []); -% jRadioSolverTypeDg = gui_component('radio', jPanelSolverType, 'br', 'DG: Discontinuous Galerkin', jGroupSolverType, '', [], []); -% switch lower(OPTIONS.SolverType) -% case 'cg', jRadioSolverTypeCg.setSelected(1); -% case 'dg', jRadioSolverTypeDg.setSelected(1); -% end -% c.gridy = 2; -% jPanelLeft.add(jPanelSolverType, c); % ==== PANEL RIGHT: FEM SOURCE MODEL ==== jPanelSrcModel = gui_river([1,1], [0,6,6,6], 'FEM source model'); @@ -164,50 +141,6 @@ c.gridy = 3; jPanelLeft.add(jPanelSrcModel, c); - % ==== PANEL RIGHT: VENANT OPTIONS ==== -% jPanelOptVen = gui_river([3,3], [0,6,6,6], 'Venant options'); -% % Number of moments -% gui_component('label', jPanelOptVen, [], 'Number of moments (1-5): ', [], '', [], []); -% jTextNbMoments = gui_component('texttime', jPanelOptVen, 'tab', '', [], '', [], []); -% gui_validate_text(jTextNbMoments, [], [], 1:5, '', 0, OPTIONS.SrcNbMoments, []); -% % Reference length -% gui_component('label', jPanelOptVen, 'br', 'Reference length (1-100): ', [], '', [], []); -% jTextRefLen = gui_component('texttime', jPanelOptVen, 'tab', '', [], '', [], []); -% gui_validate_text(jTextRefLen, [], [], 1:100, '', 0, OPTIONS.SrcRefLen, []); -% % Weighting exponent -% gui_component('label', jPanelOptVen, 'br', 'Weighting exponent (1-3): ', [], '', [], []); -% jTextWeightExp = gui_component('texttime', jPanelOptVen, 'tab', '', [], '', [], []); -% gui_validate_text(jTextWeightExp, [], [], 1:3, '', 0, OPTIONS.SrcWeightExp, []); -% % Relaxation Factor -% gui_component('label', jPanelOptVen, 'br', 'Relaxation factor exponent (3-9): ', [], '', [], []); -% jTextRelaxFactor = gui_component('texttime', jPanelOptVen, 'tab', sprintf('%e', OPTIONS.SrcRelaxFactor), [], '', [], []); -% gui_validate_text(jTextRelaxFactor, [], [], 3:9, '', 0, OPTIONS.SrcRelaxFactor, []); -% % Mixed moments -% jCheckMixedMoments = gui_component('checkbox', jPanelOptVen, 'br', 'Mixed moments', [], '', [], []); -% if (OPTIONS.SrcMixedMoments == 1) -% jCheckMixedMoments.setSelected(1); -% end -% % Restrict -% jCheckRestrict = gui_component('checkbox', jPanelOptVen, 'br', 'Restrict', [], '', [], []); -% if (OPTIONS.SrcRestrict == 1) -% jCheckRestrict.setSelected(1); -% end -% c.gridy = 1; -% jPanelRight.add(jPanelOptVen, c); - - % ==== PANEL RIGHT: SUBTRACTION OPTIONS ==== -% jPanelOptSub = gui_river([3,3], [0,6,6,6], 'Subtraction options'); -% % Number of moments -% gui_component('label', jPanelOptSub, [], 'intorderadd (1-5): ', [], '', [], []); -% jTextIntorderadd = gui_component('texttime', jPanelOptSub, 'tab', '', [], '', [], []); -% gui_validate_text(jTextIntorderadd, [], [], 0:5, '', 0, OPTIONS.SrcIntorderadd, []); -% % Number of moments -% gui_component('label', jPanelOptSub, 'br', 'intorderadd_lb (1-5): ', [], '', [], []); -% jTextIntorderadd_lb = gui_component('texttime', jPanelOptSub, 'tab', '', [], '', [], []); -% gui_validate_text(jTextIntorderadd_lb, [], [], 0:5, '', 0, OPTIONS.SrcIntorderadd_lb, []); -% c.gridy = 2; -% jPanelRight.add(jPanelOptSub, c); - % ==== PANEL RIGHT: INPUT OPTIONS ==== jPanelInput = gui_river([1,1], [0,6,6,6], 'Source space'); % Shrink source space @@ -293,21 +226,9 @@ % Create the BstPanel object that is returned by the function ctrl = struct('jCheckLayer', jCheckLayer, ... 'jTextCond', jTextCond, ... - ... 'jRadioFemTypeFit', jRadioFemTypeFit, ... - ... 'jRadioFemTypeUnfit', jRadioFemTypeUnfit, ... - ... 'jRadioSolverTypeCg', jRadioSolverTypeCg, ... - ... 'jRadioSolverTypeDg', jRadioSolverTypeDg, ... 'jRadioSrcModelVen', jRadioSrcModelVen, ... 'jRadioSrcModelSub', jRadioSrcModelSub, ... 'jRadioSrcModelPar', jRadioSrcModelPar, ... - ... 'jTextNbMoments', jTextNbMoments, ... - ... 'jTextRefLen', jTextRefLen, ... - ... 'jTextWeightExp', jTextWeightExp, ... - ... 'jTextRelaxFactor', jTextRelaxFactor, ... - ... 'jCheckMixedMoments', jCheckMixedMoments, ... - ... 'jCheckRestrict', jCheckRestrict, ... - ... 'jTextIntorderadd', jTextIntorderadd, ... - ... 'jTextIntorderadd_lb', jTextIntorderadd_lb, ... 'jTextSrcShrink', jTextSrcShrink, ... 'jCheckSrcForceInGM', jCheckSrcForceInGM, ... 'jCheckSaveTransfer', jCheckSaveTransfer, ... @@ -320,8 +241,6 @@ bstPanelNew = BstPanel(panelName, jPanelNew, ctrl); % Update comments UpdatePanel(1); - - %% ================================================================================= % === LOCAL CALLBACKS ============================================================ @@ -359,11 +278,7 @@ function UpdatePanel(isForced) ExpertMode = bst_get('ExpertMode'); % Show/hide panels jPanelRight.setVisible(ExpertMode); -% jPanelType.setVisible(ExpertMode); -% jPanelSolverType.setVisible(ExpertMode); jPanelSrcModel.setVisible(ExpertMode); -% jPanelOptVen.setVisible(ExpertMode && jRadioSrcModelVen.isSelected()); -% jPanelOptSub.setVisible(ExpertMode && jRadioSrcModelSub.isSelected()); jPanelInput.setVisible(ExpertMode); jPanelOutput.setVisible(ExpertMode); jPanelMegComputationOption.setVisible(ExpertMode); @@ -414,18 +329,6 @@ function UpdatePanel(isForced) end end s.UseTensor = ctrl.UseTensor; -% % FEM method type -% if ctrl.jRadioFemTypeFit.isSelected() -% s.FemType = 'fitted'; -% elseif ctrl.jRadioFemTypeUnfit.isSelected() -% s.FemType = 'unfitted'; -% end - % FEM solver type - % if ctrl.jRadioSolverTypeCg.isSelected() - % s.SolverType = 'cg'; - % elseif ctrl.jRadioSolverTypeDg.isSelected() - % s.SolverType = 'dg'; - % end % Source model if ctrl.jRadioSrcModelVen.isSelected() s.SrcModel = 'multipolar_venant'; @@ -434,19 +337,6 @@ function UpdatePanel(isForced) elseif ctrl.jRadioSrcModelPar.isSelected() s.SrcModel = 'partial_integration'; end - % % Venant options - % if strcmpi(s.SrcModel, 'venant') - % s.SrcNbMoments = str2double(ctrl.jTextNbMoments.getText()); - % s.SrcRefLen = str2double(ctrl.jTextRefLen.getText()); - % s.SrcWeightExp = str2double(ctrl.jTextWeightExp.getText()); - % s.SrcRelaxFactor = str2double(ctrl.jTextRelaxFactor.getText()); - % s.SrcMixedMoments = ctrl.jCheckMixedMoments.isSelected(); - % s.SrcRestrict = ctrl.jCheckRestrict.isSelected(); - % % Subtraction options - % elseif strcmpi(s.SrcModel, 'subtraction') - % s.SrcIntorderadd = str2double(ctrl.jTextIntorderadd.getText()); - % s.SrcIntorderadd_lb = str2double(ctrl.jTextIntorderadd_lb.getText()); - % end % Input options s.SrcShrink = str2double(ctrl.jTextSrcShrink.getText()); if ~isempty(ctrl.jCheckSrcForceInGM) @@ -468,12 +358,6 @@ function UpdatePanel(isForced) else s.EnableCacheMemory = 0; end - -% if ~isempty(ctrl.jCheckMegPerBlockOfSensor) -% s.MegPerBlockOfSensor = ctrl.jCheckMegPerBlockOfSensor.isSelected(); -% else -% s.MegPerBlockOfSensor = 0; -% end end @@ -530,4 +414,4 @@ function UpdatePanel(isForced) isType = ~isempty(strfind(strName, 'skin')) || ~isempty(strfind(strName, 'scalp')) || ~isempty(strfind(strName, 'head')); end end -end +end \ No newline at end of file diff --git a/toolbox/forward/panel_headmodel.m b/toolbox/forward/panel_headmodel.m index 11d3f450d7..750d03e9a0 100644 --- a/toolbox/forward/panel_headmodel.m +++ b/toolbox/forward/panel_headmodel.m @@ -80,7 +80,7 @@ jComboMethodMEG.addItem(BstListItem('os_meg', '', 'Overlapping spheres', [])); jComboMethodMEG.addItem(BstListItem('openmeeg', '', 'OpenMEEG BEM', [])); jComboMethodMEG.addItem(BstListItem('duneuro', '', 'DUNEuro FEM', [])); - jComboMethodMEG.addItem(BstListItem('duneuro_2026', '', 'DUNEuro FEM 2026', [])); + jComboMethodMEG.addItem(BstListItem('duneuro2026', '', 'DUNEuro FEM 2026', [])); jComboMethodMEG.setSelectedIndex(1); else jCheckMethodMEG = []; @@ -96,7 +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('duneuro_2026', '', 'DUNEuro FEM2026', [])); + jComboMethodEEG.addItem(BstListItem('duneuro2026', '', 'DUNEuro FEM 2026', [])); jComboMethodEEG.setSelectedIndex(1); else jCheckMethodEEG = []; @@ -111,7 +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('duneuro_2026', '', 'DUNEuro FEM2026', [])); + jComboMethodECOG.addItem(BstListItem('duneuro2026', '', 'DUNEuro FEM 2026', [])); jComboMethodECOG.setSelectedIndex(0); else jCheckMethodECOG = []; @@ -126,7 +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('duneuro_2026', '', 'DUNEuro FEM2026', [])); + jComboMethodSEEG.addItem(BstListItem('duneuro2026', '', 'DUNEuro FEM 2026', [])); jComboMethodSEEG.setSelectedIndex(0); else jCheckMethodSEEG = []; @@ -443,7 +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, 'duneuro_2026', 'DUNEuro FEM2026'); + sMethod.Comment = strrep(sMethod.Comment, 'duneuro2026', 'DUNEuro FEM 2026'); % Grid type if strcmpi(sMethod.HeadModelType, 'volume') sMethod.Comment = [sMethod.Comment ' (volume)']; @@ -453,7 +453,7 @@ function UpdateComment(varargin) end isOpenMEEG = any(strcmpi(allMethods, 'openmeeg')); isDuneuro = any(strcmpi(allMethods, 'duneuro')); - isDuneuro2026 = any(strcmpi(allMethods, 'duneuro_2026')); + isDuneuro2026 = any(strcmpi(allMethods, 'duneuro2026')); % Get protocol description ProtocolInfo = bst_get('ProtocolInfo'); @@ -695,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; From 54ddd51a042fb575ff798b4006906cf1a6629915 Mon Sep 17 00:00:00 2001 From: tmedani Date: Wed, 11 Feb 2026 08:52:20 -0800 Subject: [PATCH 04/11] Updates and bug fixes Updates and bug fixes --- .../{bst_duneuro_2026.m => bst_duneuro2026.m} | 32 ++++++++++--------- toolbox/forward/bst_headmodeler.m | 6 ++-- 2 files changed, 20 insertions(+), 18 deletions(-) rename toolbox/forward/{bst_duneuro_2026.m => bst_duneuro2026.m} (96%) diff --git a/toolbox/forward/bst_duneuro_2026.m b/toolbox/forward/bst_duneuro2026.m similarity index 96% rename from toolbox/forward/bst_duneuro_2026.m rename to toolbox/forward/bst_duneuro2026.m index e0eb10bb19..8b62ecf399 100644 --- a/toolbox/forward/bst_duneuro_2026.m +++ b/toolbox/forward/bst_duneuro2026.m @@ -1,4 +1,4 @@ -function [Gain, errMsg] = bst_duneuro_2026(cfg) +function [Gain, errMsg] = bst_duneuro2026(cfg) % BST_DUNEURO: Call Duneuro to compute a FEM solution for Brainstorm. % % USAGE: [Gain, errMsg] = bst_duneuro(cfg) @@ -42,10 +42,10 @@ %% ===== SENSORS ===== % Select modality for DUNEuro -isEeg = strcmpi(cfg.EEGMethod, 'duneuro_2026') && ~isempty(cfg.iEeg); -isMeg = strcmpi(cfg.MEGMethod, 'duneuro_2026') && ~isempty(cfg.iMeg); -isEcog = strcmpi(cfg.ECOGMethod, 'duneuro_2026') && ~isempty(cfg.iEcog); -isSeeg = strcmpi(cfg.SEEGMethod, 'duneuro_2026') && ~isempty(cfg.iSeeg); +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) @@ -72,7 +72,7 @@ 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(iInteg)]; + MegChannels = [MegChannels; iChan, sChan.Loc(:,iInteg)', sChan.Orient(:,iInteg)', sChan.Weight]; end end % In the case where the MEG integration points are used @@ -261,10 +261,10 @@ bstdn_write_pem_electrodes(TmpDir, EegLoc', 'measurement'); end % Write the MEG sensors data -if isMeg +if isMeg || isMeeg % Write coil data - coil_to_channel_transform = eye(length(MegChannels)); - coil_to_channel_transform = coil_to_channel_transform.* MegChannels(:,end); + % 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); end @@ -288,7 +288,7 @@ transfer_matrix_config.do_eeg = 'True'; end transfer_matrix_config.residual_reduction = '1e-16'; -transfer_matrix_config.nr_threads = '-1'; % can be used as user parameters +transfer_matrix_config.nr_threads = '-1'; % can be used as user parameters ==> highlighted as super parameters % Check with Malte is there is an optimised number without overwhelming the % user computer. @@ -307,8 +307,10 @@ leadfield_config.do_meg = 'True'; leadfield_config.do_eeg = 'True'; end +% set final hard code value leadfield_config.eeg_scaling = '1e0'; % check with Malte if those value are optimised -leadfield_config.meg_scaling = '1e5'; +leadfield_config.meg_scaling = '1e5'; % Malte to check and get final value for MKSA system. +% permeability: leadfield_config.sourcemodel = cfg.SrcModel2026; % [select from the interface: 'multipolar_venant', 'local_subtraction', 'partial_integration'] leadfield_config.nr_threads = '-1'; % same as above @@ -439,9 +441,9 @@ 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 + % 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"); @@ -449,7 +451,7 @@ 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'); + h5write(io_file_path, "/measurement/sensors/magnetometers/coil_to_channel_transform", coil_to_channel_transform); % if all goes well, return 1 iOk = 1; diff --git a/toolbox/forward/bst_headmodeler.m b/toolbox/forward/bst_headmodeler.m index 927f13f7c6..cd54bde308 100644 --- a/toolbox/forward/bst_headmodeler.m +++ b/toolbox/forward/bst_headmodeler.m @@ -505,15 +505,15 @@ %% ===== COMPUTE: DUNEURO ===== if ismember('duneuro', {OPTIONS.MEGMethod, OPTIONS.EEGMethod, OPTIONS.ECOGMethod, OPTIONS.SEEGMethod}) ... - || ismember('duneuro_2026', {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 if ismember('duneuro', {OPTIONS.MEGMethod, OPTIONS.EEGMethod, OPTIONS.ECOGMethod, OPTIONS.SEEGMethod}) [Gain_dn, errMessage] = bst_duneuro(OPTIONS); - elseif ismember('duneuro_2026', {OPTIONS.MEGMethod, OPTIONS.EEGMethod, OPTIONS.ECOGMethod, OPTIONS.SEEGMethod}) - [Gain_dn, errMessage] = bst_duneuro_2026(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'}; From 84f46d1ea481e39cd6ca378e8ec05304d14acccc Mon Sep 17 00:00:00 2001 From: tmedani Date: Thu, 12 Feb 2026 14:26:19 -0800 Subject: [PATCH 05/11] Add thread detection as UI for users Replace hardcoded thread settings with a configurable NbOfThread: bst_duneuro2026 now uses cfg.NbOfThread for transfer_matrix and leadfield nr_threads instead of '-1'. Add automatic thread detection in duneuro_defaults (GetThreads) and expose cfgDef.NbOfThread and NbOfThreadMax defaults. Update the duneuro2026 panel to add a "Number of Thread" solver option, wire the new jTextNbOfThread control into the panel state (visibility, validation, and GetPanelContents). Includes minor layout/whitespace adjustments and a safe fallback for single-threaded systems. --- toolbox/forward/bst_duneuro2026.m | 5 +-- toolbox/forward/duneuro_defaults.m | 18 +++++++---- toolbox/forward/panel_duneuro2026.m | 50 ++++++++++++++++++----------- 3 files changed, 45 insertions(+), 28 deletions(-) diff --git a/toolbox/forward/bst_duneuro2026.m b/toolbox/forward/bst_duneuro2026.m index 8b62ecf399..a110c268c8 100644 --- a/toolbox/forward/bst_duneuro2026.m +++ b/toolbox/forward/bst_duneuro2026.m @@ -288,7 +288,8 @@ transfer_matrix_config.do_eeg = 'True'; end transfer_matrix_config.residual_reduction = '1e-16'; -transfer_matrix_config.nr_threads = '-1'; % can be used as user parameters ==> highlighted as super parameters +transfer_matrix_config.nr_threads = cfg.NbOfThread; +% can be used as user parameters ==> highlighted as super parameters % Check with Malte is there is an optimised number without overwhelming the % user computer. @@ -313,7 +314,7 @@ % permeability: leadfield_config.sourcemodel = cfg.SrcModel2026; % [select from the interface: 'multipolar_venant', 'local_subtraction', 'partial_integration'] -leadfield_config.nr_threads = '-1'; % same as above +leadfield_config.nr_threads = cfg.NbOfThread; % same as above %% ===== RUN DUNEURO ====== bst_progress('text', 'DUNEuro: Computing leadfield...'); % disp(['DUNEURO> System call: ' callStr]); diff --git a/toolbox/forward/duneuro_defaults.m b/toolbox/forward/duneuro_defaults.m index 7b4051624f..5d06b34445 100644 --- a/toolbox/forward/duneuro_defaults.m +++ b/toolbox/forward/duneuro_defaults.m @@ -58,6 +58,7 @@ cfgDef.SolSubstractMean = false; % [solution.solver] cfgDef.SolSolverReduction = 1e-10; +[cfgDef.NbOfThread, cfgDef.NbOfThreadMax] = GetThreads(); % [solution.source_model] cfgDef.SrcModel = 'venant'; % partial_integration, venant, subtraction cfgDef.SrcModel2026 = 'multipolar_venant'; % 'multipolar_venant', 'local_subtraction', 'partial_integration' @@ -76,7 +77,6 @@ cfgDef.BstMegTransferFile = 'meg_transfer.dat'; cfgDef.BstEegLfFile = 'eeg_lf.dat'; cfgDef.BstMegLfFile = 'meg_lf.dat'; - % [MEG computation Options] cfgDef.UseIntegrationPoint = 1; cfgDef.EnableCacheMemory = 0; @@ -88,9 +88,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 ]; @@ -135,7 +132,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 index fc3ffed8c2..2c18f6bc77 100644 --- a/toolbox/forward/panel_duneuro2026.m +++ b/toolbox/forward/panel_duneuro2026.m @@ -141,23 +141,33 @@ 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: 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 @@ -230,6 +240,7 @@ 'jRadioSrcModelSub', jRadioSrcModelSub, ... 'jRadioSrcModelPar', jRadioSrcModelPar, ... 'jTextSrcShrink', jTextSrcShrink, ... + 'jTextNbOfThread', jTextNbOfThread, ... 'jCheckSrcForceInGM', jCheckSrcForceInGM, ... 'jCheckSaveTransfer', jCheckSaveTransfer, ... 'jCheckUseIntegrationPoint', jCheckUseIntegrationPoint,... @@ -281,6 +292,7 @@ function UpdatePanel(isForced) jPanelSrcModel.setVisible(ExpertMode); jPanelInput.setVisible(ExpertMode); jPanelOutput.setVisible(ExpertMode); + jPanelOptSolver.setVisible(ExpertMode); jPanelMegComputationOption.setVisible(ExpertMode); % Update expert button if ExpertMode @@ -339,6 +351,7 @@ function UpdatePanel(isForced) 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 @@ -393,8 +406,7 @@ function UpdatePanel(isForced) end -%% ===== DETECTION FUNCTION ===== -% Check the type of a layer based on its name + % 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); From c10053401ac013724b42d508755d9c71cc538b3a Mon Sep 17 00:00:00 2001 From: tmedani Date: Sat, 14 Feb 2026 21:03:17 -0800 Subject: [PATCH 06/11] Update Configuration and add coil mapping Replace hard-coded DUNEuro parameters with cfg-driven values and add a coil->channel transform builder. --- toolbox/forward/bst_duneuro2026.m | 64 +++++++++++++++++++++++------- toolbox/forward/duneuro_defaults.m | 14 ++++++- 2 files changed, 62 insertions(+), 16 deletions(-) diff --git a/toolbox/forward/bst_duneuro2026.m b/toolbox/forward/bst_duneuro2026.m index a110c268c8..a770513ad5 100644 --- a/toolbox/forward/bst_duneuro2026.m +++ b/toolbox/forward/bst_duneuro2026.m @@ -36,7 +36,7 @@ % - via docker % - via podman % - via apptainer -runner = 'docker'; +runner = cfg.containerRunner; % Get DUNEuro container executable bst_plugin('SetProgressLogo', 'duneuro'); @@ -267,6 +267,7 @@ coil_to_channel_transform = MegChannels(:,8:end); dnbst_write_magnetometers(TmpDir, MegChannels(:,2:4), MegChannels(:,5:7), coil_to_channel_transform); end +coil_to_channel_transform = build_coil_to_channel_transform_matrix(MegChannels) %% ===== CONDUCTIVITY MODEL ===== dnbst_write_volume_conductor(TmpDir, FemMat, cfg.FemCond); @@ -287,11 +288,8 @@ transfer_matrix_config.do_meg = 'True'; transfer_matrix_config.do_eeg = 'True'; end -transfer_matrix_config.residual_reduction = '1e-16'; -transfer_matrix_config.nr_threads = cfg.NbOfThread; -% can be used as user parameters ==> highlighted as super parameters -% Check with Malte is there is an optimised number without overwhelming the -% user computer. +transfer_matrix_config.residual_reduction = cfg.residual_reduction; +transfer_matrix_config.nr_threads = num2str(cfg.NbOfThread); %% ===== LEADFIELD MATRIX CONFIGURATION ===== leadfield_config = []; @@ -309,12 +307,10 @@ leadfield_config.do_eeg = 'True'; end % set final hard code value -leadfield_config.eeg_scaling = '1e0'; % check with Malte if those value are optimised -leadfield_config.meg_scaling = '1e5'; % Malte to check and get final value for MKSA system. -% permeability: - +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 = cfg.NbOfThread; % same as above +leadfield_config.nr_threads = num2str(cfg.NbOfThread); % same as above %% ===== RUN DUNEURO ====== bst_progress('text', 'DUNEuro: Computing leadfield...'); % disp(['DUNEURO> System call: ' callStr]); @@ -448,16 +444,47 @@ 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"); + % 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); + % 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 coil_to_channel_transform = build_coil_to_channel_transform_matrix(MegChannels) + % Extract channel indices and weights + chan_idx = MegChannels(:,1); + w = MegChannels(:,end); + % Find unique channels + channels = unique(chan_idx); + nb_chan = length(channels); + % Count coils per channel (should be 4) + nCoil = sum(chan_idx == channels(1)); + % Total size + N = nb_chan * nCoil; + % Preallocate sparse matrix + Wbig = sparse(N, N); + % Build matrix + for k = 1:nb_chan + % Get rows corresponding to this channel + rows_k = find(chan_idx == channels(k)); + % Extract weights for this channel + weights_k = w(rows_k); + % Define block columns + cols = (k-1)*nCoil + (1:nCoil); + % Place weights in row k + Wbig(k, cols) = weights_k(:)'; + end + + coil_to_channel_transform = Wbig; + % not sure about this + h5create(io_file_path, "/measurement/sensors/magnetometers/coil_to_channel_transform", [nb_chan, nb_chan], Datatype="double"); + h5write(io_file_path, "/measurement/sensors/magnetometers/coil_to_channel_transform", coil_to_channel_transform); +end + function [iOk, errMsg] = bstdn_write_pem_electrodes(duneuro_io_dir, electrode_positions, electrode_type_flag) iOk = 0; errMsg = ''; @@ -541,7 +568,16 @@ 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 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; diff --git a/toolbox/forward/duneuro_defaults.m b/toolbox/forward/duneuro_defaults.m index 5d06b34445..c2f3596460 100644 --- a/toolbox/forward/duneuro_defaults.m +++ b/toolbox/forward/duneuro_defaults.m @@ -58,10 +58,8 @@ cfgDef.SolSubstractMean = false; % [solution.solver] cfgDef.SolSolverReduction = 1e-10; -[cfgDef.NbOfThread, cfgDef.NbOfThreadMax] = GetThreads(); % [solution.source_model] cfgDef.SrcModel = 'venant'; % partial_integration, venant, subtraction -cfgDef.SrcModel2026 = 'multipolar_venant'; % 'multipolar_venant', 'local_subtraction', 'partial_integration' cfgDef.SrcIntorderadd = 0; cfgDef.SrcIntorderadd_lb = 2; cfgDef.SrcNbMoments = 3; @@ -81,6 +79,18 @@ 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; From be54cfbbe295881421d1f33aa493e1cd8f9ee322 Mon Sep 17 00:00:00 2001 From: tmedani Date: Sun, 15 Feb 2026 11:04:50 -0800 Subject: [PATCH 07/11] add support for anisotropic conductivity Add options to process and write the anisotropic conductivity --- toolbox/forward/bst_duneuro2026.m | 19 ++++++++++++++++++- .../process/functions/process_fem_tensors.m | 4 +--- 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/toolbox/forward/bst_duneuro2026.m b/toolbox/forward/bst_duneuro2026.m index a770513ad5..89114af580 100644 --- a/toolbox/forward/bst_duneuro2026.m +++ b/toolbox/forward/bst_duneuro2026.m @@ -267,11 +267,27 @@ coil_to_channel_transform = MegChannels(:,8:end); dnbst_write_magnetometers(TmpDir, MegChannels(:,2:4), MegChannels(:,5:7), coil_to_channel_transform); end -coil_to_channel_transform = build_coil_to_channel_transform_matrix(MegChannels) +coil_to_channel_transform = build_coil_to_channel_transform_matrix(MegChannels); %% ===== CONDUCTIVITY MODEL ===== +% Isotropic without tensor +if ~cfg.UseTensor dnbst_write_volume_conductor(TmpDir, FemMat, cfg.FemCond); % NOTE: online viewer of hdf5 files : https://myhdf5.hdfgroup.org/ +else % With tensor (isotropic or anisotropic) + % Transformation matrix and tensor mapping on each direction + CondTensor = zeros(length(FemMat.Elements),6) ; + for ind =1 : length(FemMat.Elements) + temp0 = reshape(FemMat.Tensors(ind,:),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 + CondTensor(ind,:) = [temp(1) temp(5) temp(9) temp(4) temp(8) temp(7)]; % this is the right order + end + % write the tensors + % TODO: Need to double check if it is correct with Malte. + dnbst_write_volume_conductor(TmpDir, FemMat, cfg.FemCond); +end %% ===== TRANSFER MATRIX CONFIGURATION ===== transfer_matrix_config = []; @@ -558,6 +574,7 @@ % 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 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 From 39b53b7888fd42d64a914d66476202a7cfa28fd6 Mon Sep 17 00:00:00 2001 From: tmedani Date: Sun, 15 Feb 2026 11:53:21 -0800 Subject: [PATCH 08/11] update the tensor interface update and Optimize the tensor interface --- toolbox/forward/bst_duneuro2026.m | 37 ++++++++++++++++--------------- 1 file changed, 19 insertions(+), 18 deletions(-) diff --git a/toolbox/forward/bst_duneuro2026.m b/toolbox/forward/bst_duneuro2026.m index 89114af580..539783cb2d 100644 --- a/toolbox/forward/bst_duneuro2026.m +++ b/toolbox/forward/bst_duneuro2026.m @@ -261,33 +261,40 @@ bstdn_write_pem_electrodes(TmpDir, EegLoc', 'measurement'); end % Write the MEG sensors data -if isMeg || isMeeg +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); + coil_to_channel_transform = build_coil_to_channel_transform_matrix(MegChannels); + end -coil_to_channel_transform = build_coil_to_channel_transform_matrix(MegChannels); %% ===== CONDUCTIVITY MODEL ===== % Isotropic without tensor if ~cfg.UseTensor -dnbst_write_volume_conductor(TmpDir, FemMat, cfg.FemCond); + % Create tensor per tissue format from the isoconductivity + tensors = zeros(length(isoCond), 3, 3); + for iTissue = 1 : length(FemMat.TissueLabels) + tensors(iTissue,1,1) = Cond(iTissue); + tensors(iTissue,2,2) = Cond(iTissue); + tensors(iTissue,3,3) = Cond(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 - CondTensor = zeros(length(FemMat.Elements),6) ; - for ind =1 : length(FemMat.Elements) - temp0 = reshape(FemMat.Tensors(ind,:),3,[]); + 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 - CondTensor(ind,:) = [temp(1) temp(5) temp(9) temp(4) temp(8) temp(7)]; % this is the right order + tensors(iElem,:,:) = temp; end % write the tensors % TODO: Need to double check if it is correct with Malte. - dnbst_write_volume_conductor(TmpDir, FemMat, cfg.FemCond); end +dnbst_write_volume_conductor(TmpDir, FemMat, tensors); %% ===== TRANSFER MATRIX CONFIGURATION ===== transfer_matrix_config = []; @@ -377,7 +384,7 @@ %% ================================================================================= % === SUPPORT FUNCTIONS ========================================================= % ================================================================================= -function [iOk, errMsg] = dnbst_write_volume_conductor(duneuro_io_dir, FemMat, isoCond) +function [iOk, errMsg] = dnbst_write_volume_conductor(duneuro_io_dir, FemMat, Tensors) iOk = 0; errMsg = ''; nodes = FemMat.Vertices; elements = FemMat.Elements - 1; @@ -388,7 +395,7 @@ nr_nodes = size(nodes, 1); nr_elements = size(elements, 1); nr_labels = size(labels, 1); - nr_unique_tensors = length(isoCond); + nr_unique_tensors = length(Cond); dim = 3; nr_nodes_per_tetrahedron = 4; @@ -401,13 +408,7 @@ 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 - % Create tensor format from the isoconductivity - tensors = zeros(length(isoCond), 3, 3); - for iTissue = 1 : length(FemMat.TissueLabels) - tensors(iTissue,1,1) = isoCond(iTissue); - tensors(iTissue,2,2) = isoCond(iTissue); - tensors(iTissue,3,3) = isoCond(iTissue); - end + if (size(tensors, 2) ~= dim) || (size(tensors, 3) ~= dim) errMsg = 'The shape of the tensors array must be (K, 3, 3)'; return; end @@ -423,7 +424,7 @@ 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])); + 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'); From 9812eb58116696320a34dff46066ef2a696ebbe8 Mon Sep 17 00:00:00 2001 From: tmedani Date: Sun, 15 Feb 2026 12:31:07 -0800 Subject: [PATCH 09/11] add comments --- toolbox/forward/bst_duneuro2026.m | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/toolbox/forward/bst_duneuro2026.m b/toolbox/forward/bst_duneuro2026.m index 539783cb2d..0481ffe4e8 100644 --- a/toolbox/forward/bst_duneuro2026.m +++ b/toolbox/forward/bst_duneuro2026.m @@ -2,6 +2,7 @@ % 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: @@ -293,6 +294,9 @@ 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); From 90111052f135044b2959cdde6c6dd9222414b7e4 Mon Sep 17 00:00:00 2001 From: tmedani Date: Thu, 19 Feb 2026 16:39:29 -0800 Subject: [PATCH 10/11] Update bst_duneuro2026.m Minor typos corrections and adding options for chunk of dipoles --- toolbox/forward/bst_duneuro2026.m | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/toolbox/forward/bst_duneuro2026.m b/toolbox/forward/bst_duneuro2026.m index 0481ffe4e8..026875aed9 100644 --- a/toolbox/forward/bst_duneuro2026.m +++ b/toolbox/forward/bst_duneuro2026.m @@ -39,6 +39,8 @@ % - 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 ===== @@ -275,11 +277,11 @@ % Isotropic without tensor if ~cfg.UseTensor % Create tensor per tissue format from the isoconductivity - tensors = zeros(length(isoCond), 3, 3); + tensors = zeros(length(cfg.FemCond), 3, 3); for iTissue = 1 : length(FemMat.TissueLabels) - tensors(iTissue,1,1) = Cond(iTissue); - tensors(iTissue,2,2) = Cond(iTissue); - tensors(iTissue,3,3) = Cond(iTissue); + 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) @@ -335,9 +337,17 @@ 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.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]); @@ -388,7 +398,7 @@ %% ================================================================================= % === SUPPORT FUNCTIONS ========================================================= % ================================================================================= -function [iOk, errMsg] = dnbst_write_volume_conductor(duneuro_io_dir, FemMat, Tensors) +function [iOk, errMsg] = dnbst_write_volume_conductor(duneuro_io_dir, FemMat, tensors) iOk = 0; errMsg = ''; nodes = FemMat.Vertices; elements = FemMat.Elements - 1; @@ -399,7 +409,7 @@ nr_nodes = size(nodes, 1); nr_elements = size(elements, 1); nr_labels = size(labels, 1); - nr_unique_tensors = length(Cond); + nr_unique_tensors = size(tensors,1); dim = 3; nr_nodes_per_tetrahedron = 4; @@ -428,7 +438,7 @@ 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])); + 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'); From 776874a38cc182a8fc530680f496e1e0349fbf53 Mon Sep 17 00:00:00 2001 From: tmedani Date: Tue, 24 Mar 2026 08:58:39 -0700 Subject: [PATCH 11/11] Update bst_duneuro2026.m Update the IP definition --- toolbox/forward/bst_duneuro2026.m | 44 +++++++++++++++---------------- 1 file changed, 21 insertions(+), 23 deletions(-) diff --git a/toolbox/forward/bst_duneuro2026.m b/toolbox/forward/bst_duneuro2026.m index 026875aed9..f823f1b92a 100644 --- a/toolbox/forward/bst_duneuro2026.m +++ b/toolbox/forward/bst_duneuro2026.m @@ -269,8 +269,7 @@ % 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); - coil_to_channel_transform = build_coil_to_channel_transform_matrix(MegChannels); - + build_coil_to_channel_transform_matrix(MegChannels); end %% ===== CONDUCTIVITY MODEL ===== @@ -485,35 +484,34 @@ iOk = 1; end -function coil_to_channel_transform = build_coil_to_channel_transform_matrix(MegChannels) +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); - % Count coils per channel (should be 4) - nCoil = sum(chan_idx == channels(1)); - % Total size - N = nb_chan * nCoil; - % Preallocate sparse matrix - Wbig = sparse(N, N); - % Build matrix - for k = 1:nb_chan - % Get rows corresponding to this channel - rows_k = find(chan_idx == channels(k)); - % Extract weights for this channel - weights_k = w(rows_k); - % Define block columns - cols = (k-1)*nCoil + (1:nCoil); - % Place weights in row k - Wbig(k, cols) = weights_k(:)'; - end - coil_to_channel_transform = Wbig; + % 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_chan], Datatype="double"); - h5write(io_file_path, "/measurement/sensors/magnetometers/coil_to_channel_transform", coil_to_channel_transform); + 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)