coregistration / alignment off? [EEG data]

[SophieHerbst]

Hello!

I am processing an EEG dataset using the pipeline (version: latest main), with dig points + individual MRIs.
I did manual coregistration and subsequently updated the anatomical landmarks in my T1-json files:

# compute the new landmarks using the transformation
mri_landmarks = mne_bids.get_anat_landmarks(
            image=bids_path_t1,
            info=info,
            trans=trans,
            fs_subject=subject_fs,
            fs_subjects_dir=subjects_dir)

print(mri_landmarks.dig)
# [<DigPoint |        LPA : (14.9, 144.1, 110.8)      : MRI voxel frame>, <DigPoint |     Nasion : (93.7, 249.4, 161.0)      : MRI voxel frame>, <DigPoint  |        RPA : (179.8, 147.0, 117.0)     : MRI voxel frame>]

# Update landmarks for existing file
 mne_bids.update_anat_landmarks(
            bids_path=bids_path_t1,
            landmarks=mri_landmarks,
            on_missing="ignore")

# Plot (I played with the parameters)
fig = plot_alignment(
            info,
            trans,
            subject=subject_fs,
            dig=True,
            eeg=["original", "projected"],
            meg=[],
            coord_frame="auto",
            subjects_dir=subjects_dir,
            show_channel_names=True,
        )

From the json file: sub-PC009_ses-01_run-1_T1w.json
"AnatomicalLandmarkCoordinates": {
"LPA": [
14.94783071409175,
144.09823694342606,
110.7827822744363
],
"NAS": [
93.72576660795733,
249.35544621088334,
161.04340943202396
],
"RPA": [
179.78230137925442,
147.03877621527408,
117.03505161340044
]
}

From the trans-file:
<Transform | MRI (surface RAS)->head>
[[ 0.99954885 0.03000787 0.00128422 -0.00128606]
[-0.02784465 0.90976769 0.41418278 -0.01740106]
[ 0.0112604 -0.41403168 0.91019285 0.02344617]
[ 0. 0. 0. 1. ]]

Yet in the report, as well as a manually created alignment plot the dig-points and fiducials look off.
Also, this output worries me:
Average distance from 62 digitized points to head: 19.34 mm

Any idea what to check? Thanks in advance!

[larsoner]

I did manual coregistration and subsequently updated the anatomical landmarks in my T1-json files:

How did you do manual coreg? It looks wrong in plot_alignment. I would make sure your plot_alignment looks good before doing anything with BIDSification. (On the plus side, the plot_alignment call looks wrong in possibly the same way as the pipeline report, so the BIDSification step is probably okay.)

I would also add to your plot_alignmnet call a mri_fiducials=True so you can see where the MRI fiducials are as well (diamond shapes).

[SophieHerbst]

Thanks Eric, I redid the coregistration using the gui and now it looks fine. Must have messed smth. up before.

[larsoner]

Okay great! I'll close in that case

[SophieHerbst]

sorry, I sent my message too early. my plots look fine now but the ones in the report still look weird. I reran the pipeline with --no-cache

[SophieHerbst]

I get the following warnings:

In my config, I set

eog_channels = ['HEOG', 'VEOG']
eeg_template_montage = "standard_1005" 

Also, I don't know if it is relevant, but the dataset has 2 sessions per participant, with eeg and mri for each session.

/Users/sh254795/Documents/PROJECTS/DELTA_NEW/.venv/bin/mne_bids_pipeline:10: RuntimeWarning: No coordinate information found for channels ['VEOG', 'HEOG']. Setting channel types to misc. To avoid this warning, set channel types explicitly.
sys.exit(main())
/Users/sh254795/Documents/PROJECTS/DELTA_NEW/.venv/bin/mne_bids_pipeline:10: RuntimeWarning: Not setting positions of 2 misc channels found in montage:
['VEOG', 'HEOG']
Consider setting the channel types to be of EEG/sEEG/ECoG/DBS/fNIRS using inst.set_channel_types before calling inst.set_montage, or omit these channels when creating your montage.
sys.exit(main())

/Users/sh254795/Documents/PROJECTS/DELTA_NEW/.venv/bin/mne_bids_pipeline:10: RuntimeWarning: Other is not an MNE-Python coordinate frame for EEG data and so will be set to 'unknown'

