[FieldTrip] please help! sourcemodel and grid don't match

Hassan Aleem ha438 at georgetown.edu
Wed Mar 30 03:19:01 CEST 2016


Hi all,

When I plot the sourcemodel and the headmodel, the sourcemodel is many
magnitudes larger than the headmodel. Not sure if it has anything to do
with my code so attaching that below. Also attaching the jpgs.

%% read the parameters

%% Parameter 2.2

mrifile='Siemens_AF8-MPRAGE_20141022_1.nii';
mritype='MNI';
realign='realigned';
volumereslicing='volumeresliced';
brainthresh=0.5;
scalpthresh=0.10;
scalpsmooth='scalp_smooth';
scalp_smooth.parameter=9;
scalpvertices=1000;
skullvertices=1000;
brainvertices=1000;
meshmethod='projectmesh';
headmodel_method='dipoli';
subjnum='915';
time='pre';
timetype='hassan';
timewindow=[0 .1];
gridtype='template10';
beamformer='lcmv';
sourcemodeltype= 'templatebased';
identifier='Parameter2-2';
===================================
%%read in mri

switch mritype
    case 'MNI'
        [mri]               = ft_read_mri(mrifile);
        mri.coordsys        ='MNI';
    case 'Analyze'
        [mri]               = ft_read_mri(mrifile);
        mri.coordsys        ='Analyze';
    case 'SPM'
        [mri]               = ft_read_mri(mrifile);
        mri.coordsys        = 'SPM';
end
save(strcat(mrifile,'_',identifier,'_','mri.mat'),'mri')
cfg=[];
ft_sourceplot(cfg,mri)
==============================

%% realigning
switch realign
    case 'realigned'
        cfg=[];
        cfg.method='fiducial';
        cfg.fiducial.nas    = [6.60688      6.30230     -2.94229];
        cfg.fiducial.lpa    = [6.04082      7.65872     0.350950];
        cfg.fiducial.rpa    = [4.41106      8.71481      3.50199];
        mri_realigned=ft_volumerealign(cfg,mri);
        ft_sourceplot(cfg,mri_realigned)

save(strcat(mrifile,'_',identifier,'_','mri_realigned.mat'),'mri_realigned')
    case 'notrealigned'
        %do nothing
end

=========================
%% volume reslicing
switch volumereslicing
    case 'volumeresliced'
    cfg = [];
    cfg.resolution = 1;
    mri_resliced = ft_volumereslice(cfg, mri);

save(strcat(mrifile,'_',identifier,'_','mri_resliced.mat'),'mri_resliced')
    case 'notresliced'
        %do nothing
end
ft_sourceplot(cfg,mri_resliced)
==============================

%% segmentation parameters
switch scalpsmooth
    case 'noscalpsmooth'
        cfg                 = [];
        cfg.brainthreshold  = brainthresh;
        cfg.scalpthreshold  = scalpthresh;
        cfg.downsample      = 1; %no downsampling
        cfg.output          = {'brain' 'skull' 'scalp'};
        seg                 = ft_volumesegment(cfg, mri);
    case 'scalp_smooth'
        cfg                 = [];
        cfg.brainthreshold  = brainthresh;
        cfg.scalpthreshold  = scalpthresh;
        cfg.scalpsmooth     = scalp_smooth.parameter;
        cfg.downsample      = 1; %no downsampling
        cfg.output          = {'brain' 'skull' 'scalp'};
        seg                 = ft_volumesegment(cfg, mri);
end

seg.anatomy         = mri.anatomy;
cfg                 = [];
cfg.funparameter    = 'trishells';
cfg.funcolormap     = [1 0 0;0 1 0;0 0 1];
cfg.maskparameter=cfg.funparameter;
ft_sourceplot(cfg,seg)

save(strcat(mrifile,scalpsmooth,'_',identifier,'_','seg.mat'),'seg')
=================================
%% prepare mesh

cfg                 = [];
cfg.method= meshmethod;
cfg.tissue          = {'scalp', 'skull', 'brain'};
cfg.numvertices     = [scalpvertices skullvertices brainvertices];
bnd                 = ft_prepare_mesh(cfg, seg);
% i enlarged / reduced the boundary mesh a little bit, because they were
causing errors, probably due to overlapping
bnd(1).pnt          = bnd(1).pnt.*1.001;
bnd(2).pnt          = bnd(2).pnt.*1;
bnd(3).pnt          = bnd(3).pnt.*0.999;

