Image Reconstruction

Notebook configuration

Decide for an example with a sparse probe or a high density probe for DOT. The notebook will load example data accordingly.

Also specify, if precomputed results of the photon propagation should be used and if the 3D visualizations should be interactive.

# choose between two datasets
DATASET = "fingertappingDOT" # high-density montage
#DATASET = "fingertapping"   # sparse montage

# choose a head model
HEAD_MODEL = "colin27"
# HEAD_MODEL = "icbm152"

# choose between the monte
FORWARD_MODEL = "MCX" # photon monte carlo
#FORWARD_MODEL = "NIRFASTER" # finite element method - NOTE, you must have NIRFASTer installed via runnning <$ bash install_nirfaster.sh CPU # or GPU> from a within your cedalion root directory.

# set this flag to False to actual compute the forward model results
PRECOMPUTED_FLUENCE = True

# set this flag to True to enable interactive 3D plots
INTERACTIVE_PLOTS = False



%load_ext autoreload
%autoreload 2

import pyvista as pv

if INTERACTIVE_PLOTS:
    pv.set_jupyter_backend('server')
else:
    pv.set_jupyter_backend('static')

import time
import os

import matplotlib.pyplot as p
import numpy as np
import xarray as xr
from IPython.display import Image
import warnings
from pint.errors import UnitStrippedWarning

import cedalion
import cedalion.dataclasses as cdc
import cedalion.datasets
import cedalion.sigproc.quality as quality
import cedalion.geometry.registration
import cedalion.geometry.segmentation
import cedalion.imagereco.forward_model as fw
import cedalion.imagereco.tissue_properties
import cedalion.io
import cedalion.plots
import cedalion.vis.plot_sensitivity_matrix
from cedalion.imagereco.solver import pseudo_inverse_stacked
from cedalion import units

xr.set_options(display_expand_data=False)

<xarray.core.options.set_options at 0x7f72857430d0>

Load a DOT finger-tapping dataset

For this demo we load an example finger-tapping recording through cedalion.datasets.get_fingertapping. The file contains a single NIRS element with one block of raw amplitude data.

if DATASET == "fingertappingDOT":
    rec = cedalion.datasets.get_fingertappingDOT()
elif DATASET == "fingertapping":
    rec = cedalion.datasets.get_fingertapping()
else:
    raise ValueError("unknown dataset")

The location of the probes is obtained from the snirf metadata (i.e. /nirs0/probe/)

Note that units (‘m’) are adopted and the coordinate system is named ‘digitized’.

geo3d_meas = rec.geo3d
display(geo3d_meas)

<xarray.DataArray (label: 346, digitized: 3)> Size: 8kB
[mm] -77.82 15.68 23.17 -61.91 21.23 56.49 ... 14.23 -38.28 81.95 -0.678 -37.03
Coordinates:
    type     (label) object 3kB PointType.SOURCE ... PointType.LANDMARK
  * label    (label) <U6 8kB 'S1' 'S2' 'S3' 'S4' ... 'FFT10h' 'FT10h' 'FTT10h'
Dimensions without coordinates: digitized

xarray.DataArray

• label: 346
• digitized: 3

• [mm] -77.82 15.68 23.17 -61.91 21.23 ... -38.28 81.95 -0.678 -37.03

  Magnitude	[[-77.817871 15.680613999999998 23.172269999999997] [-61.906841 21.227732000000003 56.492802] [-85.37146 -16.079958 8.900885] ... [77.521 28.883 -39.113] [80.59 14.229 -38.278] [81.95 -0.678 -37.027]]
  Units	millimeter

• Coordinates: (2)
  • type
    (label)
    object
    PointType.SOURCE ... PointType.L...

    array([<PointType.SOURCE: 1>, <PointType.SOURCE: 1>, <PointType.SOURCE: 1>,
           ..., <PointType.LANDMARK: 3>, <PointType.LANDMARK: 3>,
           <PointType.LANDMARK: 3>], dtype=object)

  • label
    (label)
    <U6
    'S1' 'S2' 'S3' ... 'FT10h' 'FTT10h'

    array(['S1', 'S2', 'S3', ..., 'FFT10h', 'FT10h', 'FTT10h'], dtype='<U6')

• Indexes: (1)
  • label
    PandasIndex

    PandasIndex(Index(['S1', 'S2', 'S3', 'S4', 'S5', 'S6', 'S7', 'S8', 'S9', 'S10',
           ...
           'POO10h', 'PO10h', 'PPO10h', 'P10h', 'TPP10h', 'AFF10h', 'F10h',
           'FFT10h', 'FT10h', 'FTT10h'],
          dtype='object', name='label', length=346))

• Attributes: (0)

cedalion.plots.plot_montage3D(rec["amp"], geo3d_meas)

The measurement list is a pandas.DataFrame that describes which source detector pairs form channels.

meas_list = rec._measurement_lists["amp"]
display(meas_list.head(5))

	sourceIndex	detectorIndex	wavelengthIndex	wavelengthActual	wavelengthEmissionActual	dataType	dataUnit	dataTypeLabel	dataTypeIndex	sourcePower	detectorGain	moduleIndex	sourceModuleIndex	detectorModuleIndex	channel	source	detector	wavelength	chromo
0	1	1	1	None	None	1	None	raw-DC	1	None	None	None	None	None	S1D1	S1	D1	760.0	None
1	1	2	1	None	None	1	None	raw-DC	1	None	None	None	None	None	S1D2	S1	D2	760.0	None
2	1	4	1	None	None	1	None	raw-DC	1	None	None	None	None	None	S1D4	S1	D4	760.0	None
3	1	5	1	None	None	1	None	raw-DC	1	None	None	None	None	None	S1D5	S1	D5	760.0	None
4	1	6	1	None	None	1	None	raw-DC	1	None	None	None	None	None	S1D6	S1	D6	760.0	None

Event/stimulus information is also stored in a pandas.DataFrame. Here events are given more descriptive names:

if DATASET == "fingertappingDOT":
   rec.stim.cd.rename_events( {
        "1": "Control",
        "2": "FTapping/Left",
        "3": "FTapping/Right",
        "4": "BallSqueezing/Left",
        "5": "BallSqueezing/Right"
    } )
elif DATASET == "fingertapping":
    rec.stim.cd.rename_events( {
        "1.0": "Control",
        "2.0": "FTapping/Left",
        "3.0": "FTapping/Right"
    } )

display(rec.stim.groupby("trial_type")[["onset"]].count())

	onset
trial_type
BallSqueezing/Left	17
BallSqueezing/Right	16
Control	65
FTapping/Left	16
FTapping/Right	16

Perform pruning, conversion to OD and bandpass filtering

(for this demo select 20 seconds after a trial starts at t=117s and transform raw amplitudes to optical density)

Perform SNR quality check and pruning and then transform CW raw amplitudes to optical density

## Prune with SNR threshold
snr_thresh = 10 # dB
snr, rec.masks["snr_mask"] = quality.snr(rec["amp"], snr_thresh)
# prune channels using the masks and the operator "all", which will keep only channels that pass all three metrics
rec["amp_pruned"], drop_list = quality.prune_ch(rec["amp"], [rec.masks["snr_mask"]], "all")

print(drop_list)

# Convert to OD
rec["od"] = cedalion.nirs.int2od(rec["amp_pruned"])

# bandpass filter the data
rec["od_freqfiltered"] = rec["od"].cd.freq_filter(fmin=0.01, fmax=0.5, butter_order=4)

[]

Calculate block averages in optical density

# segment data into epochs
epochs = rec["od_freqfiltered"].cd.to_epochs(
        rec.stim, # stimulus dataframe
        ["FTapping/Left", "FTapping/Right"],  # select events, discard the others
        before=5 * cedalion.units.s, # seconds before stimulus
        after=30 * cedalion.units.s, # seconds after stimulus
)

# calculate baseline
baseline = epochs.sel(reltime=(epochs.reltime < 0)).mean("reltime")
# subtract baseline
epochs_blcorrected = epochs - baseline

# group trials by trial_type. For each group individually average the epoch dimension
blockaverage = epochs_blcorrected.groupby("trial_type").mean("epoch")

# Plot block averages. Please ignore errors if the plot is too small in the HD case

noPlts2 = int(np.ceil(np.sqrt(len(blockaverage.channel))))
f,ax = p.subplots(noPlts2,noPlts2, figsize=(12,10))
ax = ax.flatten()
for i_ch, ch in enumerate(blockaverage.channel):
    for ls, trial_type in zip(["-", "--"], blockaverage.trial_type):
        ax[i_ch].plot(blockaverage.reltime, blockaverage.sel(wavelength=760, trial_type=trial_type, channel=ch), "r", lw=2, ls=ls)
        ax[i_ch].plot(blockaverage.reltime, blockaverage.sel(wavelength=850, trial_type=trial_type, channel=ch), "b", lw=2, ls=ls)

    ax[i_ch].grid(1)
    ax[i_ch].set_title(ch.values)
    ax[i_ch].set_ylim(-.02, .02)
    ax[i_ch].set_axis_off()
    ax[i_ch].axhline(0, c="k")
    ax[i_ch].axvline(0, c="k")

p.suptitle("760nm: r | 850nm: b | left: - | right: --")
p.tight_layout()

Load segmented MRI scan

For this example use a segmentation of the Colin27 average brain.

if HEAD_MODEL == "colin27":
    SEG_DATADIR, mask_files, landmarks_file = cedalion.datasets.get_colin27_segmentation()
elif HEAD_MODEL == "icbm152":
    SEG_DATADIR, mask_files, landmarks_file = cedalion.datasets.get_icbm152_segmentation()
else:
    raise ValueError("unknown head model")

The segmentation masks are in individual niftii files. The dict mask_files maps mask filenames relative to SEG_DATADIR to short labels. These labels describe the tissue type of the mask.

In principle the user is free to choose these labels. However, they are later used to lookup the tissue’s optical properties. So they must be map to one of the tabulated tissue types (c.f. cedalion.imagereco.tissue_properties.TISSUE_LABELS).

The variable landmarks_file holds the path to a file containing landmark positions in scanner space (RAS). This file can be created with Slicer3D.

display(SEG_DATADIR)
display(mask_files)
display(landmarks_file)

'/home/runner/.cache/cedalion/colin27_segmentation.zip.unzip/colin27_segmentation'

{'csf': 'mask_csf.nii',
 'gm': 'mask_gray.nii',
 'scalp': 'mask_skin.nii',
 'skull': 'mask_bone.nii',
 'wm': 'mask_white.nii'}

'landmarks.mrk.json'

Coordinate systems

Up to now we have geometrical data from three different coordinate reference systems (CRS):

• The optode positions are in one space CRS='digitized' and the coordinates are in meter. In our example the origin is at the head center and y-axis pointing in the superior direction. Other digitization tools can use other units or coordinate systems.

