Precomputed forward model results

We provide precomputed fluence and sensitivity files for the example datasets. These are created by this notebook and can be obtained through cedalion.data.get_precomputed_sensitivity.

The second part of the notebook visualizes the currently available sensitivities.

# 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

# set to true, to actually compute fluence and sensitivity files
RUN_COMPUTATION = False

import pyvista as pv

pv.set_jupyter_backend("static")

import os
from functools import lru_cache
from pathlib import Path

import numpy as np
import xarray as xr
import pandas as pd

import cedalion.data
import cedalion.dataclasses as cdc
import cedalion.dot as dot
import cedalion.vis.blocks as vbx
from cedalion.geometry.landmarks import normalize_landmarks_labels
from cedalion.geometry.registration import SpringICPResult, register_optodes_spring_icp
from cedalion.io.forward_model import FluenceFile, load_Adot
from cedalion.vis.anatomy import sensitivity_matrix

xr.set_options(display_expand_data=False);

def compute_fluence_mcx(geo3d_snapped_ijk, meas_list, head, output_file, log_file):
    fwm = dot.ForwardModel(
        head, geo3d_snapped_ijk, meas_list, log_file
    )

    fwm.compute_fluence_mcx(output_file)


def compute_sensitivity(geo3d_snapped_ijk, meas_list, head, fluence_fname, sensitivity_fname, log_file):
    fwm = dot.ForwardModel(
        head, geo3d_snapped_ijk, meas_list, log_file
    )

    fwm.compute_sensitivity(fluence_fname, sensitivity_fname)


def get_fnirs_dataset(dataset):
    if dataset == "fingertappingDOT":
        rec = cedalion.data.get_fingertappingDOT()
        return rec.geo3d, rec._measurement_lists["amp"], rec["amp"]

    elif dataset == "fingertapping":
        rec = cedalion.data.get_fingertapping()
        return rec.geo3d, rec._measurement_lists["amp"], rec["amp"]

    elif dataset == "nn22_resting":
        rec = cedalion.data.get_nn22_resting_state()
        return rec.geo3d, rec._measurement_lists["amp"], rec["amp"]

    elif dataset == "ninja_cap_56x144":
        geo3d, landmarks, meas_list = cedalion.data.get_ninja_cap_probe()
        geo3d = xr.concat((geo3d, landmarks), dim="label")
        geo3d = geo3d.pint.quantify("mm")
        ts = cedalion.dataclasses.empty_timeseries_from_measurement_list(meas_list)
        return geo3d, meas_list, ts

    elif dataset == "ninja_uhd_cap_164x496":
        geo3d, landmarks, meas_list = cedalion.data.get_ninja_uhd_cap_probe()
        geo3d = xr.concat((geo3d, landmarks), dim="label")
        geo3d = geo3d.pint.quantify("mm")
        ts = cedalion.dataclasses.empty_timeseries_from_measurement_list(meas_list)
        return geo3d, meas_list, ts

    elif dataset == "lumo_testdataset":
        rec = cedalion.data.get_lumo_testdataset()
        geo3d = normalize_landmarks_labels(rec.geo3d)

        ch_mask = (
            cedalion.nirs.channel_distances(rec["amp"], geo3d) < 3.5 * cedalion.units.cm
        )
        rec["amp"] = rec["amp"].sel(channel=ch_mask)
        meas_list = rec._measurement_lists["amp"]

        meas_list = pd.DataFrame(
            [
                r
                for _, r in meas_list.iterrows()
                if f"{r.source}{r.detector}" in rec["amp"].channel.values
            ]
        )

        return geo3d, meas_list, rec["amp"]

    elif dataset == "artinis_testdataset":
        rec = cedalion.data.get_artinis_testdataset()
        geo3d = normalize_landmarks_labels(rec.geo3d)
        geo3d = geo3d.pint.dequantify().pint.quantify("cm").pint.to("mm")

        # dataset lacks landmarks but is in MNI305. Adopt Colin landmarks
        head_ijk = cedalion.dot.get_standard_headmodel("colin27")
        head_ras = head_ijk.apply_transform(head_ijk.t_ijk2ras)
        colin_landmarks = head_ras.landmarks.sel(label=["Nz", "Iz", "LPA", "RPA", "Cz"])
        colin_landmarks = colin_landmarks.points.set_crs("pos")

        geo3d = xr.concat(
            (geo3d, colin_landmarks),
            dim="label",
        )

        return geo3d, rec._measurement_lists["amp"], rec["amp"]

    elif dataset == "kernel_testdataset":
        rec = cedalion.data.get_kernel_testdataset()
        geo3d = normalize_landmarks_labels(rec.geo3d)
        return geo3d, rec._measurement_lists["amp"], rec["amp"]

