01-BasicConcepts.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Core concepts\n",
    "\n",
    "## Predefined atlases, parcellations, and reference spaces\n",
    "\n",
    "The `siibra.core` module provides classes for the core concepts of brain atlases. These include \n",
    "\n",
    " - the `Atlas` class as the basic entry point for working with a brain atlas \n",
    " - the `Parcellation` class, giving access to a particular brain segregation scheme\n",
    " - the `Space` class, givin access to a particular reference space\n",
    " \n",
    "The above classes are all derived from the basic `AtlasConcept` class, and represent semantic concepts. They are special in the sense that `siibra` automatically builds a registry with predefined objects for each of them as soon as you import the package. To configure the predefined objects, siibra will retrieve configuration details from EBRAINS and some other online resources when you import the package: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import siibra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To performing this bootstrap process, an internet connection is required to import siibra the first time. However, the configuration information is then stored in a local cache folder on your system, so after the first time you will be able to import your `siibra` package without being online."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accessing predefined atlases\n",
    "\n",
    "Siibra represents the core concepts via the registries  `siibra.atlases`, `siibra.parcellations`, and `siibra.spaces`. Elements in each registry can be accessed in different ways:\n",
    "\n",
    " - You can iterate over all objects\n",
    " - An integer index gives sequential access to individual elements\n",
    " - A string index will be matched against the name or key of objects. If it does not match exactly, an inexact string matching will be used to see if a unique entry can be found.\n",
    " - Object keys can be tab-completed as attributes of the registry\n",
    "\n",
    "Let's try this out for the `siibra.atlases` registry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Which atlases are provided? We can iterate over objects in the registry.\n",
    "for atlas in siibra.atlases:\n",
    "    print(repr(atlas),\"-\",atlas.name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Access the first element in the registry\n",
    "siibra.atlases[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Access elements by their name as attributes\n",
    "# with autocompletion by most Python interpreters\n",
    "siibra.atlases.MULTILEVEL_HUMAN_ATLAS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how both alternatives provide the same object. The easiest and recommended way to access items from a registry however, is to use keywords for accessing elements. `siibra` will try to figure out the matching item, or inform you if you need to be more precise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The easiest way: Using string matching of keywords\n",
    "siibra.atlases['human']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accessing parcellations and regions\n",
    "\n",
    "While the registry `siibra.parcellations` gives access to all available parcellations, the recommended way is to access parcellations via an atlas object. Each atlas object provices a filtered registry of its supported parcellations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "atlas = siibra.atlases['human']\n",
    "list(atlas.parcellations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can get access to any parcellation object using `atlas.get_parcellation`, or of course by directly using the registry:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jubrain = atlas.get_parcellation('julich') # will return the latest version per default\n",
    "\n",
    "# let's look at some metadata\n",
    "print(\"Name:    \",jubrain.name)\n",
    "print(\"Id:      \",jubrain.id)\n",
    "print(\"Modality:\",jubrain.modality)\n",
    "print()\n",
    "print(jubrain.description)\n",
    "print()\n",
    "for p in jubrain.publications:\n",
    "    print(p['citation'])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting parcellation is a semantic object. It does not itself represent a map, but rather the defintion of the parcellation. Thus it provides\n",
    "\n",
    "- the region hierarchy\n",
    "- functions to find and access regions\n",
    "- functions to access different forms of maps\n",
    "\n",
    "Let's search a region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# search regions known by the parcellation\n",
    "for r in jubrain.find_regions('amygdala'):\n",
    "    print(r.name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you see, areas often appear three times: Julich-Brain defines them separately for the left and right hemisphere, and additionally defines a common parent region. The parent object actually represents the corresponding subtree. We can more easily access concrete regions by using `decode_region` instead of `find_regions` - it aims to resolve a single element instead of returning all possible candidates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the whole amygdala subtree\n",
    "jubrain.decode_region('amygdala')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# only VTM, a small subtree\n",
    "jubrain.decode_region('VTM')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the exact left VTM - this is not a tree, it is a leaf of the region hierarchy\n",
    "jubrain.decode_region('vtm left')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To search for regions, we need not explicitly fetch the parcellation object. The atlas object provides a similar function, however, it will return matching regions from all parcellations it knows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# search all regions known by the atlas\n",
    "for r in atlas.find_regions('amygdala'):\n",
    "    print(f\"{r.name:30.30} {r.parcellation}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accessing maps\n",
    "\n",
    "Different from the above semantic objects - atlases, spaces, parcellations, regions - a parcellation map or region map is a spatial object. In some atlases, parcellations are defined in different spaces, so `siibra` expects you to specify one, or will use a default and inform you. Let's get the maximum probability map of the Julich-Brain in the MNI152 space to see how that works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpm = atlas.get_map(space=\"mni152\",parcellation=\"julich\")\n",
    "type(mpm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Labelled parcellation maps\n",
    "\n",
    "The returned map still does not contain image data. This is because `siibra` uses a lazy strategy for loading data. To access the actual image data, we call the `fetch()` method. We will use the excellent `nilearn` library for plotting the map. It plots in the MNI152 space by default, so as long as we use this space it works easily."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn import plotting\n",
    "plotting.plot_stat_map(mpm.fetch())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now this is only the left hemisphere! This is because Julich-Brain ships the left and right hemispheres in different volumes, so corresponding regions can use the same label index and still be distinguished. In fact, if a parcellation provides multiple maps, we can iterate them using `fetchall()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for img in mpm.fetchall():\n",
    "    plotting.plot_stat_map(img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probabilistic maps\n",
    "\n",
    "Julich-Brain, like some other parcellations, is a probabilistic map. The labelled volumes above are just a simplified representation, displaying for each voxel the btain region of highes probability. We can access the much richer information of the probability maps, which provide a continuous distribution in the reference space for each parcticular region.\n",
    "\n",
    "We received the labelled volumes above because `siibra` uses labelled volumes as the default map type. To retrieve probability maps, we explictiely use `siibra.maptype.CONTINUOUS`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmaps = atlas.get_map('mni152','julich',maptype='continuous')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we can iterate over all >300 pmaps using `fetchall()`. For simplicity, we will just a random index here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pmap = pmaps.fetch(mapindex=102)\n",
    "plotting.plot_stat_map(pmap)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linking region objects with map and label indices\n",
    "\n",
    "Now you migh wonder which region this refers to. Of course, the parcellation objects in `siibra` help you to safely translate indices into regions, and vice versa:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# which region corresponds to map index 10?\n",
    "pmaps.decode_label(mapindex=102)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# vice versa, what is the index of that region?\n",
    "pmaps.decode_region('hoc5 left')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, a region can be linked to a map index (the index of the image volume) and to a label index (the color in a labelled parcellation map, which does not apply to continuous maps). So, in order to find the label index of CM left in the maximum probability map, we would do exactly the same but get a different index:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpm.decode_region(\"hoc5 left\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can verify the index by thresholding the MPM with the label index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fetch the left hemishpere maximum probability map\n",
    "mpm_l = mpm.fetch()\n",
    "\n",
    "# build a mask by thresholding the label index reported for hoc5\n",
    "import numpy as np\n",
    "from nibabel import Nifti1Image\n",
    "A = mpm_l.get_fdata()\n",
    "A[A!=6] = 0\n",
    "mask = Nifti1Image(A,mpm_l.affine)\n",
    "plotting.plot_roi(mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, in order to get a regional mask, you need not take these steps. `siibra`'s region objects can do that right away."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = atlas.get_region('hoc5 left').build_mask(\"mni152\")\n",
    "plotting.plot_roi(mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}