Merge pull request #26 from int-brain-lab/streamlines

Streamlines

Author: oliche

examples/atlas_streamlines.ipynb

@@ -0,0 +1,166 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "5823062b-e25f-45aa-96d7-65c83845af61",
+   "metadata": {},
+   "source": [
+    "# Working with cortical streamlines\n",
+    "\n",
+    "In this example we walk through various methods that use the Allen streamlines. Streamlines are the paths that connect the surface of the isocortex to the lower white matter surface while following the curvature of these surfaces."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "46e7e0e7-4295-40a5-a5b3-1bc0fc84f07e",
+   "metadata": {},
+   "source": [
+    "## Depth lookup\n",
+    "We can use the streamlines to find the depth of coordinates in 3D space relative to the surface of the cortex. Note as the streamlines are only defined in the cortex so this lookup will only work for coordinates that lie within the Isocortex of the Allen atlas volume. In this example we obtain the depths in percentage, where 0 % indicates the surface of the isocortex and 100 % the white matter boundary."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0a3fddd1-9ddd-4dca-b256-3dd7d29bad67",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from iblatlas.streamlines.utils import xyz_to_depth\n",
+    "from iblatlas.atlas import AllenAtlas\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "ba = AllenAtlas()\n",
+    "# In this example we are going to get out the coordinates of the voxels that lie within the region MOp\n",
+    "region_id = ba.regions.acronym2index('MOp')[1][0][0]\n",
+    "# MOp is in the beryl mapping so we map our label volume to this mapping\n",
+    "ba.label = ba.regions.mappings['Beryl-lr'][ba.label]\n",
+    "# These are the indices in the volume, order is AP, ML, DV\n",
+    "ixyz = np.where(ba.label == region_id)\n",
+    "# Convert these to xyz coordinates relative to bregma in order ML, AP, DV\n",
+    "xyz = ba.bc.i2xyz(np.c_[ixyz[1], ixyz[0], ixyz[2]])\n",
+    "\n",
+    "depth_per = xyz_to_depth(xyz) * 100\n",
+    "\n",
+    "# We can assign a volume\n",
+    "depth_per_vol = np.full(ba.image.shape, np.nan)\n",
+    "depth_per_vol[ixyz[0], ixyz[1], ixyz[2]] = depth_per\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ba.plot_cslice(620/1e6, volume='image', ax=ax)\n",
+    "ba.plot_cslice(620/1e6, volume='volume', region_values=depth_per_vol, ax=ax, cmap='viridis', vmin=0, vmax=100)\n",
+    "ax.set_axis_off()\n",
+    "cbar = fig.colorbar(ax.get_images()[1], ax=ax)\n",
+    "cbar.set_label('Depth from surface of cortex (%)')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c1da61f0-7281-45d8-8bd7-33be0d373302",
+   "metadata": {},
+   "source": [
+    "We can also do this look up in um relative to the surface of the cortex."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "047509db-fb6b-46f3-8cb2-817b12223bed",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "depth_um = xyz_to_depth(xyz, per=False)\n",
+    "\n",
+    "# We can assign a volume\n",
+    "depth_um_vol = np.full(ba.image.shape, np.nan)\n",
+    "depth_um_vol[ixyz[0], ixyz[1], ixyz[2]] = depth_um\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ba.plot_cslice(620/1e6, volume='image', ax=ax)\n",
+    "ba.plot_cslice(620/1e6, volume='volume', region_values=depth_um_vol, ax=ax, cmap='tab20', vmin=0, vmax=2000)\n",
+    "ax.set_axis_off()\n",
+    "cbar = fig.colorbar(ax.get_images()[1], ax=ax)\n",
+    "cbar.set_label('Depth from surface of cortex (um)')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "bdbe5c0b-d7e3-441e-82ab-8bbcbd013bdd",
+   "metadata": {},
+   "source": [
+    "## Flatmap projection\n",
+    "The streamlines can also be used to project a volume or xyz coordinates onto a flattened view of the dorsal cortex. The volume or points can be aggregated along the streamlines in different ways, e.g max, mean, median.\n",
+    "\n",
+    "In this example we will project the Allen Atlas dwi volume onto the dorsal cortex projection. We will take the maximum value along the streamline for each voxel on the surface"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7b61405d-8589-4cb8-8afc-6f8a87dc4588",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from iblatlas.streamlines.utils import project_volume_onto_flatmap\n",
+    "\n",
+    "ba = AllenAtlas(25)\n",
+    "proj, fig, ax = project_volume_onto_flatmap(ba.image, res_um=25, aggr='max', plot=True, cmap='bone')\n",
+    "ax.set_axis_off()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "702b9abd-f97d-4a4e-967a-e9967e630293",
+   "metadata": {},
+   "source": [
+    "We can also do the same with xyz coordinates (we will use the same xyz coordinates for the region MOp and assign random values to each coordinate)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "19bf31f8-8f22-4f0e-91bb-7998f9379dc6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from iblatlas.streamlines.utils import project_points_onto_flatmap, get_mask\n",
+    "\n",
+    "# Get a mask to plot our results on, options are 'image', 'annotation' or 'boundary\n",
+    "mask = get_mask('boundary')\n",
+    "\n",
+    "values = np.random.rand(xyz.shape[0])\n",
+    "\n",
+    "proj = project_points_onto_flatmap(xyz, values, aggr='mean', plot=False)\n",
+    "# Mask values with 0 with nan for display purposes\n",
+    "proj[proj == 0] = np.nan\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.imshow(mask, cmap='Greys')\n",
+    "ax.imshow(proj)\n",
+    "ax.set_axis_off()"
+   ]
+  }
+ ],
+ "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.11.9"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

iblatlas/atlas.py

@@ -915,7 +915,8 @@ def _take_remap(vol, ind, axis, mapping):
                 index = bc.xyz2i(np.array([coordinate] * 3))[axis]
             return _take(region_values, index, axis=axis)
 
-    def compute_boundaries(self, values):
+    @staticmethod
+    def compute_boundaries(values):
         """
         Compute the boundaries between regions on slice
         :param values:
@@ -1811,54 +1812,3 @@ def get_bc(res_um=10):
     bc = BrainCoordinates(nxyz=nxyz, xyz0=-bc.i2xyz(iorigin), dxyz=dxyz)
 
     return bc
-
-
-def _download_depth_files(file_name):
-    """
-    These files have been generated using this script
-    https://github.com/int-brain-lab/ibldevtools/blob/master/Mayo/flatmaps/2025-03-20_depths_from_streamlines.py
-    :param file_name:
-    :return:
-    """
-    file_path = BrainAtlas._get_cache_dir().joinpath('depths', file_name)
-    if not file_path.exists():
-        file_path.parent.mkdir(exist_ok=True, parents=True)
-        aws.s3_download_file(f'atlas/depths/{file_path.name}', file_path)
-
-    return np.load(file_path)
-
-
-def xyz_to_depth(xyz, res_um=25):
-    """
-    For a given xyz coordinates find the depth from the surface of the cortex. Note the lookup will only
-    work for xyz cooordinates that are in the Isocortex of the Allen volume. If coordinates outside of this
-    region are given then the depth is returned as nan.
-
-    Parameters
-    ----------
-    xyz : numpy.array
-        An (n, 3) array of Cartesian coordinates. The order is ML, AP, DV and coordinates should be given in meters
-        relative to bregma.
-
-    res_um : float or int
-        The resolution of the brain atlas to do the depth lookup
-
-    Returns
-    -------
-        numpy.array
-        The depths from the surface of the cortex for each cartesian coordinate. If the coordinate does not lie within
-        the Isocortex, depth value returned is nan
-    """
-
-    ind_flat = _download_depth_files(f'depths_ind_{res_um}.npy')
-    depths = _download_depth_files(f'depths_um_{res_um}.npy')
-    bc = get_bc(res_um=res_um)
-
-    ixyz = bc.xyz2i(xyz, mode='clip')
-    iravel = np.ravel_multi_index((ixyz[:, 1], ixyz[:, 0], ixyz[:, 2]), (bc.ny, bc.nx, bc.nz))
-    a, b = ismember(iravel, ind_flat)
-
-    lookup_depths = np.full(iravel.shape, np.nan, dtype=np.float32)
-    lookup_depths[a] = depths[b]
-
-    return lookup_depths

iblatlas/streamlines/__init__.py

@@ -0,0 +1,54 @@
+"""A package for working with Allen streamlines data.
+
+This package provides a way to use the Allen streamline data to
+1. Do a lookup of xyz position to depth from the surface of the cortex
+2. Project a volume or xyz positions along the streamlines to display them on a flattened dorsal cortex flatmap
+
+Streamlines are the paths that connect the surface of the isocortex to the white matter surface while following the
+curvature of these surfaces.
+
+For the xyz to depth position lookup the following files have been generated with 10, 25 and 50 um resolution
+
+     depths/
+     ├── depths_ind_{res}.npy
+     ├── depths_per_{res}.npy
+     ├── depths_um_{res}.npy
+
+- The depths_ind_{res}.npy contains the flattened index for all voxels contained within the isocortex of the Allen
+volume at the specified resolution. The ordering of the flattened index is AP, DV, ML
+- The depths_per_{res}.npy contains the depth in percentage from the surface of the cortex of each voxel in the isocortex.
+The percentage values were taken from the laplacian provided by the Allen institute. Voxels at the surface have a percentage
+ depth of 0 while those at the white matter boundary have a percentage depth of 1.
+- The depths_um_{res}.npy contains the depth in um from the surface of the cortex of each voxel in the isocortex. The depth
+values have been computed using the streamlines provided by the Allen institute. For each streamline the length of
+the streamline was computed and voxels assigned a depth based on the steps along the streamline. When multiple streamlines
+pass through a given voxel the median depth has been taken.
+
+The 25 and 50 um resolution datasets have been computed by subsampling the 10 um datasets.
+
+
+To enable the projection onto the dorsal flatmap the following files have been generated with 10, 25 and 50 um resolution
+
+     depths/
+     ├── paths_{res}.pqt
+     ├── dorsal_flatmap.npy
+     ├── dorsal_annotation.npy
+     ├── dorsal_image.npy
+
+
+- The paths_{res}.pqt files contain information about the voxels that correspond to each streamline. They have been
+derived from the dorsal_flatmap_paths_10.h5 paths provided by the Allen institute. We store the values as a pqt file and remove
+values that correspond to no streamline (paths=0) and duplicate voxels per streamline. As we have removed duplicate
+voxels per streamline this means that when taking the average or std along the streamline the behavior will differ
+from the Allen code.
+- The dorsal_flatmap.npy is a 2D numpy array with the dorsal cortex flatmap. The value in each voxel corresponds the
+streamline index that corresponds to that voxel.
+- The dorsal_annotation is 2D numpy array with the dorsal cortex flatmap. The value in each voxel corresponds to the
+region id on the surface of the cortex that corresponds to that voxel
+- The dorsal_image is a 2D numpy array with the dorsal cortex flatmap. The value in each voxel corresponds to the Allen
+image volume projected onto the flatmap display
+
+[1] Wang, Quanxin et al., “The Allen Mouse Brain Common Coordinate Framework: A 3D Reference Atlas”
+ Cell, Volume 181, Issue 4, 936 - 953.e20 May. 2020, https://doi.org/10.1016/j.cell.2020.04.007
+
+"""