{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Xarray Data Structures - an fNIRS example\n", "\n", "This exampple 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": "2025-04-14T12:53:45.045459Z", "iopub.status.busy": "2025-04-14T12:53:45.045270Z", "iopub.status.idle": "2025-04-14T12:53:46.697565Z", "shell.execute_reply": "2025-04-14T12:53:46.697044Z" } }, "outputs": [], "source": [ "import cedalion\n", "import cedalion.nirs\n", "import cedalion.xrutils\n", "import cedalion.xrutils as xrutils\n", "from cedalion.datasets 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.datasets`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Load amplitude data from the snirf file." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2025-04-14T12:53:46.700326Z", "iopub.status.busy": "2025-04-14T12:53:46.699717Z", "iopub.status.idle": "2025-04-14T12:53:46.816118Z", "shell.execute_reply": "2025-04-14T12:53:46.815660Z" } }, "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": 3, "metadata": { "execution": { "iopub.execute_input": "2025-04-14T12:53:46.818467Z", "iopub.status.busy": "2025-04-14T12:53:46.818118Z", "iopub.status.idle": "2025-04-14T12:53:46.823280Z", "shell.execute_reply": "2025-04-14T12:53:46.822818Z" } }, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "recordings" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Amplitude data" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2025-04-14T12:53:46.851441Z", "iopub.status.busy": "2025-04-14T12:53:46.850993Z", "iopub.status.idle": "2025-04-14T12:53:46.864184Z", "shell.execute_reply": "2025-04-14T12:53:46.863659Z" } }, "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": 5, "metadata": { "execution": { "iopub.execute_input": "2025-04-14T12:53:46.866287Z", "iopub.status.busy": "2025-04-14T12:53:46.865970Z", "iopub.status.idle": "2025-04-14T12:53:46.874681Z", "shell.execute_reply": "2025-04-14T12:53:46.874241Z" } }, "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": [ "Specify differential path length factors (DPF). Obtain a matrix of tabulated extinction coefficients for the wavelengths of our dataset and calculate the inverse. Cedalion offers dedicated functions for mBLL conversion ( nirs.int2od(), nirs.od2conc(), and nirs.beer-lambert() functions from the nirs subpackage) - but we do not use them here to better showcase how Xarrays work. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2025-04-14T12:53:46.893848Z", "iopub.status.busy": "2025-04-14T12:53:46.893462Z", "iopub.status.idle": "2025-04-14T12:53:46.906835Z", "shell.execute_reply": "2025-04-14T12:53:46.906311Z" } }, "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$\");" ] } ], "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.8" } }, "nbformat": 4, "nbformat_minor": 2 }