cedalion.dot package

Submodules

cedalion.dot.forward_model module

Forward model for simulating light transport in the head.

NOTE: Cedalion currently supports two ways to compute fluence: 1) via monte-carlo simulation using the MonteCarloXtreme (MCX) package, and 2) via the finite element method (FEM) using the NIRFASTer package. While MCX is automatically installed using pip, NIRFASTER has to be manually installed runnning <$ bash install_nirfaster.sh CPU # or GPU> from a within your cedalion root directory.

class cedalion.dot.forward_model.ForwardModel(

head_model: TwoSurfaceHeadModel,

geo3d: cdt.LabeledPoints,

measurement_list: DataFrame,

)

Bases: object

Forward model for simulating light transport in the head.

…

Args: head_model (TwoSurfaceHeadModel): Head model containing voxel projections to brain

and scalp surfaces.

optode_pos (cdt.LabeledPoints): Optode positions. optode_dir (xr.DataArray): Optode orientations (directions of light beams). tissue_properties (xr.DataArray): Tissue properties for each tissue type. volume (xr.DataArray): Voxelated head volume from segmentation masks. unitinmm (float): Unit of head model, optodes expressed in mm. measurement_list (pd.DataFrame): List of measurements of experiment with source,

detector, channel, and wavelength.

compute_fluence(nphoton)

Compute fluence for each channel and wavelength from photon simulation.

compute_sensitivity(fluence_all, fluence_at_optodes)

Compute sensitivity matrix from fluence.

compute_fluence_mcx(fluence_fname: str | Path, **kwargs)

Compute fluence for each channel and wavelength using MCX package.

Parameters:

• fluence_fname – the output hdf5 file to store the fluence

• kwargs – key-value pairs are passed to MCX’s configuration dict. For example nphoton (int) to control the number of photons to simulate. See https://pypi.org/project/pmcx for further options.

Returns:

Fluence in each voxel for each channel and wavelength.

Return type:

xr.DataArray

References

(Fang and Boas [FB09]) Qianqian Fang and David A. Boas, “Monte Carlo Simulation of Photon Migration in 3D Turbid Media Accelerated by Graphics Processing Units,” Optics Express, vol.17, issue 22, pp. 20178-20190 (2009).

(Yu et al. [YNPKF18]) Leiming Yu, Fanny Nina-Paravecino, David Kaeli, Qianqian Fang, “Scalable and massively parallel Monte Carlo photon transport simulations for heterogeneous computing platforms,” J. Biomed. Opt. 23(1), 010504 (2018).

(Yan and Fang [YF20]) Shijie Yan and Qianqian Fang* (2020), “Hybrid mesh and voxel based Monte Carlo algorithm for accurate and efficient photon transport modeling in complex bio-tissues,” Biomed. Opt. Express, 11(11) pp. 6262-6270. https://www.osapublishing.org/boe/abstract.cfm?uri=boe-11-11-6262

compute_fluence_nirfaster(

fluence_fname: str | Path,

meshingparam=None,

)

Compute fluence for each channel and wavelength using NIRFASTer package.

Parameters:

• fluence_fname – the output hdf5 file to store the fluence

• meshingparam (ff.utils.MeshingParam) – Parameters to be used by the CGAL mesher. Note: they should all be double

Returns: xr.DataArray: Fluence in each voxel for each channel and wavelength.

References

(Dehghani et al. [DEY+09]) Dehghani, Hamid, et al. “Near infrared optical tomography using NIRFAST: Algorithm for numerical model and image reconstruction.” Communications in numerical methods in engineering 25.6 (2009): 711-732.

compute_sensitivity(

fluence_fname: str | Path,

sensitivity_fname: str | Path,

)

Compute sensitivity matrix from fluence.

Parameters:

• fluence_fname – the input hdf5 file to store the fluence

• sensitivity_fname – the output netcdf file for the sensitivity

compute_stacked_sensitivity(

spectrum: str = 'prahl',

) → DataArray

Stack sensitivity matrices and incorporate extinction coefficients.

For image reconstruction the 3D sensitivity arrays must be transformed into invertible 2D matrices and extinction coefficients must be incorporated. This function accepts sensitivities both in image space (Adot, vertices) or in spatial basis function space (H, kernels). The dimensions get transformed from (channel, vertex) to (flat_channel, flat_vertex) or (channel, kernel) to (flat_channel, flat_kernel), respectively.

