Image Reconstruction

[1]:

# This cells setups the environment when executed in Google Colab.
try:
    import google.colab
    !curl -s https://raw.githubusercontent.com/ibs-lab/cedalion/dev/scripts/colab_setup.py -o colab_setup.py
    # Select branch with --branch "branch name" (default is "dev")
    %run colab_setup.py
except ImportError:
    pass

Notebook configuration

Decide for an example dataset 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.

[2]:

# 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

[3]:

import pyvista as pv

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

import os
import time
import warnings
from pathlib import Path
from tempfile import TemporaryDirectory

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

import cedalion
import cedalion.dataclasses as cdc
import cedalion.datasets
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.sigproc.quality as quality
import cedalion.sigproc.motion_correct as motion_correct
import cedalion.vis.plot_sensitivity_matrix
from cedalion import units
from cedalion.imagereco.solver import pseudo_inverse_stacked
from cedalion.io.forward_model import FluenceFile, load_Adot

import cedalion.xrutils as xrutils
xrutils.unit_stripping_is_error()


xr.set_options(display_expand_data=False);

[4]:

# helper function to display gifs in rendered notbooks
def display_image(fname : str):
    display(Image(data=open(fname,'rb').read(), format='png'))

Working Directory

In this notebook the output of the fluence and sensitivity calculations are stored in a temporary directory. This will be deleted when the notebook ends.

[5]:

temporary_directory = TemporaryDirectory()
tmp_dir_path = Path(temporary_directory.name)

Load a finger-tapping dataset

For this demo we load an example finger-tapping recording through either cedalion.datasets.get_fingertapping or cedalion.datasets.get_fingertappingDOT.

The snirf files of these datasets contains a single NIRS element with one block of raw amplitude data.

[6]:

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’.

[7]:

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

[8]:

# visualize the montage
cedalion.plots.plot_montage3D(rec["amp"], geo3d_meas)

../../_images/examples_head_models_40_image_reconstruction_12_0.png

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

[9]:

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.

[10]:

rec.stim

[10]:

	onset	duration	value	trial_type
0	23.855104	10.0	1.0	1
1	54.132736	10.0	1.0	1
2	84.410368	10.0	1.0	1
3	114.688000	10.0	1.0	1
4	146.112512	10.0	1.0	1
...	...	...	...	...
125	1431.535616	10.0	1.0	5
126	1526.038528	10.0	1.0	5
127	1650.819072	10.0	1.0	5
128	1805.418496	10.0	1.0	5
129	1931.116544	10.0	1.0	5

130 rows × 4 columns

For clarity, events are given more descriptive names:

[11]:

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"
    } )

# count number of trials per trial_type
display(
    rec.stim.groupby("trial_type")[["onset"]]
    .count()
    .rename({"onset": "#trials"}, axis=1)
)

	#trials
trial_type
BallSqueezing/Left	17
BallSqueezing/Right	16
Control	65
FTapping/Left	16
FTapping/Right	16

Preprocessing

Perform motion correction, conversion to optical density and bandpass filtering.

[12]:

rec["od"] = cedalion.nirs.int2od(rec["amp"])
rec["od_tddr"] = motion_correct.tddr(rec["od"])
rec["od_wavelet"] = motion_correct.wavelet(rec["od_tddr"])

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

Calculate block averages in optical density

[13]:

# segment data into epochs
epochs = rec["od_freqfiltered"].cd.to_epochs(
    rec.stim,  # stimulus dataframe
    ["FTapping/Left", "FTapping/Right"],  # select fingertapping events, discard others
    before=5 * units.s,  # seconds before stimulus
    after=30 * 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()

../../_images/examples_head_models_40_image_reconstruction_22_0.png

Load segmented MRI scan

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

[14]:

if HEAD_MODEL == "colin27":
    SEG_DATADIR, mask_files, landmarks_file = cedalion.datasets.get_colin27_segmentation()
    PARCEL_DIR = cedalion.datasets.get_colin27_parcel_file()
elif HEAD_MODEL == "icbm152":
    SEG_DATADIR, mask_files, landmarks_file = cedalion.datasets.get_icbm152_segmentation()
    PARCEL_DIR = cedalion.datasets.get_icbm152_parcel_file()
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.

[15]:

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

'/home/runner/.cache/cedalion/v25.1.0/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.

[16]:

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"),
    landmarks_ras_file=landmarks_file,
    parcel_file=PARCEL_DIR,
    brain_face_count=None,
    scalp_face_count=None
)

