cedalion.imagereco package

Submodules

cedalion.imagereco.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.imagereco.forward_model.TwoSurfaceHeadModel(
segmentation_masks: xr.DataArray,
brain: cdc.Surface,
scalp: cdc.Surface,
landmarks: cdt.LabeledPointCloud,
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.LabeledPointCloud 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] = {'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[source]

Constructor from binary masks as gained from segmented MRI scans.

Parameters:

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

  • mask_files (Dict[str, str]) – Dictionary mapping segmentation types to NIFTI filenames.

  • landmarks_ras_file (Optional[str]) – Filename of the landmarks in RAS space.

  • brain_seg_types (list[str]) – List of segmentation types to be included in the brain surface.

  • scalp_seg_types (list[str]) – List of segmentation types to be included in the scalp surface.

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

  • brain_face_count (Optional[int]) – Number of faces for the brain surface.

  • scalp_face_count (Optional[int]) – Number of faces for the scalp surface.

  • fill_holes (bool) – 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: Path | 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[source]

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

Parameters:

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

  • mask_files (dict[str, str]) – Dictionary mapping segmentation types to NIFTI filenames.

  • brain_surface_file (str) – Path to the brain surface.

  • scalp_surface_file (str) – Path to the scalp surface.

  • landmarks_ras_file (Optional[str]) – Filename of the landmarks in RAS space.

  • brain_seg_types (list[str]) – List of segmentation types to be included in the brain surface.

  • scalp_seg_types (list[str]) – List of segmentation types to be included in the scalp surface.

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

  • brain_face_count (Optional[int]) – Number of faces for the brain surface.

  • scalp_face_count (Optional[int]) – Number of faces for the scalp surface.

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

  • parcel_file – Path to parcel json file.

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.

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.LabeledPointCloud,
mode: str = 'trans_rot_isoscale',
) → cdt.LabeledPointCloud[source]

Align and snap optodes or points to the scalp surface.

Parameters:

  • points (cdt.LabeledPointCloud) – 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.LabeledPointCloud

snap_to_scalp_voxels(
points: cdt.LabeledPointCloud,
) → cdt.LabeledPointCloud[source]

Snap optodes or points to the closest scalp voxel.

Parameters:

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

Returns:

Points aligned and snapped to the closest scalp

voxel.

Return type:

cdt.LabeledPointCloud

class cedalion.imagereco.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

static compute_stacked_sensitivity(sensitivity: DataArray)[source]

Compute stacked HbO and HbR sensitivity matrices from fluence.

Parameters:

sensitivity (xr.DataArray) – Sensitivity matrix for each vertex and wavelength.

Returns:

Stacked sensitivity matrix for each channel and vertex.

Return type:

xr.DataArray

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:
cedalion.imagereco.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.imagereco.forward_model.stack_flat_vertex(array: DataArray)[source]
cedalion.imagereco.forward_model.unstack_flat_vertex(array: DataArray)[source]
cedalion.imagereco.forward_model.stack_flat_channel(array: DataArray)[source]
cedalion.imagereco.forward_model.unstack_flat_channel(array: DataArray)[source]

cedalion.imagereco.solver module

Solver for the image reconstruction problem.

cedalion.imagereco.solver.pseudo_inverse_stacked(
Adot: DataArray,
alpha: float = 0.01,
Cmeas: _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None = None,
alpha_spatial: float | None = None,
)[source]

Calculate the pseudo-inverse of a stacked sensitivity matrix.

Parameters:

  • Adot – Stacked matrix

  • alpha – Thikonov regularization parameter

  • Cmeas – Optional measurement regularization parameters. If specified this can be either a vector of size nchannel or a matrix of size nchannelxnchannel.

  • alpha_spatial – Optional spatial regularization parameter. Suggested default is 1e-3, or 1e-2 when spatial basis functions are used.

Returns:

Pseudo-inverse of the stacked matrix.

Return type:

xr.DataArray

cedalion.imagereco.tissue_properties module

Tissue properties for light transport simulation.

class cedalion.imagereco.tissue_properties.TissueType(value)[source]

Bases: Enum

SKIN = 1[source]
SKULL = 2[source]
DM = 3[source]
CSF = 4[source]
GM = 5[source]
WM = 6[source]
OTHER = 7[source]
cedalion.imagereco.tissue_properties.get_tissue_properties(
segmentation_masks: DataArray,
) → ndarray[source]

Return tissue properties for the given segmentation mask.

cedalion.imagereco.utils module

Utility functions for image reconstruction.

cedalion.imagereco.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.imagereco.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.imagereco.utils.create_mock_activation_below_point(
head_model: TwoSurfaceHeadModel,
point: Annotated[DataArray, DataArraySchema(dims='label', coords='label', 'label', 'type')],
time_length: Quantity,
sampling_rate: Quantity,
spatial_size: Quantity,
vmax: Quantity,
)[source]

Create a mock activation below a point.

Parameters:

  • head_model (cedalion.imagereco.forward_model.TwoSurfaceHeadModel) – The head model.

  • point (cdt.LabeledPointCloud) – The point below which to create the activation.

  • time_length (units.Quantity) – The length of the activation.

  • sampling_rate (units.Quantity) – The sampling rate.

  • spatial_size (units.Quantity) – The spatial size of the activation.

  • vmax (units.Quantity) – The maximum value of the activation.

Returns:

The activation.

Return type:

xr.DataArray

Module contents