fieldtrip · DanielMiklody · Mar 16, 2023 · Mar 18, 2023 · Mar 18, 2023 · Apr 28, 2023
diff --git a/external/openmeeg/openmeeg_msm.m b/external/openmeeg/openmeeg_msm.m
@@ -0,0 +1,108 @@
+function [msm] = openmeeg_msm(pos, headmodel, NAflag, source)
+
+% openmeeg_msm creates the boundary condition part b of a volume conduction
+% model of the head using the boundary element method (BEM). This function 
+% takes as input the triangulated surfaces that describe the boundaries and
+% the source position of single monopoles. It returns as output a vector of
+% boundary conditions which can be used together with a volume conduction 
+% model to compute leadfields.
+%
+% This function implements
+%   Gramfort et al. OpenMEEG: opensource software for quasistatic
+%   bioelectromagnetics. Biomedical engineering online (2010) vol. 9 (1) pp. 45
+%   http://www.biomedical-engineering-online.com/content/9/1/45
+%   doi:10.1186/1475-925X-9-45
+% and
+%   Kybic et al. Generalized head models for MEG/EEG: boundary element method
+%   beyond nested volumes. Phys. Med. Biol. (2006) vol. 51 pp. 1333-1346
+%   doi:10.1088/0031-9155/51/5/021
+%
+% This link with FieldTrip is derived from the OpenMEEG project
+% with contributions from Daniel Wong and Sarang Dalal, and uses external
+% command-line executables. See http://openmeeg.github.io/
+%
+% Use as
+%   dsm = ft_sysmat_openmeeg(pos, headmodel, sens, nonadaptive_flag)
+%
+%
+% See also FT_HEADMODEL_OPENMEEG, FT_SENSINTERP_OPENMEEG, FT_PREPARE_LEADFIELD
+
+%$Id$
+
+ft_hastoolbox('openmeeg', 1);  % add to path (if not yet on path)
+openmeeg_license;              % show the license (only once)
+prefix = om_checkombin;        % check the installation of the binaries
+if(~ispc) % if Linux/Mac, set number of threads
+    omp_num_threads = feature('numCores');
+    prefix = ['export OMP_NUM_THREADS=' num2str(omp_num_threads) ';' prefix];
+end
+
+ecog = ft_getopt(headmodel,'ecog','no');
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% this uses an implementation that was contributed by INRIA Odyssee Team
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% store the current path and change folder to the temporary one
+cwd = pwd;
+bndom = headmodel.bnd;
+
+workdir = fullfile(tempdir,['ft_om_' datestr(now,'ddmmyyHHMMSSFFF') '_' num2str(randi(1000000))]);
+mkdir(workdir);
+
+try
+    % Write the triangulations to file, named after tissue type.
+    % OpenMEEG v2.3 and up internally adjusts the convention for surface
+    % normals, but OpenMEEG v2.2 expects surface normals to point inwards;
+    % this checks and corrects if needed
+    bndfile = fullfile(workdir,strcat(headmodel.tissue,'.tri'));
+    for ii=1:length(bndom)
+        om_save_tri(bndfile{ii}, bndom(ii).pos, bndom(ii).tri);
+    end
+
+    condfile  = fullfile(workdir,'om.cond');
+    geomfile  = fullfile(workdir,'om.geom');
+    dipfile = fullfile(workdir,'dip.txt');
+    msmfile = fullfile(workdir,'msm.bin'); % extension added directly in om_assemble function call below
+
+    % write conductivity and mesh files
+    bndlabel = {};
+    for i=1:length(bndom)
+        [dum,bndlabel{i}] = fileparts(bndfile{i});
+    end
+
+    om_write_geom(geomfile,bndfile,bndlabel);
+    om_write_cond(condfile,headmodel.cond,bndlabel);
+
+    % handle dipole file
+    ndip = size(pos,1);
+    pos = [pos , ones(ndip,1)]; % save pos with each 3D orientation
+    om_save_full(pos,dipfile,'ascii');
+
+    omp_num_threads = feature('numCores');
+
+    if(strcmp(NAflag,'yes'))
+        str = ' -MSMNA';
+    else
+        str = ' -MSM';
+    end
+
+    om_status = system([prefix 'om_assemble' str ' ' geomfile ' ' condfile ' ' dipfile ' ' msmfile]);
+
+    if(om_status ~= 0) % status = 0 if successful
+        ft_error(['Aborting OpenMEEG pipeline due to above error.']);
+    end
+
+    msm = om_load_full(msmfile,'binary');
+catch
+    rethrow(lasterror)
+end
+
+try
+    rmdir(workdir,'s'); % remove workdir with intermediate files
+catch
+    disp(lasterr);
+    rmdir(workdir,'s'); % remove workdir with intermediate files
+    ft_error('an error occurred while running OpenMEEG');
+end
diff --git a/forward/ft_compute_leadfield.m b/forward/ft_compute_leadfield.m
@@ -472,7 +472,7 @@
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       lf = eeg_leadfieldb(dippos, sens.elecpos, headmodel);
 
