Skip to content

Commit

Permalink
Merge branch 'master' of github.com:fieldtrip/fieldtrip
Browse files Browse the repository at this point in the history
  • Loading branch information
@robertoostenveld
robertoostenveld committed Jul 13, 2018
2 parents b0bd5ec + ae722b2 commit db973c9
Show file tree
Hide file tree
Showing 6 changed files with 83 additions and 64 deletions.
37 changes: 24 additions & 13 deletions external/openmeeg/testOpenMEEGeeg.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,16 +5,16 @@

% Copyright (C) 2010-2017, OpenMEEG developers

addpath(cd) % Make sure current folder is in the path
addpath(pwd) % Make sure current folder is in the path

%% Set the position of the probe dipole
dippos = [0 0 70];

%% Set the radius and conductivities of each of the compartments

% 4 Layers
r = [85 88 92 100];
c = [1 1/20 1/80 1];
r = [100 92 88 85];
c = [1 1/80 1/20 1];

[rdms,mags] = run_bem_computation(r,c,dippos);

Expand All @@ -28,7 +28,7 @@
assert(norm(mags-[0.84467 0.84469 0.83887])<thr)

% 3 Layers
r = [88 92 100];
r = [100 92 88];
c = [1 1/80 1];

[rdms,mags] = run_bem_computation(r,c,dippos);
Expand All @@ -38,8 +38,8 @@
assert(norm(mags-[1.0498 1.0498 1.0207])<thr)

% 2 Layers
r = [92 100];
c = [1 1/4];
r = [100 92];
c = [1/4 1];

[rdms,mags] = run_bem_computation(r,c,dippos);
% assertElementsAlmostEqual(rdms, [0.15514 0.15514 0.1212], 'absolute', 1e-3)
Expand Down Expand Up @@ -75,24 +75,35 @@
end

%% Create a triangulated mesh, the first boundary is inside
vol = [];
bnd = [];
for ii=1:length(r)
vol.bnd(ii).pos = pos * r(ii);
vol.bnd(ii).tri = tri;
bnd(ii).pos = pos * r(ii);
bnd(ii).tri = tri;
end


%% Compute the BEM model
cfg = [];
cfg.conductivity = c;
cfg.method = 'openmeeg';

vol_bem = ft_prepare_headmodel(cfg, vol);

cfg.headmodel = vol_bem;
% vol structure only needs these elements;
% ft_prepare_headmodel should *not* be used!
vol_bem = [];
vol_bem.bnd = bnd;
vol_bem.cond = c;
vol_bem.type = 'openmeeg';

% ft_prepare_headmodel has a bug; likely the geom file does not match
% the order of the layers
%cfg.conductivity = c;
%vol_bem = ft_prepare_headmodel(cfg, bnd);
%vol_bem = rmfield(vol_bem,'mat');

cfg.vol = vol_bem;
cfg.grid.pos = dippos;
cfg.grid.unit = 'mm';
cfg.elec = sens;
% grid = ft_prepare_leadfield(cfg);
grid = ft_prepare_leadfield(cfg);

lf_openmeeg = grid.leadfield{1};
Expand Down
2 changes: 1 addition & 1 deletion forward/ft_headmodel_openmeeg.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,7 +96,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% store the current path and change folder to the temporary one
tmpfolder = cd;
tmpfolder = pwd;
bndom = headmodel.bnd;

try
Expand Down
64 changes: 36 additions & 28 deletions forward/private/leadfield_openmeeg.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,22 +101,23 @@
method = ft_getopt(vol,'method','hminv');
ecog = ft_getopt(vol,'ecog','no');

% Use basefile and basepath for saving files
if isfield(vol,'basefile')
basefile = vol.basefile;
else
basefile = tempname;
end
% Use basefile and workdir for saving files
if isfield(vol,'path')
path = vol.path;
workdir = vol.path;
cleanup_flag = false;
else
path = fullfile(tempdir,'ft-om');
workdir = fullfile(tempdir,'ft-om');
cleanup_flag = true;
end
mkdir(path);
mkdir(workdir);

sensorFile = fullfile(path, [basefile '_sensorcoords.txt']);
if isfield(vol,'basefile')
basefile = vol.basefile;
else
[~,basefile] = fileparts(tempname);
end

sensorFile = fullfile(workdir, [basefile '_sensorcoords.txt']);
if exist(sensorFile,'file')
disp('Sensor coordinate file already exists. Skipping...')
else
Expand Down Expand Up @@ -166,15 +167,15 @@
fclose(fid);
end

condFile = fullfile(path, [basefile '.cond']);
condFile = fullfile(workdir, [basefile '.cond']);
if exist(condFile,'file')
disp('Conductivity file already exists. Skipping...')
else
disp('Writing conductivity file...')
write_cond(vol,condFile);
end

geomFile = fullfile(path, [basefile '.geom']);
geomFile = fullfile(workdir, [basefile '.geom']);
if exist(geomFile,'file')
disp('Geometry descriptor file already exists. Skipping...')
else
Expand All @@ -184,7 +185,7 @@

