cedalion.dot.forward_model

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.

Functions

`apply_inv_sensitivity`(od, inv_sens)	Apply the inverted sensitivity matrix to optical density data.
`image_to_channel_space`(Adot, img[, spectrum])
`stack_flat_channel`(array)
`stack_flat_vertex`(array)
`unstack_flat_channel`(array)
`unstack_flat_vertex`(array)

Classes

ForwardModel(head_model, geo3d, measurement_list)

Forward model for simulating light transport in the head.

class cedalion.dot.forward_model.ForwardModel( head_model: TwoSurfaceHeadModel, geo3d: cdt.LabeledPointCloud, measurement_list: DataFrame, )[source]

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.LabeledPointCloud): 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)[source]: Compute fluence for each channel and wavelength from photon simulation.

compute_sensitivity(fluence_all, fluence_at_optodes)[source]: Compute sensitivity matrix from fluence.

compute_fluence_mcx(fluence_fname: str | Path, **kwargs)[source]

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, )[source]

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, )[source]

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[source]

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, )[source]

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][source]

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)[source]

cedalion.dot.forward_model.unstack_flat_vertex(array: DataArray)[source]

cedalion.dot.forward_model.stack_flat_channel(array: DataArray)[source]

cedalion.dot.forward_model.unstack_flat_channel(array: DataArray)[source]

cedalion.dot.forward_model.image_to_channel_space( Adot: DataArray, img: DataArray, spectrum: str | None = None, )[source]