• The segmentation masks are in voxel space (CRS='ijk') in which the voxel edges are aligned with the coordinate axes. Each voxel has unit edge length, i.e. coordinates are dimensionless. Axis-aligned grids are computationally efficient, which is why the photon simulation code (MCX) uses this coordinate system.

• The voxel space (CRS='ijk') is related to scanner space (CRS='ras' or CRS='aligned') in which coordinates have physical units and coordinate axes point to the (r)ight, (a)nterior and s(uperior) directions. The relation between both spaces is given through an affine transformation (e.g. t_ijk2ras). When loading the segmentation masks in Slicer3D this transformation is automatically applied. Hence, the picked landmark coordinates are exported in RAS space.

  The niftii file provides a string label for the scanner space. In this example the RAS space is called ‘aligned’ because the masks are aligned to another MRI scan.

To avoid confusion between these different coordinate systems, cedalion tries to be explicit about which CRS a given point cloud or surface is in.

The TwoSurfaceHeadModel

The photon propagation considers the complete MRI scan, in which each voxel is attributed to one tissue type with its respective optical properties. However, the image reconstruction does not intend to reconstruct absorption changes in each voxel. The inverse problem is simplified, by considering only two surfaces (scalp and brain) and reconstruct only absorption changes in voxels close to these surfaces.

The class cedalion.imagereco.forward_model.TwoSurfaceHeadModel groups together the segmentation mask, landmark positions and affine transformations as well as the scalp and brain surfaces. The brain surface is calculated by grouping together white and gray matter masks. The scalp surface encloses the whole head.

head = fw.TwoSurfaceHeadModel.from_surfaces(
    segmentation_dir=SEG_DATADIR,
    mask_files = mask_files,
    brain_surface_file= os.path.join(SEG_DATADIR, "mask_brain.obj"),
    # scalp_surface_file= os.path.join(SEG_DATADIR, "mask_scalp.obj"), # FIXME!
    landmarks_ras_file=landmarks_file,
    brain_face_count=None,
    scalp_face_count=None
)

head.segmentation_masks

<xarray.DataArray (segmentation_type: 5, i: 181, j: 217, k: 181)> Size: 36MB
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
Coordinates:
  * segmentation_type  (segmentation_type) <U5 100B 'csf' 'gm' ... 'skull' 'wm'
Dimensions without coordinates: i, j, k

xarray.DataArray

• segmentation_type: 5
• i: 181
• j: 217
• k: 181

• 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

  array([[[[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, 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, 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],
           [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, 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]]]], dtype=uint8)

• Coordinates: (1)
  • segmentation_type
    (segmentation_type)
    <U5
    'csf' 'gm' 'scalp' 'skull' 'wm'

    array(['csf', 'gm', 'scalp', 'skull', 'wm'], dtype='<U5')

• Indexes: (1)
  • segmentation_type
    PandasIndex

    PandasIndex(Index(['csf', 'gm', 'scalp', 'skull', 'wm'], dtype='object', name='segmentation_type'))

• Attributes: (0)

head.landmarks

<xarray.DataArray (label: 4, ijk: 3)> Size: 96B
[] 89.95 205.8 35.86 91.76 25.05 17.14 18.05 109.9 17.34 165.9 113.0 17.92
Coordinates:
  * label    (label) <U3 48B 'Nz' 'Iz' 'LPA' 'RPA'
    type     (label) object 32B PointType.LANDMARK ... PointType.LANDMARK
Dimensions without coordinates: ijk

xarray.DataArray

• label: 4
• ijk: 3

• [] 89.95 205.8 35.86 91.76 25.05 ... 109.9 17.34 165.9 113.0 17.92

  Magnitude	[[89.95007602503179 205.78946685798357 35.858879089397426] [91.75606727599917 25.050003051755283 17.14201354980326] [18.05000305178696 109.88059425358084 17.3430213928456] [165.94999694827212 112.96190643314736 17.9221572876216]]
  Units	dimensionless

• Coordinates: (2)
  • label
    (label)
    <U3
    'Nz' 'Iz' 'LPA' 'RPA'

    array(['Nz', 'Iz', 'LPA', 'RPA'], dtype='<U3')

  • type
    (label)
    object
    PointType.LANDMARK ... PointType...

    array([<PointType.LANDMARK: 3>, <PointType.LANDMARK: 3>,
           <PointType.LANDMARK: 3>, <PointType.LANDMARK: 3>], dtype=object)

• Indexes: (1)
  • label
    PandasIndex

    PandasIndex(Index(['Nz', 'Iz', 'LPA', 'RPA'], dtype='object', name='label'))

• Attributes: (0)

head.brain

TrimeshSurface(mesh=<trimesh.Trimesh(vertices.shape=(15002, 3), faces.shape=(29988, 3))>, crs='ijk', units=<Unit('dimensionless')>)

head.scalp

TrimeshSurface(mesh=<trimesh.Trimesh(vertices.shape=(152572, 3), faces.shape=(304080, 3))>, crs='ijk', units=<Unit('dimensionless')>)

TwoSurfaceHeadModel.from_segmentation converts everything into voxel space (CRS='ijk')

head.crs

'ijk'

The transformation matrix to translate from voxel to scanner space:

head.t_ijk2ras

