MNE-BIDS defacing

jstout211 · May 3, 2021, 4:17pm

Does the Defacing algorithm in mne_bids require freesurfer processed files for defacing?

I used a nifti file converted using AFNI to display the fiducials. The landmarks appear correct (our anatomical positions and head coils are the same - so NAS is not the anatomical nasion). When I use defacing=True, too much of the head is removed. I can import this same nifti file into freesurfer and do the same fiducial plot and defacing. The defacing algorithm works fine with the T1.mgz file as an input. Additionally if I convert the T1.mgz file to nifti, the algorithm still works, so it does not appear to be a nifti issue. As you can see the fiducial locations are the same in all images.

Code used - From mne-bids examples

>>

t1w_bids_path = write_anat(
image=t1_fname, # path to the MRI scan
bids_path=bids_path,
raw=raw, # the raw MEG data file connected to the MRI
trans=trans, # our transformation matrix
deface={‘inset’:0},
overwrite=True,
verbose=True # this will print out the sidecar file
)
anat_dir = t1w_bids_path.directory
##<<

Below are the “info” outputs from the freesurfer generated nifti and afni generated nifti. The voxel sizes appear the same, but the number of slices are not. The translation column also appears to be different.

**# Info for dataset created using >> mri_convert T1.mgz T1.nii **
mri_info T1.nii

Volume information for T1.nii
type: nii
dimensions: 256 x 256 x 256
voxel sizes: 1.000000, 1.000000, 1.000000
type: UCHAR (0)
fov: 256.000
dof: 0
xstart: -128.0, xend: 128.0
ystart: -128.0, yend: 128.0
zstart: -128.0, zend: 128.0
TR: 0.00 msec, TE: 0.00 msec, TI: 0.00 msec, flip angle: 0.00 degrees
nframes: 1
PhEncDir: UNKNOWN
FieldStrength: 0.000000
ras xform present
xform info: x_r = -1.0000, y_r = 0.0000, z_r = 0.0000, c_r = 1.0390
: x_a = 0.0000, y_a = 0.0000, z_a = 1.0000, c_a = 27.3820
: x_s = 0.0000, y_s = -1.0000, z_s = 0.0000, c_s = -13.1950
Orientation : LIA
Primary Slice Direction: coronal

voxel to ras transform:
-1.0000 0.0000 0.0000 129.0390
0.0000 0.0000 1.0000 -100.6180
0.0000 -1.0000 0.0000 114.8050
0.0000 0.0000 0.0000 1.0000

voxel-to-ras determinant -1

ras to voxel transform:
-1.0000 -0.0000 -0.0000 129.0390
-0.0000 -0.0000 -1.0000 114.8050
-0.0000 1.0000 -0.0000 100.6180
-0.0000 -0.0000 -0.0000 1.0000

##AFNI conversion to nifti. File used as input to the Freesurfer recon
mri_info APBWVFAR.nii

Volume information for APBWVFAR.nii
type: nii
dimensions: 256 x 256 x 208
voxel sizes: 1.000000, 1.000000, 1.000000
type: SHORT (4)
fov: 256.000
dof: 0
xstart: -128.0, xend: 128.0
ystart: -128.0, yend: 128.0
zstart: -104.0, zend: 104.0
TR: 0.00 msec, TE: 0.00 msec, TI: 0.00 msec, flip angle: 0.00 degrees
nframes: 1
PhEncDir: UNKNOWN
FieldStrength: 0.000000
ras xform present
xform info: x_r = -0.0000, y_r = -0.0000, z_r = 1.0000, c_r = 1.0390
: x_a = -1.0000, y_a = -0.0000, z_a = -0.0000, c_a = 27.3820
: x_s = 0.0000, y_s = -1.0000, z_s = 0.0000, c_s = -13.1950
Orientation : PIR
Primary Slice Direction: sagittal

voxel to ras transform:
-0.0000 -0.0000 1.0000 -102.9610
-1.0000 -0.0000 -0.0000 155.3820
0.0000 -1.0000 0.0000 114.8050
0.0000 0.0000 0.0000 1.0000

voxel-to-ras determinant 1