Parameters:

• sensitivity – Sensitivity matrix (Adot/H ) for each vertex/kernel and wavelength.

• spectrum – name of the extinction coefficient spectrum

Returns:

Stacked sensitivity matrix for each channel and vertex. shape = (flat_channel, flat_vertex)

static parcel_sensitivity(

Adot: DataArray,

chan_droplist: list = None,

dOD_thresh: float = 0.001,

minCh: int = 1,

dHbO: float = 10,

dHbR: float = -3,

)

Calculate a mask for parcels based on their effective cortex sensitivity.

Parcels are considered good, if a change in HbO and HbR [µM] in the parcel leads to an observable change of at least dOD in at least one wavelength of one channel. Sensitivities of all vertices in the parcel are summed up in the sensitivity matrix Adot. Bad channels in an actual measurement that are pruned can be considered by providing a boolean channel_mask, where False indicates bad channels that are dropped and not considered for parcel sensitivity. Requires headmodel with parcelation coordinates.

Parameters:

• Adot (channel, vertex, wavelength)) – Sensitivity matrix with parcel coordinate belonging to each vertex

• chan_droplist – list of channel names to be dropped from consideration of sensitivity (e.g. pruned channels due to bad signal quality)

• dOD_thresh – threshold for minimum dOD change in a channel that should be observed from a hemodynamic change in a parcel

• minCh – minimum number of channels per parcel that should see a change above dOD_thresh

• dHbO – change in HbO conc. in the parcel in [µM] used to calculate dOD

• dHbR – change in HbR conc. in the parcel in [µM] used to calculate dOD

Returns:

A tuple (parcel_dOD, parcel_mask), where parcel_dOD (channel, parcel, wavelength) contains the delta OD observed in a channel for each wavelength given the assumed dHb change in a parcel, and parcel_mask is a boolean DataArray with parcel coords from Adot that is true for parcels for which dOD_thresh is met.

Initial Contributors:

• Alexander von Lühmann | vonluehmann@tu-berlin.de | 2025

cedalion.dot.forward_model.apply_inv_sensitivity(

od: cdt.NDTimeSeries,

inv_sens: DataArray,

) → tuple[DataArray, DataArray]

Apply the inverted sensitivity matrix to optical density data.

Parameters:

• od – time series of optical density data

• inv_sens – the inverted sensitivity matrix

Returns:

Two DataArrays for the brain and scalp with the reconcstructed time series per vertex and chromophore.

cedalion.dot.forward_model.stack_flat_vertex(array: DataArray)

cedalion.dot.forward_model.unstack_flat_vertex(array: DataArray)

cedalion.dot.forward_model.stack_flat_channel(array: DataArray)

cedalion.dot.forward_model.unstack_flat_channel(array: DataArray)

cedalion.dot.forward_model.image_to_channel_space(

Adot: DataArray,

img: DataArray,

spectrum: str | None = None,

)

cedalion.dot.head_model module

class cedalion.dot.head_model.TwoSurfaceHeadModel(

segmentation_masks: xr.DataArray,

brain: cdc.Surface,

scalp: cdc.Surface,

landmarks: cdt.LabeledPoints,

t_ijk2ras: cdt.AffineTransform,

t_ras2ijk: cdt.AffineTransform,

voxel_to_vertex_brain: scipy.sparse.spmatrix,

voxel_to_vertex_scalp: scipy.sparse.spmatrix,

)

Bases: object

Head Model class to represent a segmented head.

Its main functions are reduced to work on voxel projections to scalp and cortex surfaces.

segmentation_masks

xr.DataArray Segmentation masks of the head for each tissue type.

Type:

xarray.core.dataarray.DataArray

brain

cdc.Surface Surface of the brain.

Type:

cedalion.dataclasses.geometry.Surface

scalp

cdc.Surface Surface of the scalp.

Type:

cedalion.dataclasses.geometry.Surface

landmarks

cdt.LabeledPoints Anatomical landmarks in RAS space.

Type:

xarray.core.dataarray.DataArray

t_ijk2ras

cdt.AffineTransform Affine transformation from ijk to RAS space.

