[FieldTrip] Combined EEG and MEG source reconstruction
"Jörn M. Horschig"
jm.horschig at donders.ru.nl
Mon Mar 11 09:47:24 CET 2013
Dear Mushfa,
I never tried what you did, however from what I know of the matter and
from what Robert wrote, I can see that you made some errors in your
script that I can point to. I can, however, neither guarantee that this
solves all or any of your problems nor that it's all correct what you
are doing. So, I can only tell you some minor things:
1) When defining trial, setting cfg.trialdef.eventtype to '?' should
actually only be done if you do not know the eventtype. Also note that
it might be better to use all your trials rather than a subset - better
signal-to-noise that way ;)
2) No reason to segment the brain twice - this operation should be done
only once and the segmented brain you base the two leadfield-matrices on
*must* be exactly the same. So you can save time by removing the
computation of vol2.
3) Specifying only one frequency might not be a good idea, because of
imperfect sampling.
4) 'singlesphere' is created specifically for MEG data. BEM models are
good for EEG, not sure how good they are for MEG, but the dipoli
implementation only works under unix (Linux/Mac) and openMEEG has to be
downloaded manually first to make it work.
5) As Robert said, units are not always properly controlled for using
FieldTrip. You can use ft_convert_units to make them match.
6) Please check again how 'grid' is defined, e.g. the .dim field
shouldn't be concatenated.
7) You do not only need to concatenate the grids, but also the data,
hope you did that ;) But I would really rather concatenated the two freq
structures you used to compute the sourcemodel/leadfield than to compute
it newly. If you concatenate you can be sure that the channel order is
persistent.
If you keep on having problems, I would suggest to use breakpoints or
/dbstop if error/ and then step through the function to see where the
error comes from, i.e. why there are NaNs or Infs in (lf*Cf*lf), and
start off with looking at if there are NaNs or Infs in lf or Cf already.
That's all I can help you with, cause that's were my knowledge ends ;)
Maybe someone else is wiser than I am. Good luck!
Best,
Jörn
On 3/8/2013 4:12 PM, Mushfa Yousuf wrote:
>
> Hello;
>
> I have few problems while calculating combine "MEG & EEG" source
> reconstruction and I need your help to troubleshoot the problems.
>
> I have a data from a Parkinson's disease patient who underwent DBS and
> we recorded 306 MEG channels and 64 EEG channels simultaneously, first
> I have calculated the leadfield individually for both EEG and MEG.
>
> I have defined all the calculations step by step so that I can
> describe my problems very precisely.
>
> *_CALCULATE THE LEADFIELD FOR MEG CHANNELS:-_*
>
> // I have defined the Pre-Stimulus data i.e. Stim-off data in the
> original code.
>
> _Reading raw data set for MEG channels_
>
> __
>
> cfg = [];
>
> cfg.dataset = 'ap_rest_stimon_bi_2_5ma.fif';
>
> cfg.trialdef.eventtype = '?';
>
> cfg.trialdef.triallength = 1;
>
> cfg.trialdef.poststim = 88;
>
> cfg = ft_definetrial(cfg);
>
> cfg.trl = cfg.trl(1:88,:); // Read 88 trials
>
> cfg.channel = {'MEG'}; *// In order to read raw data only for MEG
> channels*
>
> dataFIC12 = ft_preprocessing(cfg);
>
> After that I have defined the _time of window of interest_ for the MEG
> raw data set
>
> cfg.toilim = [0 0.5];
>
> dataPost = ft_redefinetrial(cfg, dataFIC12);
>
> and calculate the _Power Spectral density matrix_ for the MEG channels
> data set
>
> cfg = [];
>
> cfg.method = 'mtmfft';
>
> cfg.output = 'powandcsd';
>
> cfg.tapsmofrq = 2;
>
> cfg.foilim = [129 129]; // Frequency of interest is 129 Hz the
> stimulation frequency
>
> cfg.pad = 1;
>
> freqPost = ft_freqanalysis(cfg, dataPost);
>
> Now next step for me is to choose the method for the Head Model. I
> have noticed that Singlesphere is the common method which can be used
> for both EEG and MEG. So I have selected the singlesphere method for
> the calculation of volume conduction model.
>
> mri = ft_read_mri('ms1879865-0004-00001-000192-01.img'); // READ MRI
>
> cfg = [];
>
> cfg.write = 'no';
>
> cfg.coordsys = 'ctf';
>
> segmentedmri1 = ft_volumesegment(cfg, mri); // Calculate the segmented mri
>
> cfg = [];
>
> cfg.method = 'singlesphere'; // *METHOD*
>
> vol1 = ft_prepare_headmodel(cfg, segmentedmri1);
>
> Now next step is to calculate the LEADFIELD
>
> cfg = [];
>
> cfg.grad = freqPost.grad;
>
> cfg.vol = vol1;
>
> cfg.reducerank = 2; // MEG
>
> cfg.channel = {'MEG'};
>
> cfg.grid.resolution = 1; // I choose a resolution 1
>
> [grid1] = ft_prepare_leadfield(cfg);
>
> *disp(grid1)*
>
> * xgrid: [-7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7]*
>
> * ygrid: [-7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8]*
>
> * zgrid: [-4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11]*
>
> * dim: [15 16 16]*
>
> * pos: [3840x3 double]*
>
> * unit: 'cm'*
>
> * inside: [1x2037 double]*
>
> * outside: [1x1803 double]*
>
> * cfg: [1x1 struct]*
>
> * leadfield: {1x3840 cell}*
>
> **
>
> **
>
> ØIt has to be noticed that Inside Dipole is 2037.
>
> **
>
> *_CALCULATE THE LEADFIELD FOR EEG CHANNELS:-_*
>
> _Reading raw data set for EEG channels_
>
> __
>
> cfg = [];
>
> cfg.dataset = 'ap_rest_stimon_bi_2_5ma.fif';
>
> cfg.trialdef.eventtype = '?';
>
> cfg.trialdef.triallength = 1;
>
> cfg.trialdef.poststim = 88;
>
> cfg = ft_definetrial(cfg);
>
> cfg.trl = cfg.trl(1:88,:); // Read 88 trials
>
> cfg.channel = {'EEG'}; *// In order to read raw data only for EEG
> channels*
>
> dataFIC12 = ft_preprocessing(cfg);
>
> After that I have defined the _time of window of interest_ for the EEG
> raw data set
>
> cfg.toilim = [0 0.5];
>
> dataPost = ft_redefinetrial(cfg, dataFIC12);
>
> and calculate the _Power Spectral density matrix_ for the EEG channels
> data set
>
> cfg = [];
>
> cfg.method = 'mtmfft';
>
> cfg.output = 'powandcsd';
>
> cfg.tapsmofrq = 2;
>
> cfg.foilim = [129 129]; // Frequency of interest is 129 Hz
>
> cfg.pad = 1;
>
> freqPost = ft_freqanalysis(cfg, dataPost);
>
> _HEAD MODEL FOR EEG_
>
> mri = ft_read_mri('ms1879865-0004-00001-000192-01.img');
>
> cfg = [];
>
> cfg.write = 'no';
>
> cfg.coordsys = 'ctf';
>
> segmentedmri2 = ft_volumesegment(cfg, mri);
>
> cfg = [];
>
> cfg.method = 'singlesphere'; // METHOD
>
> vol2 = ft_prepare_headmodel(cfg, segmentedmri2);
>
> _LEADFIELD:___
>
> cfg = [];
>
> cfg.elec = freqPost.elec;
>
> cfg.vol = vol2;
>
> cfg.reducerank = 3; % 3 for EEG
>
> cfg.channel = {'EEG'};
>
> cfg.grid.resolution = 9.54; // Increase the resolution so that the
> inside dipoles for
>
> EEG should be close to the inside dipoles of MEG
>
> [grid2] = ft_prepare_leadfield(cfg);
>
> *disp(grid2)*
>
> * xgrid: [1x16 double]*
>
> * ygrid: [1x16 double]*
>
> * zgrid: [1x13 double]*
>
> * dim: [16 16 13]*
>
> * pos: [3328x3 double]*
>
> * unit: 'mm'*
>
> * inside: [1x2038 double]*
>
> * outside: [1x1290 double]*
>
> * cfg: [1x1 struct]*
>
> * leadfield: {1x3328 cell}*
>
> ØIt has to be noticed that Inside Dipole is 2038.
>
> *_Concatenate of MEG and EEG LEADFIELD:-_*
>
> Now I have the individual LEADFIELD of both MEG and EEG grid1 and
> grid2 respectively.
>
> *disp(grid1)*
>
> * xgrid: [-7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7]*
>
> * ygrid: [-7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8]*
>
> * zgrid: [-4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11]*
>
> * dim: [15 16 16]*
>
> * pos: [3840x3 double]*
>
> * unit: 'cm'*
>
> * inside: [1x2037 double]*
>
> * outside: [1x1803 double]*
>
> * cfg: [1x1 struct]*
>
> * leadfield: {1x3840 cell}*
>
> **
>
> *disp(grid2)*
>
> * xgrid: [1x16 double]*
>
> * ygrid: [1x16 double]*
>
> * zgrid: [1x13 double]*
>
> * dim: [16 16 13]*
>
> * pos: [3328x3 double]*
>
> * unit: 'mm'*
>
> * inside: [1x2038 double]*
>
> * outside: [1x1290 double]*
>
> * cfg: [1x1 struct]*
>
> * leadfield: {1x3328 cell}*
>
> **
>
> *Few things to be noticed here *
>
> **
>
> Ø*unit is different for both MEG 'cm' and EEG 'mm'*
>
> Ø*I have used the same method for calculating the Volume for both MEG
> & EEG => "Singlesphere" *
>
> Ø*Concatenating the "pos" vector of EEG and MEG and the other vectors*
>
> **
>
> *Grid.pos = cat(1,grid1.pos,grid2.pos);*
>
> *Grid.xgrid = [grid1.xgrid grid2.xgrid];*
>
> *Grid.ygrid = [grid1.ygrid grid2.ygrid];*
>
> *Grid.zgrid = [grid1.zgrid grid2.zgrid];*
>
> *Grid.inside = [grid1.inside grid2.inside];*
>
> *Grid.outside = [grid1.outside grid2.outside];*
>
> *Grid.leadfield = [grid1.leadfield grid2.leadfield];*
>
> *Grid.dim = [grid1.dim grid2.dim];*
>
> **
>
> *_Combined Source Reconstruction:-_*
>
> cfg = [];
>
> cfg.dataset = 'ap_rest_stimon_bi_2_5ma.fif';
>
> cfg.trialdef.eventtype = '?';
>
> cfg.trialdef.triallength = 1;
>
> cfg.trialdef.poststim = 88;
>
> cfg = ft_definetrial(cfg);
>
> cfg.trl = cfg.trl(1:88,:); // Read 88 trials
>
> cfg.channel = {'MEG','EEG'}; *// In order to read raw data for MEG &
> EEG channels*
>
> dataFIC12 = ft_preprocessing(cfg);
>
> After that as previously defined, I have calculated the Spectral density.
>
> Now I have to calculate the inverse solution for both EEG & MEG data
> trials.
>
> _Source Analysis: without contrasting condition:-_
>
> __
>
> cfg = [];
>
> cfg.method = 'dics';
>
> cfg.frequency = 129;
>
> cfg.grid = Grid; // Concatenated grid
>
> cfg.vol = vol1; // As vol1 = vol2 : both singlesphere
>
> cfg.powmethod = 'trace';
>
> cfg.dics.projectnoise = 'yes';
>
> cfg.dics.lambda = 0.0001;
>
> sourcePost_nocon = ft_sourceanalysis(cfg, freqPost); // *freqpost
> contains both MEG *
>
> *and EEG row data*
>
> ****
>
> I have received the following error:
>
> scanning grid 2037/4075
>
> ??? Error using ==> svd
>
> Input to SVD must not contain NaN or Inf.
>
> Error in ==> beamformer_dics>pinv at 568
>
> [U,S,V] = svd(A,0);
>
> Error in ==> beamformer_dics at 314
>
> filt = pinv(lf' * invCf * lf) * lf' * invCf; % Gross eqn. 3,
>
> use PINV/SVD to cover rank deficient leadfield
>
> Error in ==> ft_sourceanalysis at 595
>
> dip(i) = beamformer_dics(grid, sens, vol, [], squeeze(Cf(i,:,:)),
>
> optarg{:});
>
> Error in ==> code at 227
>
> sourcePost_nocon = ft_sourceanalysis(cfg, freqPost);
>
> ØIf you notice that till 2037 grid it represents the MEG grid and
> after that EEG grid is concatenated. i.e 2037 +2038 = 4075 (which also
> represent inside dipole).
>
> ØSecondly I try to use other method for the head modeling i.e
> openmeeg, dipoli and bemcp method which leads me to different errors
> which I was unable to trouble shoot.
>
> ·Openmeeg = OpenMEEG binaries are not correctly installed
>
> ·Dipoli =vol.mat file was missing after calculating the headmodel
>
> ·Bemcp =vol.mat file contains NaN fiels which leads to wrong
>
> leadfield computation.
>
> *_I apologize for the long email. I wanted to show the steps which I
> used so that it is easier to debug the problem._*
>
> *__*
>
> Thank you;
>
> Mushfa Yousuf
>
>
>
> On Fri, Feb 8, 2013 at 2:36 PM, Robert Oostenveld
> <r.oostenveld at donders.ru.nl <mailto:r.oostenveld at donders.ru.nl>> wrote:
>
> Dear Mushfa
>
> In principle it can be made to work, but you'll have to do a
> little bit of work outside fieldtrip. You would start with
> 1) ft_prepare_headmodel and ft_prepare_leadfield for the MEG part
> 2) ft_prepare_headmodel and ft_prepare_leadfield for the EEG part
> Subsequently you would (manually) concatenate the leadfields of
> the MEG and EEG. Here you have to make sure that the channel order
> is consistent with that in your combined data.
> The resulting structure with combined leadfields can be passed
> into ft_sourceanalysis as cfg.grid.
>
> However, I foresee a problem with the scaling. At this moment
> FieldTrip is not (yet) very precise in maintaining the correct
> units throughout. Consequently you might get different units for
> the EEG and MEG contribution to each dipole, causing weird
> cancelation effects. You'll just have to check carefully that the
> leadfield units and the data units consistently scaled between
> megmag, megplanar and eeg.
>
> This brings me to another known problem for the combined
> vectorview magnetometer and planar gradiometer channels: at this
> moment the forward computation computes both in T (or fT), whereas
> the data is represented as T for the magnetometers and T/m for the
> gradiometers. So the gradiometer leadfields are off by a factor
> "1/baseline", which is approximately 1/70 if I recall correctly.
> Many neuromag users use only the planar channels, and therefore
> don't notice. We are working on fixing the neuromag
> channel/gradiometer units. That is not as simple as it appears,
> because conversions of the geometrical dimensions from m to cm or
> mm also affect the gradiometer baseline. So converting the data
> from m to mm should result in the planar channel data getting 1000
> times smaller, whereas the magnetometer channels shoudl stay the
> same. See http://bugzilla.fcdonders.nl/show_bug.cgi?id=963 and
> http://bugzilla.fcdonders.nl/show_bug.cgi?id=963 if you want to
> follow this.
>
> best regards
> Robert
>
>
>
> On 28 Jan 2013, at 20:20, Mushfa Yousuf wrote:
>
>> Hello Fieldtrippers,
>>
>> I am currently doing source analysis using beamformer method in fieldtrip.
>>
>> The data i work on is a simultaneous measurement of MEG (306) with EEG
>> (128) on PD patients.
>>
>> I was successful in calculating the sources for both the modalities alone.
>>
>> I tried according the tutorial "Combined EEG and MEG source reconstruction" and was successful till the estimation
>>
>>
>> of the lead field. But, then when i want to use the created vol for the
>> source analysis i have two structures one for EEG and the second one for
>> MEG. Now, i can only define only one to estimate the sources. How to use
>>
>>
>> this for both the modalities? Any help would be appreciated.
>>
>> With regards
>> Mushfa
>> _______________________________________________
>> fieldtrip mailing list
>> fieldtrip at donders.ru.nl <mailto:fieldtrip at donders.ru.nl>
>> http://mailman.science.ru.nl/mailman/listinfo/fieldtrip
>
>
> _______________________________________________
> fieldtrip mailing list
> fieldtrip at donders.ru.nl <mailto:fieldtrip at donders.ru.nl>
> http://mailman.science.ru.nl/mailman/listinfo/fieldtrip
>
>
>
>
> _______________________________________________
> fieldtrip mailing list
> fieldtrip at donders.ru.nl
> http://mailman.science.ru.nl/mailman/listinfo/fieldtrip
--
Jörn M. Horschig
PhD Student
Donders Institute for Brain, Cognition and Behaviour
Centre for Cognitive Neuroimaging
Radboud University Nijmegen
Neuronal Oscillations Group
FieldTrip Development Team
P.O. Box 9101
NL-6500 HB Nijmegen
The Netherlands
Contact:
E-Mail: jm.horschig at donders.ru.nl
Tel: +31-(0)24-36-68493
Web: http://www.ru.nl/donders
Visiting address:
Trigon, room 2.30
Kapittelweg 29
NL-6525 EN Nijmegen
The Netherlands
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20130311/bc844fbd/attachment-0002.html>
More information about the fieldtrip
mailing list