search

mne-tools / mne-bids

MNE-BIDS is a Python package that allows you to read and write BIDS-compatible datasets with the help of MNE-Python.
https://mne.tools/mne-bids/
BSD 3-Clause "New" or "Revised" License
132 stars 87 forks source link

Read fiducials from BIDS files #210

Closed jasmainak closed 5 years ago

jasmainak commented 5 years ago

See the anatomical landmarks section here. And from this you should be able to generate trans for coregistration automatically.

Based on our discussion at OHBM

jasmainak commented 5 years ago

You might have to finish the work of @monkeyman192 in #77 to make this work though. So you can do a round trip and testing is easier ...

sappelhoff commented 5 years ago

I'll tackle this on Monday :+1: it'd be really cool to generate trans from this ... without having to store any derivative beyond the three coordinates of LPA, RPA, NAS

And directly after, I'll build this in to the somato dataset PR to bring it closer to a merge: https://github.com/mne-tools/mne-python/pull/6414

jasmainak commented 5 years ago

great!

sappelhoff commented 5 years ago

I am a bit stuck - let me present the issue:

Let's look at our subject folder (somato dataset):

└── sub-01
    ├── anat
    │   └── sub-01_T1w.nii.gz
    ├── meg
    │   ├── sub-01_coordsystem.json
    │   ├── sub-01_task-somato_channels.tsv
    │   ├── sub-01_task-somato_events.tsv
    │   ├── sub-01_task-somato_meg.fif
    │   └── sub-01_task-somato_meg.json
    └── sub-01_scans.tsv

Important files for us are

sub-01_task-somato_meg.fif

In this file we have the original physical locations (3D) of Nasion, RPA, and LPA ... as digitized by some measurement device (e.g. Polhemus)

sub-01_coordsystem.json

duplication of Nasion RPA and LPA (same as in .fif file) ... see MEG section on coordsystem.json ... the relevant field is HeadCoilCoordinates

sub-01_T1w.nii.gz

This is the MRI that we want to coregister the subject's head to. Currently, it's simply a NIfTI file ... there is no information on 3D physical locations of Nasion, LPA, and RPA.

Next steps

... then we would have:

└── sub-01
    ├── anat
    │   ├── sub-01_T1w.nii.gz
    │   └── sub-01_T1w.json
    ├── meg
    │   ├── sub-01_coordsystem.json
    │   ├── sub-01_task-somato_channels.tsv
    │   ├── sub-01_task-somato_events.tsv
    │   ├── sub-01_task-somato_meg.fif
    │   └── sub-01_task-somato_meg.json
    └── sub-01_scans.tsv

and with this, we could use

... to end up with 2 sets of each LPA, RPA, NAS ...

jasmainak commented 5 years ago

we also have AnatomicalLandmarkCoordinates in the sub-01_ccordsystem.json, no?

You could use "fit fiducials" in MNE. Another way to look at this is that we want to learn a rotation matrix. I am not quite sure why MNE-Python does not do it this way. Maybe @larsoner can shed more light on it. In any case, you have the function in MNE that finds you the trans file.

You don't need to save it, MNE almost everywhere accept the trans object ...

larsoner commented 5 years ago

It's not just a rotation matrix, but rotation + translation

sappelhoff commented 5 years ago

Okay, but one question I have is still not resolved. These are the required inputs:

    src_pts : array, shape = (n, 3)
        Points to which the transform should be applied.
    tgt_pts : array, shape = (n, 3)
        Points to which src_pts should be fitted. Each point in tgt_pts should
        correspond to the point in src_pts with the same index.

... but how do I obtain the other one (src_pts or tgt_pts) from the MRI?

larsoner commented 5 years ago

Typically you do not get trans from doing this matched point business, but instead from a manual process. It should be more accurate than just using the LPA/Nasion/RPA, even once you mark them (manually or automatically) on the MRI.

The info['dev_head_t'] is often calculated from matched HPI point locations, though.

jasmainak commented 5 years ago

You have AnatomicalLandmarkCoordinates in _coordsystem.json?

sappelhoff commented 5 years ago

You have AnatomicalLandmarkCoordinates in _coordsystem.json?