Type:

xarray.core.dataarray.DataArray

t_ras2ijk

cdt.AffineTransform Affine transformation from RAS to ijk space.

Type:

xarray.core.dataarray.DataArray

voxel_to_vertex_brain

scipy.sparse.spmatrix Mapping from voxel to brain vertices.

Type:

scipy.sparse._matrix.spmatrix

voxel_to_vertex_scalp

scipy.sparse.spmatrix Mapping from voxel to scalp vertices.

Type:

scipy.sparse._matrix.spmatrix

crs

str Coordinate reference system of the head model.

from_segmentation(cls, segmentation_dir, mask_files, landmarks_ras_file,

brain_seg_types, scalp_seg_types, smoothing, brain_face_count, scalp_face_count): Construct instance from segmentation masks in NIfTI format.

apply_transform(transform)

Apply a coordinate transformation to the head model.

save(foldername)

Save the head model to a folder.

load(foldername)

Load the head model from a folder.

align_and_snap_to_scalp(points)

Align and snap optodes or points to the scalp surface.

segmentation_masks: DataArray

brain: Surface

scalp: Surface

landmarks: Annotated[DataArray, DataArraySchema(dims='label', coords='label', 'label', 'type')]

t_ijk2ras: DataArray

t_ras2ijk: DataArray

voxel_to_vertex_brain: spmatrix

voxel_to_vertex_scalp: spmatrix

classmethod from_segmentation(

segmentation_dir: str,

mask_files: dict[str, str] = {'csf': 'csf.nii', 'gm': 'gm.nii', 'scalp': 'scalp.nii', 'skull': 'skull.nii', 'wm': 'wm.nii'},

landmarks_ras_file: str | None = None,

brain_seg_types: list[str] = ['gm', 'wm'],

scalp_seg_types: list[str] = ['scalp'],

smoothing: float = 0.5,

brain_face_count: int | None = 180000,

scalp_face_count: int | None = 60000,

fill_holes: bool = True,

) → TwoSurfaceHeadModel

Constructor from binary masks as gained from segmented MRI scans.

Parameters:

• segmentation_dir – Folder containing the segmentation masks in NIFTI format.

• mask_files – Dictionary mapping segmentation types to NIFTI filenames.

• landmarks_ras_file – Filename of the landmarks in RAS space.

• brain_seg_types – List of segmentation types to be included in the brain surface.

• scalp_seg_types – List of segmentation types to be included in the scalp surface.

• smoothing – Smoothing factor for the brain and scalp surfaces.

• brain_face_count – Number of faces for the brain surface.

• scalp_face_count – Number of faces for the scalp surface.

• fill_holes – Whether to fill holes in the segmentation masks.

classmethod from_surfaces(

segmentation_dir: str,

mask_files: dict[str, str] = {'csf': 'csf.nii', 'gm': 'gm.nii', 'scalp': 'scalp.nii', 'skull': 'skull.nii', 'wm': 'wm.nii'},

brain_surface_file: str = None,

scalp_surface_file: str = None,

landmarks_ras_file: str | None = None,

brain_seg_types: list[str] = ['gm', 'wm'],

scalp_seg_types: list[str] = ['scalp'],

smoothing: float = 0.5,

brain_face_count: int | None = 180000,

scalp_face_count: int | None = 60000,

fill_holes: bool = False,

parcel_file: Path | str | None = None,

) → TwoSurfaceHeadModel

Constructor from seg.masks, brain and head surfaces as gained from MRI scans.

Parameters:

• segmentation_dir – Folder containing the segmentation masks in NIFTI format.

• mask_files – Dictionary mapping segmentation types to NIFTI filenames.

• brain_surface_file – Path to the brain surface.

• scalp_surface_file – Path to the scalp surface.

• landmarks_ras_file – Filename of the landmarks in RAS space.

• brain_seg_types – List of segmentation types to be included in the brain surface.

• scalp_seg_types – List of segmentation types to be included in the scalp surface.

• smoothing – Smoothing factor for the brain and scalp surfaces.

• brain_face_count – Number of faces for the brain surface.

• scalp_face_count – Number of faces for the scalp surface.

• fill_holes – Whether to fill holes in the segmentation masks.

• parcel_file – path to the json file mapping vertices to parcels.