sys.exit(main())
/Users/sh254795/Documents/PROJECTS/DELTA_NEW/.venv/bin/mne_bids_pipeline:10: RuntimeWarning: There are channels without locations (n/a) that are not marked as bad: ['VEOG', 'HEOG']
sys.exit(main())
/Users/sh254795/Documents/PROJECTS/DELTA_NEW/.venv/bin/mne_bids_pipeline:10: RuntimeWarning: Not setting positions of 2 eog channels found in montage:
['VEOG', 'HEOG']
Consider setting the channel types to be of EEG/sEEG/ECoG/DBS/fNIRS using inst.set_channel_types before calling inst.set_montage, or omit these channels when creating your montage.
sys.exit(main())
/Users/sh254795/Documents/PROJECTS/DELTA_NEW/.venv/bin/mne_bids_pipeline:10: RuntimeWarning: Unable to map the following column(s) to MNE:
dose_mg_per_kg: 0.256756757
sys.exit(main())

[SophieHerbst]

I noted a difference between the transformation information that I recreate using the code in get_head_mri_trans (mne bids / read.py), and my manually made one (in the last column). Could that be the issue?

fs_t1_mgh = nib.load(str(fs_t1_path))
t1_nifti = nib.load(str(t1w_bids_path.fpath))

t1w_json_path = _find_matching_sidecar(
        bids_path=t1w_bids_path, extension=".json", on_error="ignore"
    )

t1w_json = json.loads(t1w_json_path.read_text(encoding="utf-8"))
mri_coords_dict = t1w_json.get("AnatomicalLandmarkCoordinates", dict())

    # landmarks array: rows: [LPA, NAS, RPA]; columns: [x, y, z]
    
mri_landmarks = np.full((3, 3), np.nan)
for landmark_name, coords in mri_coords_dict.items():
    if landmark_name.upper() == ("LPA").upper():
        mri_landmarks[0, :] = coords
    elif landmark_name.upper() == ("RPA").upper():
        mri_landmarks[2, :] = coords
    elif (
        landmark_name.upper() == ("NAS" ).upper()
        or landmark_name.lower() == ("nasion" + suffix).lower()
    ):
        mri_landmarks[1, :] = coords
    else:
        continue

# Convert to MGH format to access vox2ras method
t1_mgh = nib.MGHImage(t1_nifti.dataobj, t1_nifti.affine)

# convert to scanner RAS
mri_landmarks = apply_trans(t1_mgh.header.get_vox2ras(), mri_landmarks)

# convert to FreeSurfer T1 voxels (same scanner RAS as T1)
mri_landmarks = apply_trans(fs_t1_mgh.header.get_ras2vox(), mri_landmarks)

# now extract transformation matrix and put back to RAS coordinates of MRI
vox2ras_tkr = fs_t1_mgh.header.get_vox2ras_tkr()
mri_landmarks = apply_trans(vox2ras_tkr, mri_landmarks)
mri_landmarks = mri_landmarks * 1e-3

# EEG landmarks
pos = raw.get_montage().get_positions()
eeg_landmarks = np.asarray((pos["lpa"], pos["nasion"], pos["rpa"]))

trans_fitted = fit_matched_points(src_pts=eeg_landmarks, tgt_pts=mri_landmarks)
trans_from_json = mne.transforms.Transform(fro="head", to="mri", trans=trans_fitted)

[larsoner]

I noted a difference between the transformation information that I recreate using the code in get_head_mri_trans (mne bids / read.py), and my manually made one (in the last column). Could that be the issue?

Yes could be. Note that BIDS uses standard MRI coordinates (vox2ras in nibabel terms), whereas MNE uses FreeSurfer's MRI coordinates for "mri" (this is vox2ras_tkr in nibabel terms). These two transforms differ only by a translation. The difference between the two you posted is a translation, so it's probably related to that. By way of example (here rounding just to make the numbers more readable):

i>>> import nibabel as nib
>>> import mne
>>> mri_path = mne.datasets.sample.data_path() / "subjects" / "sample" / "mri" / "T1.mgz"
>>> img = nib.load(mri_path)
>>> img.header.get_vox2ras().round(1)
array([[  -1. ,    0. ,   -0. ,  122.7],
       [   0. ,    0. ,    1. , -119. ],
       [   0. ,   -1. ,    0. ,  100.7],
       [   0. ,    0. ,    0. ,    1. ]], dtype=float32)
>>> img.header.get_vox2ras_tkr().round(1)
array([[  -1.,    0.,    0.,  128.],
       [   0.,    0.,    1., -128.],
       [   0.,   -1.,    0.,  128.],
       [   0.,    0.,    0.,    1.]], dtype=float32)
>>> import numpy as np
>>> ras2ras_tkr = np.linalg.inv(img.header.get_vox2ras()) @ img.header.get_vox2ras_tkr()
>>> ras2ras_tkr.round(1)
array([[  1. ,   0. ,  -0. ,  -5.3],
       [  0. ,   1. ,   0. , -27.3],
       [ -0. ,  -0. ,   1. ,  -9. ],
       [  0. ,   0. ,   0. ,   1. ]], dtype=float32)

It looks like you're being careful to use the correct one in your code but it makes me wonder if it's the difference here...