[FieldTrip] Issue with source grand averaging

Mohd Faizal Mohd Zulkifly mfzl_mzly at yahoo.com
Sat May 2 14:51:41 CEST 2020


 Dear Eelke,
Thank you for your reply and yes you are right. Now I solved source grandaveraging issue if I include a line "source.pos= sourcemodel.grid.pos" after doing source analysis (https://mailman.science.ru.nl/pipermail/fieldtrip/2015-October/022580.html. ). I have some more questions which you might help:1) I tried to understand why do we need to include source.pos=sourcemodel.grid.pos after ft_sourceanalysis. If I excluded this line ft_sourcegrandaveraging throws an error (“Error using ft_sourcegrandaverage (line 116) the input sources vary in the field inside”) ?
2) When should I do ft_sourceinterpolation and ft_volumenormalise? Should it be on the individual subject or on the source grandaverage (average subject)? It seems that after including "source.pos= sourcemodel.grid.pos" and do  ft_sourceinterpolation and ft_volumenormalise on individual subject, no source activation anymore after ft_sourceplot of source grandaverage data. However, if I do not include the "source.pos= sourcemodel.grid.pos" , then interpolate and volume normalize individual subject, there was source activation. I am not sure this problem was due to"source.pos= sourcemodel.grid.pos" or the way I did ft_interpolation and ft_volumenormalise. Hope you can suggest something.  Here is the script:************************************************************ 
    cfg                   = [];    cfg.channel           = datall_avg.label;    cfg.elec              = elec_aligned.elec_aligned;    cfg.method            = 'lcmv';    cfg.headmodel         = headmodel.hdm;    cfg.sourcemodel.grid   = leadfield;     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 N20 latency range    ind    = find(eegsrc_act_basnrm.time>=0.016 & eegsrc_act_basnrm.time<=0.024); %at 16ms-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.pos = sourcemodel.grid.pos   %%  <<<I Include here??
 % interpolate and normalize to individual mri    %first, interpolate        cfg           = [];    cfg.parameter = 'pow';    eeg_source_interp  = ft_sourceinterpolate(cfg, source, mri);  %%< interpolate on individual     source??
    % second, normalise    cfg = [];    cfg.nonlinear     = 'no';    eeg_source_interp_norm = ft_volumenormalise(cfg, eeg_source_interp);

Thank you.Faizal















    On Friday, 1 May 2020, 04:43:20 pm GMT+2, Eelke Spaak <e.spaak at donders.ru.nl> wrote:  
 
 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

_______________________________________________
fieldtrip mailing list
https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
https://doi.org/10.1371/journal.pcbi.1002202
  
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20200502/ab68131f/attachment.htm>


More information about the fieldtrip mailing list