Returns:

An instance of the TwoSurfaceHeadModel class.

Return type:

TwoSurfaceHeadModel

property crs

Coordinate reference system of the head model.

apply_transform(

transform: cdt.AffineTransform,

) → TwoSurfaceHeadModel

Apply a coordinate transformation to the head model.

Parameters:

transform – Affine transformation matrix (4x4) to be applied.

Returns:

Transformed head model.

save(foldername: str)

Save the head model to a folder.

Parameters:

foldername (str) – Folder to save the head model into.

Returns:

None

classmethod load(foldername: str)

Load the head model from a folder.

Parameters:

foldername (str) – Folder to load the head model from.

Returns:

Loaded head model.

Return type:

TwoSurfaceHeadModel

align_and_snap_to_scalp(

points: cdt.LabeledPoints,

mode: str = 'general',

) → cdt.LabeledPoints

Align and snap optodes or points to the scalp surface.

Parameters:

• points (cdt.LabeledPoints) – Points to be aligned and snapped to the scalp surface.

• mode – method to derive the affine transform. Could be either ‘trans_rot_isoscale’ or ‘general’. See cedalion.geometry.registraion for details.

Returns:

Points aligned and snapped to the scalp surface.

Return type:

cdt.LabeledPoints

snap_to_scalp_voxels(points: cdt.LabeledPoints) → cdt.LabeledPoints

Snap optodes or points to the closest scalp voxel.

Parameters:

points (cdt.LabeledPoints) – Points to be snapped to the closest scalp voxel.

Returns:

Points aligned and snapped to the closest scalp

voxel.

Return type:

cdt.LabeledPoints

scale_to_landmarks(

target_landmarks: cdt.LabeledPoints,

) → TwoSurfaceHeadModel

scale_to_headsize(

circumference: cdt.QLength,

nz_cz_iz: cdt.QLength,

lpa_cz_rpa: cdt.QLength,

) → TwoSurfaceHeadModel

cedalion.dot.head_model.get_standard_headmodel(model: str) → TwoSurfaceHeadModel

Create a TwoSurfaceHeadmodel for common atlases.

Onc created, this function caches a head model.

Parameters:

model – either colin27 or icbm152

Returns:

The loaded head model with 1010-landmarks and parcel labels assigned.

cedalion.dot.image_recon module

Solver for the image reconstruction problem.

cedalion.dot.image_recon.REG_PAPER_MUA_SBF = {'alpha_meas': 10000.0, 'alpha_spatial': 0.01, 'apply_c_meas': True}

Optimal set of regularization parameters according to an optimization study for a ball squeezing dataset. .

cedalion.dot.image_recon.SBF_GAUSSIANS_DENSE = {'mask_threshold': -2, 'sigma_brain': <Quantity(1, 'millimeter')>, 'sigma_scalp': <Quantity(5, 'millimeter')>, 'threshold_brain': <Quantity(1, 'millimeter')>, 'threshold_scalp': <Quantity(5, 'millimeter')>}

Optimal set of Gaussian SBF parameters according to an optimization study for a ball squeezing dataset. .

cedalion.dot.image_recon.SBF_GAUSSIANS_SPARSE = {'mask_threshold': -2, 'sigma_brain': <Quantity(5, 'millimeter')>, 'sigma_scalp': <Quantity(20, 'millimeter')>, 'threshold_brain': <Quantity(5, 'millimeter')>, 'threshold_scalp': <Quantity(20, 'millimeter')>}

A sparse set of gaussians SBFs.

class cedalion.dot.image_recon.SpatialBasisFunctions

Bases: ABC

Base for SBF implementations.

abstract property H: DataArray

The sensitivity in kernel space.

abstractmethod kernel_to_image_space_conc(

conc_img: DataArray,

) → DataArray

Transform an image from kernel to image space in recon. mode ‘conc’.

abstractmethod kernel_to_image_space_mua(

mua_img: DataArray,

) → DataArray

Transform an image from kernel to image space in recon. mode ‘mua’.

abstractmethod to_file(fname: Path | str)

Serialize prepared spatial basis functions into a file.

Parameters:

fname – path of the output file.

abstractmethod classmethod from_file(

fname: Path | str,

) → SpatialBasisFunctions

