[FieldTrip] Inner matrix mismatch in source interpolate

[Herring, J.D. (Jim)]

Hi Russ,

 

I ran into a similar error which was avoided by following the steps in the
extended beamformer tutorial
(http://fieldtrip.fcdonders.nl/tutorial/beamformingextended). The problem
seems to rather be that ft_sourceinterpolate wrongly detects your
functional data as being 2 dimensional (if I remember correctly) because a
.dim field is not present. I did not have the time yet to find-out if it
is a problem with the other beamformer tutorial or an actual bug.

 

In any case, if I compare the steps in the extended beamformer tutorial
with the steps you took I noticed that in the extended beamformer tutorial
a sourcemodel (ft_prepare_sourcemodel) is calculated after calculating the
headmodel, which is then used to calculate the leadfield
(ft_prepare_leadfield), followed by using this leadfield in
ft_sourceanalysis. Could you try and see if that works? 

 

Best,

 

Jim

 

From: fieldtrip-bounces at science.ru.nl
[mailto:fieldtrip-bounces at science.ru.nl] On Behalf Of Russell G Port
Sent: vrijdag 15 maart 2013 0:18
To: FieldTrip discussion list
Subject: [FieldTrip] Inner matrix mismatch in source interpolate

 

Hi Fieldtrippers,

 

When running ft source interpolate on my CTF MEG data, I get an error
where it reaches line 237; tmp    = interpmat*tmp. The error is an inner
matrix dimension mismatch. Here interpmat is  a 16777216(256x256x256) x642
matrix, (mri voxels X leaderfield interpolation matrix). tmp is a 1 * 642
matrix of power at the leaderfield coordinates. If this is just matrix
math, then I assume the problem is just that the matrix have to be mxn and
nxp. Is there any reason I cannot just use tmp=tmp' to make the inner
dimensions match?

 

Cheers Russ

 

---------- Forwarded message ----------
From: Russell G Port <russgport at gmail.com>
Date: Sat, Mar 9, 2013 at 8:38 PM
Subject: source interpolate - tmp = interpmat*tmp; inner matrix must agree
To: FieldTrip discussion list <fieldtrip at science.ru.nl>


Hi Fieldtrippers,

I have a slight issue that I was hoping someone can help with. I have some
CTF data, that I was trying to perform DICS beamforming and then source
interpolation on. I am trying to beamform to the source of gamma frequency
in the right hemisphere, comparing pre and post stimuli.  Heres the basic
script:

People can skip the first part of code, I included it as possibly useful
information for this plead 

%%%%% Data has already been preprocessed and timelocked%%%%%%

     load('processed_data_500Hz.mat')
     load('timelock_data_500Hz.mat')

%%%% following the beamforming tutorial %%%%%%%

%%%%%%%%%%%%%split and combine%%%%%%%%%%%%%  
     
    cfg        = [];                                           
    cfg.toilim = baseline;  
    cfg.channel=channel
    total_bsl   = ft_redefinetrial(cfg, data);
      
    
    cfg        = [];                                           
    cfg.toilim = interest;                       
    cfg.channel=channel
    total_interest   = ft_redefinetrial(cfg, data);
    
    cfg      = [];
    total_cmb = ft_appenddata(cfg, total_bsl, total_interest);
     
    total_cmb.trialinfo = [zeros(length(total_bsl.trial), 1);
ones(length(total_interest.trial), 1)];


    %%%%%%%%%%%%%cross spectrial density for dics and tf for common
filter%%%%%%%%%%%%%%
    
cfg = [];
cfg.method    = 'mtmfft';
cfg.output    = 'powandcsd';
cfg.tapsmofrq = 25;
cfg.foilim    = 55;
freqPre_total = ft_freqanalysis(cfg, total_bsl);

cfg = [];
cfg.method    = 'mtmfft';
cfg.output    = 'powandcsd';
cfg.tapsmofrq = 25;
cfg.foilim    = 55;
freqPost_total = ft_freqanalysis(cfg, total_interest);

cfg = [];
cfg.method    = 'mtmfft';
cfg.output    = 'powandcsd';
cfg.tapsmofrq = 25;
cfg.foilim    = 55;
freqALL_total = ft_freqanalysis(cfg, total_cmb);


%%%%%%%%%%%%%%%%%%%%%%%% calc leadfield %%%%%%%%%%%%%%%%%%%%%%%%%%%

% 
T=ft_read_mri('R128_V2.mri')

       cfg.fiducial.nas  = T.hdr.fiducial.mri.nas
    cfg.fiducial.lpa  = T.hdr.fiducial.mri.lpa
    cfg.fiducial.rpa  = T.hdr.fiducial.mri.rpa
    cfg.method         = 'fiducial'
    cfg.coordsys       = 'ctf' 
   
    [mri] = ft_volumerealign(cfg, T)
 
cfg          = [];
cfg.coordsys = 'ctf'; % the MRI is expressed in the CTF coordinate system,
see below
segmentedmri = ft_volumesegment(cfg, mri);

cfg        = [];
cfg.method = 'singleshell';
hdm        = ft_prepare_headmodel(cfg, segmentedmri);

 vol = ft_convert_units(hdm, 'cm');

%%%%%%%%%%%%%%%


cfg                 = [];
cfg.grad            = freqPost_total.grad;
cfg.vol          = vol;
cfg.reducerank      = 2;
cfg.channel         = 'MR'
[grid] = ft_prepare_leadfield(cfg);



Upto this point everything works. My volume conduction model, grads, and
grid looks good when verified visually...

I then compute the common filter and apply it to the beamforming of each
condition


%%%%%%%%%%%%%%%%%%%%% computer common filter%%%%%%%%%%%

cfg              = [];
cfg.method       = 'dics';
cfg.grad              = freqALL_total.grad;
cfg.keeptrials        = 'yes';
cfg.frequency    = freqALL_total.freq;
cfg.grid         = grid;
cfg.vol          = vol;
cfg.dics.projectnoise = 'yes';
cfg.dics.lambda       = '5%';
cfg.dics.keepfilter   = 'yes';
cfg.dics.realfilter   = 'yes';
cfg.channel         = 'MR'
sourceAll = ft_sourceanalysis(cfg, freqALL_total);

%%%%%%%%%%%%%%%%%% source for bsl + interest %%%%%%%%%%%%%

cfg.grid.filter = sourceAll.avg.filter;
cfg.grid         = grid;
sourcePre_con  = ft_sourceanalysis(cfg, freqPre_total );
sourcePost_con = ft_sourceanalysis(cfg, freqPost_total);


%%%%%%%%%% Calc difference %%%%%%%%%%%%%%%%%

sourceDiff = sourcePost_con;
sourceDiff.avg.pow = (sourcePost_con.avg.pow - sourcePre_con.avg.pow) ./
sourcePre_con.avg.pow;



Heres where it falls apart, because I try to source interpolate


%%%%%%%%%%%%%%%%%%%%% graph time %%%%%%%%%%%%%%%%%%

% slicemri = ft_volumereslice([], mri);
% % mri = ft_volumereslice([],segmentedmri) %% didnt appear to help with
the current problem

cfg            = [];
cfg.downsample = 2;
cfg.parameter  = 'avg.pow';
cfg.coordsys     = 'ctf';
cfg.channel         = 'MR'
source_diffInt  = ft_sourceinterpolate(cfg, sourceDiff , mri); &&& I also
tried segmentedmri, but no difference

It works until it reachs line 237; tmp    = interpmat*tmp; where there is
the matrixs are mismatch (interpmat is  a 16777216(256x256x256) x642
matrix, and tmp is a 1 * 642 matrix. From my basic skills I gather that
tmp is sourceDiff.avg.pow (the power difference for each leadfield point).
interpmat is the interpolation matrix return by

interpmat = interp_ungridded(functional.pos (position of the leadfield
locations), pos (mri voxels), 'projmethod', cfg.interpmethod (nearest),
'sphereradius', cfg.sphereradius) (0.05);

what I get out of this is a 16777216(256x256x256 voxels) x642 (leadfield
position) matrix. it then undergoes  interpmat(~anatomical.inside(:), :) =
0. My anatomical.inside is all voxels (all of the 256 x 256 x 256 matrix=
1), as a result of the defaults in ft_checkdata.

When I try to find the minimum value in this matrix I get (via
min(interpmat) =  All zero sparse: 1-by-642

I think what is supposed to happen is that the interpmat is a matrix that
relates meg to mri space, and as such then mutlipling with a vector that
represents the power in meg space, I can see how that would get points of
power in mri space. Though this must be wrong since the math is not
allowed.... I cannot see where I went off the beaten track. It seems as if
other people have not had this problem. 

Sorry for the bother,

Russ Port

Graduate Student @ Upenn


P.S. thanks to anyone who tries to help in advance...









 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20130315/6b19f2bb/attachment-0002.html>