03-ProbabilisticAssignment.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import siibra\n",
    "assert siibra.__version__ >= \"1.0a08\"\n",
    "from nilearn import plotting\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define two points in MNI space\n",
    "\n",
    "We start by specifying some points in the reference space. Here, we use the MNI space. We can find such points, amongst other possibilities, by clicking them in the [siibra-explorer](https://siibra-explorer.apps.hbp.eu/). `siibra` has specific data types for points and point sets, which are aware of the reference space used, and can be warped between reference spaces on demand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# These are copy-pasted from the interactive atlas viewer:\n",
    "points = siibra.PointSet(\n",
    "    (\n",
    "        [-25.65, -2.75, -33.75],\n",
    "        [-37.35, -81.05, -6.3],\n",
    "    ),\n",
    "    space='mni152', sigma_mm=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.view_markers(list(map(tuple,points)), ['red', 'cyan'], marker_size=10)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assign brain regions to the 3D points\n",
    "\n",
    "We assign the points to brain regions from the Julich-Brain cytoarchitectonic atlas, using a certain location tolerance of 5mm standard deviation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "julich_pmaps = siibra.get_map(\n",
    "    parcellation=siibra.parcellations.get('julich'),\n",
    "    space=siibra.spaces.get('mni152'),\n",
    "    maptype=siibra.MapType.STATISTICAL\n",
    ")\n",
    "assignments = julich_pmaps.assign(points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result of the assignment is a pandas dataframe, where each row describes a region assigned to one of the input structures (which in this case are the provided two points, 0 and 1). Each assignment has some scores describing the relation of the point, represented as a 3D Gaussian kernel with bandwidth corresponding to the point uncertainty, and the respectivty assigned statistical map from the parcellation:\n",
    "\n",
    "- \"Value\" refers to the average value of the statistcal map under the point kernel\n",
    "- \"Correlation\" is Pearson's correlation coefficient between the statistical map and the point kernel\n",
    "- \"IoU\" is the intersection over union between both\n",
    "- \"Contains\" is the intersection over area of the statistical map (indicating containedness of the point)\n",
    "- \"Contained\" is the intersection over area of the point (indicating containedness of the statistical map)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at two of the matched areas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assignments.sort_values(by='intersection over union', ascending=False, inplace=True)\n",
    "\n",
    "a1 = assignments[assignments['input structure']==0].iloc[0]\n",
    "a2 = assignments[assignments['input structure']==1].iloc[0]\n",
    "\n",
    "for a in [a1, a2]:\n",
    "    \n",
    "    # fetch the probability map corresponding to the assigned regions\n",
    "    statmap = julich_pmaps.fetch(region=a.region)\n",
    "\n",
    "    # plot the map, center at the point\n",
    "    view = plotting.plot_stat_map(\n",
    "        statmap, \n",
    "        title=f\"{a.region} ({a['correlation']:.2f})\", \n",
    "        cmap='viridis', \n",
    "        cut_coords=a.centroid)\n",
    "\n",
    "    # add the point\n",
    "    view.add_markers([a.centroid])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find connectivity in terms of streamlines from DTI\n",
    "\n",
    "We also investigate the connectivtiy of the first region as measured by in-vivo imaging. To do so, we select the first region in the atlas, and search for connectivity profiles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get structural connectivity matrices for  Julich Brain\n",
    "features = siibra.features.get(\n",
    "    julich_pmaps.parcellation,\n",
    "    siibra.features.connectivity.StreamlineCounts\n",
    ")\n",
    "\n",
    "# get one of the connectivity matrices from the first compound\n",
    "sc = features[0].get_element(\"000\")\n",
    "\n",
    "# get the connectivity profile with the strongest 10 connections for the first matched region\n",
    "profile = sc.get_profile(a1.region, max_rows=10)\n",
    "\n",
    "profile.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "c6134497cd5410b1b275ebc88f99d14855849379b696eb9e04cff5dd9aa5e77a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}