[17]:

head.segmentation_masks

[17]:

<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

[18]:

head.landmarks

[18]:

<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

[19]:

head.brain

[19]:

TrimeshSurface(mesh=<trimesh.Trimesh(vertices.shape=(15002, 3), faces.shape=(29988, 3))>, crs='ijk', units=<Unit('dimensionless')>, vertex_coords={'parcel': array(['VisCent_ExStr_8_LH', 'VisCent_ExStr_8_LH', 'VisCent_Striate_2_LH',
       ..., 'SalVentAttnA_ParMed_6_RH', 'ContA_PFCl_5_RH',
       'Background+FreeSurfer_Defined_Medial_Wall_RH'],
      shape=(15002,), dtype='<U44')})

[20]:

head.scalp

[20]:

TrimeshSurface(mesh=<trimesh.Trimesh(vertices.shape=(10050, 3), faces.shape=(20096, 3))>, crs='ijk', units=<Unit('dimensionless')>, vertex_coords={})

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

[21]:

head.crs

[21]:

'ijk'

The transformation matrix to translate from voxel to scanner space:

[22]:

head.t_ijk2ras

[22]:

<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

Changing between coordinate systems:

[23]:

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')>, vertex_coords={'parcel': array(['VisCent_ExStr_8_LH', 'VisCent_ExStr_8_LH', 'VisCent_Striate_2_LH',
       ..., 'SalVentAttnA_ParMed_6_RH', 'ContA_PFCl_5_RH',
       'Background+FreeSurfer_Defined_Medial_Wall_RH'],
      shape=(15002,), dtype='<U44')})

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.

[24]:

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

<xarray.DataArray (label: 346, ijk: 3)> Size: 8kB
[] 11.61 128.0 90.59 22.65 134.1 125.5 ... 172.7 138.4 33.95 173.3 127.9 33.21
Coordinates:
    type     (label) object 3kB PointType.SOURCE ... PointType.LANDMARK
  * label    (label) <U6 8kB 'S1' 'S2' 'S3' 'S4' ... 'FFT10h' 'FT10h' 'FTT10h'
Dimensions without coordinates: ijk

[25]:

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()

../../_images/examples_head_models_40_image_reconstruction_42_0.png

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.

[26]:

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).

[27]:

if PRECOMPUTED_FLUENCE:
    if FORWARD_MODEL == "MCX":
        fluence_fname = cedalion.datasets.get_precomputed_fluence(DATASET, HEAD_MODEL)
    elif FORWARD_MODEL == "NIRFASTER":
        raise NotImplementedError(
            "Currently there are no precomputed NIRFASTER results available"
        )
else:
    fluence_fname = tmp_dir_path / "fluence.h5"

    if FORWARD_MODEL == "MCX":
        fwm.compute_fluence_mcx(fluence_fname)
    elif FORWARD_MODEL == "NIRFASTER":
        fwm.compute_fluence_nirfaster(fluence_fname)

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

fluence_all is a xr.DataArray with dimensions: (‘label’, ‘wavelength’, ‘i’, ‘j’, ‘k’),

i.e. for each optode and wavelength it stores the 3D image of the computed fluence in each voxel
fluence_at_optodes is a xr.DataArray with dimensions: (‘optode1’, ‘optode2’, ‘wavelength’).

It contains the fluence directly at the position of the optodes, used for normalization purposes.

Both arrays are stored on disk in the hdf5 file at fluence_fname and should be queried through cedalion.io.forward_model.FluenceFile.

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

Plot fluence

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

[28]:

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

[29]:

time.sleep(1)

plt = pv.Plotter()

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

# fluence_file.get_fluence returns a 3D numpy array with the fluence
# for a specified source and wavelength.
with FluenceFile(fluence_fname) as fluence_file:
    f = fluence_file.get_fluence(src, wl) * fluence_file.get_fluence(det, wl)

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")
cog = cog.pint.dequantify().values
plt.camera.position = cog + [-300, 30, 100]
plt.camera.focal_point = cog
plt.camera.up = [0, 0, 1]

plt.show()

../../_images/examples_head_models_40_image_reconstruction_52_0.png

Calculate concentration changes

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

[35]:

blockaverage

[35]:

