Skip to content

Commit

Permalink
added a pythonic implementation of tck2connectom with radial search (…
Browse files Browse the repository at this point in the history 
…also see #30)
  • Loading branch information
@sina-mansour
sina-mansour committed Mar 28, 2022
1 parent ef0b9a4 commit 88a4123
Show file tree
Hide file tree
Showing 2 changed files with 428 additions and 0 deletions.
356 changes: 356 additions & 0 deletions scripts/python/notebooks/UKbiobank_python_tck2connectome.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,356 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "british-galaxy",
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"import gc\n",
"import sys\n",
"import random\n",
"import datetime\n",
"import importlib\n",
"import itertools\n",
"import numpy as np\n",
"from scipy import spatial\n",
"import scipy.sparse as sparse\n",
"import scipy.stats as stats\n",
"import pandas as pd\n",
"import nibabel as nib\n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sns\n",
"\n",
"from mayavi import mlab\n",
"\n",
"# os.chdir('/media/sina/Windows1/Users/smansourlako/Documents/Reserach/Codes/fMRI')\n",
"os.chdir('/home/sina/Documents/Research/Codes/fMRI')\n",
"\n",
"# import constants as cs\n",
"import myconstants as cs\n",
"# import importlib.util\n",
"# spec = importlib.util.spec_from_file_location('cs', '/home/sina/Documents/Research/Codes/fMRI/constants.py')\n",
"# cs = importlib.util.module_from_spec(spec)\n",
"# sys.modules['cs'] = cs\n",
"# spec.loader.exec_module(cs)\n",
"import utils\n",
"import niutils\n",
"import hcp\n"
]
},
{
"cell_type": "markdown",
"id": "sophisticated-boards",
"metadata": {},
"source": [
"---\n",
"\n",
"Let's first load all the files needed:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "desirable-ordering",
"metadata": {},
"outputs": [],
"source": [
"# load the atlas file\n",
"\n",
"atlas_file = \"/home/sina/Documents/Research/Datasets/UK_biobank/sample/1000243_2_0/dMRI/dMRI/atlases/combinations/native.dMRI_space.aparc.a2009s+Tian_Subcortex_S1_3T.nii.gz\"\n",
"# atlas_file = \"/home/sina/Documents/Research/Codes/UKB-connectomics/data/temporary/subjects/1000243_2_0/tractography/atlases/native.dMRI_space.aparc.a2009s.nii.gz\"\n",
"\n",
"atlas=nib.load(atlas_file)\n",
"\n",
"# load the tractography endpoint information\n",
"\n",
"endpoint_file = \"/home/sina/Documents/Research/Codes/UKB-connectomics/data/temporary/subjects/1000243_2_0/tractography/endpoints/tracks_10M_endpoints.npy\"\n",
"\n",
"endpoints = np.load(endpoint_file)\n"
]
},
{
"cell_type": "markdown",
"id": "bottom-modem",
"metadata": {},
"source": [
"---\n",
"\n",
"Now we need to extract appropriate voxel coordinates from the atlas:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "283ffcf6-170c-4976-8359-73c129ad857a",
"metadata": {},
"outputs": [],
"source": [
"node_indices = np.arange(np.prod(atlas.shape)).reshape(atlas.shape)\n",
"\n",
"ind_i, ind_j, ind_k = np.meshgrid(\n",
" np.arange(atlas.shape[0]),\n",
" np.arange(atlas.shape[1]),\n",
" np.arange(atlas.shape[2]), indexing='ij',\n",
")\n",
"\n",
"node_ijk = np.array([ind_i.reshape(-1), ind_j.reshape(-1), ind_k.reshape(-1),]).T\n",
"\n",
"node_xyz = nib.affines.apply_affine(atlas.affine, node_ijk)\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "778784f8-cec9-418a-a4dd-bb6894c1fc4f",
"metadata": {},
"outputs": [],
"source": [
"# only select voxels with a label greater than zero\n",
"selection_mask = (atlas.get_fdata() > 0)\n",
"\n",
"selection_indices = node_ijk[selection_mask.reshape(-1), :]\n",
"\n",
"selection_xyz = node_xyz[selection_mask.reshape(-1), :]\n",
"\n",
"selection_labels = atlas.get_fdata()[selection_mask].astype(int)\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "ed2d85c4-10a3-45ee-bfc3-f6bea14a909f",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 6.75 s, sys: 281 ms, total: 7.03 s\n",
"Wall time: 7.03 s\n"
]
}
],
"source": [
"%%time\n",
"# build a kdtree for spatial queries\n",
"kdtree = spatial.cKDTree(selection_xyz)\n",
"\n",
"# get the list of endpoints\n",
"starts = endpoints[:,0,:]\n",
"ends = endpoints[:,-1,:]\n",
"\n",
"# query for closest coordinate from selection\n",
"start_dists, start_indices = kdtree.query(starts)\n",
"end_dists, end_indices = kdtree.query(ends)\n",
"\n",
"# mask points that are further than the search radius from all selection coordinates\n",
"search_radius = 4\n",
"distance_mask = (start_dists < search_radius) & (end_dists < search_radius)\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "3bcdc681-b951-431e-a2a7-56b3d7f9826d",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 194 ms, sys: 12.2 ms, total: 206 ms\n",
"Wall time: 205 ms\n"
]
}
],
"source": [
"%%time\n",
"# now generate a connectivity matrix\n",
"\n",
"# only keep valid endpoints according to the search radius\n",
"valid_start_indices = start_indices[distance_mask]\n",
"valid_end_indices = end_indices[distance_mask]\n",
"\n",
"# number of regions/nodes\n",
"node_count = selection_labels.max()\n",
"\n",
"# generate connectivity matrix\n",
"adj = np.zeros((node_count, node_count), dtype=np.float32)\n",
"np.add.at(adj, (selection_labels[valid_start_indices] - 1, selection_labels[valid_end_indices] - 1), 1)\n",
"adj = adj + adj.T\n",
"adj[np.diag_indices_from(adj)] /= 2\n"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "4ec83bbd-8a5c-49e4-a0c3-fea3a56d2f3d",
"metadata": {},
"outputs": [],
"source": [
"adj[np.diag_indices_from(adj)] /= 2\n"
]
},
{
"cell_type": "markdown",
"id": "exact-valuable",
"metadata": {},
"source": [
"---\n",
"\n",
"Let's also compare with mrtix generate adjacency:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "coral-defensive",
"metadata": {},
"outputs": [],
"source": [
"mradj = np.loadtxt(\"/home/sina/Documents/Research/Codes/UKB-connectomics/data/temporary/subjects/1000243_2_0/tractography/connectomes/aparc.a2009s+Tian_Subcortex_S1_3T/connectome_streamline_count_10M.csv\", delimiter=',')"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "76fd8c7e-5f59-4c5e-b445-ffd2e6404255",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[7.370e+02, 1.100e+01, 0.000e+00, ..., 4.400e+01, 7.600e+02,\n",
" 4.280e+02],\n",
" [1.100e+01, 3.048e+03, 0.000e+00, ..., 1.000e+00, 1.060e+02,\n",
" 2.000e+00],\n",
" [0.000e+00, 0.000e+00, 5.647e+03, ..., 2.990e+02, 4.890e+02,\n",
" 1.170e+02],\n",
" ...,\n",
" [4.400e+01, 1.000e+00, 2.990e+02, ..., 6.200e+01, 1.404e+03,\n",
" 6.500e+01],\n",
" [7.600e+02, 1.060e+02, 4.890e+02, ..., 1.404e+03, 5.933e+03,\n",
" 2.620e+02],\n",
" [4.280e+02, 2.000e+00, 1.170e+02, ..., 6.500e+01, 2.620e+02,\n",
" 2.864e+03]])"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mradj"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "ba8032ff-f0e1-4f1e-a2e5-ecbcb988dada",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2407"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.sum(adj!=mradj)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "7881f788-aa73-4484-bd55-e5ff9644708a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"24489"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.sum(adj==mradj)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "8363cf7a-f90b-4b15-b463-9184d3396b2c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[1. , 0.99999966],\n",
" [0.99999966, 1. ]])"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.corrcoef(adj.reshape(-1), mradj.reshape(-1))"
]
},
{
"cell_type": "markdown",
"id": "altered-acceptance",
"metadata": {},
"source": [
"---\n",
"\n",
"Store output as csv:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "included-family",
"metadata": {},
"outputs": [],
"source": [
"np.savetxt('tmp.csv', adj, delimiter=',')"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.10"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Loading

0 comments on commit 88a4123

Please sign in to comment.