ras to voxel transform:
0.0000 -1.0000 0.0000 155.3820
0.0000 0.0000 -1.0000 114.8050
1.0000 0.0000 0.0000 102.9610
0.0000 0.0000 0.0000 1.0000

MNE-Python version: 0.22
mne-bids: 0.7
operating system: Ubuntu 18.04.1

alexrockhill · May 4, 2021, 3:01am

It might be that the inset parameter is too large as it is the default for a 256 x 256 x 256 image if yours is not that large. You could try passing deface=dict(inset=0)

mainakjas · May 4, 2021, 10:47am

Hi Jeff,

It seems very mysterious what is happening. Are you able to reproduce this issue with MNE sample dataset by following the same steps that you did? I’m not sure if any of the MNE-BIDS developers have AFNI to investigate further. If you can provide us with a file – perhaps converting the MNE sample data using AFNI, then we can dig deeper.

Freesurfer might be doing some magic coordinate transformations that we are not aware of. Did you make these plots with plot_anat in nilearn or Freesurfer? It might be worth doing everything in Python so we know exactly what is going on.

Mainak

jstout211 · May 4, 2021, 8:00pm

@alexrockhill - the example above actually used deface={inset:0}. I tried to use negative values as well, but that threw an error.

@mainakjas - I did some testing today. It doesn’t appear to be related to AFNI. I started from scratch and used dcm2niix for the dicom to nifti conversion. There was a similar result. All of the images displayed are using the plot_anat function and the code is pretty much the example code from (07. Save and load T1-weighted MRI scan along with anatomical landmarks in BIDS — mne_bids 0.7 documentation), except for the use of the nifti file.

Below is the same with the CAMCAN dataset (sub-CC110033) and approximate transform. Even though this is already defaced, the mne-bids takes a much larger chunk.

#CODE - mostly from mne-bids website:
import os
import os.path as op
import shutil
import glob

import numpy as np
import matplotlib.pyplot as plt

import nibabel as nib
from nilearn.plotting import plot_anat

import mne
from mne.datasets import sample
from mne.source_space import head_to_mri

from mne_bids import (write_raw_bids, BIDSPath, write_anat,
get_head_mri_trans, print_dir_tree)

data_path=’/fast/BIDS/test_for_MNE/test_mne_bids’
os.chdir(data_path)
subjid = ‘sub-CC110033’
raw_fname=f’{subjid}_ses-rest_task-rest.fif’

output_path = op.join(data_path, ‘bids_out’)
if not os.path.exists(output_path): os.mkdir(output_path)

raw = mne.io.read_raw_fif(raw_fname)
raw.info[‘line_freq’] = 50 # specify power line frequency as required by BIDS

sub = ‘03’
ses = ‘01’
task = ‘rest’
run = ‘01’
bids_path = BIDSPath(subject=sub, session=ses, task=task,
run=run, root=output_path)
if not os.path.exists(’./bids_out’): os.mkdir(’./bids_out’)
write_raw_bids(raw, bids_path)

Get the path to our MRI scan

t1_fname = op.join(data_path, f’{subjid}_T1w.nii.gz’)

Load the transformation matrix and show what it looks like

trans_fname = op.join(data_path, f’{subjid}-trans.fif’)

trans = mne.read_trans(trans_fname)
print(trans)

First create the BIDSPath object.

t1w_bids_path =
BIDSPath(subject=sub, session=ses, root=output_path, suffix=‘T1w’)

We use the write_anat function

t1w_bids_path = write_anat(
image=t1_fname, # path to the MRI scan
bids_path=t1w_bids_path,
raw=raw, # the raw MEG data file connected to the MRI
trans=trans, # our transformation matrix
verbose=True, # this will print out the sidecar file
overwrite=True
)
anat_dir = t1w_bids_path.directory

estim_trans = get_head_mri_trans(bids_path=bids_path)
pos = np.asarray((raw.info[‘dig’][0][‘r’],
raw.info[‘dig’][1][‘r’],
raw.info[‘dig’][2][‘r’]))

We now use the `head_to_mri` function from MNE-Python to convert MEG

coordinates to MRI scanner RAS space. For the conversion we use our

estimated transformation matrix and the MEG coordinates extracted from the

raw file. `subjects` and `subjects_dir` are used internally, to point to

the T1-weighted MRI file: `t1_mgh_fname`. Coordinates are is mm.

