06-CorticalPatchSampling.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0a7d492",
   "metadata": {},
   "outputs": [],
   "source": [
    "import siibra\n",
    "assert siibra.__version__ >= \"1.0a08\"\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib notebook\n",
    "import numpy as np\n",
    "from nilearn import plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02c460e6",
   "metadata": {},
   "source": [
    "## Select a probability map of a cortical brain region from the Julich-Brain atlas\n",
    "\n",
    "We start this tutorial by selecting a cortical brain region from the Julich-Brain cytoarchitectonic atlas. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01b3935d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we will search for the region by keywords\n",
    "regionspec = '4p'\n",
    "hemisphere = 'right'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7ade335",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the parcellation, its map in MNI space\n",
    "julichbrain = siibra.parcellations.get('julich 3')\n",
    "region = julichbrain.get_region(f'{regionspec} {hemisphere}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4befd21",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get access to the corresponding probabilistc maps in MNI space, \n",
    "# and fetch the one for the selected region.\n",
    "julich_pmaps = julichbrain.get_map(space='mni152', maptype='statistical')\n",
    "pmap = julich_pmaps.fetch(region=region)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1c17f3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The result is a NIfTI image object.\n",
    "# We can plot it with common neuroscience libraries; \n",
    "# here we juse nilearn's wonderful 'plotting' module.\n",
    "plotting.plot_glass_brain(pmap, cmap='viridis', title=region.name)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7d47d20",
   "metadata": {},
   "source": [
    "## Find a 1-micron section intersecting the brain region\n",
    "\n",
    "We can use region maps to query spatial features.\n",
    "Here we search for 1 micron sections of the BigBrain which intersect with the region map.\n",
    "Note that siibra automatically resovles the mismatch between the template spaces:\n",
    "It will warp the bounding boxes of BigBrain sections to MNI space for the query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8de57a5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# find 1 micron section in bigbrain space which overlaps with the region\n",
    "pmap_volume = siibra.volumes.from_nifti(pmap, space='mni152', name=region.name)\n",
    "sections = siibra.features.get(pmap_volume, siibra.features.cellular.CellbodyStainedSection)\n",
    "section = sections[len(sections)//2]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "178c2334",
   "metadata": {},
   "source": [
    "## Intersect surface meshes with the image plane\n",
    "\n",
    "We are interested in cortical structure. Here it is very hepful to project the cortical layer surfaces to the image section. This is a 3D to 2D projection, leading to sets of (possibly split) contours of the layer surfaces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6f11475",
   "metadata": {},
   "outputs": [],
   "source": [
    "# access cortical layer maps in BigBrain space\n",
    "layers = siibra.parcellations.get('cortical layers')\n",
    "layermap = layers.get_map(space='bigbrain')\n",
    "\n",
    "# Derive a 3D plan from the 1 micron section.\n",
    "# Note that the 1 micron section is a 3D image volume,\n",
    "# however its y-dimension is flat. \n",
    "# Siibra will create a plane from the smallest image dimension.\n",
    "# The plane is defined in the physical brain reference space.\n",
    "imgplane = siibra.experimental.Plane3D.from_image(section)\n",
    "print(f\"Extracted a plane in {imgplane.space} space, with normal {imgplane.normal}.\")\n",
    "\n",
    "# Intersect all meshes in the parcellation map with this plane\n",
    "# to obtain layer contours. \n",
    "# The contours are represented as ordered sets of points.\n",
    "layer_contours = {\n",
    "    l: imgplane.intersect_mesh(layermap.fetch(region=l, format='mesh'))\n",
    "    for l in layermap.regions\n",
    "}\n",
    "\n",
    "for layername, contours in layer_contours.items():\n",
    "    print(\n",
    "        f\" - {layername+':':30} {len(contours)} intersected contours.\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2346bf2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the layer contours with a lower resolution version of the image section.\n",
    "# We do this with matplotlib, using the actual pixel space of the image,\n",
    "# so the image is not properly oriented here.\n",
    "\n",
    "# fetch a 0.8mm version\n",
    "im_low_res = section.fetch(resolution_mm=0.8)\n",
    "\n",
    "plt.figure()\n",
    "plt.imshow(im_low_res.get_fdata().squeeze(), cmap='gray', vmin=0, vmax=2**16)\n",
    "plt.axis('off')\n",
    "\n",
    "# siibra's map objects typically provide colormaps, let's use them!\n",
    "layercolors = {\n",
    "    l:layermap.get_colormap().colors[layermap.get_index(l).label]\n",
    "    for l in layermap.regions\n",
    "}\n",
    "\n",
    "# plot the contours\n",
    "for layername, contour_list in layer_contours.items():\n",
    "    for contours in contour_list:\n",
    "        phys2pix = np.linalg.inv(im_low_res.affine)\n",
    "        pixels = contours.transform(phys2pix, space=None).homogeneous\n",
    "        plt.plot(pixels[:, 2], pixels[:, 0], '-', color=layercolors[layername]) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4a19325",
   "metadata": {},
   "source": [
    "## Find a closeby cortical profile in the brain region\n",
    "\n",
    "We want to find a cortical profile which is likely in the selected brain region, and also as close as possible to the image section.\n",
    "\n",
    "We do so by looking up the values of layer IV contour points in the region map, and selecting the point with highest probability.\n",
    "\n",
    "Then we find the \"column\" of corresponding surface vertices in the cortical layer maps closest to the contour point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f918ed4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# find the most probable contour point in the layer 4 surface\n",
    "l4layer = layers.get_region(f\"4 {hemisphere}\")\n",
    "\n",
    "# collect all layer IV contour points into one set\n",
    "l4points = siibra.PointSet(\n",
    "    [p for contour in layer_contours[l4layer.name] for p in contour], \n",
    "    space='bigbrain'\n",
    ")\n",
    "\n",
    "# use the probabilistic map volume to lookup their values \n",
    "probs = pmap_volume.evaluate_points(l4points)\n",
    "\n",
    "# get the highest ranked point\n",
    "l4point = l4points[probs.argmax()]\n",
    "\n",
    "# siibra has an experimental class for querying cortical profiles\n",
    "profile = siibra.experimental.CorticalProfileSampler().query(l4point)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6c97fa6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the image crop around the brain region with layers and the profile\n",
    "crop_voi = section.get_boundingbox().intersection(pmap_volume)\n",
    "crop = section.fetch(voi=crop_voi, resolution_mm=0.2)\n",
    "phys2pix = np.linalg.inv(crop.affine)\n",
    "\n",
    "# plot the image\n",
    "plt.figure()\n",
    "plt.imshow(crop.get_fdata().squeeze(), cmap='gray', vmin=0, vmax=2**16)\n",
    "\n",
    "# plot the contour segments\n",
    "for layername, contours in layer_contours.items():\n",
    "    for contour in contours:\n",
    "        for segment in contour.crop(crop_voi):\n",
    "            pixels = segment.transform(phys2pix, space=None).homogeneous\n",
    "            plt.plot(pixels[:, 2], pixels[:, 0], '-', ms=4, color=layercolors[layername])\n",
    "\n",
    "# plot the profile points\n",
    "for p in profile:\n",
    "    x, y, z = p.transform(phys2pix, space=None)\n",
    "    plt.plot(z, x, 'r.', ms=10)\n",
    "    \n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10252153",
   "metadata": {},
   "source": [
    "## Define the cortical image patch including the profile\n",
    "\n",
    "Now extract a cortical patch from the section which is centered around the layer 4 point and approximately rotated orthogonal to the brain surface.\n",
    "\n",
    "siibra's Plane3D class can extract a rotated patch definition from a set of points, which are projected onto the plane.\n",
    "The rotation is derived from the principal axes of the projected points.\n",
    "The patch can then extract image data in the new orientation from the original image volume.\n",
    "\n",
    "Since this function is not aware of the special ordering of the profile points however, the patch might be extracted upside down. \n",
    "We can detect this using a bit of a hack, by checking to which of the corner points the inner surface point is closest, and flip if required."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eaba9b97",
   "metadata": {},
   "outputs": [],
   "source": [
    "patch_canvas = imgplane.get_enclosing_patch(profile)\n",
    "i = np.argmin(\n",
    "    np.linalg.norm((patch_canvas.corners.coordinates - profile[0].coordinate), axis=1)\n",
    ")\n",
    "if i in [0, 3]:\n",
    "    patch_canvas.flip()\n",
    "    \n",
    "# now fetch the image data\n",
    "patch = patch_canvas.extract_volume(section, resolution_mm=0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de4d23b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the patch individually, and in 3D context.\n",
    "# Here we first use nilearn plotting functions with the resulting NIfTI files,\n",
    "# so the extracted patch is still shown in its anatomical orientation.\n",
    "view = plotting.plot_img(patch.fetch(), cmap='gray', display_mode='y', annotate=False, title=f\"Rotated patch\")\n",
    "template = siibra.get_template(\"bigbrain\")\n",
    "bigbrain_lowres = template.fetch(resolution_mm=0.8)\n",
    "view = plotting.plot_img(patch.fetch(), bg_img=bigbrain_lowres, cmap='gray', title=\"Rotated patch in reference space\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4af5faee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# to verify that the underlying pixels have been resampled\n",
    "# into local cortical orientation, we display the underlying image\n",
    "# using matplotlib.\n",
    "\n",
    "plt.figure()\n",
    "\n",
    "# access the underlying image data\n",
    "patchimg = patch.fetch()\n",
    "phys2vox = np.linalg.inv(patchimg.affine)\n",
    "patchdata = patchimg.get_fdata().squeeze()\n",
    "\n",
    "# plot the pure image array\n",
    "plt.imshow(patchdata, cmap='gray', vmin=0, vmax=2**16)\n",
    "plt.axis('off')\n",
    "\n",
    "\n",
    "# overlay fragments of intersected layer surface contours.\n",
    "for layername, contours in layer_contours.items():\n",
    "    for contour in contours:\n",
    "        for s in contour.crop(patch.get_boundingbox()):\n",
    "            pixels = s.transform(phys2vox, space=None).homogeneous\n",
    "            plt.plot(pixels[:, 2], pixels[:, 0], '-', ms=4, color=layercolors[layername])\n",
    "\n",
    "# let's show the title of the extracted patch for fun...\n",
    "import textwrap\n",
    "plt.title(\"\\n\".join(textwrap.wrap(patch.name, 50)))\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53afcc33",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}