disp('Writing OpenMEEG mesh files...')

write_mesh(vol,path,basefile);
write_mesh(vol,workdir,basefile);

disp('Validating mesh...')
[om_status om_msg] = system([fullfile(OPENMEEG_PATH, 'om_check_geom'), ' -g ', geomFile])
Expand All @@ -196,7 +197,7 @@
chunks = ceil(size(voxels,1)/VOXCHUNKSIZE);
dipFile = cell(chunks,1);
for ii = 1:chunks
dipFile{ii} = fullfile(path, [basefile '_voxels' num2str(ii) om_ext]);
dipFile{ii} = fullfile(workdir, [basefile '_voxels' num2str(ii) om_ext]);
if exist(dipFile{ii},'file')
fprintf('\t%s already exists. Skipping...\n', dipFile{ii});
else
Expand All @@ -207,7 +208,7 @@
end


hmFile = fullfile(path, [basefile '_hm' om_ext]);
hmFile = fullfile(workdir, [basefile '_hm' om_ext]);
if exist(hmFile,'file')
disp('Head matrix already exists. Skipping...')
else
Expand All @@ -219,7 +220,7 @@
end

if strcmp(method,'hminv')
hminvFile = fullfile(path, [basefile '_hminv' om_ext]);
hminvFile = fullfile(workdir, [basefile '_hminv' om_ext]);
if exist(hminvFile,'file')
disp('Inverse head matrix already exists. Skipping...');
else
Expand All @@ -237,7 +238,7 @@

dsmFile = cell(chunks,1);
for ii = 1:chunks
dsmFile{ii} = fullfile(path, [basefile '_dsm' num2str(ii) om_ext]);
dsmFile{ii} = fullfile(workdir, [basefile '_dsm' num2str(ii) om_ext]);
if exist(dsmFile{ii},'file')
fprintf('\t%s already exists. Skipping...\n', dsmFile{ii});
else
Expand All @@ -255,14 +256,14 @@

if ft_senstype(sens, 'eeg')
if strcmp(ecog,'yes')
ohmicFile = fullfile(path, [basefile '_h2ecogm']);
ohmicFile = fullfile(workdir, [basefile '_h2ecogm']);
cmd = '-h2ecogm';
else
ohmicFile = fullfile(path, [basefile '_h2em' om_ext]);
ohmicFile = fullfile(workdir, [basefile '_h2em' om_ext]);
cmd = '-h2em';
end
else
ohmicFile = fullfile(path, [basefile '_h2mm' om_ext]);
ohmicFile = fullfile(workdir, [basefile '_h2mm' om_ext]);
cmd = '-h2mm';
end
if exist(ohmicFile,'file')
Expand All @@ -279,7 +280,7 @@
disp('Contribution of all sources to the MEG sensors')
scFile = cell(chunks,1);
for ii = 1:chunks
scFile{ii} = fullfile(path, [basefile '_ds2mm' num2str(ii) om_ext]);
scFile{ii} = fullfile(workdir, [basefile '_ds2mm' num2str(ii) om_ext]);
if exist(scFile{ii},'file')
fprintf('\t%s already exists. Skipping...\n',scFile{ii})
else
Expand All @@ -295,9 +296,9 @@
bemFile = cell(chunks,1);
for ii = 1:chunks
if ft_senstype(sens, 'eeg')
bemFile{ii} = fullfile(path, [basefile '_eeggain' num2str(ii) om_ext]);
bemFile{ii} = fullfile(workdir, [basefile '_eeggain' num2str(ii) om_ext]);
else
bemFile{ii} = fullfile(path, [basefile '_meggain' num2str(ii) om_ext]);
bemFile{ii} = fullfile(workdir, [basefile '_meggain' num2str(ii) om_ext]);
end

if exist(bemFile{ii},'file')
Expand All @@ -324,14 +325,21 @@
end

% Import lead field/potential
[g, voxels_in] = import_gain(path, basefile, ft_senstype(sens, 'eeg'));
if (voxels_in ~= voxels) && (nargout == 1); ft_warning('Imported voxels from OpenMEEG process not the same as function input.'); end;
[g, voxels_in] = import_gain(workdir, basefile, ft_senstype(sens, 'eeg'));
if find(voxels_in ~= voxels) & (nargout == 1)
ft_warning('Imported voxels from OpenMEEG process not the same as function input.');
end

% multiplication with sens.tra matrix should happen in ft_prepare_leadfield
% instead:
sens.tra = eye(size(g,1))
lp = sens.tra*g; % Mchannels x (3 orientations x Nvoxels)

%lp = g; % Mchannels x (3 orientations x Nvoxels)

% Cleanup
if cleanup_flag
rmdir(basepath,'s')
rmdir(workdir,'s')
end
if (isunix && ~ismac)
setenv('LD_LIBRARY_PATH',ldLibraryPath0);
Expand All @@ -351,14 +359,14 @@ function write_cond(vol,filename)
fclose(fid);


