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, log_file: Path | None = None, )[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.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)[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 via the adjoint Monte Carlo method.

The sensitivity matrix (Jacobian) maps absorption changes in each voxel/vertex to changes in the measured optical density. It is computed using the adjoint Monte Carlo approach (Boas and Dale [BD05], Yao et al. [YIF16]): the fluence from source and detector optodes is multiplied element-wise and projected onto the head surface.

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

Stack chromo and vertex dimensions into a single flat_vertex dim.

Parameters:: array – DataArray with "chromo" and "vertex" dimensions.
Returns:: DataArray with a new stacked "flat_vertex" dimension.
Raises:: ValueError – If array is missing either the "chromo" or "vertex" dimension.

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

Unstack the flat_vertex dimension back into chromo and vertex.

Parameters:: array – DataArray with a "flat_vertex" multi-index dimension.
Returns:: DataArray with separate "chromo" and "vertex" dimensions.

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

Stack wavelength and channel dims into a single flat_channel dim.

Parameters:: array – DataArray with "wavelength" and "channel" dimensions.
Returns:: DataArray with a new stacked "flat_channel" dimension.
Raises:: ValueError – If array is missing either the "wavelength" or "channel" dimension.

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

Unstack the flat_channel dimension back into wavelength and channel.

Parameters:: array – DataArray with a "flat_channel" multi-index dimension.
Returns:: DataArray with separate "wavelength" and "channel" dimensions.

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

Project an image-space quantity into channel (measurement) space via Adot.

Performs the forward projection y = Adot @ img where the shared vertex dimension is contracted. If img is in concentration units, the chromophore-stacked sensitivity matrix is used; if it is in absorption units (1/length), the raw Adot is used directly.

Parameters:

Adot – Sensitivity matrix with a "vertex" dimension.
img – Image DataArray with a "vertex" dimension. Must be quantified in either concentration ("[concentration]") or absorption ("[1/length]") units.
spectrum – Extinction coefficient spectrum (e.g. "prahl"). Required when img has concentration units.

Returns:

DataArray in channel space with "wavelength" and "channel" dimensions.

Raises:

ValueError – If img has incompatible units or spectrum is None when concentration units are detected.

cedalion.dot.head_model module

Two-surface head model and standard atlas loading for DOT forward modelling.

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

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

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

Type:: xarray.core.dataarray.DataArray

brain[source]

cdc.Surface Surface of the brain.

Type:: cedalion.dataclasses.geometry.Surface

scalp[source]

cdc.Surface Surface of the scalp.

Type:: cedalion.dataclasses.geometry.Surface

landmarks[source]

cdt.LabeledPoints Anatomical landmarks in RAS space.

Type:: xarray.core.dataarray.DataArray

t_ijk2ras[source]

cdt.AffineTransform Affine transformation from ijk to RAS space.

Type:: xarray.core.dataarray.DataArray

t_ras2ijk[source]

cdt.AffineTransform Affine transformation from RAS to ijk space.

Type:: xarray.core.dataarray.DataArray

voxel_to_vertex_brain[source]

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

Type:: scipy.sparse._matrix.spmatrix

voxel_to_vertex_scalp[source]

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

Type:: scipy.sparse._matrix.spmatrix

crs[source]: 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)[source]: Apply a coordinate transformation to the head model.

save(foldername)[source]: Save the head model to a folder.

load(foldername)[source]: Load the head model from a folder.

align_and_snap_to_scalp(points)[source]: Align and snap optodes or points to the scalp surface.

segmentation_masks: DataArray[source]

brain: Surface[source]

scalp: Surface[source]

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

t_ijk2ras: DataArray[source]

t_ras2ijk: DataArray[source]

voxel_to_vertex_brain: spmatrix[source]

voxel_to_vertex_scalp: spmatrix[source]

classmethod from_segmentation( segmentation_dir: str, mask_files: dict[str, str] | None = None, landmarks_ras_file: str | None = None, brain_seg_types: list[str] | None = None, scalp_seg_types: list[str] | None = None, smoothing: float = 0.0, brain_face_count: int | None = 180000, scalp_face_count: int | None = 60000, fill_holes: bool = True, ) → TwoSurfaceHeadModel[source]

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] | None = None, brain_surface_file: str = None, scalp_surface_file: str = None, landmarks_ras_file: str | None = None, brain_seg_types: list[str] | None = None, scalp_seg_types: list[str] | None = None, smoothing: float = 0.0, brain_face_count: int | None = None, scalp_face_count: int | None = None, fill_holes: bool = False, coordinates_file: Path | str | None = None, voxel_to_vertex_mapping_file_brain: Path | None = None, ) → TwoSurfaceHeadModel[source]

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.
coordinates_file – path to file with additional vertex coordinates.
voxel_to_vertex_mapping_file_brain – precomputed voxel to vertex mapping