-    case 'openmeeg'
+    case {'openmeeg','openmeegHArtMuT_dipole'}
       ft_hastoolbox('openmeeg', 1);
 
       dsm         = ft_getopt(varargin, 'dsm');
@@ -483,7 +483,113 @@
         dsm            = ft_sysmat_openmeeg(dippos, headmodel, sens, nonadaptive);
       end
       lf               = ds2sens + h2sens*headmodel.mat*dsm;
-
+
+    case {'openmeegHArtMuT_tripole'}
+            ft_hastoolbox('openmeeg', 1);
+            nonadaptive = ft_getopt(varargin, 'nonadaptive');
+
+
+            if isfield(headmodel,'locs')
+                locs = headmodel.triplocs;
+            else
+                locs=[-1.6 0 3.2]; %adapt units?
+            end
+            if isfield(headmodel,'amps')
+                amps = headmodel.tripamps;
+            else
+                amps=[-1 2 -1];
+            end
+
+            [h2sens,ds2sens] = ft_sensinterp_openmeeg(dippos, headmodel, sens);
+
+            brainsources=surface_inside(dippos,headmodel.bnd(end).pos,headmodel.bnd(end).tri)==1;
+
+            if isfield(headmodel,'eyes')
+                lm=headmodel.eyes.lm;
+                rm=headmodel.eyes.rm;
+                rad=headmodel.eyes.rad;
+            else
+                lm=[-33.7, 58.5, -36.4];
+                rm=[35.2, 59.5, -37.3];
+                rad = 12;
+            end
+            eyesourcesl = sqrt(sum(bsxfun(@minus,dippos,lm).^2,2))<=rad;
+            eyesourcesr = sqrt(sum(bsxfun(@minus,dippos,rm).^2,2))<=rad;
+            musclesources = ~brainsources & ~(eyesourcesl|eyesourcesr);
+            lf=nan(size(h2sens,1),size(dippos,1)*3);
+            brainsources2=reshape(repmat(brainsources,1,3)',[],1);
+            musclesources2=reshape(repmat(musclesources,1,3)',[],1);
+            eyesourcesl2=reshape(repmat(eyesourcesl,1,3)',[],1);
+            eyesourcesr2=reshape(repmat(eyesourcesr,1,3)',[],1);
+            if sum(brainsour
8000
ces)>0
+                dsm            = ft_sysmat_openmeeg(dippos(brainsources,:), headmodel, sens, nonadaptive);
+                lf(:,brainsources2)      = h2sens*headmodel.mat*dsm;
+            end
+            if sum(eyesourcesl|eyesourcesr)>0
+                %mirror each position to the other eye to get symmetric
+                %dipoles
+                eyemirr_l = dippos(eyesourcesl,:)-lm+rm;
+                dsml = ft_sysmat_openmeeg(dippos(eyesourcesl,:), headmodel, sens, nonadaptive);
+                dsmlm = ft_sysmat_openmeeg(eyemirr_l, headmodel, sens, nonadaptive);
+                eyemirr_r = dippos(eyesourcesr,:)-rm+lm;
+                dsmr = ft_sysmat_openmeeg(dippos(eyesourcesr,:), headmodel, sens, nonadaptive);
+                dsmrm = ft_sysmat_openmeeg(eyemirr_r, headmodel, sens, nonadaptive);
+
+                lf(:,eyesourcesl2)      = h2sens*headmodel.mat*(dsml+dsmlm);
+                lf(:,eyesourcesr2)      = h2sens*headmodel.mat*(dsmr+dsmrm);
+            end
+            if sum(musclesources)>0
+                trippos=dippos(musclesources,:);
+                trippos=repmat(trippos,1,1,3);
+                trippos=permute(trippos,[1 3 2]);
+                tripposx=trippos;
+                tripposx(:,:,1)=bsxfun(@plus,locs,tripposx(:,:,1));
+                tripposy=trippos;
+                tripposy(:,:,2)=bsxfun(@plus,locs,tripposy(:,:,2));
+                tripposz=trippos;
+                tripposz(:,:,3)=bsxfun(@plus,locs,tripposz(:,:,3));
+                failsx=~is_inside(tripposx(:,3,:),headmodel.bnd(1))|~is_inside(tripposx(:,1,:),headmodel.bnd(1))|...
+                    is_inside(tripposx(:,3,:),headmodel.bnd(2))|is_inside(tripposx(:,1,:),headmodel.bnd(2));
+                failsy=~is_inside(tripposy(:,3,:),headmodel.bnd(1))|~is_inside(tripposy(:,1,:),headmodel.bnd(1))|...
+                    is_inside(tripposy(:,3,:),headmodel.bnd(2))|is_inside(tripposy(:,1,:),headmodel.bnd(2));
+                failsz=~is_inside(tripposz(:,3,:),headmodel.bnd(1))|~is_inside(tripposz(:,1,:),headmodel.bnd(1))|...
+                    is_inside(tripposz(:,3,:),headmodel.bnd(2))|is_inside(tripposz(:,1,:),headmodel.bnd(2));
+                tripmod=locs;
+                while   any(failsx)
+                    tripmod=tripmod/2;
+                    tripposx(failsx,:,1)=bsxfun(@plus,tripmod,trippos(failsx,:,1));
+                    failsx=~is_inside(tripposx(:,3,:),headmodel.bnd(1))|~is_inside(tripposx(:,1,:),headmodel.bnd(1))|...
+                        is_inside(tripposx(:,3,:),headmodel.bnd(2))|is_inside(tripposx(:,1,:),headmodel.bnd(2));
+                end
+                tripmod=locs;
+                while   any(failsy)
+                    tripmod=tripmod/2;
+                    tripposy(failsy,:,2)=bsxfun(@plus,tripmod,trippos(failsy,:,2));
+                    failsy=~is_inside(tripposy(:,3,:),headmodel.bnd(1))|~is_inside(tripposy(:,1,:),headmodel.bnd(1))|...
+                        is_inside(tripposy(:,3,:),headmodel.bnd(2))|is_inside(tripposy(:,1,:),headmodel.bnd(2));                   
+                end
+                tripmod=locs;
+                while   any(failsz)
+                    tripmod=tripmod/2;
+                    tripposz(failsz,:,3)=bsxfun(@plus,tripmod,trippos(failsz,:,3));
+                    failsz=~is_inside(tripposz(:,3,:),headmodel.bnd(1))|~is_inside(tripposz(:,1,:),headmodel.bnd(1))|...
+                        is_inside(tripposz(:,3,:),headmodel.bnd(2))|is_inside(tripposz(:,1,:),headmodel.bnd(2));
+                end
+                msm              = openmeeg_msm(reshape(tripposx(:,2,:),size(tripposx,1),size(tripposx,3)), headmodel, sens, nonadaptive);
+                msm1_1           = openmeeg_msm(reshape(tripposx(:,1,:),size(tripposx,1),size(tripposx,3)), headmodel, sens, nonadaptive);
+                msm1_2           = openmeeg_msm(reshape(tripposx(:,3,:),size(tripposx,1),size(tripposx,3)), headmodel, sens, nonadaptive);
+                msm2_1           = openmeeg_msm(reshape(tripposy(:,1,:),size(tripposx,1),size(tripposx,3)), headmodel, sens, nonadaptive);
+                msm2_2           = openmeeg_msm(reshape(tripposy(:,3,:),size(tripposx,1),size(tripposx,3)), headmodel, sens, nonadaptive);
+                msm3_1           = openmeeg_msm(reshape(tripposz(:,1,:),size(tripposx,1),size(tripposx,3)), headmodel, sens, nonadaptive);
+                msm3_2           = openmeeg_msm(reshape(tripposz(:,3,:),size(tripposx,1),size(tripposx,3)), headmodel, sens, nonadaptive);
+                tripoles=nan(size(msm,1),size(msm,2),3);
+                tripoles(:,:,1)=amps(1)*msm1_1+amps(2)*msm+amps(3)*msm1_2;
+                tripoles(:,:,2)=amps(1)*msm2_1+amps(2)*msm+amps(3)*msm2_2;
+                tripoles(:,:,3)=amps(1)*msm3_1+amps(2)*msm+amps(3)*msm3_2;
+                tripoles=reshape(permute(tripoles,[1 3 2]),size(tripoles,1),[]);
+                lf(:,musclesources2)     = h2sens*headmodel.mat*tripoles;
+            end
+
     case {'infinite_currentdipole' 'infinite'}
       lf = eeg_infinite_dipole(dippos, sens.elecpos, headmodel);
 

diff --git a/forward/ft_inside_headmodel.m b/forward/ft_inside_headmodel.m
@@ -119,6 +119,11 @@
     [pos, tri] = headsurface(headmodel, [], 'inwardshift', inwardshift, 'surface', 'brain');
     inside = surface_inside(dippos, pos, tri);
 
+  case {'openmeegHArtMuT_dipole', 'openmeegHArtMuT_tripole'}
+    % this is a model with a realistic shape described by a triangulated boundary
+    [pos, tri] = headsurface(headmodel, [], 'inwardshift', inwardshift, 'surface', 'skin');
+    inside = surface_inside(dippos, pos, tri);
+
   case {'simbio', 'duneuro'}
     % this is a model with hexaheders or tetraheders
     if isfield(headmodel, 'tet')

diff --git a/forward/ft_prepare_vol_sens.m b/forward/ft_prepare_vol_sens.m
@@ -523,7 +523,7 @@
           headmodel.mat = headmodel.mat - repmat(avg, size(headmodel.mat,1), 1);
         end
       end
-    case  'openmeeg'
+    case  {'openmeeg','openmeegHArtMuT_tripole','openmeegHArtMuT_dipole'}
       % don't do anything, h2em or h2mm generated later in ft_prepare_leadfield
 
     case 'fns'

diff --git a/ft_prepare_headmodel.m b/ft_prepare_headmodel.m
@@ -286,7 +286,7 @@
     end
     headmodel = ft_headmodel_asa(cfg.headmodel);
 
-  case {'bemcp' 'dipoli' 'openmeeg'}
+  case {'bemcp' 'dipoli' 'openmeeg' 'openmeegHArtMuT_dipole' 'openmeegHArtMuT_tripole'}
     % the low-level functions all need a mesh
     if isfield(data, 'pos') && isfield(data, 'tri')
       if ~isfield(cfg, 'numvertices') || isempty(cfg.numvertices) || isequal(cfg.numvertices, arrayfun(@(x) size(x.pos, 1), data))