<xarray.DataArray (trial_type: 2, channel: 100, wavelength: 2, reltime: 154)> Size: 493kB
[] -0.001023 -0.001108 -0.001151 -0.001155 ... 0.0009149 0.001069 0.001239
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'

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

[36]:

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

Convert concentration changes into micromolar.

[37]:

dC_brain = dC_brain.pint.quantify().pint.to("uM").pint.dequantify()
dC_scalp = dC_scalp.pint.quantify().pint.to("uM").pint.dequantify()

Visualizing image reconstruction results

In the following, different cedalion plot functions will be showcased to visualize concentration changes on the brain and scalp surfaces.

Channel Space

The function cedalion.plots.scalp_plot_gif allows to create an animated gif of channel-space OD changes projected on the scalp.

Here, it is used to show the time course of the blockaveraged OD changes.

[38]:

from cedalion.plots import scalp_plot_gif

# 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="OD 850 nm",
)

[39]:

display_image(f"{filename_scalp}.gif")

../../_images/examples_head_models_40_image_reconstruction_72_0.png

Image Space

Single-View Animations of Activitations on the Brain

The function cedalion.plots.image_recon_view allows to create an animated gif of image-space concentration changes projected on the brain.

[40]:

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 / uM',
    filename=filename_view,
    SAVE=True,
    time_range=(-5,30,0.5)*units.s,
    fps=6,
    geo3d_plot = geo3d_plot,
    wdw_size = (1024, 768)
)

[41]:

display_image(f"{filename_view}.gif")

../../_images/examples_head_models_40_image_reconstruction_75_0.png

Alternatively, we can just select a single time point and plot activity as a still image at that time. Note the different file suffix (.png).

[42]:

# selects the nearest time sample at t=4s in X_ts
X_ts_plot = X_ts.sel(time=4, method="nearest") # note: sel does not accept quantified units

filename_view = 'image_recon_view_still'

image_recon_view(
    X_ts_plot,  # 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 / uM',
    filename=filename_view,
    SAVE=True,
    time_range=(-5,30,0.5)*units.s,
    fps=6,
    geo3d_plot = geo3d_plot,
    wdw_size = (1024, 768)
)

../../_images/examples_head_models_40_image_reconstruction_77_0.png

[43]:

display_image(f"{filename_view}.png")

../../_images/examples_head_models_40_image_reconstruction_78_0.png

Multi-View Animations of Activitations on the Brain

The function cedalion.plots.image_recon_multi_view shows the activity on the brain from all angles as still image or animated across time:

[44]:

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 / uM',
    filename=filename_multiview,
    SAVE=True,
    time_range=(-5,30,0.5)*units.s,
    fps=6,
    geo3d_plot = None, #  geo3d_plot
    wdw_size = (1024, 768)
)

[45]:

display_image(f"{filename_multiview}.gif")

../../_images/examples_head_models_40_image_reconstruction_81_0.png

Multi-View Animations of Activitations on the Scalp

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

[46]:

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 / uM',
    filename=filename_multiview_scalp,
    SAVE=True,
    time_range=(-5,30,0.5)*units.s,
    fps=6,
    geo3d_plot = geo3d_plot,
    wdw_size = (1024, 768)
)

[47]:

display_image(f"{filename_multiview_scalp}.gif")

../../_images/examples_head_models_40_image_reconstruction_84_0.png

Parcel Space

The Schaefer Atlas [SKG+18] as implemented in Cedalion provides nearly 600 labels for different regions of the brain. Each vertex of the brain surface has its correspondng parcel label assigned as a coordinate.

[48]:

head.brain.vertices

[48]:

<xarray.DataArray (label: 15002, ijk: 3)> Size: 360kB
[] 77.59 20.68 74.43 84.24 20.3 70.42 ... 130.4 152.2 100.9 96.35 105.7 87.96
Coordinates:
  * label    (label) int64 120kB 0 1 2 3 4 5 ... 14997 14998 14999 15000 15001
  * parcel   (label) <U44 3MB 'VisCent_ExStr_8_LH' ... 'Background+FreeSurfer...
Dimensions without coordinates: ijk

Using parcel labels, vertices belonging to the same brain region can be easily grouped together with the DataArray.groupby function.

To obtain the average hemodynamic response in a parcel, the baseline-subtraced concentration changes of all vertices in a parcel are averaged. As baseline the first sample is used.

[49]:

# subtract baseline
baseline =  dC_brain.isel(reltime=0) # first sample along reltime dimension
dC_brain_blsubtracted = dC_brain - baseline

# average over parcels
avg_HbO = dC_brain_blsubtracted.sel(chromo="HbO").groupby('parcel').mean()

avg_HbR = dC_brain_blsubtracted.sel(chromo="HbR").groupby('parcel').mean()

display(dC_brain_blsubtracted.rename("dC_brain_blsubtracted"))
display(avg_HbO.rename("avg_HbO"))

The montage in this dataset covers only parts of the head. Consequently, many brain regions lack significant signal coverage due to the absence of optodes.

To focus on relevant regions, a subset of parcels from the somatosensory and motor regions in both hemispheres is selected.

[50]:

selected_parcels = [
    "SomMotA_1_LH", "SomMotA_3_LH", "SomMotA_4_LH", "SomMotA_5_LH", "SomMotA_9_LH", "SomMotA_10_LH",
    "SomMotA_1_RH", "SomMotA_2_RH", "SomMotA_3_RH", "SomMotA_4_RH", "SomMotA_6_RH", "SomMotA_7_RH"
]

The following plot visualizes the montage and the selected parcels:

[51]:

# map parcel labels to colors
parcel_colors = {
    parcel: p.cm.jet(i / (len(selected_parcels) - 1))
    for i, parcel in enumerate(selected_parcels)
}

# assign colors to vertices
b = cdc.VTKSurface.from_trimeshsurface(head.brain)
b = pv.wrap(b.mesh)
b["parcels"] = np.asarray([
    parcel_colors.get(parcel, (0.8, 0.8, 0.8, .3))
    for parcel in head.brain.vertices.parcel.values
])

plt = pv.Plotter()

plt.add_mesh(
    b,
    scalars="parcels",
    rgb=True,
    smooth_shading=False
)

cedalion.plots.plot_labeled_points(plt, geo3d_plot)

legends = [a for a in parcel_colors.items()]
plt.add_legend(labels= legends, face='o', size=(0.3,0.3))

cog = head.brain.vertices.mean("label")
cog = cog.pint.dequantify().values

plt.camera.position = cog + [0,0,400]
plt.camera.focal_point = cog
plt.camera.up = [0,1,0]
plt.reset_camera()

plt.show()

../../_images/examples_head_models_40_image_reconstruction_92_0.png

Plot averaged time traces in each parcel for the ‘FTapping/Right’ and ‘FTapping/Left’ conditions:

[52]:

f,ax = p.subplots(2,6, figsize=(20,5))
ax = ax.flatten()
for i_par, par in enumerate(selected_parcels):
    ax[i_par].plot(avg_HbO.sel(parcel = par, trial_type = "FTapping/Right").reltime, avg_HbO.sel(parcel = par, trial_type = "FTapping/Right").values, "r", lw=2, ls='-')
    ax[i_par].plot(avg_HbR.sel(parcel = par, trial_type = "FTapping/Right").reltime, avg_HbR.sel(parcel = par, trial_type = "FTapping/Right").values, "b", lw=2, ls='-')

    ax[i_par].grid(1)
    ax[i_par].set_title(par, color=parcel_colors[par])
    ax[i_par].set_ylim(-.05, .2)

p.suptitle("Parcellations: HbO: r | HbR: b | FTapping/Right", y=1)
p.tight_layout()

f,ax = p.subplots(2,6, figsize=(20,5))
ax = ax.flatten()
for i_par, par in enumerate(selected_parcels):
    ax[i_par].plot(avg_HbO.sel(parcel = par, trial_type = "FTapping/Left").reltime, avg_HbO.sel(parcel = par, trial_type = "FTapping/Left").values, "r", lw=2, ls='-')
    ax[i_par].plot(avg_HbR.sel(parcel = par, trial_type = "FTapping/Left").reltime, avg_HbR.sel(parcel = par, trial_type = "FTapping/Left").values, "b", lw=2, ls='-')

    ax[i_par].grid(1)
    ax[i_par].set_title(par, color=parcel_colors[par])
    ax[i_par].set_ylim(-.05, .2)

p.suptitle("Parcellations: HbO: r | HbR: b | FTapping/Left", y=1)
p.tight_layout()

../../_images/examples_head_models_40_image_reconstruction_94_0.png

../../_images/examples_head_models_40_image_reconstruction_94_1.png

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

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