S7: Data Augmentation

This example notebook illustrates the functionality in cedalion.sim.synthetic_hrf and cedalion.sim.synthetic_artifact to create simulated datasets with added activations and motion artifacts.

It has two parts:

Adding synthetic activations
Adding synthetic motion artifacts

[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

[2]:

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

[3]:

import matplotlib.pyplot as plt
import numpy as np
import pyvista as pv
import xarray as xr

import cedalion
import cedalion.dataclasses as cdc
import cedalion.data
import cedalion.geometry.landmarks as cd_landmarks
import cedalion.dot as dot
import cedalion.models.glm as glm
import cedalion.nirs

import cedalion.vis.anatomy
import cedalion.vis.blocks as vbx

import cedalion.sigproc.quality as quality
import cedalion.sim.synthetic_hrf as synhrf
import cedalion.sim.synthetic_artifact as synart
import cedalion.xrutils as xrutils
from cedalion import units
from IPython.display import Image

xr.set_options(display_expand_data=False)

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

[4]:

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

Adding Synthetic Activations

Loading and preprocessing the dataset

This notebook uses a high-density, whole head resting state dataset recorded with a NinjaNIRS 22.

[5]:

rec = cedalion.data.get_nn22_resting_state()

geo3d = rec.geo3d
meas_list = rec._measurement_lists["amp"]

amp = rec["amp"]
amp = amp.pint.dequantify().pint.quantify("V")

display(amp)

<xarray.DataArray (channel: 567, wavelength: 2, time: 3309)> Size: 30MB
[V] 0.0238 0.02385 0.02375 0.02361 0.02364 ... 0.2887 0.2892 0.2897 0.291 0.2923
Coordinates:
  * time        (time) float64 26kB 0.0 0.1112 0.2225 ... 367.8 367.9 368.0
    samples     (time) int64 26kB 0 1 2 3 4 5 ... 3303 3304 3305 3306 3307 3308
  * channel     (channel) object 5kB 'S10D87' 'S10D94' ... 'S47D29' 'S5D137'
    source      (channel) object 5kB 'S10' 'S10' 'S10' ... 'S56' 'S47' 'S5'
    detector    (channel) object 5kB 'D87' 'D94' 'D89' ... 'D129' 'D29' 'D137'
  * wavelength  (wavelength) float64 16B 760.0 850.0
Attributes:
    data_type_group:  unprocessed raw

[6]:

cedalion.vis.anatomy.plot_montage3D(rec["amp"], geo3d)

../../_images/examples_tutorial_7_data_augmentation_8_0.png

[7]:

# Select channels which have at least a signal-to-noise ratio of 10
snr_thresh = 10  # the SNR (std/mean) of a channel.
snr, snr_mask = quality.snr(rec["amp"], snr_thresh)
amp_selected, masked_channels = xrutils.apply_mask(
    rec["amp"], snr_mask, "drop", "channel"
)

print(f"Removed {len(masked_channels)} channels because of low SNR.")

Removed 22 channels because of low SNR.

[8]:

# Calculate optical density
od = cedalion.nirs.cw.int2od(amp_selected)

Construct headmodel

We load the the Colin27 headmodel, since we need the geometry for image reconstruction.

[9]:

head_ijk = dot.get_standard_headmodel("colin27")

# transform coordinates to a RAS coordinate system
head_ras = head_ijk.apply_transform(head_ijk.t_ijk2ras)

[10]:

display(head_ijk.brain)
display(head_ras.brain)

TrimeshSurface(faces: 29988 vertices: 15002 crs: ijk units: dimensionless vertex_coords: ['parcel'])

TrimeshSurface(faces: 29988 vertices: 15002 crs: mni units: millimeter vertex_coords: ['parcel'])

[11]:

head_ras.landmarks

[11]:

<xarray.DataArray (label: 73, mni: 3)> Size: 2kB
[mm] 0.0045 -118.6 -23.08 -86.08 -19.99 ... -100.4 37.02 49.81 -99.34 21.47
Coordinates:
  * label    (label) <U3 876B 'Iz' 'LPA' 'RPA' 'Nz' ... 'PO1' 'PO2' 'PO4' 'PO6'
    type     (label) object 584B PointType.LANDMARK ... PointType.LANDMARK
Dimensions without coordinates: mni

We want to build the synthetic HRFs at C3 and C4 (green dots in the image below)

[12]:

plt_pv = pv.Plotter()
vbx.plot_surface(plt_pv, head_ras.brain, color="#d3a6a1")
vbx.plot_surface(plt_pv, head_ras.scalp, opacity=0.1)
vbx.plot_labeled_points(
    plt_pv, head_ras.landmarks.sel(label=["C3", "C4"]), show_labels=True
)
vbx.camera_at_cog(plt_pv, head_ras.brain, (-400,500,400), fit_scene=True)

plt_pv.show()

../../_images/examples_tutorial_7_data_augmentation_17_0.png

Build spatial activation pattern on brain surface for landmarks C3 and C4

The function build_spatial_activation is used to place a spatial activation pattern on the brain surface. The activation pattern is a Gaussian function of the geodesic distance to a seed vertex. Hence, the size of the activation is determined by the standard deviation of this Gaussian, specified by the parameter spatial_scale. The peak intensity in HbO is determined by the parameter intensity_scale. The intensity of HbR activation is specified relative to the HbO peak intensity. So if the HbO pattern describes an increase in Hbo then providing a negative factor smaller than 1 yields a decrease in HbR with smaller amplitude. The seed vertex (integer) can be selected as the closest vertex to a given landmark:

[13]:

# obtain the closest vertices to C3 and C4
_, c3_seed = head_ras.brain.mesh.kdtree.query(
    head_ras.landmarks.sel(label="C3").pint.magnitude
)
_, c4_seed = head_ras.brain.mesh.kdtree.query(
    head_ras.landmarks.sel(label="C4").pint.magnitude
)

[14]:

# create the spatial activation
spatial_act = synhrf.build_spatial_activation(
    head_ras.brain,
    c3_seed,
    spatial_scale=2 * cedalion.units.cm,
    intensity_scale=1 * units.micromolar,
    hbr_scale=-0.4,
)

The resulting DataArray contains an activation value for each vertex and chromophore on the brain surface.

[15]:

display(spatial_act)

<xarray.DataArray (vertex: 15002, chromo: 2)> Size: 240kB
[M] 1.115e-28 -4.461e-29 2.003e-29 -8.013e-30 5.23e-30 ... 0.0 -0.0 0.0 -0.0
Coordinates:
  * chromo   (chromo) <U3 24B 'HbO' 'HbR'
Dimensions without coordinates: vertex

[16]:

f,ax = plt.subplots(1,2,figsize=(10,5))
cedalion.vis.anatomy.plot_brain_in_axes(
        od,
        head_ras.landmarks,
        spatial_act.sel(chromo="HbO").pint.to("uM"),
        head_ras.brain,
        ax[0],
        camera_pos="C3",
        cmap="RdBu_r",
        vmin=-1,
        vmax=+1,
        cb_label=r"$\Delta$ HbO / µM",
        title=None,
    )
ax[0].set_title("HbO")
cedalion.vis.anatomy.plot_brain_in_axes(
        od,
        head_ras.landmarks,
        spatial_act.sel(chromo="HbR").pint.to("uM"),
        head_ras.brain,
        ax[1],
        camera_pos="C3",
        cmap="RdBu_r",
        vmin=-1,
        vmax=+1,
        cb_label=r"$\Delta$ HbR / µM",
        title=None,
    )
ax[1].set_title("HbR");

../../_images/examples_tutorial_7_data_augmentation_23_0.png

The following plot illustrates the effects of the spatial_scale and intensity_scale parameters:

[17]:

f, ax = plt.subplots(2, 3, figsize=(9, 6))
for i, spatial_scale in enumerate([0.5 * units.cm, 2 * units.cm, 3 * units.cm]):
    spatial_act = synhrf.build_spatial_activation(
        head_ras.brain,
        c3_seed,
        spatial_scale=spatial_scale,
        intensity_scale=1 * units.micromolar,
        hbr_scale=-0.4,
    )

    cedalion.vis.anatomy.plot_brain_in_axes(
        od,
        head_ras.landmarks,
        spatial_act.sel(chromo="HbO").pint.to("uM"),
        head_ras.brain,
        ax[0, i],
        camera_pos="C3",
        cmap="RdBu_r",
        vmin=-1,
        vmax=+1,
        cb_label=r"$\Delta$ HbO / µM",
        title=None,
    )
    ax[0, i].set_title(f"spatial_scale: {spatial_scale.magnitude} cm")

for i, intensity_scale in enumerate(
    [
        0.5 * units.micromolar,
        1.0 * units.micromolar,
        2.0 * units.micromolar,
    ]
):
    spatial_act = synhrf.build_spatial_activation(
        head_ras.brain,
        c3_seed,
        spatial_scale=2 * units.cm,
        intensity_scale=intensity_scale,
        hbr_scale=-0.4,
    )

    cedalion.vis.anatomy.plot_brain_in_axes(
        od,
        head_ras.landmarks,
        spatial_act.sel(chromo="HbO").pint.to("uM"),
        head_ras.brain,
        ax[1, i],
        camera_pos="C3",
        cmap="RdBu_r",
        vmin=-2,
        vmax=+2,
        cb_label=r"$\Delta$ HbO / µM",
        title=None,
    )
    ax[1, i].set_title(f"intensity_scale: {intensity_scale.magnitude} µM")

f.tight_layout()

../../_images/examples_tutorial_7_data_augmentation_25_0.png

For this example notebook two activations are placed below C3 and C4:

[18]:

spatial_act_c3 = synhrf.build_spatial_activation(
    head_ras.brain,
    c3_seed,
    spatial_scale=2 * cedalion.units.cm,
    intensity_scale=1 * units.micromolar,
    hbr_scale=-0.4,
)
spatial_act_c4 = synhrf.build_spatial_activation(
    head_ras.brain,
    c4_seed,
    spatial_scale=2 * cedalion.units.cm,
    intensity_scale=1 * units.micromolar,
    hbr_scale=-0.4,
)

We concatenate the two images for C3 and C4 along dimension trial_type to get a single DataArray with the spatial information for both landmarks.

[19]:

# concatenate the two spatial activations along a new dimension
spatial_imgs = xr.concat(
    [spatial_act_c3, spatial_act_c4], dim="trial_type"
).assign_coords(trial_type=["Stim C3", "Stim C4"])
spatial_imgs

[19]:

<xarray.DataArray (trial_type: 2, vertex: 15002, chromo: 2)> Size: 480kB
[M] 1.115e-28 -4.461e-29 2.003e-29 ... -2.887e-14 1.427e-21 -5.706e-22
Coordinates:
  * chromo      (chromo) <U3 24B 'HbO' 'HbR'
  * trial_type  (trial_type) <U7 56B 'Stim C3' 'Stim C4'
Dimensions without coordinates: vertex

Plots of spatial patterns

Using the helper function cedalion.vis.anatomy.plot_brain_in_axes, the created activations on the brain surface below C3 and C4 are plotted:

[20]:

f, ax = plt.subplots(1,2, figsize=(10,5))

cedalion.vis.anatomy.plot_brain_in_axes(
    od,
    head_ras.landmarks,
    spatial_imgs.sel(trial_type="Stim C3", chromo="HbO").pint.to("uM"),
    head_ras.brain,
    ax[0],
    camera_pos="C3",
    cmap="RdBu_r",
    vmin=-1,
    vmax=+1,
    cb_label=r"$\Delta$ HbO / µM",
    title="C3 Activation",
)

cedalion.vis.anatomy.plot_brain_in_axes(
    od,
    head_ras.landmarks,
    spatial_imgs.sel(trial_type="Stim C4", chromo="HbO").pint.to("uM"),
    head_ras.brain,
    ax[1],
    camera_pos="C4",
    cmap="RdBu_r",
    vmin=-1,
    vmax=+1,
    cb_label=r"$\Delta$ HbO / µM",
    title="C4 Activation",
)

f.tight_layout()

../../_images/examples_tutorial_7_data_augmentation_31_0.png

Image Reconstruction

We load the precomputed Adot matrix to be able to map from image to channel space. (For details see S5 image reconstruction tutorial).

[21]:

Adot = cedalion.data.get_precomputed_sensitivity("nn22_resting", "colin27")

[22]:

# we only consider brain vertices, not scalp and drop pruned channels
Adot_brain = Adot.sel(vertex=Adot.is_brain, channel=od.channel)
Adot_brain

The forward model and image reconstruction transform time series in image space into channel space, and vice versa. In the present case, we need to translate between optical-density time series of different wavelengths in channel space and concentration time series of different chromophores in image space. This is facilitated by the function dot.forward_model.image_to_channel_space:

[23]:

# propagate concentration changes in the brain to optical density
# changes in channel space

spatial_chan = dot.forward_model.image_to_channel_space(
    Adot_brain, spatial_imgs, spectrum="prahl"
)
spatial_chan

[23]:

<xarray.DataArray (trial_type: 2, wavelength: 2, channel: 545)> Size: 17kB
[] -1.583e-12 -6.242e-11 -3.406e-11 -2.385e-13 ... 4.92e-06 7.374e-08 5.261e-08
Coordinates:
  * wavelength  (wavelength) float64 16B 760.0 850.0
  * channel     (channel) object 4kB 'S10D112' 'S10D133' ... 'S9D94' 'S9D96'
  * trial_type  (trial_type) <U7 56B 'Stim C3' 'Stim C4'
    source      (channel) object 4kB 'S10' 'S10' 'S10' 'S10' ... 'S9' 'S9' 'S9'
    detector    (channel) object 4kB 'D112' 'D133' 'D135' ... 'D92' 'D94' 'D96'

Show the spatial activation in channel space with a scalp plot.

[24]:

fig, ax = plt.subplots(1, 2)
# adjust plot size
fig.set_size_inches(12, 6)
cedalion.vis.anatomy.scalp_plot(
    od,
    rec.geo3d,
    spatial_chan.sel(trial_type="Stim C3", wavelength=850),
    ax[0],
    cmap="YlOrRd",
    title="850nm, activation under C3",
    vmin=spatial_chan.min().item(),
    vmax=spatial_chan.max().item(),
    cb_label="max peak amplitude",
)
cedalion.vis.anatomy.scalp_plot(
    od,
    rec.geo3d,
    spatial_chan.sel(trial_type="Stim C4", wavelength=850),
    ax[1],
    cmap="YlOrRd",
    title="850nm, activation under C4",
    vmin=spatial_chan.min().item(),
    vmax=spatial_chan.max().item(),
    cb_label="Max peak amplitude",
)
plt.show()

../../_images/examples_tutorial_7_data_augmentation_38_0.png

Get top 5 channels for each trial type where synthetic activation is highest

[25]:

roi_chans_c3 = spatial_chan.channel[
    spatial_chan.sel(trial_type="Stim C3").max("wavelength").argsort()[-5:].values
].values
roi_chans_c4 = spatial_chan.channel[
    spatial_chan.sel(trial_type="Stim C4").max("wavelength").argsort()[-5:].values
].values

Concentration Scale

The activations were simulataed in image space with a peak concentration change of 1 µM in one vertex. The change in optical density in one channel reflects concentration changes in the ensemble of vertices that this channel is sensitive to.

When applying the Beer-Lambert-transformation in channel space, a change in concentration is calculated for each channel. However, the scales of these concentration changes and the concentration changes in single vertices are not the same.

Here, a correction factor is calculated to scale the activation in channel space to 1uM.

[26]:

dpf = xr.DataArray(
    [6, 6],
    dims="wavelength",
    coords={"wavelength": rec["amp"].wavelength},
)

# add time axis with one time point so we can convert to conc
spatial_chan_w_time = spatial_chan.expand_dims("time")
spatial_chan_w_time = spatial_chan_w_time.assign_coords(time=[0])
spatial_chan_w_time.time.attrs["units"] = "second"
display(spatial_chan_w_time)
spatial_chan_conc = cedalion.nirs.cw.od2conc(
    spatial_chan_w_time, geo3d, dpf, spectrum="prahl"
)

<xarray.DataArray (time: 1, trial_type: 2, wavelength: 2, channel: 545)> Size: 17kB
[] -1.583e-12 -6.242e-11 -3.406e-11 -2.385e-13 ... 4.92e-06 7.374e-08 5.261e-08
Coordinates:
  * wavelength  (wavelength) float64 16B 760.0 850.0
  * channel     (channel) object 4kB 'S10D112' 'S10D133' ... 'S9D94' 'S9D96'
  * trial_type  (trial_type) <U7 56B 'Stim C3' 'Stim C4'
    source      (channel) object 4kB 'S10' 'S10' 'S10' 'S10' ... 'S9' 'S9' 'S9'
    detector    (channel) object 4kB 'D112' 'D133' 'D135' ... 'D92' 'D94' 'D96'
  * time        (time) int64 8B 0

[27]:

# rescale so that synthetic hrfs add 1 micromolar at peak.
rescale_factor = (1* units.micromolar / spatial_chan_conc.max()).item()

display(rescale_factor)
spatial_chan *= rescale_factor

9.25210059207612 dimensionless

HRFs in channel space

So far the notebook focused on the spatial extent of the activation. To build the temporal HRF model we use the same functionality that generates hrf regressors for the GLM.

First we select a basis function, which defines the temporal shape of the HRF.

[28]:

basis_fct = glm.Gamma(tau=0 * units.s, sigma=3 * units.s, T=3 * units.s)

[29]:

od.time

A Stim DataFrame, which contains the onset, duration and amplitude of the synthetic HRFs, is created.

[30]:

stim_df = synhrf.build_stim_df(
    max_time=od.time.values[-1] * units.seconds,
    trial_types=["Stim C3", "Stim C4"],
    min_interval=10 * units.seconds,
    max_interval=20 * units.seconds,
    min_stim_dur = 10 * units.seconds,
    max_stim_dur = 10 * units.seconds,
    min_stim_value = 1.0,
    max_stim_value = 1.0,
    order="random",
)

[31]:

stim_df.head()

[31]:

	onset	duration	value	trial_type
0	14.09	10.0	1.0	Stim C4
1	41.90	10.0	1.0	Stim C4
2	69.49	10.0	1.0	Stim C4
3	89.99	10.0	1.0	Stim C4
4	114.77	10.0	1.0	Stim C3

We can now use our stim dataframe, basis function, and spatial information to create the synthetic HRF timeseries

[32]:

syn_ts = synhrf.build_synthetic_hrf_timeseries(od, stim_df, basis_fct, spatial_chan)

We get a synthetic HRF timeseries for each channel, trial_type and chromo / wavelength

[33]:

syn_ts

[33]:

<xarray.DataArray (trial_type: 2, wavelength: 2, channel: 545, time: 3309)> Size: 58MB
[] -0.0 -0.0 -0.0 -0.0 -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:
  * wavelength  (wavelength) float64 16B 760.0 850.0
  * channel     (channel) object 4kB 'S10D112' 'S10D133' ... 'S9D94' 'S9D96'
  * trial_type  (trial_type) <U7 56B 'Stim C3' 'Stim C4'
    source      (channel) object 4kB 'S10' 'S10' 'S10' 'S10' ... 'S9' 'S9' 'S9'
    detector    (channel) object 4kB 'D112' 'D133' 'D135' ... 'D92' 'D94' 'D96'
  * time        (time) float64 26kB 0.0 0.1112 0.2225 ... 367.8 367.9 368.0

Sum the synthetic timeseries over trial_type dimension, so it has the same shape as the resting state data

[34]:

syn_ts_sum = syn_ts.sum(dim='trial_type')
syn_ts_sum

[34]:

<xarray.DataArray (wavelength: 2, channel: 545, time: 3309)> Size: 29MB
[] 0.0 0.0 0.0 0.0 0.0 0.0 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:
  * wavelength  (wavelength) float64 16B 760.0 850.0
  * channel     (channel) object 4kB 'S10D112' 'S10D133' ... 'S9D94' 'S9D96'
    source      (channel) object 4kB 'S10' 'S10' 'S10' 'S10' ... 'S9' 'S9' 'S9'
    detector    (channel) object 4kB 'D112' 'D133' 'D135' ... 'D92' 'D94' 'D96'
  * time        (time) float64 26kB 0.0 0.1112 0.2225 ... 367.8 367.9 368.0

Adding HRFs to measured data

Here, the simulated activations are combined with physiological noise by adding the synthetic HRFs to the resting state dataset:

[35]:

# set to false to process only the synth. HRF in the following
ADD_RESTING_STATE = True

if ADD_RESTING_STATE:
    od_w_hrf = od + syn_ts_sum
else:
    od_w_hrf = syn_ts_sum
    od_w_hrf = od_w_hrf.assign_coords(
        {
            "samples": ("time", od.samples.values),
            "time": ("time", od.time.data),
        }
    )
    od_w_hrf = od_w_hrf.pint.quantify({"time": "s"})

Recover the HRFs again

In the following, the added activations should be extracted from the simulated dataset again. To this end, the data is frequency filtered and block averages are calculated.

[36]:

od_w_hrf_filtered = od_w_hrf.cd.freq_filter(fmin=0.02, fmax=0.5, butter_order=4)

[37]:

od_w_hrf_filtered

[37]:

<xarray.DataArray (channel: 545, wavelength: 2, time: 3309)> Size: 29MB
[] 0.002073 0.002634 0.00317 0.003663 ... 0.006999 0.006033 0.004973 0.00387
Coordinates:
  * time        (time) float64 26kB 0.0 0.1112 0.2225 ... 367.8 367.9 368.0
  * channel     (channel) object 4kB 'S10D87' 'S10D94' ... 'S47D29' 'S5D137'
  * wavelength  (wavelength) float64 16B 760.0 850.0
    samples     (time) int64 26kB 0 1 2 3 4 5 ... 3303 3304 3305 3306 3307 3308
    source      (channel) object 4kB 'S10' 'S10' 'S10' ... 'S56' 'S47' 'S5'
    detector    (channel) object 4kB 'D87' 'D94' 'D89' ... 'D129' 'D29' 'D137'

[38]:

epochs = od_w_hrf_filtered.cd.to_epochs(
    stim_df,  # stimulus dataframe
    ["Stim C3", "Stim C4"],  # select events
    before=5 * units.seconds,  # seconds before stimulus
    after=20 * units.seconds,  # 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")

n_roi = roi_chans_c3.size
# show results
f, ax = plt.subplots(2, n_roi, figsize=(16, 8))
ax = ax.flatten()
for i_ch, ch in enumerate(roi_chans_c3):
    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)
    ax[i_ch].set_ylim(-0.05, 0.05)

for i_ch, ch in enumerate(roi_chans_c4):
    for ls, trial_type in zip(["-", "--"], blockaverage.trial_type):
        ax[i_ch + n_roi].plot(
            blockaverage.reltime,
            blockaverage.sel(wavelength=760, trial_type=trial_type, channel=ch),
            "r",
            lw=2,
            ls=ls,
        )
        ax[i_ch + n_roi].plot(
            blockaverage.reltime,
            blockaverage.sel(wavelength=850, trial_type=trial_type, channel=ch),
            "b",
            lw=2,
            ls=ls,
        )
    ax[i_ch + n_roi].grid(1)
    ax[i_ch + n_roi].set_title(ch)
    ax[i_ch + n_roi].set_ylim(-0.05, 0.05)

plt.suptitle(
    "Blockaverage for channels most sensitive to C3 (top) and C4 (bottom): 760nm: r | 850nm: b | C3: - | C4: --"
)
plt.tight_layout()
plt.show()

../../_images/examples_tutorial_7_data_augmentation_62_0.png

Map block average back to brain surface

We map our extracted block averages back to the brain surface to visualize the recovered HRFs activation for Stim C3. We can compare it to the synthetic HRF image we created earlier.

[39]:

sbf = dot.GaussianSpatialBasisFunctions(head_ras, Adot, **dot.SBF_GAUSSIANS_DENSE, verbose=False)

[40]:

recon = dot.ImageRecon(
    Adot,
    recon_mode="mua2conc",
    alpha_meas=0.01,
    brain_only=True,
    spatial_basis_functions=sbf
)

blockaverage_img = recon.reconstruct(blockaverage)
blockaverage_img

[40]:

<xarray.DataArray (chromo: 2, vertex: 15002, trial_type: 2, reltime: 226)> Size: 108MB
[µM] -0.0886 -0.07553 -0.06249 -0.04967 -0.03719 ... 0.0 0.0 0.0 0.0 0.0
Coordinates:
  * chromo      (chromo) <U3 24B 'HbO' 'HbR'
  * reltime     (reltime) float64 2kB -4.995 -4.884 -4.773 ... 19.76 19.87 19.98
  * trial_type  (trial_type) object 16B 'Stim C3' 'Stim C4'
    is_brain    (vertex) bool 15kB True True True True ... True True True True
    parcel      (vertex) object 120kB 'VisCent_ExStr_8_LH' ... 'Background+Fr...
Dimensions without coordinates: vertex

[41]:

# plot HbO time trace of left and right brain hemisphere during FTapping/Right

for view in ["left_hemi", "right_hemi"]:
    trial_type = "Stim C3"
    gif_fname = "Ftapping-right" + "_HbO_" + view + ".gif"

    hbo = blockaverage_img.sel(chromo="HbO", trial_type=trial_type).pint.dequantify()
    hbo_brain = hbo.transpose("vertex", "reltime")

    ntimes = hbo.sizes["reltime"]

    b = cdc.VTKSurface.from_trimeshsurface(head_ras.brain)
    b = pv.wrap(b.mesh)
    b["reco_hbo"] = hbo_brain[:, 0] - hbo_brain[:, 0]

    p = pv.Plotter()

    p.add_mesh(
        b,
        scalars="reco_hbo",
        cmap="seismic",  # 'gist_earth_r',
        clim=(-2.5, 2.5),
        scalar_bar_args={"title": "HbO / µM"},
        smooth_shading=True,
    )

    tl = lambda tt: f"{trial_type} HbO rel. time: {tt:.3f} s"
    time_label = p.add_text(tl(0))

    cog = head_ras.brain.vertices.mean("label")
    cog = cog.pint.dequantify().values
    if view == "left_hemi":
        p.camera.position = cog + [-400, 0, 0]
    else:
        p.camera.position = cog + [400, 0, 0]
    p.camera.focal_point = cog
    p.camera.up = [0, 0, 1]
    p.reset_camera()

    p.open_gif(gif_fname)

    for i in range(0, ntimes, 3):
        b["reco_hbo"] = hbo_brain[:, i] - hbo_brain[:, 0]
        time_label.set_text("upper_left", tl(hbo_brain.reltime[i]))

        p.write_frame()

    p.close()

[42]:

display_image("Ftapping-right_HbO_left_hemi.gif")
display_image("Ftapping-right_HbO_right_hemi.gif")

../../_images/examples_tutorial_7_data_augmentation_68_0.png

../../_images/examples_tutorial_7_data_augmentation_68_1.png

[43]:

for c, channel in [("C3", roi_chans_c3[-1]), ("C4", roi_chans_c4[-1])]:
    trial_type = f"Stim {c}"

    f,ax = plt.subplots(2,4, figsize=(14,6))
    vmin,vmax = -2.5,2.5

    cedalion.vis.anatomy.plot_brain_in_axes(
        od,
        head_ras.landmarks,
        spatial_imgs.sel(trial_type=trial_type, chromo="HbO").pint.to("uM"),
        head_ras.brain,
        ax[0,0],
        camera_pos=c,
        cmap="RdBu_r",
        vmin=vmin,
        vmax=vmax,
        cb_label=r"$\Delta$ HbO / µM",
        #title=f"{c} Activation",
    )
    ax[0,0].set_title(f"true activation at {c}")

    trels = [0., 10., 15.]
    for i_wl, wl in enumerate(blockaverage.wavelength.values):
        tmp = blockaverage.sel(
            trial_type=trial_type, channel=channel, wavelength=wl
        )
        ax[1,0].plot(tmp.reltime, tmp, label=f"{wl} nm")
    for trel in trels:
        ax[1,0].axvline(trel, c="k", ls=":")
    ax[1,0].legend(loc="upper left")
    ax[1,0].set_xlabel("$t_{rel}$")
    ax[1,0].set_ylabel("OD")
    ax[1,0].set_title(channel)

    for i_col, trel in enumerate(trels, start=1):
        for i_row, cam_pos in enumerate(["C3", "C4"]):
            cedalion.vis.anatomy.plot_brain_in_axes(
                od,
                head_ras.landmarks,
                blockaverage_img.sel(chromo="HbO", trial_type=trial_type)
                .sel(reltime=trel, method="nearest"),
                head_ras.brain,
                ax[i_row, i_col],
                camera_pos=cam_pos,
                cmap="RdBu_r",
                vmin=vmin,
                vmax=vmax,
                cb_label=r"$\Delta$ HbO / µM",
                #title="$t_{rel}=" + f"{trel:.1f} s$",
            )
            ax[i_row, i_col].set_title(f"view at {cam_pos} | " + "$t_{rel}=" + f"{trel:.1f} s$")

../../_images/examples_tutorial_7_data_augmentation_69_0.png

../../_images/examples_tutorial_7_data_augmentation_69_1.png

Adding Motion Artifacts

The second part of this notebook demonstrates how to add spike and baseline-shift artifacts to existing time series.

First, we’ll load some example data.

[44]:

rec = cedalion.data.get_fingertappingDOT()
rec["od"] = cedalion.nirs.cw.int2od(rec["amp"])

f, ax = plt.subplots(1, 1, figsize=(12, 4))
ax.plot(
    rec["amp"].time,
    rec["amp"].sel(channel="S1D2", wavelength="760"),
    "r-",
    label="760nm",
)
ax.plot(
    rec["amp"].time,
    rec["amp"].sel(channel="S1D2", wavelength="850"),
    "g-",
    label="850nm",
)

plt.legend()
ax.set_xlabel("time / s")
ax.set_ylabel("intensity / V")

display(rec["od"])

<xarray.DataArray (channel: 100, wavelength: 2, time: 8794)> Size: 14MB
[] 0.02263 0.02322 0.01368 0.005969 0.01692 ... 0.01558 0.02343 0.02274 0.02761
Coordinates:
  * time        (time) float64 70kB 0.0 0.2294 0.4588 ... 2.017e+03 2.017e+03
    samples     (time) int64 70kB 0 1 2 3 4 5 ... 8788 8789 8790 8791 8792 8793
  * 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

../../_images/examples_tutorial_7_data_augmentation_72_1.png

Artifact Generation

Artifacts are generated by the functions gen_bl_shift and gen_spike.

These function take as arguments:

the time axis of time series
onset time
duration

The amplitude of the generated artifacts is set to 1.

[45]:

time = rec["amp"].time

sample_bl_shift = synart.gen_bl_shift(time, 1000) # baseline shift a t = 1000s
sample_spike = synart.gen_spike(time, 1500, 3) # spike artifact at t=1500s for 3s

display(sample_bl_shift)

fig, ax =  plt.subplots(1, 1, figsize=(12,2))
ax.plot(time, sample_bl_shift, "r-", label="bl_shift")
ax.plot(time, sample_spike, "g:", label="spike")
ax.set_xlabel('Time / s')
ax.set_ylabel('Amp')
ax.legend()

plt.tight_layout()
plt.show()

<xarray.DataArray 'time' (time: 8794)> Size: 70kB
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
Coordinates:
  * time     (time) float64 70kB 0.0 0.2294 0.4588 ... 2.017e+03 2.017e+03

../../_images/examples_tutorial_7_data_augmentation_74_1.png

Controlling Artifact Timing

Artifacts can be placed using a timing data frame with columns ‘onset’, ‘duration’, ‘trial_type’, ‘value’, and ‘channel’, i.e. it resembles the stimulus data frame with the added ‘channel’ column.

We can use the function add_event_timing to create such timing data frames. The function allows precise control over each event.

The function sel_chans_by_opt allows us to select a list of channels by way of a list of optode labels. This reflects the fact that motion artifacts usually stem from the motion of a specific optode or set of optodes, which in turn affects all related channels.

We can also use the functions random_events_num and random_events_perc to add random events to the data frame, specifying either the number of events or the percentage of the time series duration, respectively.

[46]:

# Create a list of events in the format (onset, duration)
events = [(1000, 1), (2000, 1)]

# Creates a new timing dataframe with the specified events.
# Setting channel to None indicates that the artifact applies to all channels.
timing_amp = synart.add_event_timing(events, 'bl_shift', None)

# Select channels by optode
chans = synart.sel_chans_by_opt(["S1"], rec["od"])

# Add random events to the timing dataframe
timing_od = synart.random_events_perc(time, 0.01, ["spike"], chans)

display(timing_amp)
display(timing_od)

	onset	duration	trial_type	value	channel
0	1000	1	bl_shift	1	None
1	2000	1	bl_shift	1	None

	onset	duration	trial_type	value	channel
0	1016.148457	0.333921	spike	1	[S1D1, S1D2, S1D4, S1D5, S1D6, S1D8]
1	911.217854	0.320027	spike	1	[S1D1, S1D2, S1D4, S1D5, S1D6, S1D8]
2	1297.378913	0.185465	spike	1	[S1D1, S1D2, S1D4, S1D5, S1D6, S1D8]
3	1519.040435	0.283268	spike	1	[S1D1, S1D2, S1D4, S1D5, S1D6, S1D8]
4	418.380688	0.255613	spike	1	[S1D1, S1D2, S1D4, S1D5, S1D6, S1D8]
...	...	...	...	...	...
71	971.715712	0.230447	spike	1	[S1D1, S1D2, S1D4, S1D5, S1D6, S1D8]
72	461.129298	0.223120	spike	1	[S1D1, S1D2, S1D4, S1D5, S1D6, S1D8]
73	909.346069	0.352874	spike	1	[S1D1, S1D2, S1D4, S1D5, S1D6, S1D8]
74	70.182194	0.369101	spike	1	[S1D1, S1D2, S1D4, S1D5, S1D6, S1D8]
75	59.235884	0.313672	spike	1	[S1D1, S1D2, S1D4, S1D5, S1D6, S1D8]

76 rows × 5 columns

Adding Artifacts to Data

The function add_artifacts scales artifacts and adds them to time-series data.

The artifact generation functions (see above) are passed as a dictionary. Keys correspond to entries in the column ‘trial_type’ of the timing data frame, i.e. each event specified in the timing data frame is generated using the function artifacts[trial_type]. If mode is ‘manual’, artifacts are scaled directly by the scale parameter, otherwise artifacts are automatically scaled by a parameter alpha which is calculated using a sliding window approach.

Finally, artifacts to be added to optical-density time series can be automatically scaled to yield specific concentration amplitudes with the function add_chromo_artifacts_2_od.

[47]:

artifacts = {"spike": synart.gen_spike, "bl_shift": synart.gen_bl_shift}

# Add baseline shifts to the amp data
rec["amp2"] = synart.add_artifacts(rec["amp"], timing_amp, artifacts)

# Convert the amp data to optical density
rec["od2"] = cedalion.nirs.cw.int2od(rec["amp2"])

dpf = xr.DataArray(
    [6, 6],
    dims="wavelength",
    coords={"wavelength": rec["amp"].wavelength},
)

# add spikes to od based on conc amplitudes
rec["od2"] = synart.add_chromo_artifacts_2_od(
    rec["od2"], timing_od, artifacts, rec.geo3d, dpf, 1.5
)

# Plot the OD data
channels = rec["od"].channel.values[0:6]
fig, axes = plt.subplots(len(channels), 1, figsize=(12, len(channels) * 2))
if len(channels) == 1:
    axes = [axes]
for i, channel in enumerate(channels):
    ax = axes[i]
    ax.plot(
        rec["od2"].time,
        rec["od2"].sel(channel=channel, wavelength="850"),
        "g-",
        label="850nm + artifacts",
    )
    ax.plot(
        rec["od"].time,
        rec["od"].sel(channel=channel, wavelength="850"),
        "r-",
        label="850nm - od",
    )
    ax.set_title(f"Channel: {channel}")
    ax.set_xlabel("Time / s")
    ax.set_ylabel("OD")
    ax.legend()
plt.tight_layout()

../../_images/examples_tutorial_7_data_augmentation_78_0.png

[48]:

# Plot the data in conc

rec["conc"] = cedalion.nirs.cw.od2conc(rec["od"], rec.geo3d, dpf)
rec["conc2"] = cedalion.nirs.cw.od2conc(rec["od2"], rec.geo3d, dpf)
channels = rec["od"].channel.values[0:6]
fig, axes = plt.subplots(len(channels), 1, figsize=(12, len(channels) * 2))
if len(channels) == 1:
    axes = [axes]
for i, channel in enumerate(channels):
    ax = axes[i]
    ax.plot(
        rec["conc2"].time,
        rec["conc2"].sel(channel=channel, chromo="HbR"),
        "g-",
        label="HbR + artifacts",
    )
    ax.plot(
        rec["conc"].time,
        rec["conc"].sel(channel=channel, chromo="HbR"),
        "b-",
        label="HbR",
    )
    ax.set_title(f"Channel: {channel}")
    ax.set_xlabel("Time / s")
    ax.set_ylabel("conc")
    ax.legend()
plt.tight_layout()
plt.show()

../../_images/examples_tutorial_7_data_augmentation_79_0.png

Magnitude	[[1.1153618078602912e-28 -4.461447231441165e-29] [2.0032289232306357e-29 -8.012915692922544e-30] [5.2295070523093204e-30 -2.0918028209237284e-30] ... [0.0 -0.0] [0.0 -0.0] [0.0 -0.0]]
Units	molar