How to do convert volumetric atlas to surface using nilearn?

Software Support
, ,
1

Summary of what happened:

I would like to convert a volumetric atlas to a surface based atlas using the new functionalities from nilearn 0.11.1. Even though I use interpolation='nearest' the array after volume-to-surface conversion contains invalid ids (numbers must range from 0-48)

Command used (and if a helper script was used, a link to the helper script or the command generated):

import numpy as np
from nilearn.image import load_img
from nilearn.surface import SurfaceImage
from nilearn import datasets
from nilearn.plotting import view_surf

# Load a volumetric atlas (example: Harvard-Oxford cortical atlas)
atlas = datasets.fetch_atlas_harvard_oxford('cort-maxprob-thr25-1mm')
atlas_img = load_img(atlas.maps)

# sanity check print number of unique labels
atlas_img_data = atlas_img.get_fdata()
uniques,counts = np.unique(atlas_img_data,return_counts=True)
print(f"Volumetric atlas has {len(counts)} unique ids")

# Load the fsaverage mesh(s)
fsaverage_meshes = datasets.load_fsaverage()

# Question: How to convert volumetric to surface with only valid atlas
# region ids? 
atlas_img_surf = SurfaceImage.from_volume(
    mesh=fsaverage_meshes["pial"],
    volume_img=atlas_img,
    interpolation='nearest'
)

# check if it worked
for hemi, data in atlas_img_surf.data.parts.items():
    uniques,counts = np.unique(data,return_counts=True)
    print(f"{hemi} hemi of surface atlas has {len(counts)} unique ids")

EDIT: The conversion must only be “good enough” as I am only interested in visualization (I don’t want to use the converted atlas for analyses).

2

n_samples=1 seems to go in the right direction?

3

Thx for raising. It seems that the nearest neighbor values are recombined at some point, so that we get non-integer values. This is a bug in _nearest_voxel_sampling.
Will open an issue.
Bertrand

4

@bthirion : Yes, but isn’t this simply because of averaging? It seems that even though user chooses nearest neighbor interpolation (which just replaces missing values with existing area ids), the function still averages in the next step which will “mistakenly” produce these non-existing values. I guess the whole function was built primarily having statistical images in mind. But you’re the expert here if you’re saying this is a bug, I won’t disagree :wink:

n_samples int or None, default=None

How many samples are drawn around each vertex and averaged. If None, use a reasonable default for the chosen sampling strategy (20 for ‘ball’ or 10 for ‘line’). For performance reasons, if using kind =”ball”, choose n_samples in [10, 20, 40, 80, 160] (defaults to 20 if None is passed), because cached positions are available.

5

Agreed, but this behavior is indeed problematic.

1 Like