Skip to content

Instantly share code, notes, and snippets.

@margulies

margulies/calculate_distance_matrix.py

Last active March 17, 2020 19:13
Show Gist options
  • Save margulies/0829b621c1d7f74bf3d04fb6dae2fc98 to your computer and use it in GitHub Desktop.
Save margulies/0829b621c1d7f74bf3d04fb6dae2fc98 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
calculate_distance_matrix.py
import numpy as np
import nibabel as nib
from surfdist.load import load_freesurfer_label, get_freesurfer_label
from surfdist.utils import surf_keep_cortex, translate_src
import gdist
import os
base_dir = '/Applications/freesurfer/subjects/'
surf = nib.freesurfer.read_geometry(os.path.join(base_dir, 'fsaverage5/surf/lh.pial'))
cort = np.sort(nib.freesurfer.read_label(os.path.join(base_dir, 'fsaverage5/label/lh.cortex.label')))
cortex_vertices, cortex_triangles = surf_keep_cortex(surf, cort)
# remove exceptions from label list (some hacky things were done here. May not be necessary...):
exceptions = ['Unknown', 'Medial_wall']
label_list = get_freesurfer_label(os.path.join(base_dir, 'fsaverage5/label/lh.aparc.a2009s.annot'), verbose = False)
rs = np.where([str(a)[2:-1] not in exceptions for a in label_list])[0]
rois = [str(label_list[r])[2:-1] for r in rs]
# get indices of rois
rois_trans = []
for i in range(len(rois)):
source_nodes = load_freesurfer_label(os.path.join(base_dir, 'fsaverage5/label/lh.aparc.a2009s.annot'), rois[i])
trans_source_nodes = translate_src(source_nodes, cort)
rois_trans.append(trans_source_nodes)
# parallel jobs version:
from joblib import Parallel, delayed
from itertools import combinations
import multiprocessing
roi_combos = np.asarray(list(combinations(range(len(rois)),2)))
distance_matrix = np.zeros((len(rois),len(rois)),dtype=np.float)
def calc_roi_dist(rois_to_run, distance_matrix): # rn is roi_number
result = gdist.compute_gdist(cortex_vertices,cortex_triangles,
source_indices=rois_trans[rois_to_run[0]],
target_indices=rois_trans[rois_to_run[1]])
print(rois_to_run, np.min(result))
distance_matrix[rois_to_run[0]][rois_to_run[1]] = np.min(result)
distance_matrix[rois_to_run[1]][rois_to_run[0]] = distance_matrix[rois_to_run[0]][rois_to_run[1]]
return distance_matrix
Parallel(n_jobs=multiprocessing.cpu_count(), prefer="threads")(delayed(calc_roi_dist)(i,distance_matrix) for i in roi_combos)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment