[FieldTrip] fieldtrip Digest, Vol 124, Issue 9

[Schoffelen, J.M. (Jan Mathijs)]

Hi John,

> On 8 Mar 2021, at 13:39, RICHARDS, JOHN <RICHARDS at mailbox.sc.edu> wrote:
> 
> Jan-Mathijs
> 
> Thanks.  I don't think your answer helps the problem.
> 
> Here is some of the code
> load('S9401_sourcemodel_ft_grid_tetra.mat');
> gridtetra=grid;
> load('S9401_sourcemodel_1mm_ft_grid.mat');
> grid1mm_hex=grid;
> 
> Grid-tetra was done using the tetra mesh, tetra of the gm voxels
>     xgrid: [1×132 double]
>     ygrid: [1×154 double]
>     zgrid: [1×131 double]
>       dim: [132 154 131]
>       pos: [2662968×3 double]
>      unit: 'mm'
>    inside: [2662968×1 logical]
>       cfg: [1×1 struct]
> 
> (I am not sure if this was strict ft_mesh tetra.   It uses similar code, but the tetra has only been added recently.  We may have put the dim in ourselves)
> 


Hmmm, it seems I have been spending my time on something that you did not want me to spend time on. I hate it when this happens :).

I am utterly confused now. When you say ‘we may have put the dim in ourselves’, does this mean that you have DIY’ing? If that’s the case, how then can we solve this on our end? 

If this is not the case, is it then the case that what you call gridtetra is still a well-formed 3D grid? However, in that case, you have sent me down the rabbit hole, by using the word ’tetra’ in your terminology. 

Note that a ‘dim’ field can only be added without penalty (or unexpected behavior) if the grid.pos is the implicit result of a call to ndgrid(xgrid,ygrid,zgrid). In other words, the topology and the order of the dipole positions should be as if it’s scanning a regular 3D grid in an orderly way. 

> 
> grid1mm was a 1 mm hex of the gm voxels
>      pos: [2683142×3 double]
>     xgrid: [1×131 double]
>     ygrid: [1×154 double]
>     zgrid: [1×133 double]
>       dim: [131 154 133]
>    inside: [2683142×1 logical]
>      unit: 'mm'
>       cfg: [1×1 struct]
> 
> %the "pow" was found with a eLORETA inverse filter using a segmented head model and the tetra grid
> 
> 
> %This is the sourceout structure for the input source, "functional"
> sourceout.avg.inside=gridtetra.inside;
> sourceout.avg.pos=gridtetra.pos;
> gridinside=gridtetra.inside;
> nslices=1;
> sourceout.avg.pow = nan(size(gridtetra.pos,1), nslices);
> sourceout.avg.pow(gridinside,:)=pow;  
> sourceout.time=1;
> sourceout.inside = gridtetra.inside;
> sourceout.pos = gridtetra.pos;
> sourceout.method = ['average'];
> sourceout.dim=gridtetra.dim;

I don’t understand why all this handcrafting would be needed here.

> cfg =[];
> cfg.parameter = 'pow';
> cfg.interpmethod= 'spline';
> 
> 
> interpout=ft_sourceinterpolate(cfg,sourceout,grid1mm);
> %this took 37 s
> 
> The interpout has the format of the grid1mm
>         time: 1
>          dim: [131 154 133]
>    transform: [4×4 double]
>         unit: 'mm'
>       inside: [131×154×133 logical]
>          pow: [2683142×1 double]
>          pos: [2683142×3 double]
>          cfg: [1×1 struct]
> 
> 
> these are the inside sum from the tetra and 1mm grids
> disp(sum(gridtetra.inside))
>      25936
> disp(sum(grid1mm.inside))
>      855617
>>> disp(sum(interpout.inside(:)))
>       25911

So, what’s now wrong with the above?


> %take off the dim from sourceout, use "nearest"
> sourceout=rmfield(sourceout,'dim');
> cfg.interpmethod='nearest';
> interpout=ft_sourceinterpolate(cfg,sourceout,grid1mm);
> %this takes a lot longer than the 'spline' above  this takes 3 m, the first one took 37 s
> disp(sum(interpout.inside(:)))
>      855617
> It works, SORT OF, since the inside now matches the inside of the grid1mm
> HOWEVER, the interpout.pow is mostly NaN, there are only non-NaN on fewer than 855617, 15756
>>> powinside=interpout.pow(interpout.inside);
>>> disp(size(powinside));
>      855617           1
> 
>>> b=find(isnan(powinside));
>>> disp(size(b));
>      839861           1
>>> b=find(~isnan(powinside));
>>> disp(size(b));
>       15756           1
> 
> I have no clue where then 15756 came from.  They might be exact overlap of the gridtetra and grid1mm.  They plot on the MRI volume as sparse, similar to the gridtetra.

Again, what exactly is the problem? To me it seems that the following is happening: since the nearest interpolation is an all-or-none phenomenon, it may have happened that a bunch of the to-be-interpolated-upon dipole positions (in grid1mm) were actually closer to a (tetra grid) dipole position that was originally marked as a non-inside dipole position. Thus those dipoles receive an interpolated value of ’NaN’, because the outside dipoles were explicitly created (in sourceout) as NaN. I think that this is an explanation that you have already provided yourself.


> Taking out the inside
> sourceout=rmfield(sourceout,'inside');
> sourceout.avg=rmfield(sourceout.avg,'inside');
> 
> cfg.interpmethod='nearest';
> interpout=ft_sourceinterpolate(cfg,sourceout,grid1mm);
> %this one took 220 s
> %same output, lots of NaN, 15756 non NaN
>>> powinside=interpout.pow(interpout.inside);
>>> b=find(isnan(powinside));
>>> disp(size(b));
>      839861           1
> 
>>> b=find(~isnan(powinside));
>>> disp(size(b));
>       15756           1
> 
> 
> %add the dim back, so now has dim but no inside
> sourceout.dim=gridtetra.dim;
> %try with spline
> cfg.interpmethod='spline';
> interpout=ft_sourceinterpolate(cfg,sourceout,grid1mm);
> 
> %this does not work.  Can't use "regular grid" (with dim) without inside,

The above yields the same result, because the NaNs are just in the data by construction.

Long story short, I still don’t know what is the exact problem that should be solved here. 

Best wishes,
Jan-Mathijs