[FieldTrip] Issue with source grand averaging

Wed Apr 29 13:14:48 CEST 2020

Dear all,

I am doing source analysis based on subject-specific MRI andcurrently having an issue to do source grandaveraging. I tried to use differentversion of fieldtrips (i.e. Fieldtrip-20181122, Fieldtrip-20191127 andFieltrip-20200109), and it still ended up with the same error: “Error usingft_sourcegrandaverage (line 116) the input sources vary in the field inside”. Ialso 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 knowmy mistakes in the scripts below. It is a long script and I was thinking itmight be helpful to provide a full script for others to help. Hope to hear fromyou guys soon. Thank you.

***********************************Forwardmodel************************************************************

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 EEGspace

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: Sensoradjustment           %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))};    % samelabels 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 inelec            

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: Sourcemodel       

% % 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 gridthat has just been created

% inverse-warp the template grid to subject specificcoordinates

cfg = [];

cfg.warpmni = 'yes';

cfg.template = template;

cfg.nonlinear = 'yes';

cfg.mri =  mri; %mri - not working properly/ showingscalp

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 thecenter (NAI tutorial)

cfg.normalizeparam      = 0.5;

cfg.grid.unit          = 'mm';

leadfield              = ft_prepare_leadfield(cfg);

%save leadfield

 % the data consists of fewer channels than theprecomputed

    % leadfields, the following chunk of codetakes 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 tosegemented mri transformation info

    trns = segmentedmri.transform /mri.transform;

    clear segmentedmri

    % reslice mri

    mri.coordsys = 'ctf';

    mri = ft_volumereslice([], mri);

    % transform mri to coordinates ofleadfield 

    mri = ft_transform_geometry(trns,mri);

    % Bring data as in time domain

    % Cut Trial data in baseline and evokedtime segments

    cfg = [];

    cfg.latency = [-1.0 0]; %baseline        % choose baseline as i didwith 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' createsaveraged 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 ofcurrent 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 spatialfilter  

    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 andpost stimulus data

   cfg.grid.filter       = eegsrc_all.avg.filter;

    % Compute source for baseline andactivty 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 onsetactivity 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 momentswithin the N1 latency range

    ind    =find(eegsrc_act_basnrm.time>=0.1 & eegsrc_act_basnrm.time<=0.5); %at18ms-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 ofmemory

    % Save source

    % interpolate and normalize toindividual 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 andnormalised data 

%% 

      % load sourcereconstructed data

%     keep subjects path_data nsubjpath_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 _sourcedata?

        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ße40

37075 Göttingen

Germany

Tel: +49 55139-8457

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20200429/55638b10/attachment.htm>