mri_pos = head_to_mri(pos=pos,
subject=f’{subjid}’,
mri_head_t=estim_trans,
subjects_dir=op.join(data_path, ‘SUBJECTS_DIR’))

Our MRI written to BIDS, we got `anat_dir` from our `write_anat` function

t1_nii_fname = op.join(anat_dir, f’sub-{sub}_ses-01_T1w.nii.gz’)

Plot it

fig, axs = plt.subplots(3, 1, figsize=(7, 7), facecolor=‘k’)
for point_idx, label in enumerate((‘LPA’, ‘NAS’, ‘RPA’)):
plot_anat(t1_nii_fname, axes=axs[point_idx],
cut_coords=mri_pos[point_idx, :],
title=label, vmax=360)
plt.show()

t1w_bids_path = write_anat(
image=t1_fname, # path to the MRI scan
bids_path=bids_path,
raw=raw, # the raw MEG data file connected to the MRI
trans=trans, # our transformation matrix
deface={‘inset’:0},
overwrite=True,
verbose=True # this will print out the sidecar file
)
anat_dir = t1w_bids_path.directory

Our MRI written to BIDS, we got `anat_dir` from our `write_anat` function

t1_nii_fname = op.join(anat_dir, f’sub-{sub}_ses-01_T1w.nii.gz’)

Plot it

fig, ax = plt.subplots()
plot_anat(t1_nii_fname, axes=ax, title=‘Defaced’, vmax=360)
plt.show()

alexrockhill · May 4, 2021, 8:40pm

Could you try finding the fiducials in scanner RAS by hand (e.g. using freeview) and then following option 2 on the example? That wouldn’t require any transform at all.

https://mne.tools/mne-bids/stable/auto_examples/convert_mri_and_trans.html#option-2-use-manual-landmarks-coordinates-in-scanner-ras-for-flash-image

alexrockhill · May 4, 2021, 8:41pm

i.e. fill this in with your own anatomical landmark data

ras_landmarks = \
    np.array([[-74.53102838, 19.62854953, -52.2888194],
              [-1.89454315, 103.69850925, 4.97120376],
              [72.01200673, 21.09274883, -57.53678375]]) / 1e3  # mm -> m

(the order is LPA, nasion, RPA by row)

jstout211 · May 4, 2021, 9:43pm

This does appear to work. I didn’t see it as an option, since I had an error earlier in the page. I started with a defaced image, so I had to increase the inset to prove that its working. I can recreate on the other dataset as well.

alexrockhill · May 4, 2021, 9:57pm

Nice, so then I would guess the trans file may be incorrect. If you want to continue to troubleshoot so that you can use points you found previously using a digitization system you can try working though applying the trans as in the example. If not you could just do this option and call it a day.

jstout211 · May 5, 2021, 3:47pm

I still think the transform based defacing is affected by the matrix size somehow. I didn’t readily find an example dataset in MNE python that used a t1.nii (versus a T1.mgz).

I ran through the provided example dataset on the website using: raw_fname = op.join(data_path, ‘MEG’, ‘sample’, ‘sample_audvis_raw.fif’). I then created a nifti image from the T1.mgz file. I also created a cropped version (reducing the dimension from the back end of the matrix to keep the offsets the same).

import nibabel as nb
mri = nb.load(t1_mgh_fname)
dat = mri.get_fdata().astype(np.int16)
nifti_out = nb.Nifti1Image(dat, mri.affine)
nifti_out.to_filename(’./t1w.nii’)

dat2 = dat[:-40,:-80,:-20]
nifti_cropped_out = nb.Nifti1Image(dat2, mri.affine)
nifti_cropped_out.to_filename(’./croppedT1.nii’)

All 3 images seen below (T1.mgz, T1.nii, cropped.nii) are in aligment and have the same RAS coordinates and have the same underlying affine matrix.

The identified LPA, NAS, RPA have the same reported RAS indices on the display. These are the same indices as reported on the website.

When the landmarks are written out, there is a difference in the output indices. Since the underlying images all have the same affine and voxel size, I think the transform to this space should be the same. So the output Anatomical landmarks in voxels should be the same.

Let me know if I am doing something wrong here. I have the code changes to replicate this on the example dataset, but it won’t allow me to upload a python file.