Load prepared spatial basis functions from a file.

Parameters:

fname – path of the file to read from.

Returns:

Loaded spatial basis functions instance.

Return type:

SpatialBasisFunctions

class cedalion.dot.image_recon.OriginalGaussianSpatialBasisFunctions(

head_model: TwoSurfaceHeadModel,

Adot: DataArray,

threshold_brain: Annotated[Quantity, '[length]'],

threshold_scalp: Annotated[Quantity, '[length]'],

sigma_brain: Annotated[Quantity, '[length]'],

sigma_scalp: Annotated[Quantity, '[length]'],

mask_threshold: float,

)

Bases: object

kernel_to_image_space_mua(

X: ndarray,

) → ndarray

Convert kernel space reconstructions to image space for mua.

Parameters:

X – Reconstruction values in kernel space. shape (kernel, …)

Returns:

Reconstruction values in image space.

Return type:

np.ndarray

kernel_to_image_space_conc(X) → ndarray

Convert kernel space reconstructions to image space for concentration.

Parameters:

X – Reconstruction values in kernel space.

Returns:

Reconstruction values in image space with HbO/HbR split.

Return type:

np.ndarray

to_file(fname: Path | str)

Serialize prepared Gaussian spatial basis functions to HDF5 file.

Parameters:

fname – path of the output file.

classmethod from_file(

fname: Path | str,

) → OriginalGaussianSpatialBasisFunctions

Load prepared Gaussian spatial basis functions from HDF5 group.

Parameters:

fname – path of the file to read from.

Returns:

Loaded instance.

Return type:

GaussianSpatialBasisFunctions

class cedalion.dot.image_recon.GaussianSpatialBasisFunctions(

head_model: TwoSurfaceHeadModel,

Adot: DataArray,

threshold_brain: Annotated[Quantity, '[length]'],

threshold_scalp: Annotated[Quantity, '[length]'],

sigma_brain: Annotated[Quantity, '[length]'],

sigma_scalp: Annotated[Quantity, '[length]'],

mask_threshold: float,

)

Bases: SpatialBasisFunctions

Gaussian Spatial Basis Functions.

Note: This implementation differs from the original one by a factor 2pi in the denominator of the gaussians.

Parameters:

• head_model – a TwoSurfaceHeadModel with brain and scalp surfaces

• Adot – the sensitivity matrix

• threshold_brain – the distance between kernel centers on the brain

• threshold_scalp – the distance between kernel centers on scalp

• sigma_brain – the width of the gaussians on the brain

• sigma_scalp – the width of the gaussians on the scalp

• mask_threshold – log10(sensitivity) threshold for vertices to be considered

property H

The sensitivity in kernel space.

kernel_to_image_space_mua(

X: ndarray,

) → ndarray

Convert kernel space reconstructions to image space for mua.

Parameters:

X – Reconstruction values in kernel space. shape (kernel, …)

Returns:

Reconstruction values in image space.

Return type:

np.ndarray

kernel_to_image_space_conc(X) → ndarray

Convert kernel space reconstructions to image space for concentration.

Parameters:

X – Reconstruction values in kernel space.

Returns:

Reconstruction values in image space with HbO/HbR split.

Return type:

np.ndarray

to_file(fname: Path | str)

Serialize prepared Gaussian spatial basis functions to HDF5 file.

Parameters:

fname – path of the output file.

classmethod from_file(

fname: Path | str,

) → GaussianSpatialBasisFunctions

Load prepared Gaussian spatial basis functions from HDF5 group.

Parameters:

fname – path of the file to read from.

Returns:

Loaded instance.

Return type:

GaussianSpatialBasisFunctions

class cedalion.dot.image_recon.ImageRecon(

Adot,

*,

alpha_meas: float = 0.001,

alpha_spatial: float | None = None,

apply_c_meas: bool = False,

recon_mode: Literal['conc', 'mua', 'mua2conc'] = 'mua',

brain_only: bool = False,

spatial_basis_functions: SpatialBasisFunctions | None = None,

)

Bases: object

Implements image reconstruction methods for diffuse optical tomography.

Parameters:

• Adot – the sensitivity matrix

