Precompute 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.datasets.get_precomputed_sensitivity.
The second part of the notebook visualizes the currently available sensitivities.
[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 to true, to actually compute fluence and sensitivity files
RUN_COMPUTATION = False
[3]:
import pyvista as pv
pv.set_jupyter_backend('static')
from functools import lru_cache
import os
from pathlib import Path
import xarray as xr
import cedalion.datasets
import cedalion.imagereco.forward_model as fw
from cedalion.io.forward_model import FluenceFile, load_Adot
import cedalion.vis.plot_sensitivity_matrix
from cedalion.geometry.landmarks import LandmarksBuilder1010
xr.set_options(display_expand_data=False);
[4]:
def compute_fluence_mcx(geo3d_snapped_ijk, meas_list, head, output_file):
fwm = cedalion.imagereco.forward_model.ForwardModel(
head, geo3d_snapped_ijk, meas_list
)
fwm.compute_fluence_mcx(output_file)
def compute_sensitivity(geo3d_snapped_ijk, meas_list, head, fluence_fname, sensitivity_fname):
fwm = cedalion.imagereco.forward_model.ForwardModel(
head, geo3d_snapped_ijk, meas_list
)
fwm.compute_sensitivity(fluence_fname, sensitivity_fname)
@lru_cache
def get_headmodel(model : str):
if model == "colin27":
SEG_DATADIR, mask_files, landmarks_file = (
cedalion.datasets.get_colin27_segmentation()
)
PARCEL_FILE = cedalion.datasets.get_colin27_parcel_file()
head_ijk = 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_FILE,
brain_face_count=None,
scalp_face_count=None,
)
head_ras = head_ijk.apply_transform(head_ijk.t_ijk2ras)
builder = LandmarksBuilder1010(head_ras.scalp, head_ras.landmarks)
head_ras.landmarks = builder.build()
head_ijk = head_ras.apply_transform(head_ras.t_ras2ijk)
elif model == "icbm152":
SEG_DATADIR, mask_files, landmarks_file = (
cedalion.datasets.get_icbm152_segmentation()
)
PARCEL_FILE = cedalion.datasets.get_icbm152_parcel_file()
head_ijk = 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_FILE,
brain_face_count=None,
scalp_face_count=None,
)
head_ras = head_ijk.apply_transform(head_ijk.t_ijk2ras)
builder = LandmarksBuilder1010(head_ras.scalp, head_ras.landmarks)
head_ras.landmarks = builder.build()
head_ijk = head_ras.apply_transform(head_ras.t_ras2ijk)
else:
raise ValueError("unknown head model")
return head_ijk
def get_fnirs_dataset(dataset):
if dataset == "fingertappingDOT":
rec = cedalion.datasets.get_fingertappingDOT()
return rec.geo3d, rec._measurement_lists["amp"]
elif dataset == "fingertapping":
rec = cedalion.datasets.get_fingertapping()
return rec.geo3d, rec._measurement_lists["amp"]
elif dataset == "nn22_resting":
rec = cedalion.datasets.get_nn22_resting_state()
return rec.geo3d, rec._measurement_lists["amp"]
elif dataset == "ninja_cap_56x144":
geo3d, landmarks, meas_list = cedalion.datasets.get_ninja_cap_probe()
geo3d = xr.concat((geo3d, landmarks), dim="label")
geo3d = geo3d.pint.quantify("mm")
return geo3d, meas_list
elif dataset == "ninja_uhd_cap_164x496":
geo3d, landmarks, meas_list = cedalion.datasets.get_ninja_uhd_cap_probe()
geo3d = xr.concat((geo3d, landmarks), dim="label")
geo3d = geo3d.pint.quantify("mm")
return geo3d, meas_list
def snap(head, geo3d, dataset):
if dataset in ["fingertapping", "fingertappingDOT"]: # right handed
geo3d_snapped_ijk = head.align_and_snap_to_scalp(geo3d)
elif dataset in ["nn22_resting", "ninja_cap_56x144", "ninja_uhd_cap_164x496"]: # left handed
geo3d_snapped_ijk = head.align_and_snap_to_scalp(geo3d, mode="general")
return geo3d_snapped_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 = get_headmodel(headmodel)
geo3d, meas_list = get_fnirs_dataset(dataset)
geo3d_snapped_ijk = snap(head, geo3d, dataset)
plotter = cedalion.vis.plot_sensitivity_matrix.Main(
sensitivity=Adot,
brain_surface=head.brain,
head_surface=head.scalp,
labeled_points=geo3d_snapped_ijk,
)
plotter.plot(high_th=0, low_th=-3)
plotter.plt.show()
def compute_fluence_and_sensitivity(dataset : str, headmodel : str):
geo3d, meas_list = get_fnirs_dataset(dataset)
head = get_headmodel(headmodel)
geo3d_snapped_ijk = snap(head, geo3d, dataset)
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, fluence_fname)
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, fluence_fname, sensitivity_fname)
else:
print(f"sensitivity for {dataset} {headmodel} exists")
[5]:
if RUN_COMPUTATION:
for dataset in ["fingertapping", "fingertappingDOT", "nn22_resting", "ninja_cap_56x144", "ninja_uhd_cap_164x496"]:
for headmodel in ["colin27", "icbm152"]:
compute_fluence_and_sensitivity(dataset, headmodel)
plot_sensitivity(dataset, headmodel)
Visualize precomputed sensitivities
[6]:
for dataset in ["fingertapping", "fingertappingDOT", "nn22_resting", "ninja_cap_56x144", "ninja_uhd_cap_164x496"]:
for headmodel in ["colin27", "icbm152"]:
display(f"{dataset=} {headmodel=}")
Adot = cedalion.datasets.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/v25.1.0/sensitivity_fingertapping_colin27.nc' to '/home/runner/.cache/cedalion/v25.1.0'.
/home/runner/work/cedalion/cedalion/src/cedalion/geometry/landmarks.py:242: UserWarning: WIP: distance calculation around ears
warnings.warn("WIP: distance calculation around ears")
"dataset='fingertapping' headmodel='icbm152'"
Downloading file 'sensitivity_fingertapping_icbm152.nc' from 'https://doc.ibs.tu-berlin.de/cedalion/datasets/v25.1.0/sensitivity_fingertapping_icbm152.nc' to '/home/runner/.cache/cedalion/v25.1.0'.
/home/runner/work/cedalion/cedalion/src/cedalion/geometry/landmarks.py:242: UserWarning: WIP: distance calculation around ears
warnings.warn("WIP: distance calculation around ears")
"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/v25.1.0/sensitivity_nn22_resting_icbm152.nc' to '/home/runner/.cache/cedalion/v25.1.0'.
"dataset='ninja_cap_56x144' headmodel='colin27'"
Downloading file 'sensitivity_ninja_cap_56x144_colin27.nc' from 'https://doc.ibs.tu-berlin.de/cedalion/datasets/v25.1.0/sensitivity_ninja_cap_56x144_colin27.nc' to '/home/runner/.cache/cedalion/v25.1.0'.
"dataset='ninja_cap_56x144' headmodel='icbm152'"
Downloading file 'sensitivity_ninja_cap_56x144_icbm152.nc' from 'https://doc.ibs.tu-berlin.de/cedalion/datasets/v25.1.0/sensitivity_ninja_cap_56x144_icbm152.nc' to '/home/runner/.cache/cedalion/v25.1.0'.
"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/v25.1.0/sensitivity_ninja_uhd_cap_164x496_colin27.nc' to '/home/runner/.cache/cedalion/v25.1.0'.
"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/v25.1.0/sensitivity_ninja_uhd_cap_164x496_icbm152.nc' to '/home/runner/.cache/cedalion/v25.1.0'.
[ ]: