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

Sreenivasan R. Nadar, Ph.D. sreenivasan.r.nadar at gmail.com
Wed Mar 30 16:39:26 CEST 2016


Hi Hasan,

The problem is probably related to dimensions (units: mm, cm) and same
could be changed by using the command: headmodel =
ft_convert_units(headmodel, 'cm').

Good luck!

Vasan

On Tue, Mar 29, 2016 at 9:19 PM, Hassan Aleem <ha438 at georgetown.edu> wrote:

> 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);
>
> _______________________________________________
> fieldtrip mailing list
> fieldtrip at donders.ru.nl
> http://mailman.science.ru.nl/mailman/listinfo/fieldtrip
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20160330/55018eab/attachment.html>


More information about the fieldtrip mailing list