Returns:

An instance of the TwoSurfaceHeadModel class.

Return type:

TwoSurfaceHeadModel

property crs[source]: Coordinate reference system of the head model.

apply_transform( transform: cdt.AffineTransform, ) → TwoSurfaceHeadModel[source]

Apply a coordinate transformation to the head model.

Parameters:: transform – Affine transformation matrix (4x4) to be applied.
Returns:: Transformed head model.

copy(copy_segmasks=False)[source]

Create a copy of this head model.

Parameters:: copy_segmasks – if False, the new instance gets a reference to the existing segmentation masks instead of a full copy.

save(foldername: str)[source]

Save the head model to a folder.

Parameters:: foldername (str) – Folder to save the head model into.
Returns:: None

classmethod load(foldername: str)[source]

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

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

align_and_relax_to_scalp( points: cdt.LabeledPoints, channels: list[tuple[str, str]] | cdt.NDTimeSeries | pd.DataFrame, nominal_distances: dict[tuple[str, str], float] | None = None, n_iter: int = 400, k_spring: float = 1.0, k_anchor: float = 10.0, step_size: float = 0.1, convergence_tol: float = 0.01, initial_align_mode: str = 'general', ) → tuple[cdt.LabeledPoints, 'SpringICPResult'][source]

Align and project optodes onto the scalp using spring-relaxation ICP.

Unlike align_and_snap_to_scalp(), which projects every optode independently to its nearest scalp point, this method preserves the probe geometry by coupling source–detector pairs through Hooke’s-law springs. The algorithm alternates between a spring-force step — attracting each pair toward its nominal channel distance — and an ICP projection step that keeps every optode on the scalp surface. Anatomical landmarks shared between the probe and the head model are enforced as strong anchor springs. The approach is based on the spring-relaxation registration introduced in AtlasViewer (Aasted et al. [AYucelC+15]).

Two-phase registration:

Initial alignment — a global transform (translation + rotation + optional scale) fitted to matched landmark pairs brings points into the scalp coordinate system.
Spring-relaxation ICP — optodes are iteratively refined: channel springs resist geometry distortion; anchor springs enforce anatomical positions; every optode is projected back onto the scalp mesh after each force step.

Parameters:

points – Probe coordinates (sources, detectors and landmarks).
channels – Channel definitions — accepted forms are a NDTimeSeries DataArray (coords source and detector), a measurement-list pandas.DataFrame (columns source, detector, wavelength), or a list of (source_label, detector_label) tuples.
nominal_distances – Optional pre-specified rest lengths for each channel spring, keyed by (src_label, det_label). When None, distances are measured from the initially aligned positions.
n_iter – Maximum number of spring-relaxation iterations.
k_spring – Spring constant for channel springs (Hooke’s law).
k_anchor – Spring constant for landmark-anchor springs. Should exceed k_spring to enforce anatomical constraints strongly.
step_size – Fraction of the net force vector applied per iteration. Reduce if the relaxation is numerically unstable.
convergence_tol – Early-stop threshold: halt when the maximum surface-projection displacement across all optodes falls below this value (in scalp units, typically mm).
initial_align_mode – Global alignment method applied before relaxation. One of "general" (12-DOF full affine), "trans_rot_isoscale" (7-DOF), "trans_rot" (6-DOF), or "identity" (unit conversion only).

Returns:

(registered_points, details) where registered_points is a LabeledPoints DataArray with the final optode positions on the scalp, and details is a SpringICPResult containing convergence diagnostics and per-channel quality metrics.

Return type:

tuple

References

Paper & Code: Aasted et al. [AYucelC+15]

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, 'lambda_R_conc': 1e-06}[source]: Optimal set of regularization parameters according to an optimization study for a ball squeezing dataset. Carlton et al. [CAK+26].

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')>}[source]: Optimal set of Gaussian SBF parameters according to an optimization study for a ball squeezing dataset. Carlton et al. [CAK+26].

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')>}[source]: A sparse set of gaussians SBFs.

cedalion.dot.image_recon.estimate_alpha_meas(C_meas, K=0.01)[source]

Implements a heuristic for choosing alpha_meas.

The strength of the regularization is determined by the relative scale of the C and R regularization matrices, which is encoded in the ratio:

K = median(eig( lambda_meas C )) / max(eig( lambda_R A R A’))