alexrockhill · May 6, 2021, 5:17am

Sorry there was a lot of info there, is there a problem, looks pretty reasonable to me.

If you have short enough code you can put it between “```” like so

import mne

mne.io.read_raw_fif('myfile-raw.fif')

mainakjas · May 6, 2021, 1:30pm

Hi Jeff,

You can share your code via gist: gist.github.com/

Ideally, if you can share a script that demos the problem in < 30 lines and runs out of the box using the sample data, we can dig deeper. Thanks a lot!

Mainak

jstout211 · May 6, 2021, 1:39pm

Hi @alexrockhill and @mainakjas - I initially thought this was a problem with the defacing. It appears that this is calculating the wrong LPA,RPA,NAS when using the transformation matrix in a dataset that doesn’t have the same matrix size as the freesurfer image. I am assuming this is why the defacing becomes over-cropped. I am using the example data and code to show that this isn’t related to our data/transform. Here is the code - a minor modification from the website code: Cropped_image_mne_bids.py · GitHub

alexrockhill · May 10, 2021, 5:07pm

Would you like any help figuring out the transforms still or did you get it? I imagine it would be nice to use those previous markers. Maybe scipy.ndimage.zoom would help change the image sizes, but I would guess that there is a better way to do it without that. If it’s in RAS, it’s probably zero-centered so you may just be able to scale in each dimension appropriately.

jstout211 · May 10, 2021, 10:09pm

@alexrockhill - You are correct - the code is zero-centering the matrix. Everything is the same between the images until calling _mri_landmarks_to_mri_voxels. This then calls the vox2ras_tkr, which calculates the transform from the dimensions of the marix. The ns = self._structarr[‘dims’][:3] * ds / 2.0 line causes the change. Since the matrix dimensions are smaller there is a difference in the transform, which causes the difference in the output locations.

def get_vox2ras_tkr(self):
    """ Get the vox2ras-tkr transform. See "Torig" here:
            https://surfer.nmr.mgh.harvard.edu/fswiki/CoordinateSystems
    """
    ds = self._structarr['delta']
    ns = self._structarr['dims'][:3] * ds / 2.0
    v2rtkr = np.array([[-ds[0], 0, 0, ns[0]],
                       [0, 0, ds[2], -ns[2]],
                       [0, -ds[1], 0, ns[1]],
                       [0, 0, 0, 1]], dtype=np.float32)
    return v2rtkr

Thanks for the help.

Topic		Replies	Views
adding montage and fiducials in bids format 💬 Support & Discussions mne-bids , montage	4	477	June 19, 2023
Mne_analysis Digest, Vol 120, Issue 39 ✉️ Mailing List Archive (read-only) list-archive	2	152	January 31, 2018
Multi-session BIDS-compatible MEG data set -- anatomical scan 💬 Support & Discussions meg , mne-bids	4	466	February 4, 2021
Accurate head surface for picking fiducials ✉️ Mailing List Archive (read-only) list-archive	2	134	December 15, 2008
MNE-BIDS 0.10 is out! 📢 Announcements mne-bids	0	325	March 14, 2022

MNE-BIDS defacing

Code used - From mne-bids examples

>>

Get the path to our MRI scan

Load the transformation matrix and show what it looks like

First create the BIDSPath object.

We use the write_anat function

We now use the head_to_mri function from MNE-Python to convert MEG

coordinates to MRI scanner RAS space. For the conversion we use our

estimated transformation matrix and the MEG coordinates extracted from the

raw file. subjects and subjects_dir are used internally, to point to

the T1-weighted MRI file: t1_mgh_fname. Coordinates are is mm.

Our MRI written to BIDS, we got anat_dir from our write_anat function

Plot it

Our MRI written to BIDS, we got anat_dir from our write_anat function

Plot it

Related topics

We now use the `head_to_mri` function from MNE-Python to convert MEG

raw file. `subjects` and `subjects_dir` are used internally, to point to

the T1-weighted MRI file: `t1_mgh_fname`. Coordinates are is mm.

Our MRI written to BIDS, we got `anat_dir` from our `write_anat` function

Our MRI written to BIDS, we got `anat_dir` from our `write_anat` function