[FieldTrip] fieldtrip Digest, Vol 124, Issue 9

[RICHARDS, JOHN]

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)


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;

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

%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.

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,

Thanks, John



Message: 4
Date: Mon, 8 Mar 2021 08:45:55 +0000
From: "Schoffelen, J.M. (Jan Mathijs)" <jan.schoffelen at donders.ru.nl>
To: FieldTrip discussion list <fieldtrip at science.ru.nl>
Subject: Re: [FieldTrip] ft_sourceinterpolate
Message-ID: <44234160-BAFB-4F26-9CE4-96DBDDAF9FB6 at donders.ru.nl>
Content-Type: text/plain; charset="utf-8"

Dear John and Stefania,

I am not sure whether I understand what the problem is, and whether there is a problem to begin with.

I tried to simulate (part of) your use case with the below snippet of code, and to me, it seems to behave as expected.

>From what I read (both on, and between the lines) in your message (not backed up by example data from your specific use case) your issue seems to be that the interpolation of the boolean inside field of the tetra mesh onto the regular volumetric (hexahedral) grid does not behave as expected. Indeed, this could be due to the fact that the interpolation from less well-structured point clouds (i.e. the tetrahedral mesh) is done using the ’nearest’ method, which is an all-or-none phenomenon when the functional values are either 0 or 1.

What happens if you'd, rather than relying on the ’inside’ definition of the mesh, explicitly define the ‘inside’ compartment in the mri volume (as I did just before the second call to ft_sourceinterpolate)? It seems that this results in a volume with functional values that do not ‘bleed into’ the outside compartment.

Best wishes,
Jan-Mathijs





% create a volumetric sphere, at '1 mm'
mri = [];
mri.dim = [75 75 75];
mri.transform = [eye(3) ones(3,1).*-38; 0 0 0 1];



[x,y,z] = ndgrid(-37:37,-37:37,-37:37);



pos     = [x(:) y(:) z(:)];
mri.brain = (x./37).^2 + (y./37).^2 + (z./37).^2 < 1;



% downsample the sphere to '3 mm'
cfg = [];
cfg.downsample = 3;
mri2 = ft_volumedownsample(cfg, mri);



% create a 'low res' tetrahedral mesh from the sphere cfg = []; cfg.method = 'tetrahedral'; mesh = ft_prepare_mesh(cfg, mri2);



% add some 'functional' data to the mesh, this should be approximately % equal to reflect the distance to the origin.
mesh.pow = sqrt(mesh.pos(:,1).^2+mesh.pos(:,2).^2+mesh.pos(:,3).^2);



% what happens now if we try to interpolate?
mri.anatomy = double(mri.brain);



cfg = [];
cfg.parameter = 'pow';
mesh_int = ft_sourceinterpolate(cfg, mesh, mri);



cfg = [];
cfg.funparameter = 'pow';
ft_sourceplot(cfg, mesh_int);
% conclusion: without an explicitly defined 'inside' compartment the % interpolated mesh also obtains values in 'voxels' that were not in the % boolean input volume: this makes sense, because the sourcemodel (i.e. the % tetra mesh cannot make a distinction between inside/outside, because it % cannot be defined similarly to a 3D (or hex) grid, which is regular, and % has a dim (to define the reshaping)



% what happens now if we try to interpolate while adding an inside to the % mri?
mri.inside = mri.brain;



cfg = [];
cfg.parameter = 'pow';
mesh_int2 = ft_sourceinterpolate(cfg, mesh, mri);



cfg = [];
cfg.funparameter = 'pow';
ft_sourceplot(cfg, mesh_int2);


***********************************************
John E. Richards
Carolina Distinguished Professor
Department of Psychology
University of South Carolina
Columbia, SC  29208
Dept Phone: 803 777 2079
Fax: 803 777 9558
Email: richards-john at sc.edu
https://jerlab.sc.edu
*************************************************