In past analyses K was about 0.01 .. 0.1. With this a heuristic for choosing alpha_meas can be formed:

alpha_meas = K / median(eig(C_meas))

Parameters:

C_meas – diagonal values of C_meas matrix
K – relative scale of C and R regularization matrices

class cedalion.dot.image_recon.SpatialBasisFunctions[source]

Bases: ABC

Base for SBF implementations.

abstract property H: DataArray[source]: The sensitivity in kernel space.

abstractmethod kernel_to_image_space_conc( conc_img: DataArray, ) → DataArray[source]: Transform an image from kernel to image space in recon. mode ‘conc’.

abstractmethod kernel_to_image_space_mua( mua_img: DataArray, ) → DataArray[source]: Transform an image from kernel to image space in recon. mode ‘mua’.

abstractmethod to_file(fname: Path | str)[source]

Serialize prepared spatial basis functions into a file.

Parameters:: fname – path of the output file.

abstractmethod classmethod from_file( fname: Path | str, ) → SpatialBasisFunctions[source]

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

Bases: object

kernel_to_image_space_mua( X: ndarray, ) → ndarray[source]

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

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

Serialize prepared Gaussian spatial basis functions to HDF5 file.

Parameters:: fname – path of the output file.

classmethod from_file( fname: Path | str, ) → OriginalGaussianSpatialBasisFunctions[source]

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, verbose: bool = True, )[source]

Bases: SpatialBasisFunctions

Gaussian Spatial Basis Functions.

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
verbose – controls visibility of status messages and progress bar

property H[source]: The sensitivity in kernel space.

kernel_to_image_space_mua( X: ndarray, ) → ndarray[source]

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

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

Serialize prepared Gaussian spatial basis functions to HDF5 file.

Parameters:: fname – path of the output file.

classmethod from_file( fname: Path | str, ) → GaussianSpatialBasisFunctions[source]

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

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.
lambda_R_conc – regularization parameter that sets the expected magnitude of the image covariance.
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')][source]

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

Compute image noise/variance estimates.

Parameters:: c_meas – Measurement covariance matrix.
Returns:: Image noise estimates.
Return type:: xr.DataArray

get_image_noise_posterior( c_meas: DataArray | None = None, )[source]

Compute posterior variance of reconstructed images.

Calculates the diagonal of the posterior covariance matrix: Cov(X|y) = R - R * A^T @ (F + λ*C)^(-1) @ A * R where R is the spatial prior. Returns only the diagonal (variance at each vertex).

Parameters:: c_meas – Measurement covariance matrix.
Returns:: Posterior variance of reconstructed images.
Return type:: xr.DataArray

compute_lambda_R_indirect()[source]

Compute wavelength-specific prior scaling parameter for indirect recon.

Scales lambda_R to ensure consistency between direct (chromophore space) and indirect (wavelength space) methods. Uses extinction coefficients to relate chromophore regularization strength to OD regularization.

Returns:: xr.DataArray Wavelength-specific parameter with dimension (wavelength,). Scaled to match direct method’s effective regularization strength.
Return type:: lambda_R_indirect

cedalion.dot.tissue_properties module

Tissue properties for light transport simulation.

class cedalion.dot.tissue_properties.TissueType(value)[source]

Bases: Enum

Canonical tissue-type labels used to look up optical properties.

SKIN = 1[source]

SKULL = 2[source]

DM = 3[source]

CSF = 4[source]

GM = 5[source]

WM = 6[source]

OTHER = 7[source]

cedalion.dot.tissue_properties.get_tissue_properties( segmentation_masks: DataArray, wavelengths: list, ) → ndarray[source]

Assemble a tissue-property array for Monte Carlo light transport simulation.

For each tissue type present in segmentation_masks the absorption, scattering, anisotropy, and refraction coefficients are looked up from the module-level dictionaries and stored in the output array. Index 0 is reserved for the background (vacuum).

Parameters:

segmentation_masks – xr.DataArray with dimension "segmentation_type" whose integer values encode tissue identity.
wavelengths – List of wavelengths for which properties are required. Currently the properties are wavelength-independent (FIXME).

Returns:

NumPy array of shape (n_tissues + 1, 4, n_wavelengths) where axis 1 encodes [absorption, scattering, anisotropy, refraction].

Raises:

ValueError – If a segmentation type string is not in TISSUE_LABELS.

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

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

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

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

cedalion.dot package

Submodules

cedalion.dot.forward_model module

cedalion.dot.head_model module

cedalion.dot.image_recon module

cedalion.dot.tissue_properties module

cedalion.dot.utils module

Module contents