forward_model_eeg_pilot.m
% https://www.fieldtriptoolbox.org/workshop/oslo2019/forward_modeling/
close all, clear, clc
%% Initialize FieldTrip
addpath('path\to\fieldtrip-20201117')
ft_defaults
%% Read MRI
fdir='path\to\jpilot';
fname = 'jpilot.nii';
mri = ft_read_mri(fullfile(fdir, fname));
%% Plot original MRI
cfg = [];
ft_sourceplot(cfg, mri);
print('-dpng', fullfile(fdir,'1_original_mri.png'))
%% MRI Realignment using landmarks
cfg = [];
cfg.coordsys = 'neuromag';
% cfg.method = 'interactive';
cfg.method = 'fiducial';
cfg.fiducial.nas = [127 46 54];
cfg.fiducial.lpa = [208 139 40];
cfg.fiducial.rpa = [ 48 142 40];
mri_aligned_fiducials = ft_volumerealign(cfg, mri);
%% Re-slice MRI
cfg = [];
cfg.resolution = 1;
mri_resliced = ft_volumereslice(cfg, mri_aligned_fiducials);
save(fullfile(fdir,'mri_resliced.mat'), 'mri_resliced')
%% Plot resliced MRI
cfg = [];
ft_sourceplot(cfg, mri_resliced);
print('-dpng', fullfile(fdir,'2_mri_aligned_resliced.png'))
%% Segment the brain
cfg = [];
cfg.output = {'brain' 'skull' 'scalp'};
mri_segmented_orig = ft_volumesegment(cfg, mri_resliced);
save(fullfile(fdir,'mri_segmented_orig.mat'), 'mri_segmented_orig')
%% Ensure clearance between adjacent compartments
se = strel('sphere', 4);
binary_brain = mri_segmented_orig.brain;
binary_skull = mri_segmented_orig.skull | mri_segmented_orig.brain;
binary_scalp = mri_segmented_orig.scalp | mri_segmented_orig.skull | mri_segmented_orig.brain;
skull_eroded = imerode(binary_scalp, se) & binary_skull;
brain_eroded = imerode(skull_eroded, se) & binary_brain;
mri_segmented = mri_segmented_orig;
mri_segmented.scalp = binary_scalp & (~skull_eroded);
mri_segmented.skull = skull_eroded & (~brain_eroded);
mri_segmented.brain = brain_eroded;
save(fullfile(fdir,'mri_segmented.mat'), 'mri_segmented')
%% Troubleshoot meshes
mri_segmented_binary = mri_segmented; % make a copy
binary_brain = mri_segmented.brain;
binary_skull = mri_segmented.skull | mri_segmented.brain;
binary_scalp = mri_segmented.scalp | mri_segmented.skull | mri_segmented.brain;
mri_segmented_binary.combined = binary_scalp + binary_skull + binary_brain;
cfg = [];
cfg.funparameter = 'combined';
cfg.location = [-76 -27 32];
ft_sourceplot(cfg, mri_segmented_binary);
%% Creating meshes
cfg = [];
cfg.method = 'projectmesh';
cfg.tissue = 'brain';
cfg.numvertices = 2400;
mesh_brain = ft_prepare_mesh(cfg, mri_segmented);
mesh_brain = ft_convert_units(mesh_brain, 'm'); % Use SI Units
cfg = [];
cfg.method = 'projectmesh';
cfg.tissue = 'skull';
cfg.numvertices = 1600;
mesh_skull = ft_prepare_mesh(cfg, mri_segmented);
mesh_skull = ft_convert_units(mesh_skull, 'm'); % Use SI Units
cfg = [];
cfg.method = 'projectmesh';
cfg.tissue = 'scalp';
cfg.numvertices = 800;
mesh_scalp = ft_prepare_mesh(cfg, mri_segmented);
mesh_scalp = ft_convert_units(mesh_scalp, 'm'); % Use SI Units
mesh_eeg = [mesh_brain mesh_skull mesh_scalp];
%% Plot mesh
figure
ft_plot_mesh(mesh_brain, 'edgecolor', 'none', 'facecolor', 'r')
ft_plot_mesh(mesh_skull, 'edgecolor', 'none', 'facecolor', 'g')
ft_plot_mesh(mesh_scalp, 'edgecolor', 'none', 'facecolor', 'b')
alpha 0.3
view(132, 14)
print('-dpng', fullfile(fdir,'3_meshes.png'))
%% Head models (component 1)
cfg = [];
cfg.method = 'openmeeg';
cfg.conductivity = [1 1/80 1] *(1/3); % S/m
headmodel = ft_prepare_headmodel(cfg, mesh_eeg);
headmodel = ft_convert_units(headmodel, 'm'); % Use SI Units
save(fullfile(fdir,'headmodel.mat'), 'headmodel')
%% Plot head model
figure
ft_plot_headmodel(headmodel, 'facealpha', 0.5)
view(90, 0)
print('-dpng', fullfile(fdir,'4_headmodel.png'))
%% Load electrodes
elec = ft_read_sens(fullfile('path\to\data', 'standard_1020.elc'), 'senstype', 'eeg');
elec = ft_convert_units(elec, 'm'); % Use SI units
%% Plot electrodes
figure
hold on
ft_plot_sens(elec, 'elecsize', 40);
ft_plot_headmodel(headmodel, 'facealpha', 0.5);
view(90, 0)
%% Realign electrodes
% scalp indices differ for the two types of head models
scalp_index = 1;
cfg = [];
cfg.method = 'project'; % onto scalp surface
cfg.headshape = headmodel.bnd(scalp_index); % scalp surface
elec_realigned = ft_electroderealign(cfg, elec);
save(fullfile(fdir,'elec_realigned.mat'), 'elec_realigned')
%% Plot realigned electrodes
figure
hold on
ft_plot_sens(elec_realigned, 'elecsize', 40);
ft_plot_headmodel(headmodel, 'facealpha', 0.5);
view(90, 0)
print('-dpng', fullfile(fdir,'5_elec_headmodel_correct.png'))
%% Create source model
cfg = [];
cfg.headmodel = headmodel; % used to estimate extent of grid
cfg.resolution = 0.01; % a source per 0.01 m -> 1 cm
cfg.inwardshift = 0.005; % moving sources 5 mm inwards from the skull, ...
% since BEM models may be unstable here
sourcemodel = ft_prepare_sourcemodel(cfg);
save(fullfile(fdir,'sourcemodel.mat'), 'sourcemodel')
%% Plot source model
figure
hold on
ft_plot_mesh(sourcemodel, 'vertexsize', 20);
ft_plot_headmodel(headmodel, 'facealpha', 0.5)
view(90, 0)
print('-dpng', fullfile(fdir,'6_sourcemodel.png'))
%% Highlight sources inside brain
inside = sourcemodel;
outside = sourcemodel;
inside.pos = sourcemodel.pos(sourcemodel.inside, :);
outside.pos = sourcemodel.pos(~sourcemodel.inside, :);
figure
hold on
ft_plot_mesh(inside, 'vertexsize', 20, 'vertexcolor', 'red');
ft_plot_mesh(outside, 'vertexsize', 20)
ft_plot_headmodel(headmodel, 'facealpha', 0.1)
view(125, 10)
print('-dpng', fullfile(fdir,'7_sourcemodel_inside_outside.png'))
%% Estimate lead field
cfg = [];
cfg.sourcemodel = sourcemodel; %% where are the sources?
cfg.headmodel = headmodel; %% how do currents spread?
cfg.elec = elec_realigned; %% where are the sensors?
% how do sources and sensors connect?
sourcemodel_and_leadfield = ft_prepare_leadfield(cfg);
save(fullfile(fdir,'sourcemodel_and_leadfield.mat'), 'sourcemodel_and_leadfield')
%% Plot lead field
figure('units', 'normalized', 'outerposition', [0 0 0.5 0.5])
source_index = 1200; %% a superficial sources
sensory_dipole_current = 100e-9; % Am (realistic)
n_sensors = length(elec_realigned.label);
inside_sources = find(sourcemodel_and_leadfield.inside);
inside_index = inside_sources(source_index);
lead = sourcemodel_and_leadfield.leadfield{inside_index};
xs = zeros(1, n_sensors);
ys = zeros(1, n_sensors);
zs = zeros(1, n_sensors);
voltages = zeros(1, n_sensors);
titles = {'Lead field (x)' 'Lead field (y)' 'Lead field (z)'};
% get the xyz and norm
for sensor_index = 1:n_sensors
this_x = lead(sensor_index, 1);
this_y = lead(sensor_index, 2);
this_z = lead(sensor_index, 3);
this_norm = norm(lead(sensor_index, :));
xs(sensor_index) = this_x * sensory_dipole_current;
ys(sensor_index) = this_y * sensory_dipole_current;
zs(sensor_index) = this_z * sensory_dipole_current;
voltages(sensor_index) = this_norm * sensory_dipole_current;
end
% plot xyz
axes = {xs ys zs};
for axis_index = 1:3
this_axis = axes{axis_index};
subplot(1, 3, axis_index)
hold on
ft_plot_topo3d(elec_realigned.chanpos, this_axis, 'facealpha', 0.8)
if strcmp(headmodel.type, 'dipoli')
caxis([-10e-6, 10e-6])
end
c = colorbar('location', 'southoutside');
c.Label.String = 'Lead field (V)';
axis tight
ft_plot_mesh(mesh_brain, 'facealpha', 0.10);
ft_plot_sens(elec_realigned, 'elecsize', 20);
title(titles{axis_index})
plot3(sourcemodel_and_leadfield.pos(inside_index, 1), ...
sourcemodel_and_leadfield.pos(inside_index, 2), ...
sourcemodel_and_leadfield.pos(inside_index, 3), 'bo', ...
'markersize', 20, 'markerfacecolor', 'r')
end
print('-dpng', fullfile(fdir,'8_leadfield_xyz.png'))
% plot norm
figure('units', 'normalized', 'outerposition', [0 0 0.5 0.85])
hold on
ft_plot_topo3d(elec_realigned.chanpos, voltages, 'facealpha', 0.8)
if strcmp(headmodel.type, 'dipoli')
caxis([0, 10e-6])
end
c = colorbar('location', 'eastoutside');
c.Label.String = 'Lead field (V)';
axis tight
ft_plot_mesh(mesh_brain, 'facealpha', 0.10);
ft_plot_sens(elec_realigned, 'elecsize', 20);
title('Leadfield magnitude')
plot3(sourcemodel_and_leadfield.pos(inside_index, 1), ...
sourcemodel_and_leadfield.pos(inside_index, 2), ...
sourcemodel_and_leadfield.pos(inside_index, 3), 'bo', ...
'markersize', 20, 'markerfacecolor', 'r')
view(-90, 0)
print('-dpng', fullfile(fdir,'9_leadfield_norm.png'))
%% Plot lead field vectors
% source_index 300 is central, 1200 is superficial
figure
hold on
ft_plot_sens(elec_realigned, 'elecsize', 40);
ft_plot_headmodel(headmodel, 'facealpha', 0.5);
view(90, 0)
source_index = 1200;
sensory_dipole_current = 1e-3; % Am (unrealistic)
n_sensors = length(elec_realigned.label);
inside_sources = find(sourcemodel_and_leadfield.inside);
inside_index = inside_sources(source_index);
lead = sourcemodel_and_leadfield.leadfield{inside_index};
for sensor_index = 1:n_sensors
this_pos = elec_realigned.chanpos(sensor_index, :);
this_lead = lead(sensor_index, :);
this_voltage = this_lead * sensory_dipole_current;
quiver3(this_pos(1), this_pos(2), this_pos(3), ...
this_voltage(1), this_voltage(2), this_voltage(3), 'r', ...
'linewidth', 3, 'markersize', 100);
end
plot3(sourcemodel_and_leadfield.pos(inside_index, 1), ...
sourcemodel_and_leadfield.pos(inside_index, 2), ...
sourcemodel_and_leadfield.pos(inside_index, 3), 'o', ...
'markersize', 60, 'markerfacecolor', 'r')
print('-dpng', fullfile(fdir,'10_leadfield_vectors.png'))
source_reconstruction_mne_pilot.m
% https://www.fieldtriptoolbox.org/tutorial/minimumnormestimate/
close all, clear all, clc
% set(0,'DefaultFigureWindowStyle','docked')
%% Initialize FieldTrip
addpath('path\to\fieldtrip-20201117')
ft_defaults
%% Load electrode positions
fdir='path\to\jpilot';
load(fullfile(fdir, 'elec_realigned.mat'))
%% Read EEG data
data = load(fullfile(fdir, 'jpilot_ME_200Hz.mat'));
%% Convert to raw data format for FT
dataRaw = data.rawEEG;
dataRaw.label = elec_realigned.label(1:64);
%% Preprocess
cfg = [];
cfg.channel = {'all'};
cfg.demean = 'yes';
cfg.baselinewindow = [-2 0];
cfg.lpfilter = 'yes';
cfg.lpfreq = 40;
data_all = ft_preprocessing(cfg, dataRaw);
%% Split data by class
cfg = [];
cfg.trials = data_all.trialinfo == 0;
data_left = ft_redefinetrial(cfg, data_all);
cfg = [];
cfg.trials = data_all.trialinfo == 1;
data_right = ft_redefinetrial(cfg, data_all);
%% Averaging and noise-covariance estimation
cfg = [];
cfg.covariance = 'yes';
cfg.covariancewindow = [-inf 0]; % calculate covariance on timepoints before zero-time point in trials
tlckL = ft_timelockanalysis(cfg, data_left);
tlckR = ft_timelockanalysis(cfg, data_right);
save(fullfile(fdir,'tlck.mat'), 'tlckL', 'tlckR')
%% Load forward solution
load(fullfile(fdir, 'sourcemodel_and_leadfield.mat'))
load(fullfile(fdir, 'headmodel.mat'))
%% Inverse solution
cfg = [];
cfg.sourcemodel.leadfield = sourcemodel_and_leadfield;
% cfg.sourcemodel = sourcemodel;
cfg.headmodel = headmodel;
cfg.elec = elec_realigned;
cfg.method = 'mne';
cfg.mne.prewhiten = 'yes';
cfg.mne.lambda = 3;
cfg.mne.scalesourcecov = 'yes';
sourceL = ft_sourceanalysis(cfg, tlckL);
sourceR = ft_sourceanalysis(cfg, tlckR);
save(fullfile(fdir,'source.mat'), 'sourceL', 'sourceR');
%% Visualization
m = sourceL.avg.pow(:,100); % plotting the result at 500 ms after the zero time-point
ft_plot_mesh(sourceL, 'vertexcolor', m);
view([180 0]); h = light; set(h, 'position', [0 1 0.2]); lighting gouraud; material dull
%% Visualization - source movie
cfg = [];
cfg.projectmom = 'yes';
sdL = ft_sourcedescriptives(cfg, sourceL);
sdR = ft_sourcedescriptives(cfg, sourceR);
sdDIFF = sdR;
sdDIFF.avg.pow = sdR.avg.pow - sdL.avg.pow;
% save sd sdDIFF;
cfg = [];
cfg.funparameter = 'pow';
ft_sourcemovie(cfg,sdDIFF);
Output of source_reconstruction_mne_pilot.m
Warning: the data does not contain a trial definition
In fixsampleinfo at line 88
In ft_datatype_raw at line 149
In ft_checkdata at line 267
In ft_preprocessing at line 289
Warning: reconstructing sampleinfo by assuming that the trials are consecutive segments of a continuous recording
> In fixsampleinfo (line 102)
In ft_datatype_raw (line 149)
In ft_checkdata (line 267)
In ft_preprocessing (line 289)
In source_reconstruction_mne_pilot (line 29)
the call to "ft_selectdata" took 1 seconds
preprocessing
preprocessing trial 177 from 177
the call to "ft_preprocessing" took 4 seconds
the input is raw data with 64 channels and 177 trials
selecting 91 trials
the call to "ft_selectdata" took 0 seconds
the call to "ft_redefinetrial" took 0 seconds
the input is raw data with 64 channels and 177 trials
selecting 86 trials
the call to "ft_selectdata" took 0 seconds
the call to "ft_redefinetrial" took 0 seconds
the input is raw data with 64 channels and 91 trials
the call to "ft_selectdata" took 0 seconds
the call to "ft_selectdata" took 0 seconds
the call to "ft_timelockanalysis" took 1 seconds
the input is raw data with 64 channels and 86 trials
the call to "ft_selectdata" took 0 seconds
the call to "ft_selectdata" took 0 seconds
the call to "ft_timelockanalysis" took 1 seconds
the input is timelock data with 64 channels and 1901 timebins
the call to "ft_selectdata" took 0 seconds
using precomputed leadfields
Warning: the labels are missing for the precomputed leadfields, assuming that they were computed with the same channel
selection
Index exceeds the number of array elements (2).
Error in getdimord (line 468)
if isequalwithoutnans(datsiz, [npos nchan nori]) || (isequal(datsiz([1 3]), [npos nori]) && isinf(nchan))
Error in ft_selectdata (line 188)
dimord{i} = getdimord(varargin{1}, datfield{i});
Error in ft_sourceanalysis (line 435)
sourcemodel = ft_selectdata(tmpcfg, sourcemodel);
Error in source_reconstruction_mne_pilot (line 67)
sourceL = ft_sourceanalysis(cfg, tlckL);