enh: read bioimage suite mgrid

[smg211]

No description provided.

[StolkArjen]

Hey Robert/JM,

This is a new functionality, for reading in electrode coordinates highlighted in BioImage Suite (and stored in so-called '.mgrid' files). There are a couple of idiosyncrasies. First, Bioimage Suite does not save out the electrode positions in voxel-coordinates (ijk) but in mm's (xyz). We solve this by converting from xyz2ijk in ft_read_sens, after reading the mgrid file with read_bioimage_mgrid (the user needs to specify the transformation matrix of the corresponding CT/MR to correct for any voxel scaling).

Second, the same transformation matrix is used to warp from ijk2ras. However, the issue here is that we do not know what RAS coordinate system exactly this entails (e.g. tal, mni, etc.)*. Accordingly, we decided to have ft_read_sens spit out sens.coordsys = 'ras' (without 'transform' as input, sens.coordsys = 'xyz'). Would this cause any issues later on (since 'ras' is not a ft coordsys), or should we alternatively output is as 'tal' anyway? At some point the brain scan is fused with for instance the Freesurfer MR, and the electrodes are warped along to end up in the same space as the Freesurfer mesh. In other words, the two scans that are fused need to be in Talairach space anyway. What are your thoughts on this?

*All we could find about Bioimage coordinates systems:
pages 39, 43, and 44 of http://bioimagesuite.yale.edu/manual/guide/bioimagesuite_manual_tcm501-95522_tcm501-284-32.pdf

[StolkArjen]

Or should we label it 'spm', see ft_determine_coordsys:

% interactively determine origin
origin = ' ';
while ~any(strcmp(origin, {'a', 'i', 'n'}))
origin = input('Is the origin of the coordinate system at the a(nterior commissure), i(nterauricular), n(ot a landmark)? ', 's');
end

if origin=='a' && strcmp(orientation, 'ras')
coordsys = 'spm';
elseif origin=='i' && strcmp(orientation, 'als')
coordsys = 'ctf';
elseif origin=='i' && strcmp(orientation, 'ras')
coordsys = 'neuromag'; % also used for itab
else
% just use the orientation
coordsys = orientation;
end

[robertoostenveld]

Why should these idiosyncrasies be represented in ft_read_sens? Most electrode (or MEG gradiometer) file formats are in xyz Cartesian mm space. So there seems to be no difference there.

If you want to convert from some form of head coordinates to voxel coordinates, you would in general do

elec_vox = ft_transform_sens(inv(mri.transform), elec_head);
and subsequently apply other transformations from vox to other coordinates in the same way. Coordinate handling is confusing, that is why it should remain explicit to the FT user.

I prefer general solutions (which may require extra documentation on the wiki) over specific ones.

[StolkArjen]

We would agree if this transformation didnt require an additional scaling step using the same transformation matrix (or using the xsize, ysize, zsize variables in th hdr), see ft read sens. That scaling step still needs to go somewhere. Second, what coordsys would you give the (un)transformed elecs? 'xyz' and 'spm'?

…

On Feb 21, 2017, at 2:37 AM, Robert Oostenveld ***@***.***> wrote: Why should these idiosyncrasies be represented in ft_read_sens? Most electrode (or MEG gradiometer) file formats are in xyz Cartesian mm space. So there seems to be no difference there. If you want to convert from some form of head coordinates to voxel coordinates, you would in general do elec_vox = ft_transform_sens(inv(mri.transform), elec_head); and subsequently apply other transformations from vox to other coordinates in the same way. Coordinate handling is confusing, that is why it should remain explicit to the FT user. I prefer general solutions (which may require extra documentation on the wiki) over specific ones. — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or mute the thread.

[robertoostenveld]

ras~=tal
ras~=mni
ras~=neuromag

All three on the right are examples of ras, but all are different. If the coordsys is not properly defined, you can just leave it out. When needed, ft_determine_coordsys will pop up and ask for it.

What is the order of the transformations? The electrode file content is in mm? Is the scaling that is specified in the electrode header done before or after the inverse transformation (translationrotationscaling) that comes from the MRI?

It might help if you think of the units of the 3x3 upper left section of the transform. Usually that is in mm-per-voxels, i.e. if you multiply it with voxels, you get mm. If you take the inverse, it is voxels-per-mm.

[StolkArjen]

The scaling (xyz2ijk) is done before the (regular, ijk2ras) transformation.

