diff --git a/docs/examples/Makefile b/docs/examples/Makefile index 1b2c97bc..aa70fb11 100644 --- a/docs/examples/Makefile +++ b/docs/examples/Makefile @@ -29,6 +29,7 @@ EXAMPLE_NOTEBOOKS = getting_started_io/00_test_installation.ipynb \ head_models/46_precompute_fluence.ipynb \ head_models/48_headmodel_landmarks_verification.ipynb \ head_models/51_spring_relaxation_registration.ipynb \ + head_models/52_mni_atlas_labels_aal3_brodmann.ipynb \ augmentation/61_synthetic_artifacts_example.ipynb \ augmentation/62_synthetic_hrfs_example.ipynb \ physio/71_ampd_heartbeat.ipynb \ diff --git a/examples/head_models/52_mni_atlas_labels_aal3_brodmann.ipynb b/examples/head_models/52_mni_atlas_labels_aal3_brodmann.ipynb new file mode 100644 index 00000000..dd2b1a49 --- /dev/null +++ b/examples/head_models/52_mni_atlas_labels_aal3_brodmann.ipynb @@ -0,0 +1,858 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Assigning AAL3 and Brodmann Atlas Labels to Brain Surfaces\n", + "\n", + "This notebook assigns AAL3 and Brodmann atlas labels to the Colin27 and ICBM152 brain surfaces, visualizes the labels, and summarizes their overlap with the Schaefer2018 parcels.\n", + "\n", + "## Overview\n", + "\n", + "Neuroimaging atlases — parcellations that divide the cerebral cortex into anatomically\n", + "or functionally defined regions — are a cornerstone of group-level fNIRS and DOT\n", + "analysis. Atlases allow researchers to:\n", + "\n", + "* Report results in a standardised anatomical vocabulary that is comparable across\n", + " studies and laboratories.\n", + "* Aggregate channel or vertex signals into region-of-interest (ROI) time series,\n", + " reducing the multiple-comparisons burden.\n", + "* Cross-reference DOT activations with the large fMRI literature, where the same\n", + " atlas labels are widely used.\n", + "\n", + "### MNI space as a shared coordinate system\n", + "\n", + "Most brain atlases are defined and distributed in **MNI (Montreal Neurological\n", + "Institute) space**: a standardised brain coordinate frame derived from averaging many\n", + "individual MRI scans into a common template.\n", + "Every standard head model in Cedalion — Colin27\n", + "(Holmes et al., 1998)\n", + "and ICBM152\n", + "(Fonov et al., 2011) — carries pre-computed MNI152\n", + "coordinates for every brain surface vertex. This means that **any** volumetric atlas\n", + "provided as a NIfTI file in MNI space can be transferred to the brain surface without\n", + "any additional registration step: the shared coordinate system acts as the bridge\n", + "between the volumetric atlas and the surface model.\n", + "\n", + "### The two MNI spaces: MNI152 and MNI305\n", + "\n", + "Two MNI variants are in common use and it is important to keep them distinct:\n", + "\n", + "| Space | Template | Typical users |\n", + "|---|---|---|\n", + "| **MNI152** | Average of 152 adult brains (ICBM152 template) | FSL, recent SPM, most modern atlases |\n", + "| **MNI305** | Average of 305 brains; the original MNI reference | FreeSurfer, older SPM, some legacy atlases |\n", + "\n", + "MNI152 and MNI305 differ by a small affine shift (a few millimetres in some regions)\n", + "and are **not** interchangeable. Cedalion stores a pre-computed affine transform\n", + "(`cedalion.dot.utils.mni305_to_mni152`) and applies it automatically when a NIfTI is\n", + "declared as `voxel_label_crs='mni305'`, so you do not need to handle the conversion\n", + "manually. Both spaces are fully supported; you simply specify which one your NIfTI\n", + "file uses.\n", + "\n", + "### What this notebook demonstrates\n", + "\n", + "This notebook shows the general workflow for applying **any** user-supplied parcellation\n", + "scheme — provided as a NIfTI label volume in either MNI152 or MNI305 space — to the\n", + "Cedalion head models. The workflow is illustrated with two widely used atlases:\n", + "\n", + "* **AAL3** (Automated Anatomical Labelling atlas, version 3)\n", + " (Rolls et al., 2020) — 170 macro-anatomical regions\n", + " covering the entire cerebral cortex and subcortex.\n", + "* **Brodmann areas** (Mai–Majtanik atlas)\n", + " (Mai & Majtanik, 2017) — the classical cytoarchitectonic\n", + " map, available as a volumetric MNI label file.\n", + "\n", + "Both atlases are bundled with Cedalion and loaded via `cedalion.data.get_atlas_files()`.\n", + "The same code applies unchanged to any other NIfTI atlas you supply.\n", + "\n", + "**For more detail, see also:**\n", + "- [43a_head_models_overview.ipynb](43a_head_models_overview.ipynb) — introduction to\n", + " `TwoSurfaceHeadModel` and the Schaefer2018 parcellation bundled with the standard models\n", + "- [43b_individualized_head_models.ipynb](43b_individualized_head_models.ipynb) — building an\n", + " individualized head model from a subject's own MRI using FreeSurfer, which yields\n", + " surface-native parcel labels with sharper boundaries than the volumetric approach shown here" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "import json\n", + "\n", + "import numpy as np\n", + "import pyvista as pv\n", + "from scipy.spatial import KDTree\n", + "\n", + "import cedalion\n", + "import cedalion.dot\n", + "from cedalion.vis.anatomy import get_vertex_colors_from_coord, plot_brain_views_grid\n", + "\n", + "pv.set_jupyter_backend(\"static\")" + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "## Loading Colin27 and ICBM152 Head Models\n", + "\n", + "We load both standard atlas head models. Each is stored in voxel space (`crs='ijk'`);\n", + "`assign_parcels_via_mni_coords` works in MNI152 space internally so the coordinate\n", + "system of the loaded head model does not need to match that of the atlas.\n", + "\n", + "Both models carry pre-computed `mni152_r/a/s` vertex coordinates on the brain surface,\n", + "which is what makes the atlas transfer possible. The inflated cortex surfaces are loaded\n", + "for visualization: inflation removes sulci and gyri so that buried cortex becomes\n", + "visible, making parcel boundaries much easier to inspect." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "colin_ijk = cedalion.dot.get_standard_headmodel(\"colin27\")\n", + "colin_inflated = cedalion.dot.get_inflated_cortex_surface(\"colin27\")\n", + "\n", + "icbm_ijk = cedalion.dot.get_standard_headmodel(\"icbm152\")\n", + "icbm_inflated = cedalion.dot.get_inflated_cortex_surface(\"icbm152\")" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "## Loading Atlases\n", + "\n", + "Atlases distributed for use with standard MNI templates typically come in two files:\n", + "\n", + "1. A **NIfTI volume** (`.nii` or `.nii.gz`) in which every voxel inside a labelled\n", + " region carries a positive integer, and background voxels carry 0 (or another\n", + " reserved value).\n", + "2. A **label mapping** (`.json`, `.csv`, or `.txt`) that translates those integers to\n", + " human-readable region names.\n", + "\n", + "`cedalion.data.get_atlas_files(name)` returns both paths for the bundled atlases.\n", + "For a custom atlas you can pass any NIfTI path you have downloaded or created.\n", + "\n", + "### AAL3 — Automated Anatomical Labelling atlas (version 3)\n", + "\n", + "AAL3 (Rolls et al., 2020) divides the cerebral cortex\n", + "and subcortical structures into 170 macro-anatomical regions based on sulcal and gyral\n", + "landmarks visible in the MNI152 template. It is one of the most widely used atlases in\n", + "the fMRI and fNIRS literature and provides a common vocabulary for reporting results.\n", + "The NIfTI file bundled with Cedalion is defined in **MNI152** space." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "aal3_voxel_label_niftii, aal3_labels_json = cedalion.data.get_atlas_files(\"aal3\")\n", + "\n", + "# dictionary to map numeric voxel labels in nifti to string labels\n", + "with aal3_labels_json.open(\"r\") as fin:\n", + " aal3_num2label = json.load(fin)\n", + " aal3_num2label = {i[\"index\"] : i[\"name\"] for i in aal3_num2label[\"labels\"]}\n", + "\n", + "aal3_num2label" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "### Brodmann Areas — Mai–Majtanik Atlas\n", + "\n", + "Brodmann areas are the classical cytoarchitectonic map of the human cortex, originally\n", + "defined by Korbinian Brodmann in 1909 from histological sections and still widely used\n", + "to communicate the location of activations in the neuroscience literature. The\n", + "Mai–Majtanik atlas (Mai & Majtanik, 2017) provides\n", + "these areas as a digitised volumetric label map registered to **MNI152** space, making\n", + "it directly compatible with the Cedalion head models. Areas such as BA44/45\n", + "(Broca's area) or BA17 (primary visual cortex) are convenient reference landmarks when\n", + "linking fNIRS results to the broader neuroimaging literature." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "brodmann_voxel_label_niftii, brodmann_labels_json = cedalion.data.get_atlas_files(\"brodmann\")\n", + "\n", + "# dictionary to map numeric voxel labels in nifti to string labels\n", + "with brodmann_labels_json.open(\"r\") as fin:\n", + " brodmann_num2label = json.load(fin)\n", + " brodmann_num2label = {i[\"index\"] : i[\"name\"] for i in brodmann_num2label[\"labels\"]}\n", + "\n", + "brodmann_num2label" + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "## Brain Vertex Coordinates Before Atlas Assignment\n", + "\n", + "The `brain.vertices` xarray `DataArray` already carries several vertex coordinates set\n", + "when the head model was built from FreeSurfer outputs:\n", + "\n", + "| Coordinate | Meaning |\n", + "|---|---|\n", + "| `parcel` | Schaefer2018 parcel label (Schaefer et al., 2018) |\n", + "| `fsaverage_vertex` | Corresponding vertex index in the FreeSurfer `fsaverage` template |\n", + "| `mni152_r/a/s` | MNI152 RAS coordinates of the vertex in mm |\n", + "\n", + "Calling `assign_parcels_via_mni_coords` will **add** a new named coordinate to this\n", + "`DataArray` without removing or modifying any existing one. Multiple atlas labels can\n", + "therefore coexist on the same vertex, which is what we exploit below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "colin_ijk.brain.vertices" + ] + }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "## Assigning Atlas Labels to Brain Vertices\n", + "\n", + "`TwoSurfaceHeadModel.assign_parcels_via_mni_coords` transfers labels from a volumetric\n", + "NIfTI atlas to the brain surface by exploiting the MNI152 coordinates shared by both.\n", + "The algorithm works as follows:\n", + "\n", + "1. **Coordinate alignment**: The NIfTI affine transform is used to compute the MNI\n", + " coordinate of every voxel centre. If the atlas is in MNI305 space\n", + " (`voxel_label_crs='mni305'`), the pre-computed affine\n", + " `cedalion.dot.utils.mni305_to_mni152` is applied first to bring voxel coordinates\n", + " into the same MNI152 frame as the head model vertices.\n", + "\n", + "2. **Nearest-labelled-voxel search**: A KD-tree is built from the voxel centres. For\n", + " each brain surface vertex, all voxels within a ball of radius `mni_eps` mm in MNI152\n", + " space are retrieved. Background voxels (label 0 or whichever integer maps to\n", + " `background_label`) are excluded, and the closest *labelled* voxel is selected.\n", + "\n", + "3. **Label assignment**: The integer label of the winning voxel is translated to a\n", + " string via `label_mapping` and stored as a new vertex coordinate. If no labelled\n", + " voxel falls within `mni_eps` mm, the vertex receives `background_label`.\n", + "\n", + "**Key parameters:**\n", + "\n", + "| Parameter | Effect |\n", + "|---|---|\n", + "| `coordinate_label` | Name of the new vertex coordinate (e.g. `'parcel_aal3'`) |\n", + "| `label_mapping` | `{int: str}` dict mapping numeric voxel labels to region names |\n", + "| `voxel_label_niftii` | Path to the NIfTI atlas file |\n", + "| `voxel_label_crs` | MNI variant of the NIfTI: `'mni152'` (default) or `'mni305'` |\n", + "| `mni_eps` | Search radius in mm (default `5`). Increase if many vertices receive `'Background'` |\n", + "\n", + "Each call returns a **new** head model object; the original is not modified\n", + "(immutable-style API). We chain two calls to assign both atlases in sequence.\n", + "\n", + "Assigning AAL3 and Brodmann labels to the Colin27 head:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "# Colin27\n", + "\n", + "colin_ijk_labeled = colin_ijk.assign_parcels_via_mni_coords(\n", + " coordinate_label=\"parcel_aal3\",\n", + " label_mapping=aal3_num2label,\n", + " voxel_label_niftii=aal3_voxel_label_niftii,\n", + " voxel_label_crs=\"mni152\",\n", + " mni_eps=5\n", + ")\n", + "\n", + "colin_ijk_labeled = colin_ijk_labeled.assign_parcels_via_mni_coords(\n", + " coordinate_label=\"parcel_brodmann\",\n", + " label_mapping=brodmann_num2label,\n", + " voxel_label_niftii=brodmann_voxel_label_niftii,\n", + " voxel_label_crs=\"mni152\",\n", + " mni_eps=5\n", + ")\n", + "\n", + "\n", + "colin_ijk_labeled.brain.vertices" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "Because ICBM152 and Colin27 share the same MNI152 coordinate frame, the *identical*\n", + "NIfTI atlas files and label dictionaries work unchanged for both head models. This\n", + "is the practical benefit of the shared coordinate system: you write the atlas assignment\n", + "code once and it applies to any head model that has MNI152 vertex coordinates.\n", + "\n", + "Assigning AAL3 and Brodmann labels to the ICBM152 head:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "# ICBM-152\n", + "\n", + "icbm_ijk_labeled = icbm_ijk.assign_parcels_via_mni_coords(\n", + " coordinate_label=\"parcel_aal3\",\n", + " label_mapping=aal3_num2label,\n", + " voxel_label_niftii=aal3_voxel_label_niftii,\n", + " voxel_label_crs=\"mni152\",\n", + " mni_eps=5\n", + ")\n", + "\n", + "icbm_ijk_labeled = icbm_ijk_labeled.assign_parcels_via_mni_coords(\n", + " coordinate_label=\"parcel_brodmann\",\n", + " label_mapping=brodmann_num2label,\n", + " voxel_label_niftii=brodmann_voxel_label_niftii,\n", + " voxel_label_crs=\"mni152\",\n", + " mni_eps=5\n", + ")\n", + "\n", + "\n", + "icbm_ijk_labeled.brain.vertices" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "## Plotting Parcellation Schemes\n", + "\n", + "Cedalion's surface plotting functions accept a list of per-vertex colors (RGB tuples\n", + "or matplotlib color specs), one entry per vertex in vertex order. The helper\n", + "`cedalion.vis.anatomy.get_vertex_colors_from_coord` builds this list from a named\n", + "vertex coordinate and a `color_mapping` that translates string labels to colors.\n", + "\n", + "Three `color_mapping` modes are supported:\n", + "\n", + "| `color_mapping` value | Behaviour |\n", + "|---|---|\n", + "| `dict` (label → color) | Each label is looked up; vertices not in the dict get `default_color` |\n", + "| `None` | A random but deterministic color is generated for every unique label |\n", + "| A single color string (e.g. `'r'`) | Every vertex whose label is in `labels` gets that color; the rest get `default_color` |\n", + "\n", + "We first demonstrate with the **Schaefer2018** parcellation that is bundled with the\n", + "head model and has an official color map:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "# example: Schaefer parcel colors\n", + "schaefer_color_dict = cedalion.data.get_colin27_headmodel_files().load_parcel_colors()\n", + "display(schaefer_color_dict)" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "`get_vertex_colors_from_coord` iterates over every brain surface vertex, reads its\n", + "label from the named coordinate, and returns the corresponding color from the mapping.\n", + "Vertices whose label is not found in the mapping receive `default_color` (grey by\n", + "default). The result is a list in the same vertex order as `brain.vertices`, ready to\n", + "pass directly to any of the `cedalion.vis` plotting functions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17", + "metadata": {}, + "outputs": [], + "source": [ + "vertex_colors = get_vertex_colors_from_coord(colin_ijk_labeled.brain, \"parcel\", color_mapping=schaefer_color_dict)\n", + "print(f\"The brain surface has {colin_ijk_labeled.brain.nvertices} vertices.\")\n", + "print(f\"The list vertex_colors has {len(vertex_colors)} entries. These are the first entries:\")\n", + "vertex_colors[:4]" + ] + }, + { + "cell_type": "markdown", + "id": "18", + "metadata": {}, + "source": [ + "`plot_brain_views_grid` renders the brain surface from five standard viewpoints\n", + "(superior, left, right, anterior, posterior) in a single figure.\n", + "We pass the **inflated** surface here — which has the same vertex order as the pial\n", + "surface — so sulcal cortex is unfolded and every parcel is visible at once.\n", + "\n", + "The Schaefer2018 parcellation on the inflated Colin27 cortex:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19", + "metadata": {}, + "outputs": [], + "source": [ + "plot_brain_views_grid(colin_inflated, vertex_colors, reset_camera=True)" + ] + }, + { + "cell_type": "markdown", + "id": "20", + "metadata": {}, + "source": [ + "### AAL3 Labels on Colin27\n", + "\n", + "We now visualize the freshly assigned AAL3 labels. Because AAL3 does not ship with a\n", + "canonical per-region color map, we pass `color_mapping=None` to let\n", + "`get_vertex_colors_from_coord` generate colors automatically. The algorithm uses a\n", + "deterministic golden-ratio hue spacing so colors are consistent across calls.\n", + "\n", + "We show both the **pial** (folded) and the **inflated** surface. The inflated view\n", + "makes buried sulcal cortex visible and allows you to judge where parcel boundaries\n", + "fall relative to major anatomical landmarks. Notice that compared to the Schaefer2018\n", + "map above, some AAL3 borders appear slightly less crisp — this is the expected\n", + "consequence of the volumetric-to-surface transfer (discussed in detail at the end\n", + "of this notebook)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21", + "metadata": {}, + "outputs": [], + "source": [ + "vertex_colors = get_vertex_colors_from_coord(colin_ijk_labeled.brain, \"parcel_aal3\", color_mapping=None, default_color=\"magenta\")\n", + "\n", + "plot_brain_views_grid(colin_ijk_labeled.brain, vertex_colors)\n", + "plot_brain_views_grid(colin_inflated, vertex_colors, reset_camera=True)" + ] + }, + { + "cell_type": "markdown", + "id": "22", + "metadata": {}, + "source": [ + "### Brodmann Labels on Colin27\n", + "\n", + "The same workflow applied to the Brodmann atlas. Brodmann areas are fewer and larger\n", + "than AAL3 regions. On the inflated surface the classical areas are clearly recognisable:\n", + "BA4/6 along the central sulcus (motor/premotor), BA17/18/19 in the occipital lobe\n", + "(visual cortex), and BA44/45 in the left inferior frontal gyrus (Broca's area)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23", + "metadata": {}, + "outputs": [], + "source": [ + "vertex_colors = get_vertex_colors_from_coord(colin_ijk_labeled.brain, \"parcel_brodmann\", color_mapping=None)\n", + "\n", + "plot_brain_views_grid(colin_ijk_labeled.brain, vertex_colors)\n", + "plot_brain_views_grid(colin_inflated, vertex_colors, reset_camera=True)" + ] + }, + { + "cell_type": "markdown", + "id": "24", + "metadata": {}, + "source": [ + "### AAL3 Labels on ICBM152\n", + "\n", + "The same atlas applied to the ICBM152 head model. Any differences from the Colin27\n", + "result reflect genuine anatomical differences between the two templates (single-subject\n", + "average vs. group average) rather than any difference in the atlas or the transfer\n", + "procedure." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "25", + "metadata": {}, + "outputs": [], + "source": [ + "vertex_colors = get_vertex_colors_from_coord(icbm_ijk_labeled.brain, \"parcel_aal3\", color_mapping=None)\n", + "\n", + "plot_brain_views_grid(icbm_ijk_labeled.brain, vertex_colors)\n", + "plot_brain_views_grid(icbm_inflated, vertex_colors, reset_camera=True)" + ] + }, + { + "cell_type": "markdown", + "id": "26", + "metadata": {}, + "source": [ + "### Brodmann Labels on ICBM152" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "27", + "metadata": {}, + "outputs": [], + "source": [ + "vertex_colors = get_vertex_colors_from_coord(icbm_ijk_labeled.brain, \"parcel_brodmann\", color_mapping=None)\n", + "\n", + "plot_brain_views_grid(icbm_ijk_labeled.brain, vertex_colors)\n", + "plot_brain_views_grid(icbm_inflated, vertex_colors, reset_camera=True)" + ] + }, + { + "cell_type": "markdown", + "id": "28", + "metadata": {}, + "source": [ + "## Customizing Colors for ROI Highlighting\n", + "\n", + "In practice you will often want to highlight a small set of regions of interest rather\n", + "than visualize the entire atlas at once. `get_vertex_colors_from_coord` supports two\n", + "convenience patterns for this." + ] + }, + { + "cell_type": "markdown", + "id": "29", + "metadata": {}, + "source": [ + "**Pattern 1 — selected parcels with distinct colors:**\n", + "\n", + "Pass a partial `color_mapping` dict containing only the regions you care about. All\n", + "other vertices receive `default_color` (grey). Use this when the highlighted regions\n", + "should be distinguished from one another by color — for example when comparing two\n", + "adjacent AAL3 frontal regions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "30", + "metadata": { + "tags": [ + "nbsphinx-thumbnail" + ] + }, + "outputs": [], + "source": [ + "# select only two parcels and specify colors (any matplotlib color spec works)\n", + "# parcels not contained in the mapping get the default color\n", + "aal3_color_mapping = {\"Frontal_Inf_Oper_R\" : \"r\", \"Frontal_Inf_Tri_R\" : \"g\"}\n", + "\n", + "vertex_colors = get_vertex_colors_from_coord(colin_ijk_labeled.brain, \"parcel_aal3\", aal3_color_mapping)\n", + "plot_brain_views_grid(colin_ijk_labeled.brain, vertex_colors)\n" + ] + }, + { + "cell_type": "markdown", + "id": "31", + "metadata": {}, + "source": [ + "**Pattern 2 — highlight a set of parcels in a single color:**\n", + "\n", + "Pass a single color string as `color_mapping` and provide a `labels` list. Every\n", + "vertex whose coordinate matches a label in the list gets the specified color; all\n", + "other vertices get `default_color`. This gives the clearest signal-vs-background\n", + "contrast when you want to show, for example, bilateral Broca's area (BA44 + BA45)\n", + "against a neutral grey background." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "32", + "metadata": {}, + "outputs": [], + "source": [ + "# specify only a single color, that will be used for all parcels. Then select only a subset of parcels to color.\n", + "\n", + "vertex_colors = get_vertex_colors_from_coord(\n", + " colin_ijk_labeled.brain,\n", + " \"parcel_brodmann\",\n", + " color_mapping=\"r\",\n", + " labels=[\"right_BA44\", \"right_BA45\"],\n", + ")\n", + "plot_brain_views_grid(colin_ijk_labeled.brain, vertex_colors)" + ] + }, + { + "cell_type": "markdown", + "id": "33", + "metadata": {}, + "source": [ + "## Spot-Check: Verifying Label Assignments at Known MNI Coordinates\n", + "\n", + "As a sanity check we verify that the assigned labels are anatomically plausible at\n", + "three well-known MNI152 locations. We find the nearest brain surface vertex in MNI152\n", + "space and read off its Brodmann label.\n", + "\n", + "| MNI coordinate | Expected region |\n", + "|---|---|\n", + "| `[-10, -90, 0]` | Occipital pole / calcarine sulcus (BA17/18) |\n", + "| `[-40, -20, 50]` | Left somatomotor cortex (BA4/6) |\n", + "| `[ 40, 20, 30]` | Right prefrontal cortex (BA44/45/46) |" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "34", + "metadata": {}, + "outputs": [], + "source": [ + "tests = np.asarray([\n", + " [-10, -90, 0],\n", + " [-40, -20, 50],\n", + " [40, 20, 30],\n", + "], dtype=float)\n", + "\n", + "vertex_coords = icbm_ijk_labeled.get_brain_mni152_coords().pint.dequantify().values\n", + "vertex_labels = icbm_ijk_labeled.brain.vertices.coords[\"parcel_brodmann\"].values\n", + "_, nearest = KDTree(vertex_coords).query(tests)\n", + "\n", + "for mni, vertex_idx in zip(tests, nearest):\n", + " print(f\"MNI {mni.tolist()} -> nearest surface label {vertex_labels[vertex_idx]}\")" + ] + }, + { + "cell_type": "markdown", + "id": "35", + "metadata": {}, + "source": [ + "## Parcel-Level Summary Tables\n", + "\n", + "`TwoSurfaceHeadModel.parcel_summary_from_vertex_coordinate` converts the per-vertex\n", + "atlas labels into a compact table with **one row per Schaefer2018 parcel** (the\n", + "parcellation bundled with the head model). For each Schaefer parcel it reports:\n", + "\n", + "| Column | Meaning |\n", + "|---|---|\n", + "| `atlas_label` | Dominant atlas region covering that parcel (plurality vote over vertices) |\n", + "| `matching_vertices` | Number of vertices whose atlas label equals `atlas_label` |\n", + "| `parcel_vertices` | Total number of vertices in the Schaefer parcel |\n", + "| `fraction_of_parcel` | Coverage fraction (`matching_vertices / parcel_vertices`) |\n", + "| `mni152_r/a/s` | MNI152 centroid of a representative vertex from that parcel |\n", + "\n", + "Medial-wall and background parcels are excluded by default because they do not\n", + "correspond to cortical tissue. The function returns a `pandas.DataFrame` that can be\n", + "saved with `.to_csv(..., sep='\\t')` for use in downstream analyses or supplementary\n", + "tables in manuscripts." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "36", + "metadata": {}, + "outputs": [], + "source": [ + "colin_ijk_labeled.parcel_summary_from_vertex_coordinate(\"parcel_aal3\", \"colin27\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "37", + "metadata": {}, + "outputs": [], + "source": [ + "colin_ijk_labeled.parcel_summary_from_vertex_coordinate(\"parcel_brodmann\", \"colin27\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "38", + "metadata": {}, + "outputs": [], + "source": [ + "icbm_ijk_labeled.parcel_summary_from_vertex_coordinate(\"parcel_aal3\", \"icbm152\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "39", + "metadata": {}, + "outputs": [], + "source": [ + "icbm_ijk_labeled.parcel_summary_from_vertex_coordinate(\"parcel_brodmann\", \"icbm152\")" + ] + }, + { + "cell_type": "markdown", + "id": "40", + "metadata": {}, + "source": [ + "### Optional: Save Summary Tables\n", + "\n", + "The summary tables above are ordinary `pandas.DataFrame` objects. To write them\n", + "to disk, set `SAVE_TABLES = True` in the next cell. Files are saved in the\n", + "`atlas_parcel_summaries/` folder relative to the notebook kernel's current\n", + "working directory.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "41", + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "\n", + "SAVE_TABLES = False\n", + "OUTPUT_DIR = Path(\"atlas_parcel_summaries\")\n", + "\n", + "if SAVE_TABLES:\n", + " OUTPUT_DIR.mkdir(parents=True, exist_ok=True)\n", + " summaries = {\n", + " \"aal3_parcel_summary_colin27.tsv\": colin_ijk_labeled.parcel_summary_from_vertex_coordinate(\"parcel_aal3\", \"colin27\"),\n", + " \"brodmann_parcel_summary_colin27.tsv\": colin_ijk_labeled.parcel_summary_from_vertex_coordinate(\"parcel_brodmann\", \"colin27\"),\n", + " \"aal3_parcel_summary_icbm152.tsv\": icbm_ijk_labeled.parcel_summary_from_vertex_coordinate(\"parcel_aal3\", \"icbm152\"),\n", + " \"brodmann_parcel_summary_icbm152.tsv\": icbm_ijk_labeled.parcel_summary_from_vertex_coordinate(\"parcel_brodmann\", \"icbm152\"),\n", + " }\n", + " for filename, summary in summaries.items():\n", + " summary.to_csv(OUTPUT_DIR / filename, sep=\"\\t\", index=False)\n", + " print(f\"Saved {len(summaries)} summary tables to {OUTPUT_DIR.resolve()}\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "42", + "metadata": {}, + "source": [ + "## Limitations of the Volumetric MNI Approach\n", + "\n", + "### Why parcel borders appear fuzzy on the inflated surface\n", + "\n", + "The volumetric-to-surface label transfer is inherently approximate. The label of a\n", + "brain surface vertex is determined by the *nearest labelled voxel* in 3-D MNI152 space.\n", + "Because sulcal folds can bring two cortically distant regions to within a few\n", + "millimetres of each other in 3-D, vertices near a sulcal wall may be assigned the\n", + "label of the anatomically adjacent — but cortically distant — region instead of their\n", + "own.\n", + "\n", + "This artefact is most visible on the **inflated surface**: inflation changes vertex\n", + "positions in 3-D while leaving labels fixed. Vertices that sit at the fundus (bottom)\n", + "of a deep sulcus can therefore show \"bleeding\" — incorrect labels — at the inflated\n", + "parcel border, even though the label assignment on the folded pial surface appears\n", + "correct.\n", + "\n", + "The `mni_eps` search radius controls the trade-off: a *smaller* radius reduces\n", + "cross-sulcal contamination but increases the number of vertices that receive\n", + "`'Background'` (no labelled voxel found within reach). The default of 5 mm is a\n", + "practical compromise for atlases with 1 mm isotropic voxels.\n", + "\n", + "### When to use FreeSurfer-based labeling instead\n", + "\n", + "For applications where parcel boundary precision matters — for example when comparing\n", + "fine-grained DOT image reconstruction results to specific cortical regions, or when\n", + "building a classifier that relies on exact parcel membership — **surface-native\n", + "labeling via FreeSurfer** is the more accurate alternative.\n", + "\n", + "FreeSurfer (Fischl, 2012) assigns parcels directly\n", + "on the cortical surface mesh by registering each vertex to the `fsaverage` template,\n", + "where parcel boundaries are defined as vertex-level annotations. Because the boundary\n", + "is expressed in the same surface space as the vertex, no voxel-lookup approximation is\n", + "introduced and region borders are exact, even on the inflated surface.\n", + "\n", + "Cedalion provides a complete FreeSurfer-based pipeline for individualized head models\n", + "in\n", + "[43b_individualized_head_models.ipynb](43b_individualized_head_models.ipynb), which\n", + "walks through FreeSurfer cortical reconstruction, registration to `fsaverage`, and\n", + "surface-native parcel annotation. The standard Colin27 and ICBM152 models were built\n", + "with exactly this pipeline; their bundled `parcel` coordinate (Schaefer2018) is a\n", + "surface-native annotation — which is why Schaefer borders look sharper than AAL3 or\n", + "Brodmann borders on the same inflated surface.\n", + "\n", + "**Practical guidance:**\n", + "\n", + "| Use case | Recommended approach |\n", + "|---|---|\n", + "| Quick atlas comparison, ROI definitions, group-level reporting | Volumetric MNI lookup (this notebook) |\n", + "| High-precision boundary analysis, individual MRI available | FreeSurfer surface-native labeling |\n", + "| Bridging a new atlas to an existing surface-native parcellation | Volumetric MNI lookup, compare overlap with Schaefer via `parcel_summary_from_vertex_coordinate` |" + ] + }, + { + "cell_type": "markdown", + "id": "43", + "metadata": {}, + "source": [ + "## References" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "44", + "metadata": {}, + "outputs": [], + "source": [ + "cedalion.bib.dump_to_notebook()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "cedalion_260507", + "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.15" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/cedalion/bibliography/references.bib b/src/cedalion/bibliography/references.bib index 75c83b59..a29ea4a8 100644 --- a/src/cedalion/bibliography/references.bib +++ b/src/cedalion/bibliography/references.bib @@ -594,4 +594,25 @@ @article{Luke2025 doi = {10.1038/s41597-024-04136-9}, } +@article{Rolls2020, + title = {Automated anatomical labelling atlas 3}, + journal = {NeuroImage}, + volume = {206}, + pages = {116189}, + year = {2020}, + issn = {1053-8119}, + doi = {https://doi.org/10.1016/j.neuroimage.2019.116189}, + author = {Edmund T. Rolls and Chu-Chung Huang and Ching-Po Lin and Jianfeng Feng and Marc Joliot}, +} + +@book{Mai2017, + title = {Human Brain in Standard {{MNI}} Space: Structure and Function: A Comprehensive Pocket Atlas}, + shorttitle = {Human Brain in Standard {{MNI}} Space}, + author = {Mai, J{\"u}rgen K. and Majtanik, Milan}, + year = 2017, + publisher = {Academic Press, an imprint of Elsevier}, + address = {London ; San Diego, CA}, + isbn = {978-0-12-811275-5}, + lccn = {RC386.6.M34 M347 2017}, +} diff --git a/src/cedalion/data/__init__.py b/src/cedalion/data/__init__.py index cec8d8c2..3aff55cd 100644 --- a/src/cedalion/data/__init__.py +++ b/src/cedalion/data/__init__.py @@ -101,6 +101,37 @@ def get(fname: str | Path) -> Path: return files("cedalion.data").joinpath(fname) +@dataclass +class AtlasFiles: + """Paths for a bundled atlas volume and label lookup.""" + + basedir: Path + nifti: Path + labels_json: Path + + +def get_aal3_atlas_files() -> AtlasFiles: + """Return paths to the bundled AAL3 atlas volume and label lookup.""" + + atlas_dir = get("atlases/aal3") + return AtlasFiles( + basedir=atlas_dir, + nifti=atlas_dir / "ROI_MNI_V7_1mm.nii", + labels_json=atlas_dir / "aal_labels.json", + ) + + +def get_brodmann_atlas_files() -> AtlasFiles: + """Return paths to the bundled Brodmann atlas volume and label lookup.""" + + atlas_dir = get("atlases/brodmann") + return AtlasFiles( + basedir=atlas_dir, + nifti=atlas_dir / "Brodmann_Mai_Matajnik.nii", + labels_json=atlas_dir / "brodmann_labels.json", + ) + + def get_ninja_cap_probe(): """Load the fullhead Ninja NIRS cap probe.""" @@ -544,10 +575,12 @@ def get_atlas_files(atlas : str) -> tuple[Path, Path]: fnames = DATASETS.fetch("atlas_aal3.zip", processor=pooch.Unzip()) fname_volume = [Path(i) for i in fnames if i.endswith(".nii")][0] fname_labels = [Path(i) for i in fnames if i.endswith(".json")][0] + cite("Rolls2020") elif atlas == "brodmann": - fnames = DATASETS.fetch("atlas_aal3.zip", processor=pooch.Unzip()) + fnames = DATASETS.fetch("atlas_brodmann.zip", processor=pooch.Unzip()) fname_volume = [Path(i) for i in fnames if i.endswith(".nii")][0] fname_labels = [Path(i) for i in fnames if i.endswith(".json")][0] + cite("Mai2017") else: raise ValueError(f"atlas must be one of {AVAILABLE_ATLASES}") diff --git a/src/cedalion/data/atlases/aal3/ROI_MNI_V7_1mm.nii b/src/cedalion/data/atlases/aal3/ROI_MNI_V7_1mm.nii new file mode 100644 index 00000000..cf9cba4b Binary files /dev/null and b/src/cedalion/data/atlases/aal3/ROI_MNI_V7_1mm.nii differ diff --git a/src/cedalion/data/atlases/aal3/aal_labels.json b/src/cedalion/data/atlases/aal3/aal_labels.json new file mode 100644 index 00000000..474655e6 --- /dev/null +++ b/src/cedalion/data/atlases/aal3/aal_labels.json @@ -0,0 +1,672 @@ +{ + "atlas_name": "AAL3", + "space": "MNI152", + "description": "AAL3 atlas warped to ICBM152 space using SPM coregistration", + "source_nifti": "w_ROI_MNI_V7_1mm.nii", + "labels": [ + { + "index": 1, + "name": "Precentral_L" + }, + { + "index": 2, + "name": "Precentral_R" + }, + { + "index": 3, + "name": "Frontal_Sup_2_L" + }, + { + "index": 4, + "name": "Frontal_Sup_2_R" + }, + { + "index": 5, + "name": "Frontal_Mid_2_L" + }, + { + "index": 6, + "name": "Frontal_Mid_2_R" + }, + { + "index": 7, + "name": "Frontal_Inf_Oper_L" + }, + { + "index": 8, + "name": "Frontal_Inf_Oper_R" + }, + { + "index": 9, + "name": "Frontal_Inf_Tri_L" + }, + { + "index": 10, + "name": "Frontal_Inf_Tri_R" + }, + { + "index": 11, + "name": "Frontal_Inf_Orb_2_L" + }, + { + "index": 12, + "name": "Frontal_Inf_Orb_2_R" + }, + { + "index": 13, + "name": "Rolandic_Oper_L" + }, + { + "index": 14, + "name": "Rolandic_Oper_R" + }, + { + "index": 15, + "name": "Supp_Motor_Area_L" + }, + { + "index": 16, + "name": "Supp_Motor_Area_R" + }, + { + "index": 17, + "name": "Olfactory_L" + }, + { + "index": 18, + "name": "Olfactory_R" + }, + { + "index": 19, + "name": "Frontal_Sup_Medial_L" + }, + { + "index": 20, + "name": "Frontal_Sup_Medial_R" + }, + { + "index": 21, + "name": "Frontal_Med_Orb_L" + }, + { + "index": 22, + "name": "Frontal_Med_Orb_R" + }, + { + "index": 23, + "name": "Rectus_L" + }, + { + "index": 24, + "name": "Rectus_R" + }, + { + "index": 25, + "name": "OFCmed_L" + }, + { + "index": 26, + "name": "OFCmed_R" + }, + { + "index": 27, + "name": "OFCant_L" + }, + { + "index": 28, + "name": "OFCant_R" + }, + { + "index": 29, + "name": "OFCpost_L" + }, + { + "index": 30, + "name": "OFCpost_R" + }, + { + "index": 31, + "name": "OFClat_L" + }, + { + "index": 32, + "name": "OFClat_R" + }, + { + "index": 33, + "name": "Insula_L" + }, + { + "index": 34, + "name": "Insula_R" + }, + { + "index": 37, + "name": "Cingulate_Mid_L" + }, + { + "index": 38, + "name": "Cingulate_Mid_R" + }, + { + "index": 39, + "name": "Cingulate_Post_L" + }, + { + "index": 40, + "name": "Cingulate_Post_R" + }, + { + "index": 41, + "name": "Hippocampus_L" + }, + { + "index": 42, + "name": "Hippocampus_R" + }, + { + "index": 43, + "name": "ParaHippocampal_L" + }, + { + "index": 44, + "name": "ParaHippocampal_R" + }, + { + "index": 45, + "name": "Amygdala_L" + }, + { + "index": 46, + "name": "Amygdala_R" + }, + { + "index": 47, + "name": "Calcarine_L" + }, + { + "index": 48, + "name": "Calcarine_R" + }, + { + "index": 49, + "name": "Cuneus_L" + }, + { + "index": 50, + "name": "Cuneus_R" + }, + { + "index": 51, + "name": "Lingual_L" + }, + { + "index": 52, + "name": "Lingual_R" + }, + { + "index": 53, + "name": "Occipital_Sup_L" + }, + { + "index": 54, + "name": "Occipital_Sup_R" + }, + { + "index": 55, + "name": "Occipital_Mid_L" + }, + { + "index": 56, + "name": "Occipital_Mid_R" + }, + { + "index": 57, + "name": "Occipital_Inf_L" + }, + { + "index": 58, + "name": "Occipital_Inf_R" + }, + { + "index": 59, + "name": "Fusiform_L" + }, + { + "index": 60, + "name": "Fusiform_R" + }, + { + "index": 61, + "name": "Postcentral_L" + }, + { + "index": 62, + "name": "Postcentral_R" + }, + { + "index": 63, + "name": "Parietal_Sup_L" + }, + { + "index": 64, + "name": "Parietal_Sup_R" + }, + { + "index": 65, + "name": "Parietal_Inf_L" + }, + { + "index": 66, + "name": "Parietal_Inf_R" + }, + { + "index": 67, + "name": "SupraMarginal_L" + }, + { + "index": 68, + "name": "SupraMarginal_R" + }, + { + "index": 69, + "name": "Angular_L" + }, + { + "index": 70, + "name": "Angular_R" + }, + { + "index": 71, + "name": "Precuneus_L" + }, + { + "index": 72, + "name": "Precuneus_R" + }, + { + "index": 73, + "name": "Paracentral_Lobule_L" + }, + { + "index": 74, + "name": "Paracentral_Lobule_R" + }, + { + "index": 75, + "name": "Caudate_L" + }, + { + "index": 76, + "name": "Caudate_R" + }, + { + "index": 77, + "name": "Putamen_L" + }, + { + "index": 78, + "name": "Putamen_R" + }, + { + "index": 79, + "name": "Pallidum_L" + }, + { + "index": 80, + "name": "Pallidum_R" + }, + { + "index": 83, + "name": "Heschl_L" + }, + { + "index": 84, + "name": "Heschl_R" + }, + { + "index": 85, + "name": "Temporal_Sup_L" + }, + { + "index": 86, + "name": "Temporal_Sup_R" + }, + { + "index": 87, + "name": "Temporal_Pole_Sup_L" + }, + { + "index": 88, + "name": "Temporal_Pole_Sup_R" + }, + { + "index": 89, + "name": "Temporal_Mid_L" + }, + { + "index": 90, + "name": "Temporal_Mid_R" + }, + { + "index": 91, + "name": "Temporal_Pole_Mid_L" + }, + { + "index": 92, + "name": "Temporal_Pole_Mid_R" + }, + { + "index": 93, + "name": "Temporal_Inf_L" + }, + { + "index": 94, + "name": "Temporal_Inf_R" + }, + { + "index": 95, + "name": "Cerebelum_Crus1_L" + }, + { + "index": 96, + "name": "Cerebelum_Crus1_R" + }, + { + "index": 97, + "name": "Cerebelum_Crus2_L" + }, + { + "index": 98, + "name": "Cerebelum_Crus2_R" + }, + { + "index": 99, + "name": "Cerebelum_3_L" + }, + { + "index": 100, + "name": "Cerebelum_3_R" + }, + { + "index": 101, + "name": "Cerebelum_4_5_L" + }, + { + "index": 102, + "name": "Cerebelum_4_5_R" + }, + { + "index": 103, + "name": "Cerebelum_6_L" + }, + { + "index": 104, + "name": "Cerebelum_6_R" + }, + { + "index": 105, + "name": "Cerebelum_7b_L" + }, + { + "index": 106, + "name": "Cerebelum_7b_R" + }, + { + "index": 107, + "name": "Cerebelum_8_L" + }, + { + "index": 108, + "name": "Cerebelum_8_R" + }, + { + "index": 109, + "name": "Cerebelum_9_L" + }, + { + "index": 110, + "name": "Cerebelum_9_R" + }, + { + "index": 111, + "name": "Cerebelum_10_L" + }, + { + "index": 112, + "name": "Cerebellum_10_R" + }, + { + "index": 113, + "name": "Vermis_1_2" + }, + { + "index": 114, + "name": "Vermis_3" + }, + { + "index": 115, + "name": "Vermis_4_5" + }, + { + "index": 116, + "name": "Vermis_6" + }, + { + "index": 117, + "name": "Vermis_7" + }, + { + "index": 118, + "name": "Vermis_8" + }, + { + "index": 119, + "name": "Vermis_9" + }, + { + "index": 120, + "name": "Vermis_10" + }, + { + "index": 121, + "name": "Thal_AV_L" + }, + { + "index": 122, + "name": "Thal_AV_R" + }, + { + "index": 123, + "name": "Thal_LP_L" + }, + { + "index": 124, + "name": "Thal_LP_R" + }, + { + "index": 125, + "name": "Thal_VA_L" + }, + { + "index": 126, + "name": "Thal_VA_R" + }, + { + "index": 127, + "name": "Thal_VL_L" + }, + { + "index": 128, + "name": "Thal_VL_R" + }, + { + "index": 129, + "name": "Thal_VPL_L" + }, + { + "index": 130, + "name": "Thal_VPL_R" + }, + { + "index": 131, + "name": "Thal_IL_L" + }, + { + "index": 132, + "name": "Thal_IL_R" + }, + { + "index": 133, + "name": "Thal_Re_L" + }, + { + "index": 134, + "name": "Thal_Re_R" + }, + { + "index": 135, + "name": "Thal_MDm_L" + }, + { + "index": 136, + "name": "Thal_MDm_R" + }, + { + "index": 137, + "name": "Thal_MDl_L" + }, + { + "index": 138, + "name": "Thal_MDl_R" + }, + { + "index": 139, + "name": "Thal_LGN_L" + }, + { + "index": 140, + "name": "Thal_LGN_R" + }, + { + "index": 141, + "name": "Thal_MGN_L" + }, + { + "index": 142, + "name": "Thal_MGN_R" + }, + { + "index": 143, + "name": "Thal_PuI_L" + }, + { + "index": 144, + "name": "Thal_PuI_R" + }, + { + "index": 145, + "name": "Thal_PuM_L" + }, + { + "index": 146, + "name": "Thal_PuM_R" + }, + { + "index": 147, + "name": "Thal_PuA_L" + }, + { + "index": 148, + "name": "Thal_PuA_R" + }, + { + "index": 149, + "name": "Thal_PuL_L" + }, + { + "index": 150, + "name": "Thal_PuL_R" + }, + { + "index": 151, + "name": "ACC_sub_L" + }, + { + "index": 152, + "name": "ACC_sub_R" + }, + { + "index": 153, + "name": "ACC_pre_L" + }, + { + "index": 154, + "name": "ACC_pre_R" + }, + { + "index": 155, + "name": "ACC_sup_L" + }, + { + "index": 156, + "name": "ACC_sup_R" + }, + { + "index": 157, + "name": "Vent_Str_L" + }, + { + "index": 158, + "name": "Vent_Str_R" + }, + { + "index": 159, + "name": "VTA_L" + }, + { + "index": 160, + "name": "VTA_R" + }, + { + "index": 161, + "name": "SN_pc_L" + }, + { + "index": 162, + "name": "SN_pc_R" + }, + { + "index": 163, + "name": "SN_pr_L" + }, + { + "index": 164, + "name": "SN_pr_R" + }, + { + "index": 165, + "name": "Red_N_L" + }, + { + "index": 166, + "name": "Red_N_R" + }, + { + "index": 167, + "name": "LC_L" + }, + { + "index": 168, + "name": "LC_R" + }, + { + "index": 169, + "name": "Raphe_D" + }, + { + "index": 170, + "name": "Raphe_M" + } + ] +} \ No newline at end of file diff --git a/src/cedalion/data/atlases/brodmann/Brodmann_Mai_Matajnik.nii b/src/cedalion/data/atlases/brodmann/Brodmann_Mai_Matajnik.nii new file mode 100644 index 00000000..339092c5 Binary files /dev/null and b/src/cedalion/data/atlases/brodmann/Brodmann_Mai_Matajnik.nii differ diff --git a/src/cedalion/data/atlases/brodmann/brodmann_labels.json b/src/cedalion/data/atlases/brodmann/brodmann_labels.json new file mode 100644 index 00000000..aeeaca52 --- /dev/null +++ b/src/cedalion/data/atlases/brodmann/brodmann_labels.json @@ -0,0 +1,524 @@ +{ + "atlas_name": "Brodmann_Mai_Matajnik", + "space": "ICBM152", + "description": "Brodmann Mai-Matajnik atlas used directly for surface label sampling", + "source_nifti": "Brodmann_Mai_Matajnik.nii", + "labels": [ + { + "index": 1, + "name": "right_Ent", + "hemisphere": "right", + "ba_number": 0 + }, + { + "index": 2, + "name": "right_BA20", + "hemisphere": "right", + "ba_number": 20 + }, + { + "index": 3, + "name": "right_BA38", + "hemisphere": "right", + "ba_number": 38 + }, + { + "index": 4, + "name": "right_BA21", + "hemisphere": "right", + "ba_number": 21 + }, + { + "index": 5, + "name": "right_BA36", + "hemisphere": "right", + "ba_number": 36 + }, + { + "index": 6, + "name": "right_BA11", + "hemisphere": "right", + "ba_number": 11 + }, + { + "index": 7, + "name": "right_BA37", + "hemisphere": "right", + "ba_number": 37 + }, + { + "index": 8, + "name": "right_BA25", + "hemisphere": "right", + "ba_number": 25 + }, + { + "index": 9, + "name": "right_BA12", + "hemisphere": "right", + "ba_number": 12 + }, + { + "index": 10, + "name": "right_BA19", + "hemisphere": "right", + "ba_number": 19 + }, + { + "index": 11, + "name": "right_BA47", + "hemisphere": "right", + "ba_number": 47 + }, + { + "index": 12, + "name": "right_BA22", + "hemisphere": "right", + "ba_number": 22 + }, + { + "index": 13, + "name": "right_BA18", + "hemisphere": "right", + "ba_number": 18 + }, + { + "index": 14, + "name": "right_BA10", + "hemisphere": "right", + "ba_number": 10 + }, + { + "index": 15, + "name": "right_BA17", + "hemisphere": "right", + "ba_number": 17 + }, + { + "index": 16, + "name": "right_BA32", + "hemisphere": "right", + "ba_number": 32 + }, + { + "index": 17, + "name": "right_BA24", + "hemisphere": "right", + "ba_number": 24 + }, + { + "index": 18, + "name": "right_BA46", + "hemisphere": "right", + "ba_number": 46 + }, + { + "index": 19, + "name": "right_BA45", + "hemisphere": "right", + "ba_number": 45 + }, + { + "index": 20, + "name": "right_BA33", + "hemisphere": "right", + "ba_number": 33 + }, + { + "index": 21, + "name": "right_BA44", + "hemisphere": "right", + "ba_number": 44 + }, + { + "index": 22, + "name": "right_BA6", + "hemisphere": "right", + "ba_number": 6 + }, + { + "index": 23, + "name": "right_BA23", + "hemisphere": "right", + "ba_number": 23 + }, + { + "index": 24, + "name": "right_BA42", + "hemisphere": "right", + "ba_number": 42 + }, + { + "index": 25, + "name": "right_BA43", + "hemisphere": "right", + "ba_number": 43 + }, + { + "index": 26, + "name": "right_BA30", + "hemisphere": "right", + "ba_number": 30 + }, + { + "index": 27, + "name": "right_BA41", + "hemisphere": "right", + "ba_number": 41 + }, + { + "index": 28, + "name": "right_BA29", + "hemisphere": "right", + "ba_number": 29 + }, + { + "index": 29, + "name": "right_BA26", + "hemisphere": "right", + "ba_number": 26 + }, + { + "index": 30, + "name": "right_BA31", + "hemisphere": "right", + "ba_number": 31 + }, + { + "index": 31, + "name": "right_BA4", + "hemisphere": "right", + "ba_number": 4 + }, + { + "index": 32, + "name": "right_BA1", + "hemisphere": "right", + "ba_number": 1 + }, + { + "index": 33, + "name": "right_BA9", + "hemisphere": "right", + "ba_number": 9 + }, + { + "index": 34, + "name": "right_BA39", + "hemisphere": "right", + "ba_number": 39 + }, + { + "index": 35, + "name": "right_BA2", + "hemisphere": "right", + "ba_number": 2 + }, + { + "index": 36, + "name": "right_BA3", + "hemisphere": "right", + "ba_number": 3 + }, + { + "index": 37, + "name": "right_BA40", + "hemisphere": "right", + "ba_number": 40 + }, + { + "index": 38, + "name": "right_BA7", + "hemisphere": "right", + "ba_number": 7 + }, + { + "index": 39, + "name": "right_BA8", + "hemisphere": "right", + "ba_number": 8 + }, + { + "index": 40, + "name": "right_BA5", + "hemisphere": "right", + "ba_number": 5 + }, + { + "index": 41, + "name": "right_BA52", + "hemisphere": "right", + "ba_number": 52 + }, + { + "index": 42, + "name": "right_BA35", + "hemisphere": "right", + "ba_number": 35 + }, + { + "index": 43, + "name": "right_BA34", + "hemisphere": "right", + "ba_number": 34 + }, + { + "index": 101, + "name": "left_Ent", + "hemisphere": "left", + "ba_number": 0 + }, + { + "index": 102, + "name": "left_BA20", + "hemisphere": "left", + "ba_number": 20 + }, + { + "index": 103, + "name": "left_BA38", + "hemisphere": "left", + "ba_number": 38 + }, + { + "index": 104, + "name": "left_BA21", + "hemisphere": "left", + "ba_number": 21 + }, + { + "index": 105, + "name": "left_BA36", + "hemisphere": "left", + "ba_number": 36 + }, + { + "index": 106, + "name": "left_BA11", + "hemisphere": "left", + "ba_number": 11 + }, + { + "index": 107, + "name": "left_BA37", + "hemisphere": "left", + "ba_number": 37 + }, + { + "index": 108, + "name": "left_BA25", + "hemisphere": "left", + "ba_number": 25 + }, + { + "index": 109, + "name": "left_BA12", + "hemisphere": "left", + "ba_number": 12 + }, + { + "index": 110, + "name": "left_BA19", + "hemisphere": "left", + "ba_number": 19 + }, + { + "index": 111, + "name": "left_BA47", + "hemisphere": "left", + "ba_number": 47 + }, + { + "index": 112, + "name": "left_BA22", + "hemisphere": "left", + "ba_number": 22 + }, + { + "index": 113, + "name": "left_BA18", + "hemisphere": "left", + "ba_number": 18 + }, + { + "index": 114, + "name": "left_BA10", + "hemisphere": "left", + "ba_number": 10 + }, + { + "index": 115, + "name": "left_BA17", + "hemisphere": "left", + "ba_number": 17 + }, + { + "index": 116, + "name": "left_BA32", + "hemisphere": "left", + "ba_number": 32 + }, + { + "index": 117, + "name": "left_BA24", + "hemisphere": "left", + "ba_number": 24 + }, + { + "index": 118, + "name": "left_BA46", + "hemisphere": "left", + "ba_number": 46 + }, + { + "index": 119, + "name": "left_BA45", + "hemisphere": "left", + "ba_number": 45 + }, + { + "index": 120, + "name": "left_BA33", + "hemisphere": "left", + "ba_number": 33 + }, + { + "index": 121, + "name": "left_BA44", + "hemisphere": "left", + "ba_number": 44 + }, + { + "index": 122, + "name": "left_BA6", + "hemisphere": "left", + "ba_number": 6 + }, + { + "index": 123, + "name": "left_BA23", + "hemisphere": "left", + "ba_number": 23 + }, + { + "index": 124, + "name": "left_BA42", + "hemisphere": "left", + "ba_number": 42 + }, + { + "index": 125, + "name": "left_BA43", + "hemisphere": "left", + "ba_number": 43 + }, + { + "index": 126, + "name": "left_BA30", + "hemisphere": "left", + "ba_number": 30 + }, + { + "index": 127, + "name": "left_BA41", + "hemisphere": "left", + "ba_number": 41 + }, + { + "index": 128, + "name": "left_BA29", + "hemisphere": "left", + "ba_number": 29 + }, + { + "index": 129, + "name": "left_BA26", + "hemisphere": "left", + "ba_number": 26 + }, + { + "index": 130, + "name": "left_BA31", + "hemisphere": "left", + "ba_number": 31 + }, + { + "index": 131, + "name": "left_BA4", + "hemisphere": "left", + "ba_number": 4 + }, + { + "index": 132, + "name": "left_BA1", + "hemisphere": "left", + "ba_number": 1 + }, + { + "index": 133, + "name": "left_BA9", + "hemisphere": "left", + "ba_number": 9 + }, + { + "index": 134, + "name": "left_BA39", + "hemisphere": "left", + "ba_number": 39 + }, + { + "index": 135, + "name": "left_BA2", + "hemisphere": "left", + "ba_number": 2 + }, + { + "index": 136, + "name": "left_BA3", + "hemisphere": "left", + "ba_number": 3 + }, + { + "index": 137, + "name": "left_BA40", + "hemisphere": "left", + "ba_number": 40 + }, + { + "index": 138, + "name": "left_BA7", + "hemisphere": "left", + "ba_number": 7 + }, + { + "index": 139, + "name": "left_BA8", + "hemisphere": "left", + "ba_number": 8 + }, + { + "index": 140, + "name": "left_BA5", + "hemisphere": "left", + "ba_number": 5 + }, + { + "index": 141, + "name": "left_BA52", + "hemisphere": "left", + "ba_number": 52 + }, + { + "index": 142, + "name": "left_BA35", + "hemisphere": "left", + "ba_number": 35 + }, + { + "index": 143, + "name": "left_BA34", + "hemisphere": "left", + "ba_number": 34 + } + ] +} \ No newline at end of file diff --git a/src/cedalion/dot/head_model.py b/src/cedalion/dot/head_model.py index 12ea05b2..5e18183a 100644 --- a/src/cedalion/dot/head_model.py +++ b/src/cedalion/dot/head_model.py @@ -31,7 +31,8 @@ register_general_affine, register_trans_rot_isoscale, register_optodes_spring_icp, - register_identity + register_identity, + SpringICPResult ) from cedalion.geometry.segmentation import ( surface_from_segmentation, @@ -821,7 +822,7 @@ def scale_to_landmarks( Args: target_landmarks: Target landmark positions (e.g. from a digitizer) in any CRS. Must contain the same label subset as the model's - landmarks. + landmarks. mode: method to derive the affine transform. Could be either 'trans_rot_isoscale' or 'general'. See cedalion.geometry.registraion for details. @@ -1036,6 +1037,89 @@ def assign_parcels_via_mni_coords( return hm + def parcel_summary_from_vertex_coordinate( + self, + atlas_coord, + head_model_name: str, + exclude_pattern: str | None = "Background|Medial_Wall", + ): + vertices = self.brain.vertices + coords = vertices.coords + mni152 = self.get_brain_mni152_coords().pint.dequantify().values + + df = pd.DataFrame({ + "vertex": vertices.label.values, + "parcel": coords["parcel"].values, + "atlas_label": coords[atlas_coord].values, + "mni152_r": mni152[:, 0], + "mni152_a": mni152[:, 1], + "mni152_s": mni152[:, 2], + }) + + if "fsaverage_vertex" in coords: + df["fsaverage_vertex"] = coords["fsaverage_vertex"].values + elif "fsaverage_vertex_id" in coords: + df["fsaverage_vertex"] = coords["fsaverage_vertex_id"].values + else: + df["fsaverage_vertex"] = pd.NA + + df["parcel"] = df["parcel"].astype(str) + df = df[df["parcel"] != ""] + + if exclude_pattern is not None: + excluded = df["parcel"].str.contains(exclude_pattern, case=False, na=False) + df = df[~excluded] + + counts = ( + df.groupby(["parcel", "atlas_label"], dropna=False) + .size() + .reset_index(name="matching_vertices") + ) + parcel_totals = ( + df.groupby("parcel").size().rename("parcel_vertices").reset_index() + ) + summary = counts.merge(parcel_totals, on="parcel", how="left") + summary["fraction_of_parcel"] = ( + summary["matching_vertices"] / summary["parcel_vertices"] + ) + summary = summary.sort_values( + ["parcel", "matching_vertices", "atlas_label"], + ascending=[True, False, True], + ).drop_duplicates("parcel") + + representative = df.sort_values("vertex").drop_duplicates( + ["parcel", "atlas_label"] + )[ + [ + "parcel", + "atlas_label", + "vertex", + "fsaverage_vertex", + "mni152_r", + "mni152_a", + "mni152_s", + ] + ] + summary = summary.merge( + representative, on=["parcel", "atlas_label"], how="left" + ) + summary.insert(0, "model", head_model_name) + return summary[ + [ + "model", + "parcel", + "atlas_label", + "vertex", + "fsaverage_vertex", + "mni152_r", + "mni152_a", + "mni152_s", + "matching_vertices", + "parcel_vertices", + "fraction_of_parcel", + ] + ] + @lru_cache def get_standard_headmodel(model : str) -> TwoSurfaceHeadModel: diff --git a/src/cedalion/vis/anatomy/__init__.py b/src/cedalion/vis/anatomy/__init__.py index e44e973a..705525a3 100644 --- a/src/cedalion/vis/anatomy/__init__.py +++ b/src/cedalion/vis/anatomy/__init__.py @@ -1,6 +1,11 @@ """Tools for visualizing data on brain and scalp surface representations.""" -from .brain_and_scalp import plot_brain_in_axes, plot_brain_and_scalp +from .brain_and_scalp import ( + plot_brain_in_axes, + plot_brain_and_scalp, + plot_brain_views_grid, + get_vertex_colors_from_coord, +) from .image_recon import image_recon, image_recon_multi_view, image_recon_view from .montage import plot_montage3D from .optode_selector import OptodeSelector @@ -11,6 +16,8 @@ __all__ = [ "plot_brain_in_axes", "plot_brain_and_scalp", + "plot_brain_views_grid", + "get_vertex_colors_from_coord", "image_recon", "image_recon_multi_view", "image_recon_view", diff --git a/src/cedalion/vis/anatomy/brain_and_scalp.py b/src/cedalion/vis/anatomy/brain_and_scalp.py index ae2bd046..547e8188 100644 --- a/src/cedalion/vis/anatomy/brain_and_scalp.py +++ b/src/cedalion/vis/anatomy/brain_and_scalp.py @@ -1,15 +1,17 @@ -from cedalion.dataclasses import PointType +import sys + +import matplotlib.colors +import matplotlib.pyplot as p import numpy as np import pyvista as pv -import matplotlib.colors +import xarray as xr from matplotlib.typing import ColorType -import cedalion.typing as cdt from numpy.typing import ArrayLike -import cedalion.dataclasses as cdc -import xarray as xr -import sys -import matplotlib.pyplot as p +import cedalion.dataclasses as cdc +import cedalion.typing as cdt +import cedalion.vis.blocks as vbx +from cedalion.dataclasses import PointType def plot_brain_and_scalp( @@ -113,6 +115,7 @@ def plot_brain_in_axes( bad_color: ColorType = [0.7, 0.7, 0.7], cb_label: str = "", camera_pos: ArrayLike | str | None = None, + **kwargs ): """Using pyvista render a brain, colored by a metric, and display it in MPL axes.""" @@ -133,6 +136,14 @@ def plot_brain_in_axes( brain_surface = cdc.VTKSurface.from_trimeshsurface(brain_surface) brain_surface = pv.wrap(brain_surface.mesh) + if "smooth_shading" not in kwargs: + kwargs["smooth_shading"] = True + if "split_sharp_edges" not in kwargs: + kwargs["split_sharp_edges"] = True + if "feature_angle" not in kwargs: + kwargs["feature_angle"] = 90 + + plt = pv.Plotter(off_screen=True) plt.add_mesh( @@ -141,7 +152,7 @@ def plot_brain_in_axes( cmap=cmap, clim=(vmin, vmax), scalar_bar_args={"title": cb_label}, - smooth_shading=True, + **kwargs ) if camera_pos is not None: @@ -181,3 +192,147 @@ def plot_brain_in_axes( ax.xaxis.set_ticks([]) ax.yaxis.set_ticks([]) + +def camera_for_view(center, view, distance=350): + """Return camera parameters for a named orthogonal brain view. + + Args: + center: 3-element array-like with the focal point coordinates (e.g. brain + centroid) in the same units as ``distance``. + view: One of ``"superior"``, ``"left"``, ``"right"``, ``"anterior"``, + ``"posterior"``. + distance: Distance from ``center`` to the camera position along the view + axis. Defaults to 350. + + Returns: + tuple: ``(position, focal_point, up)`` where each element is a numpy + array suitable for assignment to ``pyvista.Camera`` attributes. + """ + cameras = { + "superior": ([0, 0, distance], [0, 1, 0]), + "left": ([-distance, 0, 0], [0, 0, 1]), + "right": ([distance, 0, 0], [0, 0, 1]), + "anterior": ([0, distance, 0], [0, 0, 1]), + "posterior": ([0, -distance, 0], [0, 0, 1]), + } + position_offset, up = cameras[view] + return center + np.asarray(position_offset), center, up + + +def plot_brain_views_grid( + brain_surface, vertex_colors, window_size=(1000, 600), reset_camera=False +): + """Render the brain surface from five standard views in a grid layout. + + Displays superior, anterior, posterior, left, and right views arranged in a + 2-row PyVista plotter window. + + Args: + brain_surface: A surface object whose ``.vertices`` attribute is a + pint-aware xarray with a ``"label"`` dimension. + vertex_colors: Per-vertex color array passed to ``vbx.plot_surface``. + window_size: ``(width, height)`` in pixels for the plotter window. + Defaults to ``(1000, 600)``. + reset_camera: If ``True``, call ``plt.reset_camera()`` after setting + each view to fit the surface tightly. Defaults to ``False``. + """ + brain_center = brain_surface.vertices.pint.dequantify().mean("label").values + plt = pv.Plotter( + shape=(2, 6), + groups=( + (0, slice(0, 2)), + (0, slice(2, 4)), + (0, slice(4, 6)), + (1, slice(0, 3)), + (1, slice(3, 6)), + ), + window_size=window_size, + ) + views = ("superior", "anterior", "posterior", "left", "right") + positions = [(0, 0), (0, 2), (0, 4), (1, 0), (1, 3)] + for view, subplot in zip(views, positions): + plt.subplot(*subplot) + vbx.plot_surface( + plt, + brain_surface, + color=vertex_colors, + ) + plt.add_text(view, font_size=10) + plt.camera.position, plt.camera.focal_point, plt.camera.up = camera_for_view( + brain_center, view + ) + if reset_camera: + plt.reset_camera() + plt.subplot(1, 2) + plt.add_text("", font_size=10) + plt.show() + + +def get_vertex_colors_from_coord( + brain_surface : cdc.TrimeshSurface, + label_coord: str, + color_mapping: dict, + default_color="w", + labels: list[str] | None = None, +): + """Build a per-vertex color list from a named coordinate on the brain surface. + + Each vertex is colored according to the value of ``label_coord`` at that vertex, + looked up in ``color_mapping``. Vertices whose coordinate value is not present + in the mapping (or whose label is filtered out) receive ``default_color``. + + Args: + brain_surface: Surface object whose ``.vertices`` attribute is an xarray + DataArray with named coordinates. + label_coord: Name of the coordinate on ``brain_surface.vertices`` whose + values are used as keys into ``color_mapping``. + color_mapping: Controls how coordinate values map to colors: + * ``None`` — generate a deterministic random color per unique label. + * ``dict`` — map each label to a matplotlib-compatible color spec. + * Any other single color spec assigns the same color to every label. + default_color: Matplotlib-compatible color used for vertices whose label + is absent from the resolved mapping. Defaults to ``"w"`` (white). + labels: If provided, only labels in this list are kept in the mapping; + all other vertices fall back to ``default_color``. + + Returns: + list: One RGB tuple per vertex, in the same order as + ``brain_surface.vertices``. + """ + coords = brain_surface.vertices.coords[label_coord].values + default_color = matplotlib.colors.to_rgb(default_color) + + def normalize_colors(c): + if (isinstance(c, tuple) or isinstance(c, list)) and all( + [isinstance(v, int) for v in c] + ): + c = [k/255. for k in c] + + return matplotlib.colors.to_rgb(c) + + if color_mapping is None: + # generate random colors + rng = np.random.default_rng(43) + color_mapping = { + k: rng.uniform(0.3, 1.0, size=3).tolist() for k in sorted(set(coords)) + } + elif isinstance(color_mapping, dict): + color_mapping = { + k: normalize_colors(v) for k, v in color_mapping.items() + } + elif not isinstance(color_mapping, dict): + # all coord values get the same color + c = matplotlib.colors.to_rgb(color_mapping) + color_mapping = {k: c for k in coords} + elif callable(color_mapping): + # support any kind of mapping + raise not NotImplementedError() + else: + raise ValueError("could not interprete color_mapping") + + if labels is not None: + color_mapping = {k: v for k, v in color_mapping.items() if k in labels} + + vertex_colors = [color_mapping.get(pp, default_color) for pp in coords] + + return vertex_colors diff --git a/src/cedalion/vis/anatomy/image_recon.py b/src/cedalion/vis/anatomy/image_recon.py index 26101d3f..7afbdb0e 100644 --- a/src/cedalion/vis/anatomy/image_recon.py +++ b/src/cedalion/vis/anatomy/image_recon.py @@ -514,7 +514,7 @@ def image_recon_multi_view( ts_title = title_str if view == 'scale_bar' else None p0, surf, lab = image_recon( X_ts, head, cmap=cmap, clim=clim, view_type=view_type, - view_position=view, p0=p0, title_str=ts_title, off_screen=False, + view_position=view, p0=p0, title_str=ts_title, off_screen=(SAVE and filename is not None), plotshape=subplot_shape, iax=iax, wdw_size=wdw_size ) subplots[view] = surf