yes ... I would have that for MEG data ... but I wouldn't have that for the MRI data. UNLESS somebody took the extra effort to extract these points from the MRI (using fiducials, which are objects placed on anatomical landmarks to identify them later in an MRI or similar) and then save them in the T1w.json

BUT T1w.json is optional ... to give an impression: Of all files in bids-examples, these are the T1w.json files (my point is: only very few files have been specified, because it's merely recommended, not required)

./ds000248/acq-epi_T1w.json
./hcp_example_bids/sub-100307/anat/sub-100307_T1w.json
./eeg_ds000117/sub-13/anat/sub-13_T1w.json
./eeg_ds000117/sub-04/anat/sub-04_T1w.json
./eeg_ds000117/sub-15/anat/sub-15_T1w.json
./eeg_ds000117/sub-12/anat/sub-12_T1w.json
./eeg_ds000117/sub-16/anat/sub-16_T1w.json
./eeg_ds000117/sub-06/anat/sub-06_T1w.json
./eeg_ds000117/sub-14/anat/sub-14_T1w.json
./eeg_ds000117/sub-02/anat/sub-02_T1w.json
./eeg_ds000117/sub-11/anat/sub-11_T1w.json
./eeg_ds000117/sub-08/anat/sub-08_T1w.json
./eeg_ds000117/sub-07/anat/sub-07_T1w.json
./eeg_ds000117/sub-05/anat/sub-05_T1w.json
./eeg_ds000117/sub-01/anat/sub-01_T1w.json
./eeg_ds000117/sub-10/anat/sub-10_T1w.json
./eeg_ds000117/sub-09/anat/sub-09_T1w.json
./eeg_ds000117/sub-03/anat/sub-03_T1w.json
./ieeg_visual/sub-02/ses-01/anat/sub-02_ses-01_T1w.json
./ieeg_visual/sub-01/ses-01/anat/sub-01_ses-01_T1w.json
./eeg_rest_fmri/acq-y232_T1w.json
./eeg_rest_fmri/acq-y224_T1w.json
./ieeg_visual_multimodal/sub-som682/ses-nyu3t01/anat/sub-som682_ses-nyu3t01_T1w.json
./ds000117/sub-13/ses-mri/anat/sub-13_ses-mri_acq-mprage_T1w.json
./ds000117/sub-04/ses-mri/anat/sub-04_ses-mri_acq-mprage_T1w.json
./ds000117/sub-15/ses-mri/anat/sub-15_ses-mri_acq-mprage_T1w.json
./ds000117/sub-12/ses-mri/anat/sub-12_ses-mri_acq-mprage_T1w.json
./ds000117/sub-16/ses-mri/anat/sub-16_ses-mri_acq-mprage_T1w.json
./ds000117/sub-06/ses-mri/anat/sub-06_ses-mri_acq-mprage_T1w.json
./ds000117/sub-14/ses-mri/anat/sub-14_ses-mri_acq-mprage_T1w.json
./ds000117/sub-02/ses-mri/anat/sub-02_ses-mri_acq-mprage_T1w.json
./ds000117/sub-11/ses-mri/anat/sub-11_ses-mri_acq-mprage_T1w.json
./ds000117/sub-08/ses-mri/anat/sub-08_ses-mri_acq-mprage_T1w.json
./ds000117/sub-07/ses-mri/anat/sub-07_ses-mri_acq-mprage_T1w.json
./ds000117/sub-05/ses-mri/anat/sub-05_ses-mri_acq-mprage_T1w.json
./ds000117/sub-01/ses-mri/anat/sub-01_ses-mri_acq-mprage_T1w.json
./ds000117/sub-10/ses-mri/anat/sub-10_ses-mri_acq-mprage_T1w.json
./ds000117/sub-09/ses-mri/anat/sub-09_ses-mri_acq-mprage_T1w.json
./ds000117/acq-mprage_T1w.json
./ds000117/sub-03/ses-mri/anat/sub-03_ses-mri_acq-mprage_T1w.json
./ds000117/acq-epi_T1w.json

and of these, only the ds000117 actually has the AnatomicalLandmarkCoordinates :grimacing:

sappelhoff commented 5 years ago

Typically you do not get trans from doing this matched point business, but instead from a manual process. It should be more accurate than just using the LPA/Nasion/RPA, even once you mark them (manually or automatically) on the MRI.

Mmmh okay ... we want to automate it, because in BIDS it is very convenient to specify two sets of points (from which to calculate trans) ... but it's very underspecified how to store a trans file itself.

But if quality suffers, perhaps we are on the wrong path?

jasmainak commented 5 years ago

UNLESS somebody took the extra effort to extract these points from the MRI (using fiducials, which are objects placed on anatomical landmarks to identify them later in an MRI or similar) and then save them in the T1w.json

If the trans file exists, you could reverse engineer this and update some of the datasets? At least it's possible for the sample data, somato data etc.

only the ds000117 actually has the AnatomicalLandmarkCoordinates

which is fair because nobody actually tried to use BIDS so far for doing anything real. That's why we go with the crash-testing approach. You may need to make small PRs to BIDS-examples/ to add this.

jasmainak commented 5 years ago

But if quality suffers, perhaps we are on the wrong path?

You can use the "tol" parameter in fit_matched_points perhaps. But it's true that it may not necessarily work always. Pinging @agramfort to pitch in

agramfort commented 5 years ago

why do you want to actually fit anything from the bids files? it could contain the info from the trans file. It just stores the fid points in the MRI space from which the trans is analytically obtained. For MNE users we want a way to make this points from the trans.

do I miss something @larsoner ?

larsoner commented 5 years ago

I'm not sure what the process should be here, I'm just helping people understand coordinate frames...

jasmainak commented 5 years ago

it could contain the info from the trans file.

It doesn't as far as I know. That's why we were considering using the BIDS specification for saving transforms. Which led to the discussion at OHBM that we could fit it from the points in MRI space and MEG space.

jasmainak commented 5 years ago

This came up in the context of somato data where @sappelhoff was not sure how to store the trans files ...

sappelhoff commented 5 years ago

Hearing what @agramfort says:

could contain the info from the trans file. It just stores the fid points in the MRI space from which the trans is analytically obtained. For MNE users we want a way to make this points from the trans.

I interpret this as the following:

--> If it is possible like this, it would solve my question "How do I get the landmarks in MRI space from the .nii?" .. in short: I would get them if I have a trans and the MEG landmarks ... because in that case somebody already has manually determined the MRI landmarks

--> Remaining problem (given my interpretation is correct so far) would be that calculating trans from two set of points is imprecise? Although I don't see how this can be the case, if we obtained the MRI landmarks from the trans to begin with (it's just the transformation in the other direction again?)

