{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Xarray Data Structures - an fNIRS example\n", "\n", "This example illustrates the usage of xarray-based data structures for calculating the Beer-Lambert transformation." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2026-06-15T14:49:24.745627Z", "iopub.status.busy": "2026-06-15T14:49:24.745467Z", "iopub.status.idle": "2026-06-15T14:49:24.751384Z", "shell.execute_reply": "2026-06-15T14:49:24.750726Z" } }, "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-06-15T14:49:24.752920Z", "iopub.status.busy": "2026-06-15T14:49:24.752783Z", "iopub.status.idle": "2026-06-15T14:49:28.774918Z", "shell.execute_reply": "2026-06-15T14:49:28.774059Z" } }, "outputs": [], "source": [ "import cedalion\n", "import cedalion.nirs\n", "import cedalion.xrutils\n", "import cedalion.xrutils as xrutils\n", "from cedalion.data import get_fingertapping_snirf_path\n", "import numpy as np\n", "import xarray as xr\n", "import pint\n", "import matplotlib.pyplot as p\n", "import scipy.signal\n", "import os.path\n", "xr.set_options(display_max_rows=3, display_values_threshold=50)\n", "np.set_printoptions(precision=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Loading raw CW-NIRS data from a SNIRF file\n", "This notebook uses a finger-tapping dataset in BIDS layout provided by [Rob Luke](https://github.com/rob-luke/BIDS-NIRS-Tapping). It can can be downloaded via `cedalion.data`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Load amplitude data from the snirf file." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2026-06-15T14:49:28.777553Z", "iopub.status.busy": "2026-06-15T14:49:28.777146Z", "iopub.status.idle": "2026-06-15T14:49:29.052115Z", "shell.execute_reply": "2026-06-15T14:49:29.051296Z" } }, "outputs": [], "source": [ "path_to_snirf_file = get_fingertapping_snirf_path()\n", "\n", "recordings = cedalion.io.read_snirf(path_to_snirf_file)\n", "rec = recordings[0] # there is only one NirsElement in this snirf file...\n", "amp = rec[\"amp\"] # ... which holds amplitude data\n", "\n", "# restrict to first 60 seconds and fill in missing units\n", "amp = amp.sel(time=amp.time < 60)\n", "amp = amp.pint.dequantify().pint.quantify(\"V\")\n", "geo3d = rec.geo3d\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2026-06-15T14:49:29.054352Z", "iopub.status.busy": "2026-06-15T14:49:29.054090Z", "iopub.status.idle": "2026-06-15T14:49:29.061216Z", "shell.execute_reply": "2026-06-15T14:49:29.060475Z" } }, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "recordings" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Amplitude data" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2026-06-15T14:49:29.097609Z", "iopub.status.busy": "2026-06-15T14:49:29.097427Z", "iopub.status.idle": "2026-06-15T14:49:29.111734Z", "shell.execute_reply": "2026-06-15T14:49:29.110880Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (channel: 28, wavelength: 2, time: 469)> Size: 210kB\n",
       "<Quantity([[[0.0914 0.091  0.091  ... 0.0903 0.0902 0.0899]\n",
       "  [0.1857 0.1864 0.1837 ... 0.1849 0.185  0.1847]]\n",
       "\n",
       " [[0.2275 0.2297 0.2261 ... 0.2241 0.2243 0.2257]\n",
       "  [0.6355 0.6377 0.6298 ... 0.6223 0.6237 0.6272]]\n",
       "\n",
       " [[0.1065 0.1066 0.1053 ... 0.1065 0.1062 0.1056]\n",
       "  [0.2755 0.2762 0.2727 ... 0.2737 0.2742 0.276 ]]\n",
       "\n",
       " ...\n",
       "\n",
       " [[0.2028 0.1997 0.2005 ... 0.1998 0.2007 0.2026]\n",
       "  [0.4666 0.4554 0.4562 ... 0.4482 0.4511 0.4541]]\n",
       "\n",
       " [[0.4885 0.4802 0.4818 ... 0.5005 0.5036 0.5045]\n",
       "  [0.8458 0.826  0.826  ... 0.8386 0.8441 0.8475]]\n",
       "\n",
       " [[0.6305 0.6284 0.6287 ... 0.6373 0.638  0.6392]\n",
       "  [1.2286 1.2206 1.219  ... 1.2232 1.2259 1.2278]]], 'volt')>\n",
       "Coordinates: (3/6)\n",
       "  * time        (time) float64 4kB 0.0 0.128 0.256 0.384 ... 59.65 59.78 59.9\n",
       "    samples     (time) int64 4kB 0 1 2 3 4 5 6 7 ... 462 463 464 465 466 467 468\n",
       "    ...          ...\n",
       "  * wavelength  (wavelength) float64 16B 760.0 850.0\n",
       "Attributes:\n",
       "    data_type_group:  unprocessed raw
" ], "text/plain": [ " Size: 210kB\n", "\n", "Coordinates: (3/6)\n", " * time (time) float64 4kB 0.0 0.128 0.256 0.384 ... 59.65 59.78 59.9\n", " samples (time) int64 4kB 0 1 2 3 4 5 6 7 ... 462 463 464 465 466 467 468\n", " ... ...\n", " * wavelength (wavelength) float64 16B 760.0 850.0\n", "Attributes:\n", " data_type_group: unprocessed raw" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "display(amp.round(4))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Montage information" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `geo3d` DataArray maps labels to 3D positions, thus storing the location of optodes and landmarks." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2026-06-15T14:49:29.113532Z", "iopub.status.busy": "2026-06-15T14:49:29.113387Z", "iopub.status.idle": "2026-06-15T14:49:29.123122Z", "shell.execute_reply": "2026-06-15T14:49:29.122351Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (label: 5, pos: 3)> Size: 120B\n",
       "<Quantity([[-0.0416  0.0268  0.1299]\n",
       " [-0.0648  0.0581  0.0908]\n",
       " [-0.0376  0.0632  0.1157]\n",
       " [-0.0413 -0.0118  0.1349]\n",
       " [ 0.      0.114  -0.    ]], 'meter')>\n",
       "Coordinates:\n",
       "  * label    (label) <U6 120B 'S1' 'S2' 'D1' 'D2' 'NASION'\n",
       "    type     (label) object 40B PointType.SOURCE ... PointType.LANDMARK\n",
       "Dimensions without coordinates: pos
" ], "text/plain": [ " Size: 120B\n", "\n", "Coordinates:\n", " * label (label) \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (channel: 28)> Size: 224B\n",
       "<Quantity([0.039 0.039 0.041 0.008 0.037 0.038 0.037 0.007 0.04  0.037 0.008 0.041\n",
       " 0.034 0.008 0.039 0.039 0.041 0.008 0.037 0.037 0.037 0.008 0.04  0.037\n",
       " 0.007 0.041 0.033 0.008], 'meter')>\n",
       "Coordinates:\n",
       "  * channel   (channel) object 224B 'S1D1' 'S1D2' 'S1D3' ... 'S8D8' 'S8D16'\n",
       "    source    (channel) object 224B 'S1' 'S1' 'S1' 'S1' ... 'S7' 'S8' 'S8' 'S8'\n",
       "    detector  (channel) object 224B 'D1' 'D2' 'D3' 'D9' ... 'D7' 'D8' 'D16'
" ], "text/plain": [ " Size: 224B\n", "\n", "Coordinates:\n", " * channel (channel) object 224B 'S1D1' 'S1D2' 'S1D3' ... 'S8D8' 'S8D16'\n", " source (channel) object 224B 'S1' 'S1' 'S1' 'S1' ... 'S7' 'S8' 'S8' 'S8'\n", " detector (channel) object 224B 'D1' 'D2' 'D3' 'D9' ... 'D7' 'D8' 'D16'" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "dists = xrutils.norm(geo3d.loc[amp.source] - geo3d.loc[amp.detector], dim=\"pos\")\n", "display(dists.round(3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Beer-Lambert transformation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Beer-Lambert transformation\n", "\n", "The modified Beer-Lambert law (mBLL) relates changes in optical density (OD) at two wavelengths to changes in oxy- and deoxyhaemoglobin (HbO / HbR) concentration:\n", "\n", "$$\\Delta c = (\\mathbf{E} \\cdot \text{DPF} \\cdot d)^{-1} \\Delta\text{OD}$$\n", "\n", "where:\n", "- $\\mathbf{E}$ is the **extinction coefficient matrix** — tabulated absorptivity of HbO and HbR at each wavelength (unit: 1/(mm·mM)).\n", "- $d$ is the **source–detector distance** (mm), derived from optode positions in `geo3d`.\n", "- **DPF** (differential path length factor) corrects for the longer effective optical path of diffusely scattered photons. A value of 6 is a common approximation for adult scalp tissue in the 700–900 nm range.\n", "\n", "Below we assemble these quantities manually to illustrate xarray's matrix inversion and broadcasting. Cedalion's convenience functions `cedalion.nirs.cw.int2od()` and `cedalion.nirs.cw.od2conc()` wrap the same steps." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "execution": { "iopub.execute_input": "2026-06-15T14:49:29.145449Z", "iopub.status.busy": "2026-06-15T14:49:29.145285Z", "iopub.status.idle": "2026-06-15T14:49:29.173005Z", "shell.execute_reply": "2026-06-15T14:49:29.172294Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (wavelength: 2, chromo: 2)> Size: 32B\n",
       "<Quantity([[-0.0024  0.0037]\n",
       " [ 0.0055 -0.0021]], 'millimeter * molar')>\n",
       "Coordinates:\n",
       "  * chromo      (chromo) <U3 24B 'HbO' 'HbR'\n",
       "  * wavelength  (wavelength) float64 16B 760.0 850.0
" ], "text/plain": [ " Size: 32B\n", "\n", "Coordinates:\n", " * chromo (chromo) \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (chromo: 2, channel: 28, time: 469)> Size: 210kB\n",
       "<Quantity([[[-0.0538 -0.1835  0.1612 ... -0.0725 -0.1064 -0.0939]\n",
       "  [-0.1971 -0.1772 -0.0504 ...  0.1388  0.0943  0.0257]\n",
       "  [ 0.0052 -0.0338  0.1269 ...  0.1558  0.0805 -0.1094]\n",
       "  ...\n",
       "  [-0.4716 -0.084  -0.0781 ...  0.2843  0.1805  0.1309]\n",
       "  [-0.6011 -0.1653 -0.124  ... -0.0753 -0.1786 -0.2662]\n",
       "  [-1.1848 -0.5717 -0.392  ... -0.076  -0.2826 -0.3658]]\n",
       "\n",
       " [[-0.168  -0.0692 -0.2042 ... -0.0249  0.0083  0.0364]\n",
       "  [-0.0429 -0.1653 -0.0254 ...  0.0097  0.0158 -0.0304]\n",
       "  [-0.0471 -0.0485  0.0284 ... -0.1066 -0.0428  0.0851]\n",
       "  ...\n",
       "  [ 0.0045  0.0366 -0.0132 ... -0.1122 -0.1236 -0.2155]\n",
       "  [ 0.3211  0.3949  0.3327 ... -0.2165 -0.264  -0.2554]\n",
       "  [ 0.6814  0.6478  0.554  ... -0.4091 -0.3963 -0.487 ]]], 'micromolar')>\n",
       "Coordinates: (3/6)\n",
       "  * chromo    (chromo) <U3 24B 'HbO' 'HbR'\n",
       "  * time      (time) float64 4kB 0.0 0.128 0.256 0.384 ... 59.65 59.78 59.9\n",
       "    ...        ...\n",
       "    detector  (channel) object 224B 'D1' 'D2' 'D3' 'D9' ... 'D7' 'D8' 'D16'
" ], "text/plain": [ " Size: 210kB\n", "\n", "Coordinates: (3/6)\n", " * chromo (chromo) " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "f,ax = p.subplots(1,1, figsize=(12,4))\n", "ax.plot( conc.time, conc.sel(channel=\"S1D1\", chromo=\"HbO\").pint.to(\"micromolar\"), \"r-\", label=\"HbO\")\n", "ax.plot( conc.time, conc.sel(channel=\"S1D1\", chromo=\"HbR\").pint.to(\"micromolar\"), \"b-\", label=\"HbR\")\n", "p.legend()\n", "ax.set_xlabel(\"time / s\")\n", "ax.set_ylabel(\"$\\Delta c$ / $\\mu M$\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "execution": { "iopub.execute_input": "2026-06-15T14:49:29.368679Z", "iopub.status.busy": "2026-06-15T14:49:29.368529Z", "iopub.status.idle": "2026-06-15T14:49:29.532623Z", "shell.execute_reply": "2026-06-15T14:49:29.531967Z" } }, "outputs": [ { "data": { "text/html": [ "

Methods used

[1]Luke2021cedalion.data.get_fingertapping_snirf_pathRobert Luke and David McAlpine.\n", "fNIRS Finger Tapping Data in BIDS Format.\n", "September 2021.\n", "doi:10.5281/zenodo.5529797.
[2]Tucker2022cedalion.io.snirf.read_snirfStephen Tucker, Jay Dubb, Sreekanth Kura, Alexander von Lühmann, Robert Franke, Jörn M. Horschig, Samuel Powell, Robert Oostenveld, Michael Lührs, Édouard Delaire, Zahra M. Aghajan, Hanseok Yun, Meryem A. Yücel, Qianqian Fang, Theodore J. Huppert, Blaise deB. Frederick, Luca Pollonini, David A. Boas, and Robert Luke.\n", "Introduction to the shared near infrared spectroscopy format.\n", "Neurophotonics, 10(1):013507, 2022.\n", "doi:10.1117/1.NPh.10.1.013507.
[3]Prahl1998cedalion.nirs.common.get_extinction_coefficientsScott A. Prahl.\n", "Optical absorption of hemoglobin.\n", "Oregon Medical Laser Center, online resource, 1998.\n", "URL: https://omlc.org/spectra/hemoglobin/.
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "cedalion.bib.dump_to_notebook()" ] } ], "metadata": { "kernelspec": { "display_name": "cedalion_240924", "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": 2 }