[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-0002.html>
More information about the fieldtrip
mailing list