figure;
ft_plot_mesh(bnd(1))%scalp
figure;
ft_plot_mesh(bnd(2))%skull
figure;
ft_plot_mesh(bnd(3))%brain

%% save mesh
save(strcat(mrifile,cfg.method,num2str(scalpvertices),num2str(skullvertices),num2str(brainvertices),'_',identifier,'_','bnd.mat'),'bnd')
=============================================

%% build headmodel

switch headmodel_method
    case 'dipoli'
        cfg                 = [];
        cfg.method          = 'dipoli';
    case 'bemcp'
        cfg                 = [];
        cfg.method          = 'bemcp';
end

headmodel           = ft_prepare_headmodel(cfg, bnd);
headmodel           = ft_convert_units(headmodel,'cm');

save(strcat(mrifile,cfg.method,'_',identifier,'_','headmodel.mat'),'headmodel');

====================================
%% adding electrodes to headmodel
elec=ft_read_sens('Hydrocel_GSN_128_1.0_TRIM_mod.sfp');    %load electrode
file
% visualize head surface (scalp)
figure;
ft_plot_mesh(headmodel.bnd(1),
'edgecolor','none','facealpha',0.8,'facecolor',[0.6 0.6 0.8]);
hold on;
% plot electrodes
ft_plot_sens(elec,'style', 'sk')
=================================
%most likely they are not aligned, so for aligning electrodes
cfg           = [];
cfg.method    = 'interactive';%can try the 'fiducial' method but is very
%likely not to fix it
cfg.elec      = elec;
% nas=mri_realigned.cfg.fiducial.nas;
% lpa=mri_realigned.cfg.fiducial.lpa;
% rpa=mri_realigned.cfg.fiducial.rpa;
cfg.headshape = headmodel.bnd(1);%scalp
elec_aligned  = ft_electroderealign(cfg);%use the gui to manually align the
electrodes

save(strcat(mrifile,'_',identifier,'_','elec_aligned.mat'),'elec_aligned')
==================================
%% load EEG data
load(strcat('subj',subjnum,'_',time,'_FT_bothStim_rejChan500hz0.1hz-30hz_fromContinuous.mat'),
'dataM6w')
load(strcat('subj',subjnum,'_',time,'_FT_bothStim_rejChan500hz0.1hz-30hz_fromContinuous.mat'),
'dataM0')
load(strcat('subj',subjnum,'_',time,'_FT_bothStim_rejChan500hz0.1hz-30hz_fromContinuous.mat'),
'dataM3W')
load(strcat('subj',subjnum,'_',time,'_FT_bothStim_rejChan500hz0.1hz-30hz_fromContinuous.mat'),
'dataM3Bw')
load(strcat('subj',subjnum,'_',time,'_FT_bothStim_rejChan500hz0.1hz-30hz_fromContinuous.mat'),
'dataM3Bb')
load(strcat('subj',subjnum,'_',time,'_FT_bothStim_rejChan500hz0.1hz-30hz_fromContinuous.mat'),
'dataM6b')
load(strcat('subj',subjnum,'_',time,'_FT_bothStim_rejChan500hz0.1hz-30hz_fromContinuous.mat'),
'dataM6')
====================================
%% timelock data