<xarray.DataArray (aligned: 4, ijk: 4)> Size: 128B
[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
Dimensions without coordinates: aligned, ijk

xarray.DataArray

• aligned: 4
• ijk: 4

• [mm] 1.0 0.0 0.0 -90.0 0.0 1.0 0.0 ... 0.0 1.0 -72.0 0.0 0.0 0.0 1.0

  Magnitude	[[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]]
  Units	millimeter

• Coordinates: (0)
• Indexes: (0)
• Attributes: (0)

Changing between coordinate systems:

head_ras = head.apply_transform(head.t_ijk2ras)
display(head_ras.crs)
display(head_ras.brain)

'aligned'

TrimeshSurface(mesh=<trimesh.Trimesh(vertices.shape=(15002, 3), faces.shape=(29988, 3))>, crs='aligned', units=<Unit('millimeter')>)

Optode Registration

The optode coordinates from the recording must be aligned with the scalp surface. Currently, cedaĺion offers a simple registration method, which finds an affine transformation (scaling, rotating, translating) that matches the landmark positions of the head model and their digitized counter parts. Afterwards, optodes are snapped to the nearest vertex on the scalp.

geo3d_snapped_ijk = head.align_and_snap_to_scalp(geo3d_meas)
display(geo3d_snapped_ijk)

<xarray.DataArray (label: 346, ijk: 3)> Size: 8kB
[] 9.901 129.0 88.0 23.58 135.0 127.0 ... 173.1 140.0 35.0 172.0 126.0 34.98
Coordinates:
    type     (label) object 3kB PointType.SOURCE ... PointType.LANDMARK
  * label    (label) <U6 8kB 'S1' 'S2' 'S3' 'S4' ... 'FFT10h' 'FT10h' 'FTT10h'
Dimensions without coordinates: ijk

xarray.DataArray

• label: 346
• ijk: 3

• [] 9.901 129.0 88.0 23.58 135.0 127.0 ... 140.0 35.0 172.0 126.0 34.98

  Magnitude	[[9.90067996231771 129.00337747726593 88.00475372480513] [23.582338748341726 135.01064219022152 127.0189188518341] [7.93222233511454 104.96269831836995 69.9588207909711] ... [169.3304445080471 157.024478761799 37.25591035636486] [173.07561584082652 139.9987014715881 34.99651508735311] [172.03364292496815 125.98707966194677 34.9768180480424]]
  Units	dimensionless

• Coordinates: (2)
  • type
    (label)
    object
    PointType.SOURCE ... PointType.L...

    array([<PointType.SOURCE: 1>, <PointType.SOURCE: 1>, <PointType.SOURCE: 1>,
           ..., <PointType.LANDMARK: 3>, <PointType.LANDMARK: 3>,
           <PointType.LANDMARK: 3>], dtype=object)

  • label
    (label)
    <U6
    'S1' 'S2' 'S3' ... 'FT10h' 'FTT10h'

    array(['S1', 'S2', 'S3', ..., 'FFT10h', 'FT10h', 'FTT10h'], dtype='<U6')

• Indexes: (1)
  • label
    PandasIndex

    PandasIndex(Index(['S1', 'S2', 'S3', 'S4', 'S5', 'S6', 'S7', 'S8', 'S9', 'S10',
           ...
           'POO10h', 'PO10h', 'PPO10h', 'P10h', 'TPP10h', 'AFF10h', 'F10h',
           'FFT10h', 'FT10h', 'FTT10h'],
          dtype='object', name='label', length=346))

• Attributes: (0)

plt = pv.Plotter()
cedalion.plots.plot_surface(plt, head.brain, color="w")
cedalion.plots.plot_surface(plt, head.scalp, opacity=.1)
cedalion.plots.plot_labeled_points(plt, geo3d_snapped_ijk)
plt.show()

Simulate light propagation in tissue

cedalion.imagereco.forward_model.ForwardModel is a wrapper around pmcx. Using the data in the head model it prepares the inputs for either pmcx or NIRFASTer and offers functionality to calculate the sensitivty matrix.

fwm = cedalion.imagereco.forward_model.ForwardModel(head, geo3d_snapped_ijk, meas_list)

Run the simulation

The compute_fluence_mcx and compute_fluence_nirfaster methods simulate a light source at each optode position and calculate the fluence in each voxel. By setting RUN_PACKAGE, you can choose between the pmcx or NIRFASTer package to perform this simulation. PLEASE NOTE: if you USE_CACHED data (download the example data) be aware that the file is quite big (~2GB).

if PRECOMPUTED_FLUENCE:
    if FORWARD_MODEL == "MCX":
        fluence_all, fluence_at_optodes = cedalion.datasets.get_precomputed_fluence(DATASET, HEAD_MODEL)
    elif FORWARD_MODEL == "NIRFASTER":
        raise NotImplementedError(
            "Currently there are no precomputed NIRFASTER results available"
        )
else:
    if FORWARD_MODEL == "MCX":
        fluence_all, fluence_at_optodes = fwm.compute_fluence_mcx()
    elif FORWARD_MODEL == "NIRFASTER":
        fluence_all, fluence_at_optodes = fwm.compute_fluence_nirfaster()

The photon simulation yields the fluence in each voxel for each wavelength.

fluence_all

<xarray.DataArray '/fluence_all' (label: 46, wavelength: 2, i: 181, j: 217,
                                  k: 181)> Size: 5GB
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
Coordinates:
  * label       (label) object 368B 'S1' 'S2' 'S3' 'S4' ... 'D30' 'D31' 'D32'
    type        (label) object 368B PointType.SOURCE ... PointType.DETECTOR
  * wavelength  (wavelength) float64 16B 760.0 850.0
Dimensions without coordinates: i, j, k

xarray.DataArray

'/fluence_all'

• label: 46
• wavelength: 2
• i: 181
• j: 217
• k: 181

• 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

  array([[[[[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., 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., 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.],
            [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., 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.]]]]])

• Coordinates: (3)
  • label
    (label)
    object
    'S1' 'S2' 'S3' ... 'D31' 'D32'

    array(['S1', 'S2', 'S3', 'S4', 'S5', 'S6', 'S7', 'S8', 'S9', 'S10', 'S11',
           'S12', 'S13', 'S14', 'D1', 'D2', 'D3', 'D4', 'D5', 'D6', 'D7', 'D8',
           'D9', 'D10', 'D11', 'D12', 'D13', 'D14', 'D15', 'D16', 'D17', 'D18',
           'D19', 'D20', 'D21', 'D22', 'D23', 'D24', 'D25', 'D26', 'D27', 'D28',
           'D29', 'D30', 'D31', 'D32'], dtype=object)

  • type
    (label)
    object
    PointType.SOURCE ... PointType.D...

    array([<PointType.SOURCE: 1>, <PointType.SOURCE: 1>,
           <PointType.SOURCE: 1>, <PointType.SOURCE: 1>,
           <PointType.SOURCE: 1>, <PointType.SOURCE: 1>,
           <PointType.SOURCE: 1>, <PointType.SOURCE: 1>,
           <PointType.SOURCE: 1>, <PointType.SOURCE: 1>,
           <PointType.SOURCE: 1>, <PointType.SOURCE: 1>,
           <PointType.SOURCE: 1>, <PointType.SOURCE: 1>,
           <PointType.DETECTOR: 2>, <PointType.DETECTOR: 2>,
           <PointType.DETECTOR: 2>, <PointType.DETECTOR: 2>,
           <PointType.DETECTOR: 2>, <PointType.DETECTOR: 2>,
           <PointType.DETECTOR: 2>, <PointType.DETECTOR: 2>,
           <PointType.DETECTOR: 2>, <PointType.DETECTOR: 2>,
           <PointType.DETECTOR: 2>, <PointType.DETECTOR: 2>,
           <PointType.DETECTOR: 2>, <PointType.DETECTOR: 2>,
           <PointType.DETECTOR: 2>, <PointType.DETECTOR: 2>,
           <PointType.DETECTOR: 2>, <PointType.DETECTOR: 2>,
           <PointType.DETECTOR: 2>, <PointType.DETECTOR: 2>,
           <PointType.DETECTOR: 2>, <PointType.DETECTOR: 2>,
           <PointType.DETECTOR: 2>, <PointType.DETECTOR: 2>,
           <PointType.DETECTOR: 2>, <PointType.DETECTOR: 2>,
           <PointType.DETECTOR: 2>, <PointType.DETECTOR: 2>,
           <PointType.DETECTOR: 2>, <PointType.DETECTOR: 2>,
           <PointType.DETECTOR: 2>, <PointType.DETECTOR: 2>], dtype=object)

  • wavelength
    (wavelength)
    float64
    760.0 850.0

    array([760., 850.])

• Indexes: (2)
  • label
    PandasIndex

    PandasIndex(Index(['S1', 'S2', 'S3', 'S4', 'S5', 'S6', 'S7', 'S8', 'S9', 'S10', 'S11',
           'S12', 'S13', 'S14', 'D1', 'D2', 'D3', 'D4', 'D5', 'D6', 'D7', 'D8',
           'D9', 'D10', 'D11', 'D12', 'D13', 'D14', 'D15', 'D16', 'D17', 'D18',
           'D19', 'D20', 'D21', 'D22', 'D23', 'D24', 'D25', 'D26', 'D27', 'D28',
           'D29', 'D30', 'D31', 'D32'],
          dtype='object', name='label'))

  • wavelength
    PandasIndex

    PandasIndex(Index([760.0, 850.0], dtype='float64', name='wavelength'))

• Attributes: (0)

Also, for a each combination of two optodes, the fluence in the voxels at the optode positions is calculated.

fluence_at_optodes

<xarray.DataArray '/fluence_at_optodes' (optode1: 46, optode2: 46, wavelength: 2)> Size: 34kB
0.05344 0.05344 1.021e-06 1.021e-06 ... 2.2e-07 2.2e-07 0.02938 0.02938
Coordinates:
  * optode1     (optode1) object 368B 'S1' 'S2' 'S3' 'S4' ... 'D30' 'D31' 'D32'
  * optode2     (optode2) object 368B 'S1' 'S2' 'S3' 'S4' ... 'D30' 'D31' 'D32'
  * wavelength  (wavelength) float64 16B 760.0 850.0

xarray.DataArray

'/fluence_at_optodes'

• optode1: 46
• optode2: 46
• wavelength: 2

• 0.05344 0.05344 1.021e-06 1.021e-06 ... 2.2e-07 0.02938 0.02938

  array([[[5.34380488e-02, 5.34380488e-02],
          [1.02059960e-06, 1.02059960e-06],
          [1.22795200e-06, 1.22795200e-06],
          ...,
          [7.01381910e-18, 7.01381910e-18],
          [3.00826673e-15, 3.00826673e-15],
          [5.10897914e-17, 5.10897914e-17]],
  
         [[1.53354347e-06, 1.53354347e-06],
          [5.28522953e-02, 5.28522953e-02],
          [1.33221207e-07, 1.33221207e-07],
          ...,
          [4.80825366e-17, 4.80825366e-17],
          [3.52020371e-15, 3.52020371e-15],
          [1.73195388e-13, 1.73195388e-13]],
  
         [[9.72019507e-07, 9.72019507e-07],
          [1.01156232e-07, 1.01156232e-07],
          [6.00267164e-02, 6.00267164e-02],
          ...,
  ...
          ...,
          [4.31555510e-02, 4.31555510e-02],
          [9.23773712e-07, 9.23773712e-07],
          [1.37099896e-08, 1.37099896e-08]],
  
         [[1.83557336e-19, 1.83557336e-19],
          [1.14419148e-15, 1.14419148e-15],
          [2.37683842e-13, 2.37683842e-13],
          ...,
          [1.09036716e-06, 1.09036716e-06],
          [3.56829092e-02, 3.56829092e-02],
          [1.43890887e-07, 1.43890887e-07]],
  
         [[4.04321902e-15, 4.04321902e-15],
          [6.47715768e-13, 6.47715768e-13],
          [2.08797941e-14, 2.08797941e-14],
          ...,
          [1.80744379e-08, 1.80744379e-08],
          [2.19959546e-07, 2.19959546e-07],
          [2.93827020e-02, 2.93827020e-02]]])

• Coordinates: (3)
  • optode1
    (optode1)
    object
    'S1' 'S2' 'S3' ... 'D31' 'D32'

    array(['S1', 'S2', 'S3', 'S4', 'S5', 'S6', 'S7', 'S8', 'S9', 'S10', 'S11',
           'S12', 'S13', 'S14', 'D1', 'D2', 'D3', 'D4', 'D5', 'D6', 'D7', 'D8',
           'D9', 'D10', 'D11', 'D12', 'D13', 'D14', 'D15', 'D16', 'D17', 'D18',
           'D19', 'D20', 'D21', 'D22', 'D23', 'D24', 'D25', 'D26', 'D27', 'D28',
           'D29', 'D30', 'D31', 'D32'], dtype=object)

  • optode2
    (optode2)
    object
    'S1' 'S2' 'S3' ... 'D31' 'D32'

    array(['S1', 'S2', 'S3', 'S4', 'S5', 'S6', 'S7', 'S8', 'S9', 'S10', 'S11',
           'S12', 'S13', 'S14', 'D1', 'D2', 'D3', 'D4', 'D5', 'D6', 'D7', 'D8',
           'D9', 'D10', 'D11', 'D12', 'D13', 'D14', 'D15', 'D16', 'D17', 'D18',
           'D19', 'D20', 'D21', 'D22', 'D23', 'D24', 'D25', 'D26', 'D27', 'D28',
           'D29', 'D30', 'D31', 'D32'], dtype=object)

  • wavelength
    (wavelength)
    float64
    760.0 850.0

    array([760., 850.])

• Indexes: (3)
  • optode1
    PandasIndex

    PandasIndex(Index(['S1', 'S2', 'S3', 'S4', 'S5', 'S6', 'S7', 'S8', 'S9', 'S10', 'S11',
           'S12', 'S13', 'S14', 'D1', 'D2', 'D3', 'D4', 'D5', 'D6', 'D7', 'D8',
           'D9', 'D10', 'D11', 'D12', 'D13', 'D14', 'D15', 'D16', 'D17', 'D18',
           'D19', 'D20', 'D21', 'D22', 'D23', 'D24', 'D25', 'D26', 'D27', 'D28',
           'D29', 'D30', 'D31', 'D32'],
          dtype='object', name='optode1'))

  • optode2
    PandasIndex

    PandasIndex(Index(['S1', 'S2', 'S3', 'S4', 'S5', 'S6', 'S7', 'S8', 'S9', 'S10', 'S11',
           'S12', 'S13', 'S14', 'D1', 'D2', 'D3', 'D4', 'D5', 'D6', 'D7', 'D8',
           'D9', 'D10', 'D11', 'D12', 'D13', 'D14', 'D15', 'D16', 'D17', 'D18',
           'D19', 'D20', 'D21', 'D22', 'D23', 'D24', 'D25', 'D26', 'D27', 'D28',
           'D29', 'D30', 'D31', 'D32'],
          dtype='object', name='optode2'))

  • wavelength
    PandasIndex

    PandasIndex(Index([760.0, 850.0], dtype='float64', name='wavelength'))

• Attributes: (0)

Plot fluence

To illustrate the tissue probed by light travelling from a source to the detector two fluence profiles need to be multiplied.

# for plotting use a geo3d without the landmarks
geo3d_plot = geo3d_snapped_ijk[geo3d_snapped_ijk.type != cdc.PointType.LANDMARK]

time.sleep(1)

plt = pv.Plotter()

if DATASET == "fingertappingDOT":
    src, det, wl = "S4", "D2", 760
elif DATASET == "fingertapping":
    src, det, wl = "S2", "D3", 760
else:
    raise ValueError("unknown dataset")

f = fluence_all.loc[src, wl].values * fluence_all.loc[det, wl].values
f[f <= 0] = f[f > 0].min()
f = np.log10(f)
vf = pv.wrap(f)

plt.add_volume(
    vf,
    log_scale=False,
    cmap="plasma_r",
    clim=(-10, 0),
)
cedalion.plots.plot_surface(plt, head.brain, color="w")
cedalion.plots.plot_labeled_points(plt, geo3d_plot, show_labels=False)

cog = head.brain.vertices.mean("label").values
plt.camera.position = cog + [-300, 30, 150]
plt.camera.focal_point = cog
plt.camera.up = [0, 0, 1]

plt.show()

/home/runner/miniconda3/envs/cedalion/lib/python3.11/site-packages/xarray/core/variable.py:338: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
  data = np.asarray(data)

Calculate the sensitivity matrices

The sensitivity matrix describes the effect of an absorption change at a given surface vertex in the OD recording in a given channel and at given wavelength. The coordinate is_brain holds a mask to distinguish brain and scalp voxels.

Adot = fwm.compute_sensitivity(fluence_all, fluence_at_optodes)
Adot

<xarray.DataArray (channel: 100, vertex: 167574, wavelength: 2)> Size: 268MB
7.021e-13 7.021e-13 1.155e-13 1.155e-13 2.353e-14 ... 0.0 0.0 0.0 0.0 0.0
Coordinates:
  * channel     (channel) <U6 2kB 'S1D1' 'S1D2' 'S1D4' ... 'S14D31' 'S14D32'
  * wavelength  (wavelength) float64 16B 760.0 850.0
    is_brain    (vertex) bool 168kB True True True True ... False False False
Dimensions without coordinates: vertex
Attributes:
    units:    mm

xarray.DataArray

• channel: 100
• vertex: 167574
• wavelength: 2

• 7.021e-13 7.021e-13 1.155e-13 1.155e-13 2.353e-14 ... 0.0 0.0 0.0 0.0

  array([[[7.02129923e-13, 7.02129923e-13],
          [1.15507218e-13, 1.15507218e-13],
          [2.35300605e-14, 2.35300605e-14],
          ...,
          [0.00000000e+00, 0.00000000e+00],
          [0.00000000e+00, 0.00000000e+00],
          [0.00000000e+00, 0.00000000e+00]],
  
         [[1.07253446e-13, 1.07253446e-13],
          [3.48839915e-15, 3.48839915e-15],
          [1.24968454e-14, 1.24968454e-14],
          ...,
          [0.00000000e+00, 0.00000000e+00],
          [0.00000000e+00, 0.00000000e+00],
          [0.00000000e+00, 0.00000000e+00]],
  
         [[2.76322775e-13, 2.76322775e-13],
          [5.36885523e-14, 5.36885523e-14],
          [2.08213877e-14, 2.08213877e-14],
          ...,
  ...
          ...,
          [0.00000000e+00, 0.00000000e+00],
          [0.00000000e+00, 0.00000000e+00],
          [0.00000000e+00, 0.00000000e+00]],
  
         [[7.96864700e-12, 7.96864700e-12],
          [1.33651251e-11, 1.33651251e-11],
          [5.51314787e-12, 5.51314787e-12],
          ...,
          [2.88155947e-20, 2.88155947e-20],
          [0.00000000e+00, 0.00000000e+00],
          [0.00000000e+00, 0.00000000e+00]],
  
         [[5.81238573e-11, 5.81238573e-11],
          [5.36105544e-11, 5.36105544e-11],
          [2.43515315e-11, 2.43515315e-11],
          ...,
          [0.00000000e+00, 0.00000000e+00],
          [0.00000000e+00, 0.00000000e+00],
          [0.00000000e+00, 0.00000000e+00]]])

• Coordinates: (3)
  • channel
    (channel)
    <U6
    'S1D1' 'S1D2' ... 'S14D31' 'S14D32'

    array(['S1D1', 'S1D2', 'S1D4', 'S1D5', 'S1D6', 'S1D8', 'S2D2', 'S2D3', 'S2D5',
           'S2D6', 'S2D7', 'S2D9', 'S3D1', 'S3D4', 'S3D5', 'S3D8', 'S3D10',
           'S3D11', 'S3D14', 'S4D2', 'S4D4', 'S4D5', 'S4D6', 'S4D7', 'S4D8',
           'S4D9', 'S4D10', 'S4D11', 'S4D12', 'S4D13', 'S4D15', 'S5D3', 'S5D6',
           'S5D7', 'S5D9', 'S5D12', 'S5D13', 'S5D16', 'S6D8', 'S6D10', 'S6D11',
           'S6D12', 'S6D14', 'S6D15', 'S7D9', 'S7D11', 'S7D12', 'S7D13', 'S7D15',
           'S7D16', 'S8D17', 'S8D18', 'S8D20', 'S8D21', 'S8D22', 'S8D24', 'S9D18',
           'S9D19', 'S9D21', 'S9D22', 'S9D23', 'S9D25', 'S10D17', 'S10D20',
           'S10D21', 'S10D24', 'S10D26', 'S10D27', 'S10D30', 'S11D18', 'S11D20',
           'S11D21', 'S11D22', 'S11D23', 'S11D24', 'S11D25', 'S11D26', 'S11D27',
           'S11D28', 'S11D29', 'S11D31', 'S12D19', 'S12D22', 'S12D23', 'S12D25',
           'S12D28', 'S12D29', 'S12D32', 'S13D24', 'S13D26', 'S13D27', 'S13D28',
           'S13D30', 'S13D31', 'S14D25', 'S14D27', 'S14D28', 'S14D29', 'S14D31',
           'S14D32'], dtype='<U6')

  • wavelength
    (wavelength)
    float64
    760.0 850.0

    array([760., 850.])

  • is_brain
    (vertex)
    bool
    True True True ... False False

    array([ True,  True,  True, ..., False, False, False])

• Indexes: (2)
  • channel
    PandasIndex

    PandasIndex(Index(['S1D1', 'S1D2', 'S1D4', 'S1D5', 'S1D6', 'S1D8', 'S2D2', 'S2D3', 'S2D5',
           'S2D6', 'S2D7', 'S2D9', 'S3D1', 'S3D4', 'S3D5', 'S3D8', 'S3D10',
           'S3D11', 'S3D14', 'S4D2', 'S4D4', 'S4D5', 'S4D6', 'S4D7', 'S4D8',
           'S4D9', 'S4D10', 'S4D11', 'S4D12', 'S4D13', 'S4D15', 'S5D3', 'S5D6',
           'S5D7', 'S5D9', 'S5D12', 'S5D13', 'S5D16', 'S6D8', 'S6D10', 'S6D11',
           'S6D12', 'S6D14', 'S6D15', 'S7D9', 'S7D11', 'S7D12', 'S7D13', 'S7D15',
           'S7D16', 'S8D17', 'S8D18', 'S8D20', 'S8D21', 'S8D22', 'S8D24', 'S9D18',
           'S9D19', 'S9D21', 'S9D22', 'S9D23', 'S9D25', 'S10D17', 'S10D20',
           'S10D21', 'S10D24', 'S10D26', 'S10D27', 'S10D30', 'S11D18', 'S11D20',
           'S11D21', 'S11D22', 'S11D23', 'S11D24', 'S11D25', 'S11D26', 'S11D27',
           'S11D28', 'S11D29', 'S11D31', 'S12D19', 'S12D22', 'S12D23', 'S12D25',
           'S12D28', 'S12D29', 'S12D32', 'S13D24', 'S13D26', 'S13D27', 'S13D28',
           'S13D30', 'S13D31', 'S14D25', 'S14D27', 'S14D28', 'S14D29', 'S14D31',
           'S14D32'],
          dtype='object', name='channel'))

  • wavelength
    PandasIndex

    PandasIndex(Index([760.0, 850.0], dtype='float64', name='wavelength'))

• Attributes: (1)

  units :

  mm

Plot Sensitivity Matrix

plotter = cedalion.vis.plot_sensitivity_matrix.Main(
    sensitivity=Adot,
    brain_surface=head.brain,
    head_surface=head.scalp,
    labeled_points=geo3d_plot,
)
plotter.plot(high_th=0, low_th=-3)
plotter.plt.show()

The sensitivity Adot has shape (nchannel, nvertex, nwavelenghts). To solve the inverse problem we need a matrix that relates OD in channel space to absorption in image space. Hence, the sensitivity must include the extinction coefficients to translate between OD and concentrations. Furthermore, channels at different wavelengths must be stacked as well vertice and chromophores into new dimensions (flat_channel, flat_vertex):

\[\begin{split} \left( \begin{matrix} OD_{c_1, \lambda_1} \\ \vdots \\ OD_{c_N,\lambda_1} \\ OD_{c_1,\lambda_2} \\ \vdots \\ OD_{c_N,\lambda_2} \end{matrix}\right) = A \cdot \left( \begin{matrix} \Delta c_{v_1, HbO} \\ \vdots \\ \Delta c_{v_N, HbO} \\ \Delta c_{v_1, HbR} \\ \vdots \\ \Delta c_{v_N, HbR} \end{matrix}\right)\end{split}\]

Adot_stacked = fwm.compute_stacked_sensitivity(Adot)
Adot_stacked

<xarray.DataArray (flat_channel: 200, flat_vertex: 335148)> Size: 536MB
9.474e-11 1.559e-11 3.175e-12 2.06e-11 5.291e-10 ... 7.612e-17 0.0 0.0 0.0 0.0
Coordinates:
    is_brain    (flat_vertex) bool 335kB True True True ... False False False
    chromo      (flat_vertex) <U3 4MB 'HbO' 'HbO' 'HbO' ... 'HbR' 'HbR' 'HbR'
    vertex      (flat_vertex) int64 3MB 0 1 2 3 ... 167570 167571 167572 167573
    wavelength  (flat_channel) float64 2kB 760.0 760.0 760.0 ... 850.0 850.0
Dimensions without coordinates: flat_channel, flat_vertex
Attributes:
    units:    1 / molar

xarray.DataArray

• flat_channel: 200
• flat_vertex: 335148

• 9.474e-11 1.559e-11 3.175e-12 2.06e-11 5.291e-10 ... 0.0 0.0 0.0 0.0

  array([[9.47394342e-11, 1.55855606e-11, 3.17494604e-12, ...,
          0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
         [1.44718670e-11, 4.70694883e-13, 1.68621792e-12, ...,
          0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
         [3.72846427e-11, 7.24427616e-12, 2.80946079e-12, ...,
          0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
         ...,
         [2.31970607e-09, 1.99178096e-09, 7.19693148e-10, ...,
          0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
         [1.94127001e-09, 3.25592495e-09, 1.34307726e-09, ...,
          4.58693300e-18, 0.00000000e+00, 0.00000000e+00],
         [1.41597565e-08, 1.30602549e-08, 5.93236189e-09, ...,
          0.00000000e+00, 0.00000000e+00, 0.00000000e+00]])

• Coordinates: (4)
  • is_brain
    (flat_vertex)
    bool
    True True True ... False False

    array([ True,  True,  True, ..., False, False, False])

  • chromo
    (flat_vertex)
    <U3
    'HbO' 'HbO' 'HbO' ... 'HbR' 'HbR'

    array(['HbO', 'HbO', 'HbO', ..., 'HbR', 'HbR', 'HbR'], dtype='<U3')

  • vertex
    (flat_vertex)
    int64
    0 1 2 3 ... 167571 167572 167573

    array([     0,      1,      2, ..., 167571, 167572, 167573])

  • wavelength
    (flat_channel)
    float64
    760.0 760.0 760.0 ... 850.0 850.0

    array([760., 760., 760., 760., 760., 760., 760., 760., 760., 760., 760.,
           760., 760., 760., 760., 760., 760., 760., 760., 760., 760., 760.,
           760., 760., 760., 760., 760., 760., 760., 760., 760., 760., 760.,
           760., 760., 760., 760., 760., 760., 760., 760., 760., 760., 760.,
           760., 760., 760., 760., 760., 760., 760., 760., 760., 760., 760.,
           760., 760., 760., 760., 760., 760., 760., 760., 760., 760., 760.,
           760., 760., 760., 760., 760., 760., 760., 760., 760., 760., 760.,
           760., 760., 760., 760., 760., 760., 760., 760., 760., 760., 760.,
           760., 760., 760., 760., 760., 760., 760., 760., 760., 760., 760.,
           760., 850., 850., 850., 850., 850., 850., 850., 850., 850., 850.,
           850., 850., 850., 850., 850., 850., 850., 850., 850., 850., 850.,
           850., 850., 850., 850., 850., 850., 850., 850., 850., 850., 850.,
           850., 850., 850., 850., 850., 850., 850., 850., 850., 850., 850.,
           850., 850., 850., 850., 850., 850., 850., 850., 850., 850., 850.,
           850., 850., 850., 850., 850., 850., 850., 850., 850., 850., 850.,
           850., 850., 850., 850., 850., 850., 850., 850., 850., 850., 850.,
           850., 850., 850., 850., 850., 850., 850., 850., 850., 850., 850.,
           850., 850., 850., 850., 850., 850., 850., 850., 850., 850., 850.,
           850., 850.])

• Indexes: (0)
• Attributes: (1)

  units :

  1 / molar

Invert the sensitivity matrix

B = pseudo_inverse_stacked(Adot_stacked, alpha = 0.01, alpha_spatial = 0.001)
B

<xarray.DataArray (flat_vertex: 335148, flat_channel: 200)> Size: 536MB
1.207e-17 3.591e-16 -3.076e-16 -2.678e-17 -2.313e-15 ... 0.0 0.0 0.0 0.0 0.0
Coordinates:
    is_brain    (flat_vertex) bool 335kB True True True ... False False False
    chromo      (flat_vertex) <U3 4MB 'HbO' 'HbO' 'HbO' ... 'HbR' 'HbR' 'HbR'
    vertex      (flat_vertex) int64 3MB 0 1 2 3 ... 167570 167571 167572 167573
    wavelength  (flat_channel) float64 2kB 760.0 760.0 760.0 ... 850.0 850.0
Dimensions without coordinates: flat_vertex, flat_channel
Attributes:
    units:    molar

xarray.DataArray

• flat_vertex: 335148
• flat_channel: 200

• 1.207e-17 3.591e-16 -3.076e-16 -2.678e-17 ... 0.0 0.0 0.0 0.0

  array([[ 1.20746728e-17,  3.59079687e-16, -3.07582245e-16, ...,
          -9.11794366e-17, -8.61544299e-17,  5.59867162e-17],
         [-1.77113797e-17,  2.86412523e-17, -6.02250215e-17, ...,
          -6.86585944e-17, -1.46715687e-16, -9.33802769e-17],
         [-1.75484750e-17,  2.04114876e-17, -3.12742571e-17, ...,
          -4.18232499e-17, -6.95913296e-17,  5.00637429e-18],
         ...,
         [ 2.91225399e-33,  8.09462798e-32,  4.71465976e-32, ...,
           6.19449440e-22,  8.81070269e-22, -3.21811311e-23],
         [-7.04561755e-34, -2.35266023e-33, -1.58950009e-33, ...,
           2.81527588e-23,  6.22335791e-26, -3.57629482e-24],
         [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00, ...,
           0.00000000e+00,  0.00000000e+00,  0.00000000e+00]])

• Coordinates: (4)
  • is_brain
    (flat_vertex)
    bool
    True True True ... False False

    array([ True,  True,  True, ..., False, False, False])

  • chromo
    (flat_vertex)
    <U3
    'HbO' 'HbO' 'HbO' ... 'HbR' 'HbR'

    array(['HbO', 'HbO', 'HbO', ..., 'HbR', 'HbR', 'HbR'], dtype='<U3')

  • vertex
    (flat_vertex)
    int64
    0 1 2 3 ... 167571 167572 167573

    array([     0,      1,      2, ..., 167571, 167572, 167573])

  • wavelength
    (flat_channel)
    float64
    760.0 760.0 760.0 ... 850.0 850.0

    array([760., 760., 760., 760., 760., 760., 760., 760., 760., 760., 760.,
           760., 760., 760., 760., 760., 760., 760., 760., 760., 760., 760.,
           760., 760., 760., 760., 760., 760., 760., 760., 760., 760., 760.,
           760., 760., 760., 760., 760., 760., 760., 760., 760., 760., 760.,
           760., 760., 760., 760., 760., 760., 760., 760., 760., 760., 760.,
           760., 760., 760., 760., 760., 760., 760., 760., 760., 760., 760.,
           760., 760., 760., 760., 760., 760., 760., 760., 760., 760., 760.,
           760., 760., 760., 760., 760., 760., 760., 760., 760., 760., 760.,
           760., 760., 760., 760., 760., 760., 760., 760., 760., 760., 760.,
           760., 850., 850., 850., 850., 850., 850., 850., 850., 850., 850.,
           850., 850., 850., 850., 850., 850., 850., 850., 850., 850., 850.,
           850., 850., 850., 850., 850., 850., 850., 850., 850., 850., 850.,
           850., 850., 850., 850., 850., 850., 850., 850., 850., 850., 850.,
           850., 850., 850., 850., 850., 850., 850., 850., 850., 850., 850.,
           850., 850., 850., 850., 850., 850., 850., 850., 850., 850., 850.,
           850., 850., 850., 850., 850., 850., 850., 850., 850., 850., 850.,
           850., 850., 850., 850., 850., 850., 850., 850., 850., 850., 850.,
           850., 850., 850., 850., 850., 850., 850., 850., 850., 850., 850.,
           850., 850.])

• Indexes: (0)
• Attributes: (1)

  units :

  molar

Calculate concentration changes

• the optical density has shape (nchannel, nwavelength, time). Additional dimensions like ‘trial_type’ in this example are allowed, too.

blockaverage

<xarray.DataArray (trial_type: 2, channel: 100, wavelength: 2, reltime: 154)> Size: 493kB
[] -0.0006526 -0.001002 -0.001244 -0.00137 ... -0.0003574 -0.0002813 -0.0001787
Coordinates:
  * reltime     (reltime) float64 1kB -5.038 -4.809 -4.58 ... 29.54 29.77 30.0
  * channel     (channel) object 800B 'S1D1' 'S1D2' 'S1D4' ... 'S14D31' 'S14D32'
    source      (channel) object 800B 'S1' 'S1' 'S1' 'S1' ... 'S14' 'S14' 'S14'
    detector    (channel) object 800B 'D1' 'D2' 'D4' 'D5' ... 'D29' 'D31' 'D32'
  * wavelength  (wavelength) float64 16B 760.0 850.0
  * trial_type  (trial_type) object 16B 'FTapping/Left' 'FTapping/Right'

xarray.DataArray

• trial_type: 2
• channel: 100
• wavelength: 2
• reltime: 154

• [] -0.0006526 -0.001002 -0.001244 ... -0.0003574 -0.0002813 -0.0001787

  Magnitude

  [[[[-0.0006526380230170489 -0.0010020108066544039 -0.001244170613138736 ... 0.00773597175753698 0.00774036187294241 0.0077306068132088285] [-0.0019411336370301533 -0.0022190726530263473 -0.002360941374091604 ... 0.006982190852098116 0.006992210888685543 0.0069783666883401554]] [[-0.0014936653048647128 -0.0017100033926540181 -0.0018182842959475346 ... 0.003431072674552541 0.0034254414764264512 0.003392523672424149] [-0.0022337708587332286 -0.002458924352824458 -0.0025573564327410554 ... 0.0047592312690254894 0.004770547499013742 0.004753363507245642]] [[-0.0012800779070457622 -0.0016089362141727148 -0.0018334954144404608 ... 0.006194047009384458 0.0061685418083578915 0.006154132834262327] [-0.002588081351635877 -0.0028605868292510416 -0.003002211833825404 ... 0.006485299874527477 0.006466747949168435 0.006448422412670403]] ... [[-0.0015562915106271354 -0.0015260578823110232 -0.0014728811786925752 ... 0.0009666426628268548 0.0009958313578833653 0.001010715686953993] [-0.002804133307684544 -0.002880566564134321 -0.0028121366350981948 ... 0.0019926644841683078 0.0021575501146303294 0.002287390958488447]] [[-0.0016493700459431406 -0.0016676017524846632 -0.0016956365693358995 ... 0.0017725864399436133 0.0018474336355347044 0.0019044602928343524] [-0.0033361717383126018 -0.003499238730727136 -0.003510319382752208 ... 0.0027944708467140386 0.0029614836029754634 0.0030935651268180647]] [[-0.001788621662365475 -0.0017317768120953723 -0.0016403145223998312 ... 0.0008487486977927028 0.0009201727863585511 0.0010021225563426879] [-0.003562686563345219 -0.0036259614638027293 -0.0035215065222963674 ... 0.0019603197919563615 0.0021799394569672545 0.002370829507532506]]] [[[-0.0004902055121861838 -0.000632681111234775 -0.0007813781310756729 ... 0.00021471584735898667 0.00016889633512251728 6.19193562013988e-05] [-0.0009860093268402894 -0.0011744259799335854 -0.0012995739657523358 ... 0.0021031198705829244 0.0020892470081507555 0.002033278445597703]] [[-0.00017903831472201864 -0.00022590391905864652 -0.00032575544700045976 ... -0.00045836017389492656 -0.0006530834975963834 -0.000759602447092332] [-0.00046618469123695084 -0.0005932215633016536 -0.0007059910396767626 ... -0.0004889051934100375 -0.0006795318442888155 -0.00077088502776959]] [[-0.00024925047216145106 -0.00039400964364672155 -0.0004890063702298026 ... -0.0006491643413920785 -0.0007105232072191717 -0.0008487159395451555] [-0.0007375538620971675 -0.0008864221530964177 -0.000948567388375058 ... -0.0005204037579280172 -0.0005415911312302692 -0.0006444321779022873]] ... [[-0.0005449489707075885 -0.0006312610722945527 -0.000636945066636855 ... 0.0004803751383537052 0.0005464517856373204 0.0006253175090123594] [-0.0007023172142688516 -0.0007238928008903047 -0.0007091039298522545 ... 0.00021451052094867553 0.00030227084564584896 0.00040344145979711313]] [[-0.0004503188539522116 -0.0005462637287796702 -0.0005353039971004708 ... -0.000896055946390247 -0.0007923386743413031 -0.0006994824603150336] [-0.0005145917002247065 -0.00056092892268372 -0.0005701108903421549 ... -0.0013666897027129277 -0.0012736234521558806 -0.0011736465307203726]] [[-0.0003063654330426074 -0.0003356392496822666 -0.0003155663490409872 ... -9.330757641512287e-05 -6.355442869129639e-05 -6.996910076128036e-06] [-0.0005925984192087044 -0.0005938311925914815 -0.0005584814220281202 ... -0.00035744466085031366 -0.00028133203566686434 -0.00017870029004370645]]]]

  Units

  dimensionless

• Coordinates: (6)
  • reltime
    (reltime)
    float64
    -5.038 -4.809 -4.58 ... 29.77 30.0

    array([-5.038, -4.809, -4.58 , -4.351, -4.122, -3.893, -3.664, -3.435, -3.206,
           -2.977, -2.748, -2.519, -2.29 , -2.061, -1.832, -1.603, -1.374, -1.145,
           -0.916, -0.687, -0.458, -0.229,  0.   ,  0.229,  0.458,  0.687,  0.916,
            1.145,  1.374,  1.603,  1.832,  2.061,  2.29 ,  2.519,  2.748,  2.977,
            3.206,  3.435,  3.664,  3.893,  4.122,  4.351,  4.58 ,  4.809,  5.038,
            5.267,  5.496,  5.725,  5.954,  6.183,  6.412,  6.641,  6.87 ,  7.099,
            7.328,  7.557,  7.786,  8.015,  8.244,  8.473,  8.702,  8.931,  9.16 ,
            9.389,  9.618,  9.847, 10.076, 10.305, 10.534, 10.763, 10.992, 11.221,
           11.45 , 11.679, 11.908, 12.137, 12.366, 12.595, 12.824, 13.053, 13.282,
           13.511, 13.74 , 13.969, 14.198, 14.427, 14.656, 14.885, 15.114, 15.343,
           15.572, 15.801, 16.03 , 16.259, 16.488, 16.717, 16.946, 17.175, 17.404,
           17.633, 17.862, 18.091, 18.32 , 18.549, 18.778, 19.007, 19.236, 19.465,
           19.694, 19.923, 20.152, 20.381, 20.61 , 20.839, 21.068, 21.297, 21.526,
           21.755, 21.984, 22.213, 22.442, 22.671, 22.9  , 23.129, 23.358, 23.587,
           23.816, 24.045, 24.274, 24.503, 24.732, 24.961, 25.19 , 25.419, 25.648,
           25.877, 26.106, 26.335, 26.564, 26.793, 27.022, 27.251, 27.48 , 27.709,
           27.938, 28.167, 28.396, 28.625, 28.854, 29.083, 29.312, 29.541, 29.77 ,
           29.999])

  • channel
    (channel)
    object
    'S1D1' 'S1D2' ... 'S14D31' 'S14D32'

    array(['S1D1', 'S1D2', 'S1D4', 'S1D5', 'S1D6', 'S1D8', 'S2D2', 'S2D3', 'S2D5',
           'S2D6', 'S2D7', 'S2D9', 'S3D1', 'S3D4', 'S3D5', 'S3D8', 'S3D10',
           'S3D11', 'S3D14', 'S4D2', 'S4D4', 'S4D5', 'S4D6', 'S4D7', 'S4D8',
           'S4D9', 'S4D10', 'S4D11', 'S4D12', 'S4D13', 'S4D15', 'S5D3', 'S5D6',
           'S5D7', 'S5D9', 'S5D12', 'S5D13', 'S5D16', 'S6D8', 'S6D10', 'S6D11',
           'S6D12', 'S6D14', 'S6D15', 'S7D9', 'S7D11', 'S7D12', 'S7D13', 'S7D15',
           'S7D16', 'S8D17', 'S8D18', 'S8D20', 'S8D21', 'S8D22', 'S8D24', 'S9D18',
           'S9D19', 'S9D21', 'S9D22', 'S9D23', 'S9D25', 'S10D17', 'S10D20',
           'S10D21', 'S10D24', 'S10D26', 'S10D27', 'S10D30', 'S11D18', 'S11D20',
           'S11D21', 'S11D22', 'S11D23', 'S11D24', 'S11D25', 'S11D26', 'S11D27',
           'S11D28', 'S11D29', 'S11D31', 'S12D19', 'S12D22', 'S12D23', 'S12D25',
           'S12D28', 'S12D29', 'S12D32', 'S13D24', 'S13D26', 'S13D27', 'S13D28',
           'S13D30', 'S13D31', 'S14D25', 'S14D27', 'S14D28', 'S14D29', 'S14D31',
           'S14D32'], dtype=object)

  • source
    (channel)
    object
    'S1' 'S1' 'S1' ... 'S14' 'S14'

    array(['S1', 'S1', 'S1', 'S1', 'S1', 'S1', 'S2', 'S2', 'S2', 'S2', 'S2',
           'S2', 'S3', 'S3', 'S3', 'S3', 'S3', 'S3', 'S3', 'S4', 'S4', 'S4',
           'S4', 'S4', 'S4', 'S4', 'S4', 'S4', 'S4', 'S4', 'S4', 'S5', 'S5',
           'S5', 'S5', 'S5', 'S5', 'S5', 'S6', 'S6', 'S6', 'S6', 'S6', 'S6',
           'S7', 'S7', 'S7', 'S7', 'S7', 'S7', 'S8', 'S8', 'S8', 'S8', 'S8',
           'S8', 'S9', 'S9', 'S9', 'S9', 'S9', 'S9', 'S10', 'S10', 'S10',
           'S10', 'S10', 'S10', 'S10', 'S11', 'S11', 'S11', 'S11', 'S11',
           'S11', 'S11', 'S11', 'S11', 'S11', 'S11', 'S11', 'S12', 'S12',
           'S12', 'S12', 'S12', 'S12', 'S12', 'S13', 'S13', 'S13', 'S13',
           'S13', 'S13', 'S14', 'S14', 'S14', 'S14', 'S14', 'S14'],
          dtype=object)

  • detector
    (channel)
    object
    'D1' 'D2' 'D4' ... 'D31' 'D32'

    array(['D1', 'D2', 'D4', 'D5', 'D6', 'D8', 'D2', 'D3', 'D5', 'D6', 'D7',
           'D9', 'D1', 'D4', 'D5', 'D8', 'D10', 'D11', 'D14', 'D2', 'D4',
           'D5', 'D6', 'D7', 'D8', 'D9', 'D10', 'D11', 'D12', 'D13', 'D15',
           'D3', 'D6', 'D7', 'D9', 'D12', 'D13', 'D16', 'D8', 'D10', 'D11',
           'D12', 'D14', 'D15', 'D9', 'D11', 'D12', 'D13', 'D15', 'D16',
           'D17', 'D18', 'D20', 'D21', 'D22', 'D24', 'D18', 'D19', 'D21',
           'D22', 'D23', 'D25', 'D17', 'D20', 'D21', 'D24', 'D26', 'D27',
           'D30', 'D18', 'D20', 'D21', 'D22', 'D23', 'D24', 'D25', 'D26',
           'D27', 'D28', 'D29', 'D31', 'D19', 'D22', 'D23', 'D25', 'D28',
           'D29', 'D32', 'D24', 'D26', 'D27', 'D28', 'D30', 'D31', 'D25',
           'D27', 'D28', 'D29', 'D31', 'D32'], dtype=object)

  • wavelength
    (wavelength)
    float64
    760.0 850.0

    array([760., 850.])

  • trial_type
    (trial_type)
    object
    'FTapping/Left' 'FTapping/Right'

    array(['FTapping/Left', 'FTapping/Right'], dtype=object)

• Indexes: (4)
  • reltime
    PandasIndex

    PandasIndex(Index([             -5.038,              -4.809,               -4.58,
                        -4.351,              -4.122, -3.8930000000000002,
                        -3.664,              -3.435,              -3.206,
           -2.9770000000000003,
           ...
            27.938000000000002,              28.167,              28.396,
                        28.625,  28.854000000000003,  29.083000000000002,
                        29.312,              29.541,               29.77,
            29.999000000000002],
          dtype='float64', name='reltime', length=154))

  • channel
    PandasIndex

    PandasIndex(Index(['S1D1', 'S1D2', 'S1D4', 'S1D5', 'S1D6', 'S1D8', 'S2D2', 'S2D3', 'S2D5',
           'S2D6', 'S2D7', 'S2D9', 'S3D1', 'S3D4', 'S3D5', 'S3D8', 'S3D10',
           'S3D11', 'S3D14', 'S4D2', 'S4D4', 'S4D5', 'S4D6', 'S4D7', 'S4D8',
           'S4D9', 'S4D10', 'S4D11', 'S4D12', 'S4D13', 'S4D15', 'S5D3', 'S5D6',
           'S5D7', 'S5D9', 'S5D12', 'S5D13', 'S5D16', 'S6D8', 'S6D10', 'S6D11',
           'S6D12', 'S6D14', 'S6D15', 'S7D9', 'S7D11', 'S7D12', 'S7D13', 'S7D15',
           'S7D16', 'S8D17', 'S8D18', 'S8D20', 'S8D21', 'S8D22', 'S8D24', 'S9D18',
           'S9D19', 'S9D21', 'S9D22', 'S9D23', 'S9D25', 'S10D17', 'S10D20',
           'S10D21', 'S10D24', 'S10D26', 'S10D27', 'S10D30', 'S11D18', 'S11D20',
           'S11D21', 'S11D22', 'S11D23', 'S11D24', 'S11D25', 'S11D26', 'S11D27',
           'S11D28', 'S11D29', 'S11D31', 'S12D19', 'S12D22', 'S12D23', 'S12D25',
           'S12D28', 'S12D29', 'S12D32', 'S13D24', 'S13D26', 'S13D27', 'S13D28',
           'S13D30', 'S13D31', 'S14D25', 'S14D27', 'S14D28', 'S14D29', 'S14D31',
           'S14D32'],
          dtype='object', name='channel'))

  • wavelength
    PandasIndex

    PandasIndex(Index([760.0, 850.0], dtype='float64', name='wavelength'))

  • trial_type
    PandasIndex

    PandasIndex(Index(['FTapping/Left', 'FTapping/Right'], dtype='object', name='trial_type'))

• Attributes: (0)

To apply the inverted sensitiviy matrix, the OD wavelength and channel dimensions need to be flattened. Then the inverted sensitivity matrix can be multiplied which contracts over flat_channel and the flat_vertex dimension remains. The flat_vertex dimensions containes vertices of the scalp and the brain for both chromophores. These need to be separated again. The function fw.apply_inv_sensitiviy takes care of all of this.

dC_brain, dC_scalp = fw.apply_inv_sensitivity(blockaverage, B)
display(dC_brain)
display(dC_scalp)

<xarray.DataArray (trial_type: 2, reltime: 154, chromo: 2, vertex: 15002)> Size: 74MB
-5.844e-16 -5.904e-17 -3.294e-17 -6.016e-17 ... -2.332e-16 -7.767e-14 -5.56e-18
Coordinates:
  * chromo      (chromo) <U3 24B 'HbO' 'HbR'
  * vertex      (vertex) int64 120kB 0 1 2 3 4 ... 14997 14998 14999 15000 15001
    is_brain    (chromo, vertex) bool 30kB True True True ... True True True
  * reltime     (reltime) float64 1kB -5.038 -4.809 -4.58 ... 29.54 29.77 30.0
  * trial_type  (trial_type) object 16B 'FTapping/Left' 'FTapping/Right'
Attributes:
    units:    molar

xarray.DataArray

• trial_type: 2
• reltime: 154
• chromo: 2
• vertex: 15002

• -5.844e-16 -5.904e-17 -3.294e-17 ... -2.332e-16 -7.767e-14 -5.56e-18

  array([[[[-5.84414016e-16, -5.90433392e-17, -3.29411256e-17, ...,
            -3.24082791e-16,  1.27341057e-13,  5.59471777e-18],
           [-1.02276878e-16, -1.02533695e-17, -6.96696103e-18, ...,
             1.23375371e-16,  1.32973532e-13,  1.52323649e-17]],
  
          [[-5.87542269e-16, -5.96456029e-17, -3.33657048e-17, ...,
            -3.32817943e-16,  8.26937157e-14,  4.98788679e-18],
           [-1.20297472e-16, -1.23680134e-17, -8.22833798e-18, ...,
             9.83736260e-17,  1.62906988e-13,  1.60327974e-17]],
  
          [[-5.72158538e-16, -5.83417586e-17, -3.27345476e-17, ...,
            -3.29107711e-16,  4.71541500e-14,  4.63824068e-18],
           [-1.34774423e-16, -1.40783186e-17, -9.24147701e-18, ...,
             8.11632881e-17,  1.76049866e-13,  1.63200050e-17]],
  
          ...,
  
          [[ 5.31188521e-16,  5.46279864e-17,  3.11888603e-17, ...,
             4.32994398e-16,  1.14975840e-13,  1.36208320e-17],
           [ 2.90229233e-16,  3.05963722e-17,  1.76342609e-17, ...,
  ...
            -1.82397113e-16,  2.16448579e-13,  7.60329687e-18],
           [-3.19480158e-18, -3.17454871e-19, -1.26453718e-18, ...,
             1.34969998e-16, -1.48406263e-14,  4.59260032e-18]],
  
          ...,
  
          [[ 1.58127008e-16,  1.75689545e-17,  1.07334896e-17, ...,
            -3.32539663e-16,  1.14946449e-13, -2.57931778e-17],
           [ 1.32107961e-16,  1.76015119e-17,  1.18690375e-17, ...,
            -1.82350457e-16, -6.60077361e-14, -6.26436979e-18]],
  
          [[ 1.43176092e-16,  1.60430915e-17,  9.69888575e-18, ...,
            -3.55209994e-16,  1.05338468e-13, -2.62892272e-17],
           [ 1.26324906e-16,  1.67710424e-17,  1.12818703e-17, ...,
            -2.12331804e-16, -7.70797894e-14, -6.36028058e-18]],
  
          [[ 1.33888076e-16,  1.51317479e-17,  8.99776462e-18, ...,
            -3.75709505e-16,  1.15143090e-13, -2.56065736e-17],
           [ 1.24781755e-16,  1.63786974e-17,  1.09258176e-17, ...,
            -2.33206597e-16, -7.76718185e-14, -5.56024053e-18]]]])

• Coordinates: (5)
  • chromo
    (chromo)
    <U3
    'HbO' 'HbR'

    array(['HbO', 'HbR'], dtype='<U3')

  • vertex
    (vertex)
    int64
    0 1 2 3 ... 14998 14999 15000 15001

    array([    0,     1,     2, ..., 14999, 15000, 15001])

  • is_brain
    (chromo, vertex)
    bool
    True True True ... True True True

    array([[ True,  True,  True, ...,  True,  True,  True],
           [ True,  True,  True, ...,  True,  True,  True]])

  • reltime
    (reltime)
    float64
    -5.038 -4.809 -4.58 ... 29.77 30.0

    array([-5.038, -4.809, -4.58 , -4.351, -4.122, -3.893, -3.664, -3.435, -3.206,
           -2.977, -2.748, -2.519, -2.29 , -2.061, -1.832, -1.603, -1.374, -1.145,
           -0.916, -0.687, -0.458, -0.229,  0.   ,  0.229,  0.458,  0.687,  0.916,
            1.145,  1.374,  1.603,  1.832,  2.061,  2.29 ,  2.519,  2.748,  2.977,
            3.206,  3.435,  3.664,  3.893,  4.122,  4.351,  4.58 ,  4.809,  5.038,
            5.267,  5.496,  5.725,  5.954,  6.183,  6.412,  6.641,  6.87 ,  7.099,
            7.328,  7.557,  7.786,  8.015,  8.244,  8.473,  8.702,  8.931,  9.16 ,
            9.389,  9.618,  9.847, 10.076, 10.305, 10.534, 10.763, 10.992, 11.221,
           11.45 , 11.679, 11.908, 12.137, 12.366, 12.595, 12.824, 13.053, 13.282,
           13.511, 13.74 , 13.969, 14.198, 14.427, 14.656, 14.885, 15.114, 15.343,
           15.572, 15.801, 16.03 , 16.259, 16.488, 16.717, 16.946, 17.175, 17.404,
           17.633, 17.862, 18.091, 18.32 , 18.549, 18.778, 19.007, 19.236, 19.465,
           19.694, 19.923, 20.152, 20.381, 20.61 , 20.839, 21.068, 21.297, 21.526,
           21.755, 21.984, 22.213, 22.442, 22.671, 22.9  , 23.129, 23.358, 23.587,
           23.816, 24.045, 24.274, 24.503, 24.732, 24.961, 25.19 , 25.419, 25.648,
           25.877, 26.106, 26.335, 26.564, 26.793, 27.022, 27.251, 27.48 , 27.709,
           27.938, 28.167, 28.396, 28.625, 28.854, 29.083, 29.312, 29.541, 29.77 ,
           29.999])

  • trial_type
    (trial_type)
    object
    'FTapping/Left' 'FTapping/Right'

    array(['FTapping/Left', 'FTapping/Right'], dtype=object)

• Indexes: (4)
  • reltime
    PandasIndex

    PandasIndex(Index([             -5.038,              -4.809,               -4.58,
                        -4.351,              -4.122, -3.8930000000000002,
                        -3.664,              -3.435,              -3.206,
           -2.9770000000000003,
           ...
            27.938000000000002,              28.167,              28.396,
                        28.625,  28.854000000000003,  29.083000000000002,
                        29.312,              29.541,               29.77,
            29.999000000000002],
          dtype='float64', name='reltime', length=154))

  • trial_type
    PandasIndex

    PandasIndex(Index(['FTapping/Left', 'FTapping/Right'], dtype='object', name='trial_type'))

  • chromo
    PandasIndex

    PandasIndex(Index(['HbO', 'HbR'], dtype='object', name='chromo'))

  • vertex
    PandasIndex

    PandasIndex(Index([    0,     1,     2,     3,     4,     5,     6,     7,     8,     9,
           ...
           14992, 14993, 14994, 14995, 14996, 14997, 14998, 14999, 15000, 15001],
          dtype='int64', name='vertex', length=15002))

• Attributes: (1)

  units :

  molar

<xarray.DataArray (trial_type: 2, reltime: 154, chromo: 2, vertex: 152572)> Size: 752MB
-9.254e-24 0.0 -3.873e-23 0.0 0.0 ... -1.859e-21 0.0 -3.872e-22 -3.289e-23 0.0
Coordinates:
  * chromo      (chromo) <U3 24B 'HbO' 'HbR'
  * vertex      (vertex) int64 1MB 15002 15003 15004 ... 167571 167572 167573
    is_brain    (chromo, vertex) bool 305kB False False False ... False False
  * reltime     (reltime) float64 1kB -5.038 -4.809 -4.58 ... 29.54 29.77 30.0
  * trial_type  (trial_type) object 16B 'FTapping/Left' 'FTapping/Right'
Attributes:
    units:    molar

xarray.DataArray

• trial_type: 2
• reltime: 154
• chromo: 2
• vertex: 152572

• -9.254e-24 0.0 -3.873e-23 0.0 0.0 ... 0.0 -3.872e-22 -3.289e-23 0.0

  array([[[[-9.25381026e-24,  0.00000000e+00, -3.87260986e-23, ...,
            -2.29218703e-23,  3.84882214e-24,  0.00000000e+00],
           [-1.17422273e-24,  0.00000000e+00,  8.58704201e-25, ...,
             3.79068194e-22,  3.87395404e-23,  0.00000000e+00]],
  
          [[-8.61157853e-24,  0.00000000e+00, -3.63886627e-23, ...,
            -2.34715656e-23,  4.93242160e-24,  0.00000000e+00],
           [-2.36693736e-25,  0.00000000e+00,  5.27030439e-24, ...,
             4.04450856e-22,  4.12379787e-23,  0.00000000e+00]],
  
          [[-7.92903996e-24,  0.00000000e+00, -3.39870426e-23, ...,
            -8.96135445e-24,  6.85124300e-24,  0.00000000e+00],
           [ 6.25192807e-25,  0.00000000e+00,  8.73846839e-24, ...,
             4.24346506e-22,  4.30874095e-23,  0.00000000e+00]],
  
          ...,
  
          [[ 1.04097671e-24,  0.00000000e+00,  6.46578703e-24, ...,
             1.15537721e-21,  1.02893027e-22,  0.00000000e+00],
           [-3.13782851e-24,  0.00000000e+00, -2.17673981e-23, ...,
  ...
             1.89700171e-22,  1.51146375e-23,  0.00000000e+00],
           [ 6.21274255e-24,  0.00000000e+00,  3.23768686e-23, ...,
             2.41588634e-22,  2.55136200e-23,  0.00000000e+00]],
  
          ...,
  
          [[-1.27706636e-23,  0.00000000e+00, -5.21559053e-23, ...,
            -1.26388796e-21, -1.09224813e-22,  0.00000000e+00],
           [ 7.27702270e-24,  0.00000000e+00,  4.81363306e-23, ...,
            -5.05210470e-22, -4.52499082e-23,  0.00000000e+00]],
  
          [[-1.42053712e-23,  0.00000000e+00, -5.89112745e-23, ...,
            -1.24888861e-21, -1.07100622e-22,  0.00000000e+00],
           [ 6.88584544e-24,  0.00000000e+00,  4.66233858e-23, ...,
            -4.73229593e-22, -4.19389461e-23,  0.00000000e+00]],
  
          [[-1.56174960e-23,  0.00000000e+00, -6.53757003e-23, ...,
            -1.16676285e-21, -9.84500608e-23,  0.00000000e+00],
           [ 6.63665240e-24,  0.00000000e+00,  4.60643929e-23, ...,
            -3.87225501e-22, -3.28862736e-23,  0.00000000e+00]]]])

• Coordinates: (5)
  • chromo
    (chromo)
    <U3
    'HbO' 'HbR'

    array(['HbO', 'HbR'], dtype='<U3')

  • vertex
    (vertex)
    int64
    15002 15003 15004 ... 167572 167573

    array([ 15002,  15003,  15004, ..., 167571, 167572, 167573])

  • is_brain
    (chromo, vertex)
    bool
    False False False ... False False

    array([[False, False, False, ..., False, False, False],
           [False, False, False, ..., False, False, False]])

  • reltime
    (reltime)
    float64
    -5.038 -4.809 -4.58 ... 29.77 30.0

    array([-5.038, -4.809, -4.58 , -4.351, -4.122, -3.893, -3.664, -3.435, -3.206,
           -2.977, -2.748, -2.519, -2.29 , -2.061, -1.832, -1.603, -1.374, -1.145,
           -0.916, -0.687, -0.458, -0.229,  0.   ,  0.229,  0.458,  0.687,  0.916,
            1.145,  1.374,  1.603,  1.832,  2.061,  2.29 ,  2.519,  2.748,  2.977,
            3.206,  3.435,  3.664,  3.893,  4.122,  4.351,  4.58 ,  4.809,  5.038,
            5.267,  5.496,  5.725,  5.954,  6.183,  6.412,  6.641,  6.87 ,  7.099,
            7.328,  7.557,  7.786,  8.015,  8.244,  8.473,  8.702,  8.931,  9.16 ,
            9.389,  9.618,  9.847, 10.076, 10.305, 10.534, 10.763, 10.992, 11.221,
           11.45 , 11.679, 11.908, 12.137, 12.366, 12.595, 12.824, 13.053, 13.282,
           13.511, 13.74 , 13.969, 14.198, 14.427, 14.656, 14.885, 15.114, 15.343,
           15.572, 15.801, 16.03 , 16.259, 16.488, 16.717, 16.946, 17.175, 17.404,
           17.633, 17.862, 18.091, 18.32 , 18.549, 18.778, 19.007, 19.236, 19.465,
           19.694, 19.923, 20.152, 20.381, 20.61 , 20.839, 21.068, 21.297, 21.526,
           21.755, 21.984, 22.213, 22.442, 22.671, 22.9  , 23.129, 23.358, 23.587,
           23.816, 24.045, 24.274, 24.503, 24.732, 24.961, 25.19 , 25.419, 25.648,
           25.877, 26.106, 26.335, 26.564, 26.793, 27.022, 27.251, 27.48 , 27.709,
           27.938, 28.167, 28.396, 28.625, 28.854, 29.083, 29.312, 29.541, 29.77 ,
           29.999])

  • trial_type
    (trial_type)
    object
    'FTapping/Left' 'FTapping/Right'

    array(['FTapping/Left', 'FTapping/Right'], dtype=object)

• Indexes: (4)
  • reltime
    PandasIndex

    PandasIndex(Index([             -5.038,              -4.809,               -4.58,
                        -4.351,              -4.122, -3.8930000000000002,
                        -3.664,              -3.435,              -3.206,
           -2.9770000000000003,
           ...
            27.938000000000002,              28.167,              28.396,
                        28.625,  28.854000000000003,  29.083000000000002,
                        29.312,              29.541,               29.77,
            29.999000000000002],
          dtype='float64', name='reltime', length=154))

  • trial_type
    PandasIndex

    PandasIndex(Index(['FTapping/Left', 'FTapping/Right'], dtype='object', name='trial_type'))

  • chromo
    PandasIndex

    PandasIndex(Index(['HbO', 'HbR'], dtype='object', name='chromo'))

  • vertex
    PandasIndex

    PandasIndex(Index([ 15002,  15003,  15004,  15005,  15006,  15007,  15008,  15009,  15010,
            15011,
           ...
           167564, 167565, 167566, 167567, 167568, 167569, 167570, 167571, 167572,
           167573],
          dtype='int64', name='vertex', length=152572))

• Attributes: (1)

  units :

  molar

Plot concentration changes

Using cedalion plot functions to visualie image recon results

Using Scalp Plot Functionality to Create a Gif of Brain Activity on the Scalp

This gives us the activity in channel space across time.

from cedalion.plots import scalp_plot_gif

warnings.filterwarnings("ignore", category=UnitStrippedWarning)

# configure the plot
data_ts = blockaverage.sel(wavelength=850, trial_type="FTapping/Right")
data_ts = data_ts.rename({"reltime": "time"})
geo3d = rec.geo3d
filename_scalp = 'scalp_plot_ts'

# call plot function
scalp_plot_gif(
    data_ts,
    geo3d,
    filename = filename_scalp,
    time_range=(-5,30,0.5)*units.s,
    scl=(-0.01, 0.01),
    fps=6,
    optode_size=6,
    optode_labels=True,
    str_title='850 nm'
    )

display(Image(data=open(filename_scalp+'.gif','rb').read(), format='png'))

Using Image Recon View Functionality to Create a Gif of Activity on the Brain

This gives us activity on the brain from a single view as still image or across time

from cedalion.plots import image_recon_view

filename_view = 'image_recon_view'

X_ts = xr.concat([dC_brain.sel(trial_type="FTapping/Right"), dC_scalp.sel(trial_type="FTapping/Right")], dim="vertex")
X_ts = X_ts.rename({"reltime": "time"})
X_ts = X_ts.transpose("vertex", "chromo", "time")
X_ts = X_ts.assign_coords(is_brain=('vertex', Adot.is_brain.values))

scl = np.percentile(np.abs(X_ts.sel(chromo='HbO').values.reshape(-1)),99)
clim = (-scl,scl)

image_recon_view(
    X_ts,  # time series data; can be 2D (static) or 3D (dynamic)
    head,
    cmap='seismic',
    clim=clim,
    view_type='hbo_brain',
    view_position='left',
    title_str='HbO',
    filename=filename_view,
    SAVE=True,
    time_range=(-5,30,0.5)*units.s,
    fps=6,
    geo3d_plot = geo3d_plot,
    wdw_size = (1024, 768)
)

display(Image(data=open(filename_view+'.gif','rb').read(), format='png'))

Alternatively, we can just select on time point and plot activity as a still image at that time.

# selects the nearest time sample at t=10s in X_ts
X_ts = X_ts.sel(time=4*units.s, method="nearest")

filename_view = 'image_recon_view_still'

image_recon_view(
    X_ts,  # time series data; can be 2D (static) or 3D (dynamic)
    head,
    cmap='seismic',
    clim=clim,
    view_type='hbo_brain',
    view_position='left',
    title_str='HbO',
    filename=filename_view,
    SAVE=True,
    time_range=(-5,30,0.5)*units.s,
    fps=6,
    geo3d_plot = geo3d_plot,
    wdw_size = (1024, 768)
)

display(Image(data=open(filename_view+'.png','rb').read(), format='png'))

Using Image Recon Multi View Functionality to Create a Gif of Activity on the Brain

This gives us activity on the brain after recon from all angles as still image or across time

from cedalion.plots import image_recon_multi_view

filename_multiview = 'image_recon_multiview'

# prepare data
X_ts = xr.concat([dC_brain.sel(trial_type="FTapping/Right"), dC_scalp.sel(trial_type="FTapping/Right")], dim="vertex")
X_ts = X_ts.rename({"reltime": "time"})
X_ts = X_ts.transpose("vertex", "chromo", "time")
X_ts = X_ts.assign_coords(is_brain=('vertex', Adot.is_brain.values))

scl = np.percentile(np.abs(X_ts.sel(chromo='HbO').values.reshape(-1)),99)
clim = (-scl,scl)


image_recon_multi_view(
    X_ts,  # time series data; can be 2D (static) or 3D (dynamic)
    head,
    cmap='seismic',
    clim=clim,
    view_type='hbo_brain',
    title_str='HbO',
    filename=filename_multiview,
    SAVE=True,
    time_range=(-5,30,0.5)*units.s,
    fps=6,
    geo3d_plot = None, #  geo3d_plot
    wdw_size = (1024, 768)
)

display(Image(data=open(filename_multiview+'.gif','rb').read(), format='png'))

Using Image Recon Multi View Functionality to Create a Gif of Activity on the Scalp

This gives us activity on the scalp after recon from all angles as still image or across time

from cedalion.plots import image_recon_multi_view

filename_multiview_scalp = 'image_recon_multiview_scalp'

# prepare data
X_ts = xr.concat([dC_brain.sel(trial_type="FTapping/Right"), dC_scalp.sel(trial_type="FTapping/Right")], dim="vertex")
X_ts = X_ts.rename({"reltime": "time"})
X_ts = X_ts.transpose("vertex", "chromo", "time")
X_ts = X_ts.assign_coords(is_brain=('vertex', Adot.is_brain.values))

scl = np.percentile(np.abs(X_ts.sel(chromo='HbO').values.reshape(-1)),99)
clim = (-scl,scl)


image_recon_multi_view(
    X_ts,  # time series data; can be 2D (static) or 3D (dynamic)
    head,
    cmap='seismic',
    clim=clim,
    view_type='hbo_scalp',
    title_str='HbO',
    filename=filename_multiview_scalp,
    SAVE=True,
    time_range=(-5,30,0.5)*units.s,
    fps=6,
    geo3d_plot = geo3d_plot,
    wdw_size = (1024, 768)
)

display(Image(data=open(filename_multiview_scalp+'.gif','rb').read(), format='png'))