[FieldTrip] Issue with source grand averaging
Eelke Spaak
e.spaak at donders.ru.nl
Fri May 1 16:36:08 CEST 2020
Dear Faizal,
I notice that you are creating a sourcemodel called 'grid' in your
code, but later on in the call to ft_prepare_leadfield you are not
using that sourcemodel, instead 'cfg.grid.pos = sourcemodel.pos;',
possibly using the (non-warped) template positions. This might cause
issues. I would try to make sure you are using the correct source
model/leadfield variables throughout the code, and see whether the
problem persists.
Best,
Eelke
On Wed, 29 Apr 2020 at 13:14, Mohd Faizal Mohd Zulkifly
<mfzl_mzly at yahoo.com> wrote:
>
> Dear all,
>
> I am doing source analysis based on subject-specific MRI and currently having an issue to do source grandaveraging. I tried to use different version of fieldtrips (i.e. Fieldtrip-20181122, Fieldtrip-20191127 and Fieltrip-20200109), and it still ended up with the same error: “Error using ft_sourcegrandaverage (line 116) the input sources vary in the field inside”. I also aware that this issue was discussed previously https://mailman.science.ru.nl/pipermail/fieldtrip/2015-October/022580.html. Unfortunately, I still facing the same error after specifying source.pos = template_grid.pos. It would be grateful if some of you can help and let me know my mistakes in the scripts below. It is a long script and I was thinking it might be helpful to provide a full script for others to help. Hope to hear from you guys soon. Thank you.
>
> ***********************************Forward model************************************************************
>
> subject_data= {'A’, ‘B’};
>
> subj={'A’, ‘B’};
>
> nsubj = numel(subject_data);
>
> % MRI Segmentation
>
> % if do_mri_seg
>
> % read the *mgz files from FreeSurfer
>
> fprintf(['starting forward model for: ', subj{s},'\n '])
>
> mri = ft_read_mri([subject_mri, subj{s},'/','mri/','T1.mgz']);
>
> mri.coordsys = 'ctf';
>
> % interactively set fiducials to align MRI space to EEG space
>
> cfg = [];
>
> cfg.method = 'interactive';
>
> mri = ft_volumerealign(cfg, mri);
>
>
>
> cfg = [];
>
> cfg.output = {'brain','skull','scalp'};
>
> segmentedmri = ft_volumesegment(cfg, mri);
>
> % save segmentedmri
>
>
>
> %%
>
> % Head mesh
>
>
>
> %if do_head_mesh
>
> addpath 'XXXXXX/fieldtrip-20181122/external/spm12'
>
> segmentedmri.dim = segmentedmri.dim;
>
> segmentedmri.coordsys = segmentedmri.coordsys;
>
> segmentedmri.unit = segmentedmri.unit;
>
> cfg = [];
>
> cfg.tissue = {'brain','skull','scalp'};
>
> cfg.method = 'projectmesh';
>
> cfg.numvertices = [3000 2000 1000];
>
> mri_bnd = [];
>
> mri_bnd.transform = segmentedmri.transform;
>
> mri_bnd.brain = segmentedmri.brain;
>
>
>
> braindil = imdilate(segmentedmri.brain, strel_bol(6));
>
> mri_bnd.skull = braindil > 0.5;
>
> mri_bnd.scalp = segmentedmri.scalp;
>
> mri_bnd.dim = segmentedmri.dim;
>
> mri_bnd.coordsys = segmentedmri.coordsys;
>
> mri_bnd.unit = segmentedmri.unit;
>
>
>
> src = ft_prepare_mesh(cfg, mri_bnd);
>
> % save src
>
> %%
>
> % COMPONENT 1: Headmodel
>
> addpath 'XXX/fieldtrip-20181122/external/dipoli'
>
> cfg = [];
>
> cfg.method = 'dipoli';
>
> hdm = ft_prepare_headmodel(cfg, src);
>
> % save hdm
>
> disp(hdm)
>
>
>
> %%
>
> % COMPONENT 2: Sensor adjustment % getting electrodes in the right position
>
>
>
> %if do_sensor_alignment
>
> sens = ft_read_sens([XXX, subj{s}, '.sfp']); % GEOSCAN
>
> sens = ft_convert_units(sens, 'mm');
>
>
>
> nas = mri.cfg.fiducial.nas;
>
> lpa = mri.cfg.fiducial.lpa;
>
> rpa = mri.cfg.fiducial.rpa;
>
>
>
> % warp NAS LPA and RPA
>
> nas = ft_warp_apply(mri.transform, nas, 'homogenous');
>
> lpa = ft_warp_apply(mri.transform, lpa, 'homogenous');
>
> rpa = ft_warp_apply(mri.transform, rpa, 'homogenous');
>
> fid = [];
>
> fid.elecpos = [nas; lpa; rpa]; % ctf-coordinates of fiducials
>
> fid.chanpos = [nas; lpa; rpa]; % ctf-coordinates of fiducials
>
> fid.label = {char(sens.label(1)); ...
>
> char(sens.label(2)); ...
>
> char(sens.label(3))}; % same labels as in elec
>
>
>
> % alignment
>
> cfg = [];
>
> cfg.method = 'fiducial';
>
> cfg.target = fid; % see above
>
> cfg.elec = sens;
>
> cfg.fiducial = {char(sens.label(1)); ...
>
> char(sens.label(2)); ...
>
> char(sens.label(3))}; % labels of fiducials in fid and in elec
>
> elec_aligned = ft_electroderealign(cfg);
>
>
>
> cfg = [];
>
> cfg.method = 'project';
>
> cfg.elec = elec_aligned;
>
> cfg.headshape = src(3);
>
> elec_aligned = ft_electroderealign(cfg);
>
> % save elec_aligned
>
>
>
> %%
>
> %COMPONENT 3: Source model
>
> % % load a template grid
>
> % NOTE: the path to the template file is user-specific
>
> pathin ='XXX/matlab/toolbox/';
>
> load(fullfile(pathin, 'fieldtrip-20200109/template/sourcemodel/standard_sourcemodel3d5mm'));
>
> template = sourcemodel;
>
> clear sourcemodel;
>
>
>
> % create the subject specific grid, using the template grid that has just been created
>
> % inverse-warp the template grid to subject specific coordinates
>
> cfg = [];
>
> cfg.warpmni = 'yes';
>
> cfg.template = template;
>
> cfg.nonlinear = 'yes';
>
> cfg.mri = mri; %mri - not working properly/ showing scalp
>
> cfg.unit = 'mm';
>
> grid = ft_prepare_sourcemodel(cfg);
>
> %save grid
>
>
>
> %%
>
> % COMPONENT 4: Leadfield
>
> fname_preproc_data = fullfile(XXX,subject_data{s}]);
>
> preproc_data = load(fname_preproc_data);
>
> EEG_name = cell2mat(fieldnames(preproc_data));
>
>
>
> %
>
> cfg = [];
>
> cfg.elec = elec_aligned;
>
> cfg.channel = preproc_data.conv_data_rpost.label; % ensure that rejected electrodes are not present
>
> cfg.headmodel = hdm;
>
> cfg.grid.pos = sourcemodel.pos;
>
> cfg.normalize = 'no'; %normalize 'yes' it when don't have a contrast to avaoid activity at the center (NAI tutorial)
>
> cfg.normalizeparam = 0.5;
>
> cfg.grid.unit = 'mm';
>
> leadfield = ft_prepare_leadfield(cfg);
>
> %save leadfield
>
>
>
> % the data consists of fewer channels than the precomputed
>
> % leadfields, the following chunk of code takes care of this
>
> [a,b] = match_str(preproc_data.conv_data_rpost.label, leadfield.label);
>
> for k = 1:numel(leadfield.leadfield)
>
> if ~isempty(leadfield.leadfield{k})
>
> tmp = leadfield.leadfield{k};
>
> tmp = tmp(b,:);
>
> tmp = tmp-repmat(mean(tmp,1),[size(tmp,1) 1]); % average re-ref
>
> leadfield.leadfield{k} = tmp;
>
> end
>
> end
>
> leadfield.label = leadfield.label(b);
>
> %save leadfield
>
>
>
> *****************************************************source_analysis**************************************
>
>
>
> subject_rpost_data= {‘A’, ‘B’};
>
>
>
> subj={‘A’, ‘B’ };
>
>
>
> nsubj = numel(subj);
>
>
>
>
>
> % Loop over subjects
>
>
>
> for s=1:nsubj
>
> fprintf(['starting source analysis for: ', subj{s},'\n '])
>
>
>
>
>
> load(fname_preproc_data);
>
> load(fname_sourcemodel);
>
> load(fname_leadfield);
>
> load(fname_headmodel);
>
> load(fname_elec);
>
>
>
> mri = ft_read_mri(fname_seg_mri);
>
> mri.coordsys = 'ctf';
>
>
>
> load(fname_smri )
>
> % transform original mri to segemented mri transformation info
>
> trns = segmentedmri.transform / mri.transform;
>
> clear segmentedmri
>
> % reslice mri
>
> mri.coordsys = 'ctf';
>
> mri = ft_volumereslice([], mri);
>
> % transform mri to coordinates of leadfield
>
> mri = ft_transform_geometry(trns, mri);
>
>
>
>
>
> % Bring data as in time domain
>
> % Cut Trial data in baseline and evoked time segments
>
> cfg = [];
>
> cfg.latency = [-1.0 0]; % baseline % choose baseline as i did with sensor analysis (100ms)
>
> datbasl = ft_selectdata(cfg, preproc_data.(EEG_name));
>
> cfg = [];
>
> cfg.latency = [0 3]; % evoked activity
>
> datact = ft_selectdata(cfg,preproc_data.(EEG_name));
>
> cfg = [];
>
> cfg.latency = [-1 3];
>
> datall = ft_selectdata(cfg,preproc_data.(EEG_name));
>
>
>
> %%
>
> % Compute covariance matrix
>
>
>
> cfg = [];
>
> cfg.covariance = 'yes';
>
> cfg.keeptrials = 'no'; % 'no' creates averaged trial
>
> cfg.removemean = 'yes';
>
> cfg.channel = {'all', '-E257'};
>
>
>
> datbasl_avg = ft_timelockanalysis(cfg, datbasl); %Baseline data (in time domain)
>
> datact_avg = ft_timelockanalysis(cfg, datact); %Activity data
>
> datall_avg = ft_timelockanalysis(cfg, datall); %Data with baseline and activity
>
>
>
>
>
> %% Source analysis
>
> %
>
> % adapt leadfield channels to channel of current eeg data (199 instead
>
> % of 256)
>
> [~,b] = match_str(datbasl_avg.label, leadfield.label);
>
> for k = 1:numel(leadfield.leadfield)
>
> if ~isempty(leadfield.leadfield{k})
>
> tmp = leadfield.leadfield{k};
>
> tmp = tmp(b,:);
>
> tmp = tmp-repmat(mean(tmp,1),[size(tmp,1) 1]); % average re-ref
>
> leadfield.leadfield{k} = tmp;
>
> end
>
> end
>
> leadfield.label = leadfield.label(b);
>
>
>
>
>
> % Compute common spatial filter
>
> cfg = [];
>
> cfg.channel = datall_avg.label;
>
> cfg.elec = elec_aligned.elec_aligned;
>
> cfg.method = 'lcmv';
>
> % cfg.grid = leadfield;
>
> cfg.headmodel = headmodel.hdm;
>
> cfg.sourcemodel.grid.pos = leadfield.pos; % bei FT wird hier der output von ft_leadfield angegeben
>
> cfg.grid.resolution = 5;
>
> cfg.lcmv.keepfilter = 'yes';
>
> cfg.lcmv.fixedori = 'yes';
>
> cfg.lcmv.projectnoise = 'yes';
>
> cfg.lcmv.weightnorm = 'nai';
>
> cfg.lcmv.lambda = '5%';
>
>
>
> eegsrc_all = ft_sourceanalysis(cfg, datall_avg); %source data in all time window
>
>
>
> % apply common filters to pre and post stimulus data
>
> cfg.grid.filter = eegsrc_all.avg.filter;
>
>
>
> % Compute source for baseline and activty time window
>
> eegsrc_basl = ft_sourceanalysis(cfg, datbasl_avg); %source data for baseline time window
>
> eegsrc_act = ft_sourceanalysis(cfg, datact_avg); %source data for activity time window
>
>
>
>
>
> % contrast post stimulus onset activity with respect to baseline
>
> eegsrc_act_basnrm = eegsrc_act;
>
> eegsrc_act_basnrm.avg.pow = (eegsrc_act.avg.pow - eegsrc_basl.avg.pow)./ eegsrc_basl.avg.pow;
>
>
>
>
>
> %% Choose time of interest
>
>
>
> % average across time the dipole moments within the N1 latency range
>
> ind = find(eegsrc_act_basnrm.time>=0.1 & eegsrc_act_basnrm.time<=0.5); %at 18ms-24ms
>
> timecourse = eegsrc_act_basnrm.avg.mom(eegsrc_act_basnrm.inside);
>
> tc_toi = zeros(size(eegsrc_act_basnrm.avg.pow(eegsrc_act_basnrm.inside)));
>
>
>
> for ii = 1:length(timecourse)
>
> tc_toi(ii) = mean(abs(timecourse{ii}(ind)));
>
> end
>
>
>
> source = eegsrc_act_basnrm;
>
> source.avg.pow(eegsrc_act_basnrm.inside) = tc_toi;
>
> % source.cfg = rmfield(source.cfg, {'headmodel' 'callinfo'}); % this is removed because it takes up a lot of memory
>
>
>
> % Save source
>
>
>
> % interpolate and normalize to individual mri
>
> %first, interpolate
>
> cfg = [];
>
> cfg.parameter = 'pow';
>
> eeg_source_interp = ft_sourceinterpolate(cfg, source, mri);
>
>
>
>
>
> % second, normalise
>
> cfg = [];
>
> cfg.nonlinear = 'no';
>
> eeg_source_interp_norm = ft_volumenormalise(cfg, eeg_source_interp);
>
>
>
> % save the interpolate and normalised data
>
>
>
> %%
>
> % load source reconstructed data
>
> % keep subjects path_data nsubj path_ft
>
>
>
> allsource_rpost = cell(1,nsubj);
>
> for s = 1:nsubj
>
> datapath = strcat(fname_eeg_source_interp_norm);
>
> % datafile = fullfile(datapath, [subj{s}, '_source.mat']); % take _source data?
>
> datafile = fullfile(datapath, [subj{s}, '_eeg_source_interp_norm.mat']);
>
> load(datafile)
>
>
>
> allsource_rpost{s} = eeg_source_interp_norm;
>
> % allsource_rpost{s} = source;
>
> end
>
>
>
>
>
> % grand average
>
> cfg = [];
>
> cfg.parameter = 'pow';
>
> % cfg.keepindividual='yes';
>
> gasource_rpost = ft_sourcegrandaverage(cfg, allsource_rpost{:});
>
> %save gasource_rpost
>
>
>
>
>
>
>
>
>
> With best wishes,
>
> Faizal Zulkifly
>
>
>
> --------------------------------------------------------
>
> Mohd Faizal Mohd Zulkifly
>
> University Medical Center Göttingen
>
> University of Göttingen
>
> Clinic for Clinical Neurophysiology
>
> Robert-Koch-Straße 40
>
> 37075 Göttingen
>
> Germany
>
> Tel: +49 551 39-8457
>
>
>
>
> _______________________________________________
> fieldtrip mailing list
> https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
> https://doi.org/10.1371/journal.pcbi.1002202
More information about the fieldtrip
mailing list