From 03323abdf81a29bf0d4c67116e2ab7409f5b46d8 Mon Sep 17 00:00:00 2001 From: harmening Date: Tue, 10 Sep 2024 22:15:05 +0200 Subject: [PATCH 1/7] add missing voxels_from_segmentation function --- src/cedalion/geometry/segmentation.py | 87 ++++++++++++++++++++------- 1 file changed, 66 insertions(+), 21 deletions(-) diff --git a/src/cedalion/geometry/segmentation.py b/src/cedalion/geometry/segmentation.py index ef91962b..7d6b99e0 100644 --- a/src/cedalion/geometry/segmentation.py +++ b/src/cedalion/geometry/segmentation.py @@ -12,25 +12,66 @@ import cedalion.dataclasses as cdc +def voxels_from_segmentation( + segmentation_mask: xr.DataArray, + segmentation_types: List[str], + isovalue=0.9, + fill_holes_in_mask=False, +) -> cdc.Voxels: + """ Generate voxels from a segmentation mask. + + Parameters + ---------- + segmentation_mask : xr.DataArray + Segmentation mask. + segmentation_types : List[str] + List of segmentation types. + isovalue : float, optional + Isovalue for marching cubes, by default 0.9. + fill_holes_in_mask : bool, optional + Fill holes in the mask, by default False. + + Returns + ------- + cdc.Voxels + Voxels in voxel space. + """ + combined_mask = ( + segmentation_mask.sel(segmentation_type=segmentation_types) + .any("segmentation_type") + .values + ) + + if fill_holes_in_mask: + combined_mask = binary_fill_holes(combined_mask).astype(combined_mask.dtype) + + voxels = np.argwhere(combined_mask) + + return cdc.Voxels(voxels, "ijk", cedalion.units.Unit("1")) + + def surface_from_segmentation( segmentation_mask: xr.DataArray, segmentation_types: List[str], isovalue=0.9, fill_holes_in_mask=False, ) -> cdc.Surface: - """Create a surface from a segmentation mask. - - Args: - segmentation_mask (xr.DataArray): Segmentation mask with dimensions segmentation - type, i, j, k. - segmentation_types (List[str]): A list of segmentation types to include in the - surface. - isovalue (Float): The isovalue to use for the marching cubes algorithm. - fill_holes_in_mask (Bool): Whether to fill holes in the mask before creating the - surface. - - Returns: - A cedalion.Surface object. + """ Generate a surface from a segmentation mask. + Parameters + ---------- + segmentation_mask : xr.DataArray + Segmentation mask. + segmentation_types : List[str] + List of segmentation types. + isovalue : float, optional + Isovalue for marching cubes, by default 0.9. + fill_holes_in_mask : bool, optional + Fill holes in the mask, by default False. + + Returns + ------- + cdc.Surface + Surface in voxel space. """ combined_mask = ( segmentation_mask.sel(segmentation_type=segmentation_types) @@ -50,15 +91,19 @@ def surface_from_segmentation( def cell_coordinates(volume, flat: bool = False): - """Create a DataArray with the coordinates of the cells in a volume. - - Args: - volume (xr.DataArray): The volume to get the cell coordinates from. - flat (bool): Whether to flatten the coordinates. - - Returns: - xr.DataArray: A DataArray with the coordinates of the cells in the volume. + """Generate cell coordinates from a volume. + Parameters + ---------- + volume : np.ndarray + Volume. + flat : bool, optional + Flatten the coordinates, by default False. + Returns + ------- + xr.DataArray + Cell coordinates. """ + # coordinates in voxel space i = np.arange(volume.shape[0]) j = np.arange(volume.shape[1]) From 97e41cf0f794995d1ace592f73b38a4e655ee427 Mon Sep 17 00:00:00 2001 From: Qianqian Fang Date: Sat, 14 Sep 2024 05:59:11 -0400 Subject: [PATCH 2/7] [forward] add pmcxcl and user-defined input --- src/cedalion/imagereco/forward_model.py | 32 +++++++++++++++++-------- 1 file changed, 22 insertions(+), 10 deletions(-) diff --git a/src/cedalion/imagereco/forward_model.py b/src/cedalion/imagereco/forward_model.py index 006417d3..b331fd3f 100644 --- a/src/cedalion/imagereco/forward_model.py +++ b/src/cedalion/imagereco/forward_model.py @@ -647,7 +647,7 @@ def _get_unitinmm(self): length = xrutils.norm(pts_ras[1] - pts_ras[0], pts_ras.points.crs) return length.pint.magnitude.item() - def _get_fluence_from_mcx(self, i_optode: int, nphoton: int): + def _get_fluence_from_mcx(self, i_optode: int, **kwargs): """Run MCX simulation to get fluence for one optode. Args: @@ -658,8 +658,10 @@ def _get_fluence_from_mcx(self, i_optode: int, nphoton: int): np.ndarray: Fluence in each voxel. """ + kwargs.setdefault("nphoton", 1e8) + cfg = { - "nphoton": nphoton, + "nphoton": kwargs['nphoton'], "vol": self.volume, "tstart": 0, "tend": 5e-9, @@ -670,16 +672,24 @@ def _get_fluence_from_mcx(self, i_optode: int, nphoton: int): "issrcfrom0": 1, "isnormalized": 1, "outputtype": "fluence", - "seed": int(np.floor(np.random.rand() * 1e7)), - "issavedet": 1, + "issavedet": 0, "unitinmm": self.unitinmm, } - import pmcx - result = pmcx.run(cfg) + # merging default cfg with additional positional arguments + + cfg = { **cfg, **kwargs } + + # if pmcx fails, try pmcxcl + + if "cuda" in cfg and cfg["cuda"]: + import pmcx + result = pmcx.run(cfg) + else: + import pmcxcl + result = pmcxcl.run(cfg) fluence = result["flux"][:, :, :, 0] # there is only one time bin - fluence = fluence * cfg["tstep"] / result["stat"]["normalizer"] return fluence @@ -722,11 +732,13 @@ def _fluence_at_optodes(self, fluence, emitting_opt): return result - def compute_fluence_mcx(self, nphoton: int = 1e8): + def compute_fluence_mcx(self, **kwargs): """Compute fluence for each channel and wavelength using MCX package. Args: nphoton (int): Number of photons to simulate. + along with other pmcx/pmcxcl accepted input fields, + see https://pypi.org/project/pmcx/ Returns: xr.DataArray: Fluence in each voxel for each channel and wavelength. @@ -763,9 +775,9 @@ def compute_fluence_mcx(self, nphoton: int = 1e8): label = self.optode_pos.label.values[i_opt] print(f"simulating fluence for {label}. {i_opt+1} / {n_optodes}") - # run MCX + # run MCX or MCXCL # shape: [i,j,k] - fluence = self._get_fluence_from_mcx(i_opt, nphoton=nphoton) + fluence = self._get_fluence_from_mcx(i_opt, **kwargs) # FIXME shortcut: currently tissue props are wavelength independent -> copy for i_wl in range(n_wavelength): From 4a30e766aeb108c8b1ca185ffce69c087b0ed964 Mon Sep 17 00:00:00 2001 From: harmening Date: Tue, 10 Sep 2024 22:15:05 +0200 Subject: [PATCH 3/7] add missing voxels_from_segmentation function --- src/cedalion/geometry/segmentation.py | 87 ++++++++++++++++++++------- 1 file changed, 66 insertions(+), 21 deletions(-) diff --git a/src/cedalion/geometry/segmentation.py b/src/cedalion/geometry/segmentation.py index ef91962b..7d6b99e0 100644 --- a/src/cedalion/geometry/segmentation.py +++ b/src/cedalion/geometry/segmentation.py @@ -12,25 +12,66 @@ import cedalion.dataclasses as cdc +def voxels_from_segmentation( + segmentation_mask: xr.DataArray, + segmentation_types: List[str], + isovalue=0.9, + fill_holes_in_mask=False, +) -> cdc.Voxels: + """ Generate voxels from a segmentation mask. + + Parameters + ---------- + segmentation_mask : xr.DataArray + Segmentation mask. + segmentation_types : List[str] + List of segmentation types. + isovalue : float, optional + Isovalue for marching cubes, by default 0.9. + fill_holes_in_mask : bool, optional + Fill holes in the mask, by default False. + + Returns + ------- + cdc.Voxels + Voxels in voxel space. + """ + combined_mask = ( + segmentation_mask.sel(segmentation_type=segmentation_types) + .any("segmentation_type") + .values + ) + + if fill_holes_in_mask: + combined_mask = binary_fill_holes(combined_mask).astype(combined_mask.dtype) + + voxels = np.argwhere(combined_mask) + + return cdc.Voxels(voxels, "ijk", cedalion.units.Unit("1")) + + def surface_from_segmentation( segmentation_mask: xr.DataArray, segmentation_types: List[str], isovalue=0.9, fill_holes_in_mask=False, ) -> cdc.Surface: - """Create a surface from a segmentation mask. - - Args: - segmentation_mask (xr.DataArray): Segmentation mask with dimensions segmentation - type, i, j, k. - segmentation_types (List[str]): A list of segmentation types to include in the - surface. - isovalue (Float): The isovalue to use for the marching cubes algorithm. - fill_holes_in_mask (Bool): Whether to fill holes in the mask before creating the - surface. - - Returns: - A cedalion.Surface object. + """ Generate a surface from a segmentation mask. + Parameters + ---------- + segmentation_mask : xr.DataArray + Segmentation mask. + segmentation_types : List[str] + List of segmentation types. + isovalue : float, optional + Isovalue for marching cubes, by default 0.9. + fill_holes_in_mask : bool, optional + Fill holes in the mask, by default False. + + Returns + ------- + cdc.Surface + Surface in voxel space. """ combined_mask = ( segmentation_mask.sel(segmentation_type=segmentation_types) @@ -50,15 +91,19 @@ def surface_from_segmentation( def cell_coordinates(volume, flat: bool = False): - """Create a DataArray with the coordinates of the cells in a volume. - - Args: - volume (xr.DataArray): The volume to get the cell coordinates from. - flat (bool): Whether to flatten the coordinates. - - Returns: - xr.DataArray: A DataArray with the coordinates of the cells in the volume. + """Generate cell coordinates from a volume. + Parameters + ---------- + volume : np.ndarray + Volume. + flat : bool, optional + Flatten the coordinates, by default False. + Returns + ------- + xr.DataArray + Cell coordinates. """ + # coordinates in voxel space i = np.arange(volume.shape[0]) j = np.arange(volume.shape[1]) From 135c0cacae22c862bcdcf802eb46e78b979b5326 Mon Sep 17 00:00:00 2001 From: Eike Middell Date: Sat, 14 Sep 2024 16:52:32 +0200 Subject: [PATCH 4/7] fix typo --- src/cedalion/io/forward_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cedalion/io/forward_model.py b/src/cedalion/io/forward_model.py index f33fa578..408dae91 100644 --- a/src/cedalion/io/forward_model.py +++ b/src/cedalion/io/forward_model.py @@ -44,7 +44,7 @@ def load_Adot(fn: str): def save_fluence(fn : str, fluence_all, fluence_at_optodes): """Save forward model computation results. - This method uses a lossy compressions algorithm to reduce file size. + This method uses a lossy compression algorithm to reduce file size. """ with h5py.File(fn, "w") as f: From cb085942d31a382b058b637cd771cfcfebc02b0b Mon Sep 17 00:00:00 2001 From: Qianqian Fang Date: Sat, 14 Sep 2024 11:15:10 -0400 Subject: [PATCH 5/7] [deploy] add missing dependencies --- environment_dev.yml | 4 ++-- environment_doc.yml | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/environment_dev.yml b/environment_dev.yml index c2d64425..5e8ac284 100644 --- a/environment_dev.yml +++ b/environment_dev.yml @@ -54,6 +54,6 @@ dependencies: - pywavefront==1.3 - setuptools-scm - snirf==0.8 - - pmcx==0.3.3 + - pmcx==0.3.4 + - pmcxcl==0.2.0 - sphinxcontrib-bibtex - diff --git a/environment_doc.yml b/environment_doc.yml index 2cb5dbd7..9d3e129a 100644 --- a/environment_doc.yml +++ b/environment_doc.yml @@ -55,4 +55,5 @@ dependencies: - setuptools-scm==7.1.0 - snirf==0.7.4 - pmcx + - pmcxcl - open3d==0.16 From fc3e209e4adf05e3f2a5612d77fe9337f19f4449 Mon Sep 17 00:00:00 2001 From: harmening Date: Sat, 14 Sep 2024 16:26:55 +0100 Subject: [PATCH 6/7] add voxel class --- src/cedalion/dataclasses/geometry.py | 57 ++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) diff --git a/src/cedalion/dataclasses/geometry.py b/src/cedalion/dataclasses/geometry.py index 7728b9bc..d013827e 100644 --- a/src/cedalion/dataclasses/geometry.py +++ b/src/cedalion/dataclasses/geometry.py @@ -96,6 +96,63 @@ def snap(self, points: cdt.LabeledPointCloud): return snapped +@dataclass +class Voxels(): + """3D voxels represented by a np.array. + + Attributes: + voxels (np.ndarray): The voxels. + crs (str): The coordinate reference system of the voxels. + units (pint.Unit): The units of the voxels. + """ + voxels: np.ndarray + crs: str + units: pint.Unit + + @property + def vertices(self) -> cdt.LabeledPointCloud: + result = xr.DataArray( + self.voxels, + dims=["label", self.crs], + coords={"label": np.arange(len(self.voxels))}, + attrs={"units": self.units}, + ) + result = result.pint.quantify() + + return result + + @property + def nvertices(self) -> int: + return len(self.voxels) + + def apply_transform(self, transform: cdt.AffineTransform) -> "Voxels": + # convert to homogeneous coordinates + num, dim = self.voxels.shape + hom = np.ones((num,dim+1)) + hom[:,:3] = self.voxels + # apply transformation + hom = (transform.pint.dequantify().values.dot(hom.T)).T + # backtransformation + transformed = np.array([hom[i,:3] / hom[i,3] for i in range(hom.shape[0])]) + + new_units = self.units * transform.pint.units + new_crs = transform.dims[0] + + return Voxels(transformed, new_crs, new_units) + + def _build_kdtree(self): + self._kdtree = KDTree(self.voxels) + + def __post_init__(self): + self._kdtree = None + + @property + def kdtree(self): + if self._kdtree is None: + self._build_kdtree() + return self._kdtree + + @dataclass class TrimeshSurface(Surface): """A surface represented by a trimesh object. From 731ad5acfdb29320e772a11a5f9d02230eee45bd Mon Sep 17 00:00:00 2001 From: Eike Middell Date: Sat, 14 Sep 2024 17:37:49 +0200 Subject: [PATCH 7/7] add Voxels to cdc.__init__ --- src/cedalion/dataclasses/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/cedalion/dataclasses/__init__.py b/src/cedalion/dataclasses/__init__.py index d5863df1..5005c2fd 100644 --- a/src/cedalion/dataclasses/__init__.py +++ b/src/cedalion/dataclasses/__init__.py @@ -9,6 +9,7 @@ TrimeshSurface, VTKSurface, affine_transform_from_numpy, + Voxels, ) from .schemas import ( build_labeled_points,