{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# S1: Head models and Forward Modelling\n", "\n", "This notebook introduces how Cedalion handles head models and forward modelling for diffuse optical tomography." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2026-01-13T10:35:20.602855Z", "iopub.status.busy": "2026-01-13T10:35:20.602771Z", "iopub.status.idle": "2026-01-13T10:35:20.607670Z", "shell.execute_reply": "2026-01-13T10:35:20.607262Z" } }, "outputs": [], "source": [ "# This cells setups the environment when executed in Google Colab.\n", "try:\n", " import google.colab\n", " !curl -s https://raw.githubusercontent.com/ibs-lab/cedalion/dev/scripts/colab_setup.py -o colab_setup.py\n", " # Select branch with --branch \"branch name\" (default is \"dev\")\n", " %run colab_setup.py\n", "except ImportError:\n", " pass" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2026-01-13T10:35:20.608610Z", "iopub.status.busy": "2026-01-13T10:35:20.608537Z", "iopub.status.idle": "2026-01-13T10:35:21.911101Z", "shell.execute_reply": "2026-01-13T10:35:21.909872Z" } }, "outputs": [], "source": [ "from pathlib import Path\n", "from tempfile import TemporaryDirectory\n", "\n", "import cedalion\n", "import cedalion.data\n", "import cedalion.dataclasses as cdc\n", "import cedalion.dot\n", "import cedalion.io\n", "import cedalion.nirs\n", "import cedalion.vis.anatomy\n", "import cedalion.vis.blocks as vbx\n", "import matplotlib.pyplot as p\n", "import numpy as np\n", "import pyvista as pv\n", "import xarray as xr\n", "from cedalion import units\n", "from matplotlib.colors import ListedColormap\n", "\n", "# set pyvista's jupyter backend to 'server' for interactive 3D plots\n", "pv.set_jupyter_backend(\"static\")\n", "# pv.set_jupyter_backend('server')\n", "\n", "xr.set_options(display_expand_data=False);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Configuration \n", "\n", "The constants in the next cell control the behavior of the notebook. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2026-01-13T10:35:21.913309Z", "iopub.status.busy": "2026-01-13T10:35:21.913039Z", "iopub.status.idle": "2026-01-13T10:35:21.915502Z", "shell.execute_reply": "2026-01-13T10:35:21.915128Z" } }, "outputs": [], "source": [ "# set this to either 'colin27', 'icbm152', 'custom_from_segmentation', 'custom_from_surfaces'\n", "HEAD_MODEL = \"colin27\"\n", "\n", "# set this to either 'mcx' or 'nirfaster'\n", "FORWARD_MODEL = \"mcx\"\n", "\n", "# set this to False to actually run the forward model\n", "PRECOMPUTED_SENSITIVITY = True\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create a temporary directory as a working directory." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2026-01-13T10:35:21.916398Z", "iopub.status.busy": "2026-01-13T10:35:21.916323Z", "iopub.status.idle": "2026-01-13T10:35:21.918135Z", "shell.execute_reply": "2026-01-13T10:35:21.917785Z" } }, "outputs": [], "source": [ "\n", "working_directory = TemporaryDirectory()\n", "tmp_dir_path = Path(working_directory.name)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loading and inspecting an example dataset\n", "\n", "This notebook loads one of Cedalion's example datasets. \n", "\n", "Using a function from the `cedalion.data` package, the following cell will download a snirf file with a recording of a fingertapping experiment. \n", "\n", "The content of the snirf file is returned as a `cedalion.dataclasses.Recording` object." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2026-01-13T10:35:21.920402Z", "iopub.status.busy": "2026-01-13T10:35:21.920335Z", "iopub.status.idle": "2026-01-13T10:35:22.549331Z", "shell.execute_reply": "2026-01-13T10:35:22.547672Z" } }, "outputs": [], "source": [ "rec = cedalion.data.get_fingertappingDOT()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Recording container carries time series and related objects through the program in ordered dictionaries much like a segmented tray carrying items in separate compartments.\n", "\n", "This notebook mainly needs the probe geometry and channel definitions, detailed below." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Probe Geometry\n", "\n", "The snirf file contains the 3D probe geometry. When loaded in Cedalion, the probe geometry is represented in a `xarray.DataArray` with shape (number of labeled points x 3). \n", "\n", "It is available as an attribute of the `Recording` container:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2026-01-13T10:35:22.551792Z", "iopub.status.busy": "2026-01-13T10:35:22.551695Z", "iopub.status.idle": "2026-01-13T10:35:22.558779Z", "shell.execute_reply": "2026-01-13T10:35:22.558484Z" } }, "outputs": [ { "data": { "text/html": [ "
<xarray.DataArray (label: 346, digitized: 3)> Size: 8kB\n",
"[mm] -77.82 15.68 23.17 -61.91 21.23 56.49 ... 14.23 -38.28 81.95 -0.678 -37.03\n",
"Coordinates:\n",
" type (label) object 3kB PointType.SOURCE ... PointType.LANDMARK\n",
" * label (label) <U6 8kB 'S1' 'S2' 'S3' 'S4' ... 'FFT10h' 'FT10h' 'FTT10h'\n",
"Dimensions without coordinates: digitized| \n", " | sourceIndex | \n", "detectorIndex | \n", "wavelengthIndex | \n", "wavelengthActual | \n", "wavelengthEmissionActual | \n", "dataType | \n", "dataUnit | \n", "dataTypeLabel | \n", "dataTypeIndex | \n", "sourcePower | \n", "detectorGain | \n", "moduleIndex | \n", "sourceModuleIndex | \n", "detectorModuleIndex | \n", "channel | \n", "source | \n", "detector | \n", "wavelength | \n", "chromo | \n", "
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | \n", "1 | \n", "1 | \n", "1 | \n", "None | \n", "None | \n", "1 | \n", "None | \n", "raw-DC | \n", "1 | \n", "None | \n", "None | \n", "None | \n", "None | \n", "None | \n", "S1D1 | \n", "S1 | \n", "D1 | \n", "760.0 | \n", "None | \n", "
| 1 | \n", "1 | \n", "2 | \n", "1 | \n", "None | \n", "None | \n", "1 | \n", "None | \n", "raw-DC | \n", "1 | \n", "None | \n", "None | \n", "None | \n", "None | \n", "None | \n", "S1D2 | \n", "S1 | \n", "D2 | \n", "760.0 | \n", "None | \n", "
| 2 | \n", "1 | \n", "4 | \n", "1 | \n", "None | \n", "None | \n", "1 | \n", "None | \n", "raw-DC | \n", "1 | \n", "None | \n", "None | \n", "None | \n", "None | \n", "None | \n", "S1D4 | \n", "S1 | \n", "D4 | \n", "760.0 | \n", "None | \n", "
| 3 | \n", "1 | \n", "5 | \n", "1 | \n", "None | \n", "None | \n", "1 | \n", "None | \n", "raw-DC | \n", "1 | \n", "None | \n", "None | \n", "None | \n", "None | \n", "None | \n", "S1D5 | \n", "S1 | \n", "D5 | \n", "760.0 | \n", "None | \n", "
| 4 | \n", "1 | \n", "6 | \n", "1 | \n", "None | \n", "None | \n", "1 | \n", "None | \n", "raw-DC | \n", "1 | \n", "None | \n", "None | \n", "None | \n", "None | \n", "None | \n", "S1D6 | \n", "S1 | \n", "D6 | \n", "760.0 | \n", "None | \n", "
| ... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "
| 195 | \n", "14 | \n", "27 | \n", "2 | \n", "None | \n", "None | \n", "1 | \n", "None | \n", "raw-DC | \n", "1 | \n", "None | \n", "None | \n", "None | \n", "None | \n", "None | \n", "S14D27 | \n", "S14 | \n", "D27 | \n", "850.0 | \n", "None | \n", "
| 196 | \n", "14 | \n", "28 | \n", "2 | \n", "None | \n", "None | \n", "1 | \n", "None | \n", "raw-DC | \n", "1 | \n", "None | \n", "None | \n", "None | \n", "None | \n", "None | \n", "S14D28 | \n", "S14 | \n", "D28 | \n", "850.0 | \n", "None | \n", "
| 197 | \n", "14 | \n", "29 | \n", "2 | \n", "None | \n", "None | \n", "1 | \n", "None | \n", "raw-DC | \n", "1 | \n", "None | \n", "None | \n", "None | \n", "None | \n", "None | \n", "S14D29 | \n", "S14 | \n", "D29 | \n", "850.0 | \n", "None | \n", "
| 198 | \n", "14 | \n", "31 | \n", "2 | \n", "None | \n", "None | \n", "1 | \n", "None | \n", "raw-DC | \n", "1 | \n", "None | \n", "None | \n", "None | \n", "None | \n", "None | \n", "S14D31 | \n", "S14 | \n", "D31 | \n", "850.0 | \n", "None | \n", "
| 199 | \n", "14 | \n", "32 | \n", "2 | \n", "None | \n", "None | \n", "1 | \n", "None | \n", "raw-DC | \n", "1 | \n", "None | \n", "None | \n", "None | \n", "None | \n", "None | \n", "S14D32 | \n", "S14 | \n", "D32 | \n", "850.0 | \n", "None | \n", "
200 rows × 19 columns
\n", "" ], "text/plain": [ " sourceIndex detectorIndex wavelengthIndex wavelengthActual \\\n", "0 1 1 1 None \n", "1 1 2 1 None \n", "2 1 4 1 None \n", "3 1 5 1 None \n", "4 1 6 1 None \n", ".. ... ... ... ... \n", "195 14 27 2 None \n", "196 14 28 2 None \n", "197 14 29 2 None \n", "198 14 31 2 None \n", "199 14 32 2 None \n", "\n", " wavelengthEmissionActual dataType dataUnit dataTypeLabel dataTypeIndex \\\n", "0 None 1 None raw-DC 1 \n", "1 None 1 None raw-DC 1 \n", "2 None 1 None raw-DC 1 \n", "3 None 1 None raw-DC 1 \n", "4 None 1 None raw-DC 1 \n", ".. ... ... ... ... ... \n", "195 None 1 None raw-DC 1 \n", "196 None 1 None raw-DC 1 \n", "197 None 1 None raw-DC 1 \n", "198 None 1 None raw-DC 1 \n", "199 None 1 None raw-DC 1 \n", "\n", " sourcePower detectorGain moduleIndex sourceModuleIndex \\\n", "0 None None None None \n", "1 None None None None \n", "2 None None None None \n", "3 None None None None \n", "4 None None None None \n", ".. ... ... ... ... \n", "195 None None None None \n", "196 None None None None \n", "197 None None None None \n", "198 None None None None \n", "199 None None None None \n", "\n", " detectorModuleIndex channel source detector wavelength chromo \n", "0 None S1D1 S1 D1 760.0 None \n", "1 None S1D2 S1 D2 760.0 None \n", "2 None S1D4 S1 D4 760.0 None \n", "3 None S1D5 S1 D5 760.0 None \n", "4 None S1D6 S1 D6 760.0 None \n", ".. ... ... ... ... ... ... \n", "195 None S14D27 S14 D27 850.0 None \n", "196 None S14D28 S14 D28 850.0 None \n", "197 None S14D29 S14 D29 850.0 None \n", "198 None S14D31 S14 D31 850.0 None \n", "199 None S14D32 S14 D32 850.0 None \n", "\n", "[200 rows x 19 columns]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "meas_list = rec._measurement_lists[\"amp\"]\n", "display(meas_list)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `'amp'` time series is represented as a `xarray.DataArray` with dimensions `'channel'`, `'wavelength'` and `'time'`. Three coordinate arrays are linked to the `'channel'` dimension, specifying for each channel a string label as well as the string label and the corresponding source and detector labels (e.g. the first channel `'S1D1'` is between source `'S1'` and detector `'D1'`)." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "execution": { "iopub.execute_input": "2026-01-13T10:35:22.587009Z", "iopub.status.busy": "2026-01-13T10:35:22.586921Z", "iopub.status.idle": "2026-01-13T10:35:22.594191Z", "shell.execute_reply": "2026-01-13T10:35:22.593894Z" } }, "outputs": [ { "data": { "text/html": [ "<xarray.DataArray (channel: 100, wavelength: 2, time: 8794)> Size: 14MB\n",
"[V] 0.0874 0.08735 0.08819 0.08887 0.0879 ... 0.09108 0.09037 0.09043 0.08999\n",
"Coordinates:\n",
" * time (time) float64 70kB 0.0 0.2294 0.4588 ... 2.017e+03 2.017e+03\n",
" samples (time) int64 70kB 0 1 2 3 4 5 ... 8788 8789 8790 8791 8792 8793\n",
" * channel (channel) object 800B 'S1D1' 'S1D2' 'S1D4' ... 'S14D31' 'S14D32'\n",
" source (channel) object 800B 'S1' 'S1' 'S1' 'S1' ... 'S14' 'S14' 'S14'\n",
" detector (channel) object 800B 'D1' 'D2' 'D4' 'D5' ... 'D29' 'D31' 'D32'\n",
" * wavelength (wavelength) float64 16B 760.0 850.0\n",
"Attributes:\n",
" data_type_group: unprocessed raw<xarray.DataArray (segmentation_type: 5, i: 181, j: 217, k: 181)> Size: 36MB\n",
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
"Coordinates:\n",
" * segmentation_type (segmentation_type) <U5 100B 'csf' 'gm' ... 'skull' 'wm'\n",
"Dimensions without coordinates: i, j, k<xarray.DataArray (label: 73, ijk: 3)> Size: 2kB\n",
"[] 90.0 7.435 48.92 3.924 106.0 24.01 ... 127.4 25.63 109.0 139.8 26.66 93.47\n",
"Coordinates:\n",
" * label (label) <U3 876B 'Iz' 'LPA' 'RPA' 'Nz' ... 'PO1' 'PO2' 'PO4' 'PO6'\n",
" type (label) object 584B PointType.LANDMARK ... PointType.LANDMARK\n",
"Dimensions without coordinates: ijk<xarray.DataArray (mni: 4, ijk: 4)> Size: 128B\n",
"[mm] 1.0 0.0 0.0 -90.0 0.0 1.0 0.0 -126.0 0.0 0.0 1.0 -72.0 0.0 0.0 0.0 1.0\n",
"Dimensions without coordinates: mni, ijk"
],
"text/plain": [
"<xarray.DataArray (ijk: 4, mni: 4)> Size: 128B\n",
"[1/mm] 1.0 -1.616e-14 -8.192e-15 90.0 ... 3.253e-16 -1.74e-16 -8.828e-17 1.0\n",
"Dimensions without coordinates: ijk, mni<xarray.DataArray (channel: 100, vertex: 25052, wavelength: 2)> Size: 40MB\n",
"1.212e-17 1.212e-17 3.962e-20 3.962e-20 ... 1.96e-18 3.815e-16 3.815e-16\n",
"Coordinates:\n",
" parcel (vertex) object 200kB 'VisCent_ExStr_8_LH' ... 'scalp'\n",
" is_brain (vertex) bool 25kB True True True True ... False False False\n",
" * channel (channel) object 800B 'S1D1' 'S1D2' 'S1D4' ... 'S14D31' 'S14D32'\n",
" source (channel) object 800B 'S1' 'S1' 'S1' 'S1' ... 'S14' 'S14' 'S14'\n",
" detector (channel) object 800B 'D1' 'D2' 'D4' 'D5' ... 'D29' 'D31' 'D32'\n",
" * wavelength (wavelength) float64 16B 760.0 850.0\n",
"Dimensions without coordinates: vertex\n",
"Attributes:\n",
" units: mm