function write_geom(vol,filename,basepathfile)
function write_geom(vol,filename,basefile)
fid=fopen(filename,'w');
fprintf(fid,'# Domain Description 1.0\n');
fprintf(fid,'Interfaces %i Mesh\n',length(vol.cond));
tissues={'_scalp.tri\n'; '_skull.tri\n'; '_csf.tri\n'; '_brain.tri\n'};
if length(vol.cond)==3; tissues = tissues([1 2 4]); end;
for ii = 1:length(vol.cond)
fprintf(fid,[basepathfile tissues{ii}]);
fprintf(fid,[basefile tissues{ii}]);
end
fprintf('\n');
fprintf(fid,'Domains %i\n',length(vol.cond)+1);
Expand Down
2 changes: 1 addition & 1 deletion ft_prepare_headmodel.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -262,7 +262,7 @@
tmpcfg.tissue = cfg.tissue;
geometry = ft_prepare_mesh(tmpcfg, data);
else
ft_error('Either a segmentated MRI or data with closed triangulated mesh is required as data input for the bemcp, dipoli or openmeeg method');
ft_error('Either a segmented MRI or data with closed triangulated mesh is required as data input for the bemcp, dipoli or openmeeg method');
end

if strcmp(cfg.method, 'bemcp')
Expand Down
6 changes: 3 additions & 3 deletions private/ft_fetch_sens.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -149,17 +149,17 @@
display('reading electrodes from file ''%s''\n', cfg.elecfile);
sens = ft_read_sens(cfg.elecfile);
% only keep positions and labels in case of EEG electrodes
sens = keepfields(sens, {'elecpos', 'chanpos', 'unit', 'coordsys', 'label'});
sens = keepfields(sens, {'elecpos', 'chanpos', 'unit', 'coordsys', 'label','tra'});
elseif hascfgelec
display('using electrodes specified in the configuration\n');
sens = cfg.elec;
% only keep positions and labels in case of EEG electrodes
sens = keepfields(sens, {'elecpos', 'chanpos', 'unit', 'coordsys', 'label'});
sens = keepfields(sens, {'elecpos', 'chanpos', 'unit', 'coordsys', 'label','tra'});
elseif hasdataelec
display('using electrodes specified in the data\n');
sens = data.elec;
% only keep positions and labels in case of EEG electrodes
sens = keepfields(sens, {'elecpos', 'chanpos', 'unit', 'coordsys', 'label'});
sens = keepfields(sens, {'elecpos', 'chanpos', 'unit', 'coordsys', 'label','tra'});

elseif hasoptofile
display('reading optodes from file ''%s''\n', cfg.optofile);
Expand Down
36 changes: 18 additions & 18 deletions template/atlas/aal/ROI_MNI_V4.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,24 +88,24 @@ T2AG Temporal_Pole_Mid_L 8211
T2AD Temporal_Pole_Mid_R 8212
T3G Temporal_Inf_L 8301
T3D Temporal_Inf_R 8302
CERCRU1G Cerebelum_Crus1_L 9001
CERCRU1D Cerebelum_Crus1_R 9002
CERCRU2G Cerebelum_Crus2_L 9011
CERCRU2D Cerebelum_Crus2_R 9012
CER3G Cerebelum_3_L 9021
CER3D Cerebelum_3_R 9022
CER4_5G Cerebelum_4_5_L 9031
CER4_5D Cerebelum_4_5_R 9032
CER6G Cerebelum_6_L 9041
CER6D Cerebelum_6_R 9042
CER7BG Cerebelum_7b_L 9051
CER7BD Cerebelum_7b_R 9052
CER8G Cerebelum_8_L 9061
CER8D Cerebelum_8_R 9062
CER9G Cerebelum_9_L 9071
CER9D Cerebelum_9_R 9072
CER10G Cerebelum_10_L 9081
CER10D Cerebelum_10_R 9082
CERCRU1G Cerebellum_Crus1_L 9001
CERCRU1D Cerebellum_Crus1_R 9002
CERCRU2G Cerebellum_Crus2_L 9011
CERCRU2D Cerebellum_Crus2_R 9012
CER3G Cerebellum_3_L 9021
CER3D Cerebellum_3_R 9022
CER4_5G Cerebellum_4_5_L 9031
CER4_5D Cerebellum_4_5_R 9032
CER6G Cerebellum_6_L 9041
CER6D Cerebellum_6_R 9042
CER7BG Cerebellum_7b_L 9051
CER7BD Cerebellum_7b_R 9052
CER8G Cerebellum_8_L 9061
CER8D Cerebellum_8_R 9062
CER9G Cerebellum_9_L 9071
CER9D Cerebellum_9_R 9072
CER10G Cerebellum_10_L 9081
CER10D Cerebellum_10_R 9082
VER1_2 Vermis_1_2 9100
VER3 Vermis_3 9110
VER4_5 Vermis_4_5 9120
Expand Down

0 comments on commit db973c9

Please sign in to comment.