Visualization Examples using Cedalion Plot Functions
This notebook will be continuously extended. We un-necessarily re-import cedalion dependencies to clarify which of these are needed for the plots in each corresponding cell. PLEASE NOTE: we are in the process of re-organizing the locations for our plotting functions. This notebook will be kept up to date with the latest release to enable an easy look up.
Importing Packages
3rd party plotting packages
Most of Cedalion’s plotting functionality is based on Matplotlib and Pyvista packages
[1]:
import pyvista as pv
#pv.set_jupyter_backend('server') # this enables interactive plots
pv.set_jupyter_backend('static') # this enables static rendering
#pv.OFF_SCREEN=True
import matplotlib.pyplot as p
from IPython.display import Image
Other packages that will be useful in this notebook
Dependencies for data processing and manipulation.
[2]:
import time
import numpy as np
import xarray as xr
Load Data for Visuatization
This cell fetches a variety example datasets from the cloud for visualization. This can take a bit of time.
[3]:
import cedalion
import cedalion.datasets
import cedalion.io
import cedalion.imagereco.forward_model as fw
# Loads a high-density finger tapping fNIRS example snirf dataset into a recording container
cedalion.units
rec = cedalion.datasets.get_fingertappingDOT()
# Loads a photogrammetric example scan
fname_scan, fname_snirf, fname_montage = cedalion.datasets.get_photogrammetry_example_scan()
pscan = cedalion.io.read_einstar_obj(fname_scan)
# Loads a precalculated example fluence profile that does not belong with this recording though
fluence_all, fluence_at_optodes = cedalion.datasets.get_precomputed_fluence("fingertappingDOT", "colin27")
# Lads a segmented MRI Scan (here the Colin27 average brain) and creates a TwoSurfaceHeadModel
SEG_DATADIR, mask_files, landmarks_file = cedalion.datasets.get_colin27_segmentation()
head = fw.TwoSurfaceHeadModel.from_segmentation(
segmentation_dir=SEG_DATADIR,
mask_files = mask_files,
landmarks_ras_file=landmarks_file
)
Plotting Time Series Using Matplotlib
We are working on a nice function that abstracts most of the work away from you. Until then you will have use standard python matplotlib functions.
[4]:
import cedalion.nirs as nirs
import cedalion.plots as plots
import matplotlib.pyplot as plt
# calculate HbO/HbR from raw data
dpf = xr.DataArray(
[6, 6],
dims="wavelength",
coords={"wavelength": rec["amp"].wavelength},
)
rec["conc"] = nirs.beer_lambert(rec["amp"], rec.geo3d, dpf, spectrum = "prahl")
# rename events
rec.stim.cd.rename_events(
{"1": "rest", "2": "Tapping/Left", "3": "Tapping/Right", "4": "Squeezing/Left", "5": "Squeezing/Right"}
)
# select which time series we work with
ts = rec["conc"]
# Thanks to the xarray DataArray structure, we can easily select the data we want to plot
# plot four channels and their stim markers
f, ax = plt.subplots(4, 1, sharex=True, figsize=(12, 6))
for i, ch in enumerate(["S1D1", "S1D2", "S7D9", "S7D11"]):
ax[i].plot(ts.time, ts.sel(channel=ch, chromo="HbO"), "r-", label="HbO")
ax[i].plot(ts.time, ts.sel(channel=ch, chromo="HbR"), "b-", label="HbR")
ax[i].set_title(f"Ch. {ch}")
# add stim markers using Cedalion's plot_stim_markers function
cedalion.plots.plot_stim_markers(ax[i], rec.stim, y=1)
ax[i].set_ylabel(r"$\Delta$ c / uM")
ax[0].legend(ncol=6)
ax[3].set_label("time / s")
ax[3].set_xlim(0,100)
plt.tight_layout()
Scalp Plot
Plots a metric, e.g. channel quality or amplitude at a time point on the scalp.
[5]:
import cedalion.plots as plots
import cedalion.sigproc.quality as quality
import matplotlib.pyplot as plt
n_channels = len(rec["amp"].channel)
# Calculate channel SNR to display as a metric in the plot
snr, snr_mask = quality.snr(rec["amp"], 3)
# the plots "metric" input needs dimension (nchannels,) so we focus on the 760nm wavelength for each channel
snr_metric = snr.sel(wavelength="760").values
# Create scalp plot showing SNR in each channel
fig, ax = plt.subplots(1,1)
plots.scalp_plot(
rec["amp"],
rec.geo3d,
snr_metric,
ax,
cmap="jet", title='760nm Channel SNR Scalp Plot', vmin=0, vmax=n_channels, cb_label="SNR")
/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)
Plotting Time Series Using an interactive GUI
run_vis() from the vis.timeseries package
[6]:
from cedalion.vis import time_series
# this calls a GUI to interactively select channels from a 2D probe.
# Input is a recording container, you can choose which time series (e.g. raw, OD, concentrations) in the container to plot
# UNCOMMENT to use the GUI. Commented out for documentation autogeneration purposes.
#time_series.run_vis(rec)
Plotting an fNIRS Montage in 3D
Using plot_montage3D()
[7]:
import cedalion.plots
# use the plot_montage3D() function. It requires a recording container and a geo3d object
cedalion.plots.plot_montage3D(rec["amp"], rec.geo3d)
Plotting a Headmodel
For instance the default Colin27 or ICBM152, using plot_surface()
[8]:
import cedalion.plots
plt = pv.Plotter()
cedalion.plots.plot_surface(plt, head.brain, color="#d3a6a1")
cedalion.plots.plot_surface(plt, head.scalp, opacity=.1)
plt.show()
Adding a montage to the headmodel
For this the montage has to be registered to the headmodel’s scalp first. Then we use plot_labeled_points
[9]:
geo3d_snapped = head.align_and_snap_to_scalp(rec.geo3d)
# now we plot the head same as before...
plt = pv.Plotter()
cedalion.plots.plot_surface(plt, head.brain, color="#d3a6a1")
cedalion.plots.plot_surface(plt, head.scalp, opacity=.1)
# but use the plot_labeled_points() function to add the snapped geo3d. The flag "show_labels" can be used to show the source, detector, and landmark names
cedalion.plots.plot_labeled_points(plt, geo3d_snapped, show_labels=True)
plt.show()
We can also easily remove the EEG landmarks for a better visual…
[10]:
import cedalion.dataclasses as cdc
# keep only points that are not of type "landmark", i.e. source and detector points
geo3d_snapped = geo3d_snapped[geo3d_snapped.type != cdc.PointType.LANDMARK]
# now we plot the head same as before...
plt = pv.Plotter()
cedalion.plots.plot_surface(plt, head.brain, color="#d3a6a1")
cedalion.plots.plot_surface(plt, head.scalp, opacity=.1)
# but use the plot_labeled_points() function to add the snapped geo3d. The flag "show_labels" can be used to show the source, detector, and landmark names
cedalion.plots.plot_labeled_points(plt, geo3d_snapped, show_labels=True)
plt.show()
Surface Plot of 3D Scans
Uses the same function as for head models, plot_surface()
[11]:
import cedalion.plots
plt = pv.Plotter()
get_landmarks = cedalion.plots.plot_surface(plt, pscan, opacity=1.0)
plt.show(interactive = True)
Interactive 3D Plot to select Landmarks
using plot_surface with “pick_landmarks = True”. Here we use a photogrammetric scan, and the landmarks are indicated by green dots. Right-clicking again on an existing landmark changes the label.
[12]:
import cedalion.plots
plt = pv.Plotter()
get_landmarks = cedalion.plots.plot_surface(plt, pscan, opacity=1.0, pick_landmarks = True)
plt.show(interactive = True)
[13]:
# for documentation purposes and to enable automatically rendered example notebooks we provide the hand-picked coordinates here too.
landmark_labels = ['Nz', 'Iz', 'Cz', 'Lpa', 'Rpa']
landmark_coordinates = [np.array([14.00420712, -7.84856869, 449.77840004]),
np.array([99.09920059, 29.72154755, 620.73876117]),
np.array([161.63815139, -48.49738938, 494.91210993]),
np.array([82.8771277, 79.79500128, 498.3338802]),
np.array([15.17214095, -60.56186128, 563.29621021])]
# uncommentif you want to see your own picked results: when you are done run get_landmarks() to get the landmarks.
# landmark_coordinates, landmark_labels = get_landmarks()
display (landmark_labels)
['Nz', 'Iz', 'Cz', 'Lpa', 'Rpa']
Plot Probe Fluence / Sensitivity on Cortex
Plot the fluence between a source-detector pair, or the accumulated sensitivity profile on the cortex
Fluence between two optodes
[14]:
import cedalion.plots
# pull fluence values from the corresponding source and detector pair
f = fluence_all.loc["S12", 760].values * fluence_all.loc["D19",760].values
f[f<=0] = f[f>0].min()
f = np.log10(f)
vf = pv.wrap(f)
plt = pv.Plotter()
# plot fluence values
plt.add_volume(
vf,
log_scale=False,
cmap='plasma_r',
clim=(-10,0),
)
# add head model
cedalion.plots.plot_surface(plt, head.brain, color="#d3a6a1")
cedalion.plots.plot_surface(plt, head.scalp, opacity=.1)
cedalion.plots.plot_labeled_points(plt, geo3d_snapped, show_labels=True)
plt.show()
Plot Sensitivity Profile on Cortex
Using the calculated fluence. This will be simplified in the future to make it easier for you.
[15]:
from cedalion.vis import plot_sensitivity_matrix
import cedalion.imagereco.forward_model as fwm
# to plot sensitivity on the cortex we need a forward model
fwm = cedalion.imagereco.forward_model.ForwardModel(head, geo3d_snapped, rec._measurement_lists["amp"])
Adot = fwm.compute_sensitivity(fluence_all, fluence_at_optodes)
plotter = plot_sensitivity_matrix.Main(
sensitivity=Adot,
brain_surface=head.brain,
head_surface=head.scalp,
labeled_points=geo3d_snapped,
)
plotter.plot(high_th=0, low_th=-3)
plotter.plt.show()
Plot ImageRecon HRF Activation (from HD fNIRS/DOT) on Cortex or Scalp
Since this requires a lot of preprocessing, please use the Image Reconstruction Jupyter Example Notebook HERE. There you will be able to generate the following types of plots yourself. Here, we load saved examples generated with the corresponding function calls:
[16]:
from cedalion.plots import scalp_plot_gif
# plot activity in channels on a scalp map over time
""" 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'
) """
[16]:
" scalp_plot_gif( \n data_ts, \n geo3d, \n filename = filename_scalp, \n time_range=(-5,30,0.5)*units.s,\n scl=(-0.01, 0.01), \n fps=6, \n optode_size=6, \n optode_labels=True, \n str_title='850 nm' \n ) "
[17]:
from cedalion.plots import image_recon_view
# plot reconstructed activity (here HbO) in a 3D view over time on the brain surface
""" 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)
) """
[17]:
" image_recon_view(\n X_ts, # time series data; can be 2D (static) or 3D (dynamic)\n head,\n cmap='seismic',\n clim=clim,\n view_type='hbo_brain',\n view_position='left',\n title_str='HbO',\n filename=filename_view,\n SAVE=True,\n time_range=(-5,30,0.5)*units.s,\n fps=6,\n geo3d_plot = geo3d_plot,\n wdw_size = (1024, 768)\n) "
[18]:
# plot reconstructed activity (here HbO) in a 3D view at any single time point on the brain surface
"""
# selects the nearest time sample at t=10s in X_ts
X_ts = X_ts.sel(time=10*units.s, method="nearest")
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)
) """
[18]:
' \n # selects the nearest time sample at t=10s in X_ts\n X_ts = X_ts.sel(time=10*units.s, method="nearest") \n\n image_recon_view(\n X_ts, # time series data; can be 2D (static) or 3D (dynamic)\n head,\n cmap=\'seismic\',\n clim=clim,\n view_type=\'hbo_brain\',\n view_position=\'left\',\n title_str=\'HbO\',\n filename=filename_view,\n SAVE=True,\n time_range=(-5,30,0.5)*units.s,\n fps=6,\n geo3d_plot = geo3d_plot,\n wdw_size = (1024, 768)\n) '
[19]:
from cedalion.plots import image_recon_multi_view
# plot reconstructed activity (here HbO) in from all 3D views over time on the brain surface
""" 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)
) """
[19]:
" image_recon_multi_view(\n X_ts, # time series data; can be 2D (static) or 3D (dynamic)\n head,\n cmap='seismic',\n clim=clim,\n view_type='hbo_brain',\n title_str='HbO',\n filename=filename_multiview,\n SAVE=True,\n time_range=(-5,30,0.5)*units.s,\n fps=6,\n geo3d_plot = None # geo3d_plot,\n wdw_size = (1024, 768)\n) "
[20]:
from cedalion.plots import image_recon_multi_view
# plot reconstructed activity (here HbO) in from all 3D views over time on the scalp surface
""" 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,
SAVE=True,
time_range=(-5,30,0.5)*units.s,
fps=6,
geo3d_plot = geo3d_plot,
wdw_size = (1024, 768)
) """
[20]:
" image_recon_multi_view(\n X_ts, # time series data; can be 2D (static) or 3D (dynamic)\n head,\n cmap='seismic',\n clim=clim,\n view_type='hbo_scalp',\n title_str='HbO',\n filename=filename_multiview,\n SAVE=True,\n time_range=(-5,30,0.5)*units.s,\n fps=6,\n geo3d_plot = geo3d_plot,\n wdw_size = (1024, 768)\n) "