cfg = [];
cfg.refef= 'yes'; %added at 2-5


    same = dataM0;
    same.trial = [dataM0.trial dataM3W.trial];
    same.time = [dataM0.time dataM3W.time];

    diff = dataM6;
    diff.trial = [dataM3Bb.trial dataM6.trial];
    diff.time = [dataM3Bb.time dataM6.time];

 switch timetype
     case 'simple'
        timelock_m0 = ft_timelockanalysis(cfg, dataM0);
        timelock_m6 = ft_timelockanalysis(cfg, dataM6);
        timelock_m3w = ft_timelockanalysis(cfg, dataM3W);
        timelock_m3b = ft_timelockanalysis(cfg, dataM3Bb);
        timelock_diff = ft_timelockanalysis(cfg, diff);
        timelock_same = ft_timelockanalysis(cfg, same);
     case 'hassan'
       cfg.vartrllength = 0; %was changed to 0 from 2
       cfg.covariance = 'yes';
       cfg.covariancewindow = timewindow;
       cfg.keeptrials = 'yes';
       timelock_m0 = ft_timelockanalysis(cfg, dataM0);
       timelock_m6 = ft_timelockanalysis(cfg, dataM6);
       timelock_m3w = ft_timelockanalysis(cfg, dataM3W);
       timelock_m3b = ft_timelockanalysis(cfg, dataM3Bb);
       timelock_diff = ft_timelockanalysis(cfg, diff);
       timelock_same = ft_timelockanalysis(cfg, same);

 end

save(strcat(mrifile,timetype,time,'_',identifier,'_','timelock_datam0.mat'),'timelock_m0')
save(strcat(mrifile,timetype,time,'_',identifier,'_','timelock_datam6.mat'),'timelock_m6')
save(strcat(mrifile,timetype,time,'_',identifier,'_','timelock_datam3w.mat'),'timelock_m3w')
save(strcat(mrifile,timetype,time,'_',identifier,'_','timelock_data3b.mat'),'timelock_m3b')
save(strcat(mrifile,timetype,time,'_',identifier,'_','timelock_diff.mat'),'timelock_diff')
save(strcat(mrifile,timetype,time,'_',identifier,'_','timelock_same.mat'),'timelock_same')

============================
%% make sourcemodel

switch sourcemodeltype
    case 'templatebased'
        template=load('standard_sourcemodel3d10mm.mat');
        elec=ft_read_sens('Hydrocel_GSN_128_1.0_TRIM_mod.sfp');    %load
electrode file
        cfg = [];
        cfg.headmodel = headmodel;
        cfg.mri=mri;
        cfg.grid.warpmni = 'yes';
        cfg.grid.template = template.sourcemodel;
        cfg.grid.nonlinear = 'yes';
        cfg.resolution = 2; % resolution of the template grid in mm
        cfg.moveinward = 1; % actually uses vol mesh
        cfg.grid.unit      ='mm';
        cfg.inwardshift = 0; % needs to be expressed to work with moveinward
        cfg.elec = elec;
        % cfg.headshape=('cortex_5124.surf.gii');
        mysourcemodel = ft_prepare_sourcemodel(cfg);

    case 'individual'
        load('standard_sourcemodel3d5mm.mat');
        cfg = [];
        cfg.grid.warpmni = 'yes';
        cfg.resolution = 5; % resolution of the template grid in mm
        cfg.grid.nonlinear = 'yes'; % use non linear normalization
        cfg.mri = mri; % use subject's mri in mni coordinates
        cfg.grid.unit = 'mm';
        mysourcemodel = ft_prepare_sourcemodel(cfg);
    case 'template'
        mysourcemodel=load('standard_sourcemodel3d10mm.mat');%doensn't work
end



save(strcat(mrifile,'resolution',num2str(cfg.resolution),'_',subjnum,'_',identifier,'_','sourcemodel.mat'),'mysourcemodel')
=======================
%% Plot the single subject'e head model and grid positions

figure
hold on
ft_plot_vol(headmodel,'edgecolor','none','facealpha',0.3);
ft_plot_mesh(mysourcemodel.pos(mysourcemodel.inside,:));
========================
%% Double check that source model is aligned to the segmented volume and to
the electrodes
*HERE IS WHERE I SEE THE ERROR*
figure
hold on
plot3(mysourcemodel.pos(mysourcemodel.inside,1),mysourcemodel.pos(mysourcemodel.inside,2),mysourcemodel.pos(mysourcemodel.inside,3),'.');
ft_plot_vol(headmodel,'facecolor','brain','edgecolor','none','facealpha',0.3);
ft_plot_sens(elec_aligned,'label','label','style','y*','Markersize',10);
========================================
%% make leadfield


