How to get information of voxel/vertices label in volumetric source space (setup_volume_source_space)?

💬 Support & Discussions
1
  • MNE version: 1.1.0
  • operating system: 4.19.0-16-amd64

Below code creates a volumetric source space that located under inner skull surface. I want to know how to find out all source belong to certain ROI, that is, to know the label of each source.

surface = os.path.join(subjects_dir, sub_code, 'bem', 'inner_skull.surf')
src = mne.setup_volume_source_space(subject=sub_code, subjects_dir=subjects_dir,
                                                        surface=surface）

<SourceSpaces: [<volume, shape=(29, 37, 30), n_used=10479>] MRI (surface RAS) coords, subject ‘someone’, ~72.0 MB>. So the coordinate of each source is in surface RAS coordination. It is an individual specific space.

I have tried with surface-based space, mne.read_labels_from_annot could be used to find the label and corrsponding vertices index. However, I could not find a way to use in volmatric source space.

I also tried to use label to create volumetric source space (as below) so that I could have the information of the roi and the corresponding source index:

fname_aseg_aparc = join(subjects_dir, sub_code, 'mri', 'aparc+aseg.mgz')
label_names = mne.get_volume_labels_from_aseg( fname_aseg_aparc )#fname_aseg
label_use =[label for label in label_names if label[0:3]=='ctx']
            
surface = os.path.join(subjects_dir, sub_code, 'bem', 'inner_skull.surf')
src = mne.setup_volume_source_space(
                                                subject = sub_code , subjects_dir = subjects_dir,
                                                surface = surface, bem=None ,sphere=None,   
                                                mri = fname_aseg_aparc,
                                                pos=5.0, #default
                                                volume_label= label_use, #default not specific roi but all brain
                                                   )

The src created above is a list of (more than 20, depending on the labels you choose) volumetric source space.

The problem is it is too large ( more than 2G) , could not be saved and when I plot stc using this source space, error occrus “src should be a instance of Volmatic SourceSpaces not a list”. Besides, I am not sure whether the coordinate in aparc+aseg.mgz in in surface RAS coordination and whether it is individual specific.

I wonder does volumatic source space we get has to be converted to MNI space, and then we need to find a parcellation in MNI space to get label?

May I ask what is the suggested way to get the label of each source in volumetric space?

Thanks a lot for helping!

2

Your message did a great job framing the problem; you can either make one volumetric source space that used much less memory but does not map directly to ROIs or you can make many source spaces, each within an ROI. From what I can tell, you’re asking how to get a one-to-one mapping of source voxels to ROI labels using the first approach. The issue with that is that, depending on resolution of the spacing of the source space, you might have larger or smaller source space voxels than ROI voxels in the aseg image. They are probably not going to be aligned as well so that, even if you choose the source space spacing to be 1mm to match the aseg file, the voxels will still not be one-to-one with the aseg image. What you could do is take the bounding box of each source space voxel and find all the ROI voxels that overlap it and assign it those labels. There does not currently exist an MNE function to do this as far as I know but it wouldn’t be too hard as long as you know how to convert from voxels to scanner RAS and back to voxels for the other space. I might have time to write a code snippet soon to do that but hopefully the explaination helps.

3

Thanks so much for the quick and helpful explanation!
May I ask two following question to make sure that I understand it correctly?

  1. Is it correct that the volumatric source spaces created using ROI label ( code showed below) should be in the same coordination — subject individual RAS freesurfer space— and thus, the vertices of each ROI souce space could be extract the coordinate from src[‘rr’] to plot in one single brain, and this brain is an individual specific brain that located in RAS freeesurfer coordination. (not a standard brain)
fname_aseg_aparc = join(subjects_dir, sub_code, 'mri', 'aparc+aseg.mgz')
label_names = mne.get_volume_labels_from_aseg( fname_aseg_aparc )#fname_aseg
label_use =[label for label in label_names if label[0:3]=='ctx']
            
surface = os.path.join(subjects_dir, sub_code, 'bem', 'inner_skull.surf')
src = mne.setup_volume_source_space(
                                                subject = sub_code , subjects_dir = subjects_dir,
                                                surface = surface, bem=None ,sphere=None,   
                                                mri = fname_aseg_aparc,
                                                pos=5.0, #default
                                                volume_label= label_use, #default not specific roi but all brain
                                                   )
  1. what is the relationship between the one-whole-brain volumetric source space and the volumetric souce space created using ROI labels (aparc+aseg.mgz)?

I think the answer is that they are in same coordination, subject individual RAS freesurfer space. And the differences are they have different ways to cut the space to voxels space and that is why we could directly check the overlab between voxels in one-whole-brain source and voxels in ROI source space as you suggested above.

However, I found the answer above may be wrong. When I try to plot glassbrain, I could use source of one-whole-brain to plot stc on glassbrain. However, when I try to plot stc on glass brain using source that consisted of a list of ROI sourcses space, error occored " index 1261 out of the boundary in size 12". If two source spaces are in same coordination and under same constrain( under inner skull surface), there should not be error related to out of boundary in one space not in the other.

I also tried to create single ROI source space and plot glassbrain and get no error.
It is just confusing why when using combined ROI source space, I got error, however, when using ROI source space seperately, I do not get error.

Could you please help with these questions? Thanks so much!

4

  1. All the coordinates are in the same scanner RAS (RAS in Freesurfer) but not the same surface RAS (TkReg RAS in Freesurfer). In nibabel, this is mgh_img.header.get_vox2ras or mgh_img.header.get_vox2ras_tkr.

  2. There was a lot in there but yes, they should be in the same coordinate frame. I’m not sure about the error but maybe MNE only supports one source space at once.

Maybe this code will help:

import numpy as np
import nibabel as nib
import mne
from mne.datasets import sample

data_path = sample.data_path()
subjects_dir = data_path / 'subjects'

surface = subjects_dir / 'sample' / 'bem' / 'inner_skull.surf'
vol_src = mne.setup_volume_source_space(
    subject='sample', subjects_dir=subjects_dir, surface=surface,
    pos=10, add_interpolator=False)  # just for speed!

aseg = nib.load(subjects_dir / 'sample' / 'mri' / 'aparc+aseg.mgz')
aseg_data = np.array(aseg.dataobj)

fs_lut = mne.read_freesurfer_lut()[0]
fs_lut_r = {v:k for k, v in fs_lut.items()}

lut = -1 * np.ones(vol_src[0]['shape'], dtype=np.int64, order='F')   # look up table
lut[vol_src[0]['inuse'].astype(bool).reshape(vol_src[0]['shape'], order='F')] = \
    np.arange(sum(vol_src[0]['inuse']))

spacing = np.diag(vol_src[0]['src_mri_t']['trans'])[:3] * 1000  # m -> mm
src_vox_scan_ras_t = np.dot(vol_src[0]['mri_ras_t']['trans'],
                            vol_src[0]['src_mri_t']['trans'])
src_vox_scan_ras_t[:3] *= 1e3

labels = dict()
for i, vox in enumerate(lut.flatten(order='F')):
    if vox < 0:
        continue
    coord = np.unravel_index(i, lut.shape)
    labels[coord] = set()
    mri_coord = mne.transforms.apply_trans(
        aseg.header.get_ras2vox(),
        mne.transforms.apply_trans(src_vox_scan_ras_t, coord))
    mri_coord = mri_coord.round().astype(int)
    xyz_ranges = list()
    for i in range(3):
        my_range = np.arange(aseg.shape[i])
        my_range = my_range[(my_range >= mri_coord[i] - spacing[i] / 2) &
                          (my_range <= mri_coord[i] + spacing[i] / 2)]
        xyz_ranges.append(my_range)
    for x in xyz_ranges[0]:
        for y in xyz_ranges[1]:
            for z in xyz_ranges[2]:
                labels[coord].add(fs_lut_r[aseg_data[x, y, z]])

Well the code could be cleaned up a bit but this might be a nice function to add to MNE if you want to make a PR. Things should work very similarly for the source vertices within a label but you don’t need labels for those since you already have them…