def snap(head_ras, geo3d, meas_list):
    common_landmarks = geo3d.points.common_labels(head_ras.landmarks)

    if len(common_landmarks) > 3:
        initial_align_mode = "general"
    elif len(common_landmarks) == 3:
        initial_align_mode = "trans_rot_isoscale"
    else:
        initial_align_mode = "identity"

    # align and relax to scalp - needs channel information. here provided through NDTimeSeries
    geo3d_relaxed, details = head_ras.align_and_relax_to_scalp(
        geo3d, meas_list, initial_align_mode=initial_align_mode
    )

    geo3d_relaxed_ijk = geo3d_relaxed.points.apply_transform(head_ras.t_ras2ijk)

    return geo3d_relaxed_ijk


def plot_sensitivity(dataset, headmodel, Adot=None):
    sensitivity_fname = f"sensitivity_{dataset}_{headmodel}.nc"

    if Adot is None:
        Adot = load_Adot(sensitivity_fname)

    head_ijk = cedalion.dot.get_standard_headmodel(headmodel)
    head_ras = head_ijk.apply_transform(head_ijk.t_ijk2ras)

    geo3d, meas_list, ts = get_fnirs_dataset(dataset)
    geo3d_snapped_ijk = snap(head_ras, geo3d, meas_list)

    metric = Adot.sum("channel").isel(wavelength=0).sel(vertex=Adot.is_brain)  # mean / sum
    metric = np.log10(np.clip(metric, 1e-3,10))

    plotter = sensitivity_matrix.Main(
        sensitivity=Adot,
        brain_surface=head_ijk.brain,
        head_surface=head_ijk.scalp,
        labeled_points=geo3d_snapped_ijk,
        wavelength=Adot.wavelength[0]
    )
    plotter.plot(high_th=0, low_th=-3)
    plotter.plt.show()

    plot_brain_views_grid(head_ijk.brain, metric, window_size=(1000,600), reset_camera=False)


def compute_fluence_and_sensitivity(dataset : str, headmodel : str, log_file : Path | None):
    head_ijk = cedalion.dot.get_standard_headmodel(headmodel)
    head_ras = head_ijk.apply_transform(head_ijk.t_ijk2ras)

    geo3d, meas_list, ts = get_fnirs_dataset(dataset)

    geo3d_snapped_ijk = snap(head_ras, geo3d, meas_list)

    fluence_fname = f"fluence_{dataset}_{headmodel}.h5"
    sensitivity_fname = f"sensitivity_{dataset}_{headmodel}.nc"

    if not Path(fluence_fname).exists():
        print(f"computing fluence for {dataset} {headmodel}")
        compute_fluence_mcx(geo3d_snapped_ijk, meas_list, head_ijk, fluence_fname, log_file)
    else:
        print(f"fluence for {dataset} {headmodel} exists")
    if not Path(sensitivity_fname).exists():
        print(f"computing sensitivity for {dataset} {headmodel}")
        compute_sensitivity(geo3d_snapped_ijk, meas_list, head_ijk, fluence_fname, sensitivity_fname, log_file)
    else:
        print(f"sensitivity for {dataset} {headmodel} exists")

def camera_for_view(center, view, distance=350):
    cameras = {
        "superior": ([0, 0, distance], [0, 1, 0]),
        "left": ([-distance, 0, 0], [0, 0, 1]),
        "right": ([distance, 0, 0], [0, 0, 1]),
        "anterior": ([0, distance, 0], [0, 0, 1]),
        "posterior": ([0, -distance, 0], [0, 0, 1]),
    }
    position_offset, up = cameras[view]
    return center + np.asarray(position_offset), center, up