switch gridtype
    case 'individual'
        cfg            = [];
        cfg.elec       = elec;
        cfg.headmodel  = headmodel;
        cfg.normalize  = 'yes';
        cfg.reducerank= 3;
         cfg.grid.unit = 'cm';% same unit as above, i.e. in cm
        [grid] = ft_prepare_leadfield(cfg);
    case 'templatebased'
        load('standard_sourcemodel3d10mm.mat')
        cfg                 = [];
        cfg.elec            = elec;
        cfg.template = sourcemodel.cfg.grid;
        cfg.grid=mysourcemodel.cfg.grid;
        cfg.headmodel             = headmodel;
        cfg.reducerank      = 3; % default is 3 for EEG, 2 for MEG
        cfg.channel  = 'all';
        cfg.normalize= 'yes';
        cfg.grid.unit       = 'cm';
        cfg.grid.tight      = 'yes';
        [grid] = ft_prepare_leadfield(cfg);
    case 'template10'
        load('standard_sourcemodel3d10mm.mat')
        cfg.normalize='yes';
        [grid]=sourcemodel.cfg.grid;
end


save(strcat(mrifile,gridtype,'_',identifier,'_','leadfield.mat'),'grid')
===========================
%% Interpolate
cfg = [];
cfg.interpmethod = 'linear';
cfg.parameter = 'pow';
sourceinterp_diff = ft_sourceinterpolate(cfg, sourceanalysis_diff, mri);
sourceinterp_same = ft_sourceinterpolate(cfg, sourceanalysis_same, mri);
sourceinterp_m0 = ft_sourceinterpolate(cfg, sourceanalysis_m0, mri);
sourceinterp_m6 = ft_sourceinterpolate(cfg, sourceanalysis_m6, mri);
sourceinterp_m3w = ft_sourceinterpolate(cfg, sourceanalysis_m3w, mri);
sourceinterp_m3b = ft_sourceinterpolate(cfg, sourceanalysis_m3b, mri);


%
save(strcat(mrifile,'_',identifier,'_','sourceinterp_diff.mat'),'sourceinterp_diff')
%
save(strcat(mrifile,'_',identifier,'_','sourceinterp_m0.mat'),'sourceinterp_m0')
%
save(strcat(mrifile,'_',identifier,'_','sourceinterp_m6.mat'),'sourceinterp_m6')
%
save(strcat(mrifile,'_',identifier,'_','sourceinterp_m3w.mat'),'sourceinterp_m3w')
%
save(strcat(mrifile,'_',identifier,'_','sourceinterp_m3b.mat'),'sourceinterp_m3b')
%
save(strcat(mrifile,'_',identifier,'_','sourceinterp_same.mat'),'sourceinterp_same')
=================================

%%beamforming

cfg = [];
cfg.method            = 'slice';
% cfg.slicedim= 1;
cfg.funcolormap = 'jet';
% cfg.method        = 'ortho';
% cfg.funcolorlim   = 'zeromax';
% cfg.funcolorlim   = [0.01 0.2];
cfg.funparameter= 'pow';


test2=sourceinterp_diff;
test2.pow=((sourceinterp_diff.pow-sourceinterp_same.pow)./(sourceinterp_diff.pow+sourceinterp_same.pow));
test=sourceinterp_diff;
test.pow=sourceinterp_diff.pow-sourceinterp_same.pow;
ft_sourceplot(cfg,test)
ft_sourceplot(cfg,test2)
ft_sourceplot(cfg,sourceinterp_same);
ft_sourceplot(cfg,sourceinterp_diff);
ft_sourceplot(cfg,sourceinterp_m0);
ft_sourceplot(cfg,sourceinterp_m6);
ft_sourceplot(cfg,sourceinterp_m3w);
ft_sourceplot(cfg,sourceinterp_m3b);
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20160329/99d7d215/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: figforfilegrid.jpg
Type: image/jpeg
Size: 37427 bytes
Desc: not available
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20160329/99d7d215/attachment.jpg>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: figforemail.jpg
Type: image/jpeg
Size: 28140 bytes
Desc: not available
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20160329/99d7d215/attachment-0001.jpg>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: combinedforemail.jpg
Type: image/jpeg
Size: 41474 bytes
Desc: not available
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20160329/99d7d215/attachment-0002.jpg>


More information about the fieldtrip mailing list