agramfort commented 5 years ago

yes it's what I have in mind

let me know if you need help to know where to find these info.

info['dig'] will give you LPA, RPA and NAS in head coord system and trans gives you transformation from Head to MRI

something like this:

trans = read_trans(...) trans = _get_trans(trans, 'head', 'mri) # https://github.com/mne-tools/mne-python/blob/master//mne/transforms.py#L433

then use apply_trans to dig points to get then in MRI coord.

HTH

jasmainak commented 5 years ago

calculating trans from two set of points is imprecise

I thought about this. If you have 3 points in either coordinates and they are not collinear, then it uniquely defines the transform. They would be in a plane, so there is only one set of transforms which can take them to the other coordinate frame. If you had 2 points, then you would have problems but with 3 points you are good.

sappelhoff commented 5 years ago

Current code

The handling of raw.info['dig'] feels a bit user unfriendly, e.g., I expected to be able to do something like raw.info['dig']['Nasion'], but this errors

TypeError: list indices must be integers or slices, not str

anyhow, this is my current solution:

"""Try handling trans files."""

import os.path as op

import numpy as np
import mne
from mne.datasets import somato
from mne.transforms import _get_trans, apply_trans

data_path = somato.data_path()

raw_path = op.join(data_path, 'MEG', 'somato', 'sef_raw_sss.fif')
trans_path = op.join(data_path, 'MEG', 'somato', 'sef_raw_sss-trans.fif')

raw = mne.io.read_raw_fif(raw_path)

def _extract_landmarks(raw):
    # Assuming that LPA, Nasion, RPA are always first and in that order
    meg_landmarks = np.asarray([raw.info['dig'][i]['r'] for i in range(3)])
    return meg_landmarks

meg_landmarks = _extract_landmarks(raw)

trans = mne.read_trans(trans_path)
from_to_t, trans = _get_trans(trans)

transformed_pts = apply_trans(from_to_t, meg_landmarks)

print(meg_landmarks)
print(transformed_pts)

Output

Opening raw data file /home/stefanappelhoff/mne_data/MNE-somato-data/MEG/somato/sef_raw_sss.fif...
    Range : 237600 ... 506999 =    791.189 ...  1688.266 secs
Ready.
Current compensation grade : 0
[[-7.1441993e-02  0.0000000e+00 -1.0244548e-08]
 [ 3.7252903e-09  1.0662061e-01 -9.3132257e-10]
 [ 7.4268661e-02 -2.9802322e-08 -8.3819032e-09]]
[[-0.07086861  0.02738484  0.02025573]
 [ 0.00304659  0.12651449 -0.01412104]
 [ 0.07467662  0.0222522   0.01558245]]

Is there some nice way to check this visually?

jasmainak commented 5 years ago

not that I know of. You just need to make sure that when mne-bids reads it back you are able to retrieve the original trans. That would be your test

sappelhoff commented 5 years ago

I have a couple of questions, please see my current progress:

  1. We start with getting the meg_landmarks

    • There did not seem an easy way for me to do this. What's the smartest way? I currently have: meg_landmarks = np.asarray([raw.info['dig'][i]['r'] for i in range(3)])
    • this assumes that LPA, Nasion, RPA are always the first three digPoints and always in that order, which is not robust: there must be a better way

  2. We continue with loading the transformation matrix: trans = mne.read_trans(trans_path)

    • This trans is head to MRI: <Transform | head->MRI (surface RAS)> ... transforming points on the head to the position they would have if they were on the MRI.

  3. Then we perform an intermediate step as suggested by Alex (why though?): from_to_t, trans = _get_trans(trans, fro='mri', to='head')

    • this inverts our trans (now called from_to_t), it's now from mri to head: `<Transform | MRI (surface RAS)->head>
    • Note that I am just spelling out the default values for fro and to ... should I actually edit the default values of this private function so that it stays fro='head', to='mri'?
    • why do we change naming --> why is our actual trans called from_to_t now? Previous trans gets assigned a string value now

  4. Now we want to apply the from_to_t to the meg_landmarks to obtain mri_landmarks: mri_landmarks = apply_trans(from_to_t, meg_landmarks)

    • from_to_t is in "MRI TO HEAD" format ... it seems to me we should use a "HEAD TO MRI" format (like our original trans was before the call to _get_trans) ... because after all, we are transforming MEG landmarks (head) to an MRI ... not the other way around
    • after this step, we have our two sets of points: Meg-landmarks, and MRI-landmarks that we can nicely save in BIDS

  5. Finally, we want to get back the original trans (see step 2): trans_fitted = fit_matched_points(src_pts=meg_landmarks, tgt_pts=mri_landmarks)

    • then we can do trans_reproduced = mne.transforms.Transform(fro='head', to='mri', trans=trans_fitted)
    • ... and we compare our original trans with trans_reproduced and find only minor differences (on the order e-02 to e-08)
    • is it true that I should use src_pts so that they correspond with the fro field in the mne.transforms.Transform (and vice versa for tgt_pts and the to field)? ... i.e., head=meg_landmarks, ... mri=mri_landmarks
agramfort commented 5 years ago

for 1 use:

In [10]: raw.info['dig'][2]['ident'] Out[10]: 3 (FIFFV_POINT_RPA)

In [11]: raw.info['dig'][0]['ident'] Out[11]: 2 (FIFFV_POINT_NASION)

In [12]: raw.info['dig'][1]['ident'] Out[12]: 1 (FIFFV_POINT_LPA)

for the rest, make it work and we'll make it nice after :)

jasmainak commented 5 years ago

should I actually edit the default values of this private function so that it stays fro='head', to='mri'?

I wouldn't bother. Remember you are already in a rabbit hole. You want to have datasets to test on study template, which led to somato PR, which led to the discussion on trans and that led to this.

why do we change naming --> why is our actual trans called from_to_t now?

It's confusing as I said :) But just use the function for now as Alex suggests. Otherwise, it's combinatorial branching of PRs

+1 for make it work first. Seems like you've figured out the way forward. I'm waiting for a PR ;-)