def plot_brain_views_grid(brain_surface, vertex_colors, window_size=(1000,600), reset_camera=False, **kwargs):
    brain_center = brain_surface.vertices.pint.dequantify().mean("label").values
    plt = pv.Plotter(
        shape=(2, 6),
        groups=(
            (0, slice(0, 2)), (0, slice(2, 4)), (0, slice(4, 6)),
                     (1, slice(0, 3)), (1, slice(3, 6)),
        ),
        window_size=window_size,
    )
    views=("superior", "anterior", "posterior", "left", "right")
    positions = [(0, 0), (0, 2), (0, 4), (1, 0), (1, 3)]
    for view, subplot in zip(views, positions):
        plt.subplot(*subplot)
        vbx.plot_surface(
            plt,
            brain_surface,
            color=vertex_colors,
            **kwargs
        )
        plt.add_text(view, font_size=10)
        plt.camera.position, plt.camera.focal_point, plt.camera.up = camera_for_view(brain_center, view)
        if reset_camera:
            plt.reset_camera()
    plt.subplot(1, 2)
    plt.add_text("", font_size=10)
    plt.show()

Compute sensitivities

LOGFILE = "mcx_output.log"

DATASETS = [
    "fingertapping",
    "fingertappingDOT",
    "nn22_resting",
    "ninja_cap_56x144",
    "ninja_uhd_cap_164x496",
    "lumo_testdataset",
    "kernel_testdataset",
    "artinis_testdataset",
]

HEADMODELS = ["colin27", "icbm152"]


if RUN_COMPUTATION:
    for dataset in DATASETS:
        for headmodel in HEADMODELS:
            compute_fluence_and_sensitivity(dataset, headmodel, LOGFILE)
            plot_sensitivity(dataset, headmodel)

Visualize precomputed sensitivities

for dataset in DATASETS:
    for headmodel in HEADMODELS:
        if headmodel == "icbm152" and dataset in [
            "lumo_testdataset",
            "kernel_testdataset",
            "artinis_testdataset",
        ]:
            continue

        display(f"{dataset=} {headmodel=}")
        Adot = cedalion.data.get_precomputed_sensitivity(dataset, headmodel)
        plot_sensitivity(dataset, headmodel, Adot)

"dataset='fingertapping' headmodel='colin27'"

Downloading file 'sensitivity_fingertapping_colin27.nc' from 'https://doc.ibs.tu-berlin.de/cedalion/datasets/dev/sensitivity_fingertapping_colin27.nc' to '/home/runner/.cache/cedalion/dev'.

"dataset='fingertapping' headmodel='icbm152'"

Downloading file 'sensitivity_fingertapping_icbm152.nc' from 'https://doc.ibs.tu-berlin.de/cedalion/datasets/dev/sensitivity_fingertapping_icbm152.nc' to '/home/runner/.cache/cedalion/dev'.

"dataset='fingertappingDOT' headmodel='colin27'"

"dataset='fingertappingDOT' headmodel='icbm152'"

"dataset='nn22_resting' headmodel='colin27'"

"dataset='nn22_resting' headmodel='icbm152'"

Downloading file 'sensitivity_nn22_resting_icbm152.nc' from 'https://doc.ibs.tu-berlin.de/cedalion/datasets/dev/sensitivity_nn22_resting_icbm152.nc' to '/home/runner/.cache/cedalion/dev'.

"dataset='ninja_cap_56x144' headmodel='colin27'"

Downloading file 'sensitivity_ninja_cap_56x144_colin27.nc' from 'https://doc.ibs.tu-berlin.de/cedalion/datasets/dev/sensitivity_ninja_cap_56x144_colin27.nc' to '/home/runner/.cache/cedalion/dev'.

"dataset='ninja_cap_56x144' headmodel='icbm152'"

Downloading file 'sensitivity_ninja_cap_56x144_icbm152.nc' from 'https://doc.ibs.tu-berlin.de/cedalion/datasets/dev/sensitivity_ninja_cap_56x144_icbm152.nc' to '/home/runner/.cache/cedalion/dev'.

"dataset='ninja_uhd_cap_164x496' headmodel='colin27'"

Downloading file 'sensitivity_ninja_uhd_cap_164x496_colin27.nc' from 'https://doc.ibs.tu-berlin.de/cedalion/datasets/dev/sensitivity_ninja_uhd_cap_164x496_colin27.nc' to '/home/runner/.cache/cedalion/dev'.

"dataset='ninja_uhd_cap_164x496' headmodel='icbm152'"

Downloading file 'sensitivity_ninja_uhd_cap_164x496_icbm152.nc' from 'https://doc.ibs.tu-berlin.de/cedalion/datasets/dev/sensitivity_ninja_uhd_cap_164x496_icbm152.nc' to '/home/runner/.cache/cedalion/dev'.

