Skip to content

Commit

Permalink
[pre-commit.ci] auto fixes from pre-commit.com hooks
Browse files Browse the repository at this point in the history
for more information, see https://pre-commit.ci
  • Loading branch information
pre-commit-ci[bot] committed May 17, 2024
1 parent ff83b19 commit 887bb23
Showing 1 changed file with 48 additions and 22 deletions.
70 changes: 48 additions & 22 deletions examples/forward/source_space_custom_atlas.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
We showcase on the sample dataset how to apply the Yeo atlas during source reconstruction.
You should replace the atlas with your own atlas and your own subject.
Any atlas can be used instead of Yeo, provided each region contains a single label (ie: no probabilistic atlas).
Any atlas can be used instead of Yeo, provided each region contains a single label (ie: no probabilistic atlas).
.. warning:: This tutorial uses FSL and FreeSurfer to perform MRI coregistrations. If you use a different software, replace the coregistration function appropriately.
"""
Expand All @@ -21,19 +21,26 @@

# %%

import mne
import mne.datasets
import os

import nilearn.datasets
import os

import mne
import mne.datasets

# The atlas is in a template space. We download here as an example Yeo 2011's atlas, which is in the MNI152 1mm template space.
# Replace this part with your atlas and the template space you used.

nilearn.datasets.fetch_atlas_yeo_2011() # Download Yeo 2011
yeo_path = os.path.join(nilearn.datasets.get_data_dirs()[0], "yeo_2011", "Yeo_JNeurophysiol11_MNI152")
atlas_path = os.path.join(yeo_path, "Yeo2011_7Networks_MNI152_FreeSurferConformed1mm.nii.gz")
atlas_template_T1_path = os.path.join(yeo_path, "FSL_MNI152_FreeSurferConformed_1mm.nii.gz")
nilearn.datasets.fetch_atlas_yeo_2011() # Download Yeo 2011
yeo_path = os.path.join(
nilearn.datasets.get_data_dirs()[0], "yeo_2011", "Yeo_JNeurophysiol11_MNI152"
)
atlas_path = os.path.join(
yeo_path, "Yeo2011_7Networks_MNI152_FreeSurferConformed1mm.nii.gz"
)
atlas_template_T1_path = os.path.join(
yeo_path, "FSL_MNI152_FreeSurferConformed_1mm.nii.gz"
)

# The participant's T1 data. Here, we consider the sample dataset
# The brain should be skull stripped. After freesurfer preprocessing, you can either use brain.mgz or antsdn.brain.mgz
Expand All @@ -47,7 +54,7 @@
assert os.path.exists(atlas_template_T1_path)
assert os.path.exists(T1_participant_path)

#%%
# %%
# The first step is to put the atlas in subject space.
# We show this step with FSL and freesurfer with linear coregistration. If your atlas is already in participant space,
# you can skip this step. Coregistration is done in two steps: compute the atlas template to subject T1 transform and apply this transform
Expand All @@ -56,46 +63,63 @@
# FSL does not know how to read .mgz, so we need to convert the T1 to nifti format
# With FreeSurfer:
T1_participant_nifti = T1_participant_path.replace("mgz", "nii.gz")
os.system("mri_convert {} {}".format(T1_participant_path, T1_participant_nifti))
os.system(f"mri_convert {T1_participant_path} {T1_participant_nifti}")

# Compute template to subject anatomical transform
template_to_anat_transform_path = os.path.join(mri_path, "template_to_anat.mat")
os.system("flirt -in {} -ref {} -out {} -omat {}".format(atlas_template_T1_path, T1_participant_nifti, os.path.join(mri_path, "T1_atlas_coreg"), template_to_anat_transform_path))
os.system(
"flirt -in {} -ref {} -out {} -omat {}".format(
atlas_template_T1_path,
T1_participant_nifti,
os.path.join(mri_path, "T1_atlas_coreg"),
template_to_anat_transform_path,
)
)

# Apply the transform to the atlas
atlas_participant = os.path.join(mri_path, "yeo_atlas.nii.gz")
os.system("flirt -in {} -ref {} -out {} -applyxfm -init {} -interp nearestneighbour".format(atlas_path, T1_participant_nifti, atlas_participant, template_to_anat_transform_path))
os.system(
f"flirt -in {atlas_path} -ref {T1_participant_nifti} -out {atlas_participant} -applyxfm -init {template_to_anat_transform_path} -interp nearestneighbour"
)

# Convert resulting atlas from nifti to mgz
# The filename must finish with aseg, to indicate to MNE that it is a proper atlas segmentation.
atlas_converted = atlas_participant.replace(".nii.gz", "aseg.mgz")
os.system("mri_convert {} {}".format(atlas_participant, atlas_converted))
os.system(f"mri_convert {atlas_participant} {atlas_converted}")

#%%
# %%
# With the atlas in participant space, we're still missing one ingredient.
# We need a dictionary mapping label to region ID / value in the fMRI.
# In FreeSurfer and atlases, these typically take the form of lookup tables.
# You can also build the dictionary by hand.

from mne._freesurfer import read_freesurfer_lut
atlas_labels = read_freesurfer_lut(os.path.join(yeo_path, "Yeo2011_7Networks_ColorLUT.txt"))[0]

atlas_labels = read_freesurfer_lut(
os.path.join(yeo_path, "Yeo2011_7Networks_ColorLUT.txt")
)[0]
print(atlas_labels)

# Drop the key corresponding to outer region
del atlas_labels["NONE"]

#%%
# %%
# For the purpose of source reconstruction, let's create a volumetric source estimate and source reconstruction with e.g eLORETA.
from mne.minimum_norm import apply_inverse, make_inverse_operator

vol_src = mne.setup_volume_source_space("sample", subjects_dir=subjects_mri_dir,
surface=os.path.join(subject_mri_path, "bem", "inner_skull.surf"))
vol_src = mne.setup_volume_source_space(
"sample",
subjects_dir=subjects_mri_dir,
surface=os.path.join(subject_mri_path, "bem", "inner_skull.surf"),
)

fif_path = os.path.join(data_path, "MEG", "sample")
fname_trans = os.path.join(fif_path, "sample_audvis_raw-trans.fif")
raw_fname = os.path.join(fif_path, "sample_audvis_filt-0-40_raw.fif")

model = mne.make_bem_model(subject="sample", subjects_dir=subjects_mri_dir, ico=4, conductivity=(0.33,))
model = mne.make_bem_model(
subject="sample", subjects_dir=subjects_mri_dir, ico=4, conductivity=(0.33,)
)
bem_sol = mne.make_bem_solution(model)

info = mne.io.read_info(raw_fname)
Expand Down Expand Up @@ -155,7 +179,9 @@
verbose=True,
)

#%%
# %%
# Then, we can finally use our atlas!
label_tcs = stc.extract_label_time_course(labels=(atlas_converted,atlas_labels),src=vol_src)
label_tcs.shape
label_tcs = stc.extract_label_time_course(
labels=(atlas_converted, atlas_labels), src=vol_src
)
label_tcs.shape

0 comments on commit 887bb23

Please sign in to comment.