• recon_mode –

  select reconstruction method

  • ’conc’: directly reconstruct hemoglobin concentrations from OD

    measurements at different wavelengths

  • ’mua’: reconstruct absorption changes from OD measurements for each

    wavelength separately.

  • ’mua2conc’: reconstruct absorption changes for each wavelength separately.

    Afterwards transform these to hemoglobin concentration changes.

• brain_only – if set to true, scalp vertices in Adot are ignored and the reconstruction is constrained to brain vertices

• alpha_meas – regularization parameter to adjust the balance between image noise and spatial resolution.

• alpha_spatial – regularization parameter that controls the effective depth of the reconstruction by controlling how strongly the vertex sensitivities are rescaled. A smaller alpha_spatial will more strongly suppress activation that is reconstructed on the scalp.

• apply_c_meas – controls whether the provided measurement covariance should be used for measurement regularization.

• spatial_basis_functions – if given reconstruct in the kernel space defined by the provided spatial-basis-function implementation. The result is returned in image space.

reconstruct(

y: Annotated[DataArray, DataArraySchema(dims='time', coords='time', 'time', 'samples')],

c_meas: DataArray | None = None,

) → Annotated[DataArray, DataArraySchema(dims='time', coords='time', 'time', 'samples')]

Reconstruct images from measurement data.

Parameters:

• y – optical density time series or time point data.

• c_meas – Diagonal elements of the measurement covariance matrix (optional). dims: wavelength x channel.

Returns:

Reconstructed images.

get_image_noise(c_meas: DataArray)

Compute image noise/variance estimates.

Parameters:

c_meas – Measurement covariance matrix.

Returns:

Image noise estimates.

Return type:

xr.DataArray

get_image_noise_tstat(

time_series: Annotated[DataArray, DataArraySchema(dims='time', coords='time', 'time', 'samples')],

c_meas: DataArray | None = None,

)

Compute t-statistic images from noise estimates.

Parameters:

• time_series – Time series data for statistics computation.

• c_meas – Measurement covariance matrix (optional).

Returns:

T-statistic images.

Return type:

xr.DataArray

cedalion.dot.tissue_properties module

Tissue properties for light transport simulation.

class cedalion.dot.tissue_properties.TissueType(value)

Bases: Enum

SKIN = 1

SKULL = 2

DM = 3

CSF = 4

GM = 5

WM = 6

OTHER = 7

cedalion.dot.tissue_properties.get_tissue_properties(

segmentation_masks: DataArray,

wavelengths: list,

) → ndarray

Return tissue properties for the given segmentation mask.

cedalion.dot.utils module

Utility functions for image reconstruction.

cedalion.dot.utils.map_segmentation_mask_to_surface(

segmentation_mask: DataArray,

transform_vox2ras: DataArray,

surface: Surface,

)

Find for each voxel the closest vertex on the surface.

Parameters:

• segmentation_mask (xr.DataArray) – A binary mask of shape (segmentation_type, i, j, k).

• transform_vox2ras (xr.DataArray) – The affine transformation from voxel to RAS space.

• surface (cedalion.dataclasses.Surface) – The surface to map the voxels to.

Returns:

A sparse matrix of shape (ncells, nvertices) that maps voxels to

cells.

Return type:

coo_array

cedalion.dot.utils.normal_hrf(t, t_peak, t_std, vmax)

Create a normal hrf.

Parameters:

• t (np.ndarray) – The time points.

• t_peak (float) – The peak time.

• t_std (float) – The standard deviation.

• vmax (float) – The maximum value of the HRF.

Returns:

The HRF.

Return type:

np.ndarray

cedalion.dot.utils.create_mock_activation_below_point(

head_model: TwoSurfaceHeadModel,

point: Annotated[DataArray, DataArraySchema(dims='label', coords='label', 'label', 'type')],

time_length: Annotated[Quantity, '[time]'],

sampling_rate: Annotated[Quantity, '[frequency]'],

spatial_size: Annotated[Quantity, '[length]'],

vmax: float,

)

Create a mock activation below a point.

Parameters:

• head_model – The head model.

• point – The point below which to create the activation.

• time_length – The length of the activation.

• sampling_rate – The sampling rate.

• spatial_size – The spatial size of the activation.

• vmax – The maximum value of the activation.

Returns:

The activation.

Return type:

xr.DataArray

Module contents