"dataset='lumo_testdataset' headmodel='colin27'"

"dataset='kernel_testdataset' headmodel='colin27'"

"dataset='artinis_testdataset' headmodel='colin27'"

References

cedalion.bib.dump_to_notebook()

Methods used

[1]	Fang2009	cedalion.data.get_precomputed_sensitivity	Qianqian Fang and David A Boas. Monte carlo simulation of photon migration in 3d turbid media accelerated by graphics processing units. Optics express, 17(22):20178–20190, 2009. doi:10.1364/OE.17.020178.
[2]	Yu2018	cedalion.data.get_precomputed_sensitivity	Leiming Yu, Fanny Nina-Paravecino, David Kaeli, and Qianqian Fang. Scalable and massively parallel monte carlo photon transport simulations for heterogeneous computing platforms. Journal of biomedical optics, 23(1):010504–010504, 2018. doi:10.1117/1.JBO.23.1.010504.
[3]	Yan2020	cedalion.data.get_precomputed_sensitivity	Shijie Yan and Qianqian Fang. Hybrid mesh and voxel based monte carlo algorithm for accurate and efficient photon transport modeling in complex bio-tissues. Biomedical Optics Express, 11(11):6262–6270, 2020. doi:10.1364/BOE.409468.
[4]	Holmes1998	cedalion.data.get_colin27_headmodel_files	Colin J. Holmes, Rick Hoge, Louis Collins, Roger Woods, Arthur W. Toga, and Alan C. Evans. Enhancement of mr images using registration for signal averaging. Journal of Computer Assisted Tomography, 22(2):324–333, March 1998. doi:10.1097/00004728-199803000-00032.
[5]	Fischl2012	cedalion.data.get_colin27_headmodel_files, cedalion.data.get_icbm152_headmodel_files	Bruce Fischl. FreeSurfer. NeuroImage, 62(2):774–781, 2012. doi:10.1016/j.neuroimage.2012.01.021.
[6]	Schaefer2018	cedalion.data.get_colin27_headmodel_files, cedalion.data.get_icbm152_headmodel_files	Alexander Schaefer, Ru Kong, Evan M Gordon, Timothy O Laumann, Xi-Nian Zuo, Avram J Holmes, Simon B Eickhoff, and BT Thomas Yeo. Local-global parcellation of the human cerebral cortex from intrinsic functional connectivity mri. Cerebral cortex, 28(9):3095–3114, 2018. doi:10.1093/cercor/bhx179.
[7]	Luke2021	cedalion.data.get_fingertapping	Robert Luke and David McAlpine. fNIRS Finger Tapping Data in BIDS Format. September 2021. doi:10.5281/zenodo.5529797.
[8]	Tucker2022	cedalion.io.snirf.read_snirf	Stephen Tucker, Jay Dubb, Sreekanth Kura, Alexander von Lühmann, Robert Franke, Jörn M. Horschig, Samuel Powell, Robert Oostenveld, Michael Lührs, Édouard Delaire, Zahra M. Aghajan, Hanseok Yun, Meryem A. Yücel, Qianqian Fang, Theodore J. Huppert, Blaise deB. Frederick, Luca Pollonini, David A. Boas, and Robert Luke. Introduction to the shared near infrared spectroscopy format. Neurophotonics, 10(1):013507, 2022. doi:10.1117/1.NPh.10.1.013507.
[9]	Aasted2015	cedalion.dot.head_model.TwoSurfaceHeadModel.align_and_relax_to_scalp	Christopher M. Aasted, Meryem A. Yücel, Robert J. Cooper, Jay Dubb, Daisuke Tsuzuki, Lino Becerra, Martin P. Petkov, David Borsook, Ippeita Dan, and David A. Boas. Anatomical guidance for functional near-infrared spectroscopy: AtlasViewer tutorial. Neurophotonics, 2(2):020801, 2015. doi:10.1117/1.NPh.2.2.020801.
[10]	Fonov2011	cedalion.data.get_icbm152_headmodel_files	Vladimir Fonov, Alan C. Evans, Kelly Botteron, C. Robert Almli, Robert C. McKinstry, and D. Louis Collins. Unbiased average age-appropriate atlases for pediatric studies. NeuroImage, 54(1):313–327, January 2011. doi:10.1016/j.neuroimage.2010.07.033.