[robertoostenveld]

so the transform in the elec header is in mm-to-voxels? That also means that the output of ft_read_sens can either be head (in mm) or voxel (dimensionless), and that the transform from the MRI is not needed inside the function.

There is now already the option in ft_read_sens
% 'coordsys' = string, 'head' or 'dewar' (default = 'head')

Please add the option 'voxel'.
If coordsys=voxel, then the output is in voxels. If coordsys=head (which is the default) it is in mm.

The voxel2ras transform (where ras is again in mm) is done afterward with ft_transform_sens.

[StolkArjen]

This is a weird situation maybe: bioimage suite saves out the electrode positions in xyz, which is a scaled version (in case the scan resolution is not 1) of the ijk coordinates. We believe that it takes that scaling from the CT/MR (we figured through reverse engineering). However, the mgrid itself doesn't contain any information about it: it's meant to be read in with bioimage suite and not any other software. So to get from 'head' (xyz) to 'voxel' (ijk), we multiply the electrode positions with .hdr.xsize etc. Then, it can be further transformed to ras head coordinates.

[StolkArjen]

Oh, the issue is: where do we do the 'head' (xyz) to 'voxel' (ijk) transformation (using hdr.xsize, hdr.ysize, hdr.zsize)?

[StolkArjen]

.. and the transform in the ct is vox2ras (there's no xyz2ijk transformation matrix - probably bioimage suite makes this on-the-fly when reading in the mgrid file, perhaps using the hdr.xsize etc).

[StolkArjen]

A proposed order (which could be an 'example script', or maybe a 'mgrid2ft.m'?):

% first read mgrid
sens = ft_read_sens(mgridfile)
uses: sens = read_bioimage_mgrid(mgridfile)

% then convert xyz coordinates (in mm) to ijk coordinates
sens.elecpos(:,1) = round(sens.elecpos(:,1)/abs(mri.transform(1,1))); % or use mri.hdr.xsize here
sens.elecpos(:,2) = round(sens.elecpos(:,2)/abs(mri.transform(2,2)));
sens.elecpos(:,3) = round(sens.elecpos(:,3)/abs(mri.transform(3,3)));

% bioimage suite indexes the first voxel as [0,0,0]
sens.elecpos = sens.elecpos+1;

% then warp ijk coordinates to ras coordinates
sens.elecpos = ft_transform_sens(mri.transform, sens.elecpos);
sens.chanpos = sens.elecpos;

[robertoostenveld]

I stil fail to see why you don't want to use ft_transform_sens and other existing FT functionality.

I propose
elec_head = ft_read_sens(mgridfile, 'coordsys', 'head')
or
elec_voxel = ft_read_sens(mgridfile, 'coordsys', 'voxel')

There is no external exposure of the transformation between electrodes and voxels.

And you can do
mri_ras = ft_read_mri(mrifile);
or
mri_unknown = ft_read_mri(mrifile)
cfg.coordsys = 'ras'
mri_ras = ft_volumerealign(cfg, mri_unknown);
and then
elec_ras = ft_transform_sens(elec_voxel, mri.transform);

You should only worry about representations, not about conversions. Please leave those to FT. Furthermore, please allow for mixing-and-matching. We support 10 different MRI formats and 10 different electrode formats, allowing for 100 combinations. Don't assume that bioimage suite is anything special, it is probably part of the 100 that we already support. And if it is not, we just add it to the 10 electrode formats, allowing for 11x10=110 combinations.

[robertoostenveld]

please add a test script to the PR, and send the data used in that test script. I expect an electrode file and a MRI file as data.

[StolkArjen]

elec_voxel = ft_read_sens(mgridfile, 'coordsys', 'voxel')
.. is impossible without any information from the MR/CT: the mgridfile comes in xyz ('head')

[robertoostenveld]

ok, if the electrode file does not contain information relevant to the transformation or scaling, then the output is fine as it is and ft_read_sens cannot and should not output voxel coordinates. No coordinate transformations should be done inside the function.

I pulled and cleaned up the code and will merge.

[robertoostenveld]

a test script and data would be nice to have.

[smg211]

Hi Robert,

I just uploaded a test script, mgrid file and corresponding mri to the "Pipeline" folder (in /read_mgrid_test) Arjen shared with you on Dropbox. I tried to attach them here, but github wouldn't support those filetypes in the comments (I'm new to github, so my apologies if there is a better solution that I haven't figured out yet).

Forgive me if any of this information is redundant or not relevant: You are correct that the mgrid file from BioImage Suite (BIS) does not contain the information relevant to transformation or scaling. Inside BIS, the mgrid file cannot be opened without already having an mri/ct open; within the BIS electrode editor GUI, the user has to open an mgrid file and an mri/ct in order to see any of the electrodes described in the mgrid file. Within that GUI, the mri/ct is shown in voxel (ijk) coordinates, and the electrode positions (corresponding to the mgrid file) are listed (in a separate window, where the mgrid file is loaded) in xyz (mm) coordinates. Through reverse engineering, we determined that BIS uses the voxel dimensions stored in mri.hdr.xsize (.ysize and .zsize), to scale the electrode positions from xyz (in the mgrid file) to voxel coordinates, so they can be viewed on the mri.

In addition to no transformation or scaling information in the mgrid file, it does not appear that there is a transformation matrix within the mri file that can transform the electrode positions in the mgrid file from xyz2ijk, or else the most logical solution would be to explicitly apply that transformation after reading in the mgrid file with ft_read_sens (and a subsequent transformation to ras, or any other coordinate system if desired). That's why we initially decided to do the scaling within ft_read_sens, but I agree that's not the best solution. Perhaps, it may just be best to explicitly describe (in a tutorial on the ft website) how to scale the positions from the xyz coordinates in an mgrid file to voxel coordinates. Let me know what you think.

Thanks,
Sandon

[StolkArjen]

Hey Robert,
Per the testing of the function, there's a test_ft_read_sens that seems to read the sens information from the major MEG filetypes. Would it fit there, or do we create a test_read_bioimage_mgrid function? In case of the former, could you copy over Sandon's .mgrid from the dropbox to /testdata, and add it to the 'datainfo' variable that test_ft_read_sens uses?

[robertoostenveld]

I had misplaced the test data and had to search for it. I found it, thanks.

While searching I also found a write_bioimage_mgrid.m function in my mail box as attachment. I just added that to FieldTrip under the general interface of the (new) ft_write_sens function.

mac011> git commit -a
[master 70ce59c] ENH - implemented function for writing electrode structures to disk. So far only for BioImage Suite .mgrid
3 files changed, 220 insertions(+), 1 deletion(-)
create mode 100644 fileio/ft_write_sens.m
create mode 100644 fileio/private/write_bioimage_mgrid.m

mac011> git push upstream master
Counting objects: 7, done.
Delta compression using up to 4 threads.
Compressing objects: 100% (7/7), done.
Writing objects: 100% (7/7), 2.89 KiB | 0 bytes/s, done.
Total 7 (delta 4), reused 0 (delta 0)
remote: Resolving deltas: 100% (4/4), completed with 4 local objects.
To github.com:fieldtrip/fieldtrip.git
2172eb0..70ce59c master -> master

[robertoostenveld]

I looked at the data,

>> ft_plot_sens(elec)
>> grid on

what are the two electrodes in the lower left corner?

[robertoostenveld]

I used this

mri = ft_read_mri('IR29_Bext.nii.gz');
% elec = read_bioimage_mgrid('IR29_grid.mgrid');
elec = ft_read_sens('IR29_grid.mgrid');

%% Visualize

figure
ft_plot_sens(elec);
ft_plot_axes(elec);

% this shows an intersection at [0 0 0]
ft_determine_coordsys(mri, 'interactive', false);

% plot it at the middle of the volume
location = ft_warp_apply(mri.transform, mri.dim/2);
figure
ft_plot_ortho(mri.anatomy, 'transform', mri.transform, 'location', location, 'style', 'intersect');
ft_plot_axes(mri);

The electrodes make sense, the MRI does not. Units are not correct, and the coordinate system is also not consistent with anything I am familiar with.

[robertoostenveld]

I can imagine the BIS software to be designed such that all visualization is done in voxel coordinates. In FieldTrip that would correspond to

figure
ft_plot_ortho(mri.anatomy, 'transform', eye(4), 'location', mri.dim/2, 'style', 'intersect');
elec_ijk = ft_transform_sens(inv(mri.transform), elec);
ft_plot_axes(struct('unit', 'mm'))

but that figure does not make sense either

[robertoostenveld]

Although the code is merged with master, the displacement has not been solved yet. I have opened an issue at http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=3286 and assigned it to @StolkArjen. Please follow up there.