diff --git a/simpa/core/device_digital_twins/detection_geometries/curved_array.py b/simpa/core/device_digital_twins/detection_geometries/curved_array.py index c530e4f7..7aadd3fc 100644 --- a/simpa/core/device_digital_twins/detection_geometries/curved_array.py +++ b/simpa/core/device_digital_twins/detection_geometries/curved_array.py @@ -4,8 +4,7 @@ import numpy as np from simpa.core.device_digital_twins import DetectionGeometryBase -from simpa.utils import Tags - +from simpa.utils import Tags, Settings class CurvedArrayDetectionGeometry(DetectionGeometryBase): """ @@ -82,7 +81,10 @@ def check_settings_prerequisites(self, global_settings) -> bool: return False return True - def update_settings_for_use_of_model_based_volume_creator(self, global_settings): + def update_settings_for_use_of_model_based_volume_creator(self, global_settings) -> Settings: + pass + + def update_settings_for_use_of_segmentation_based_volume_creator(self, global_settings) -> Settings: pass def get_detector_element_positions_base_mm(self) -> np.ndarray: diff --git a/simpa/core/device_digital_twins/detection_geometries/linear_array.py b/simpa/core/device_digital_twins/detection_geometries/linear_array.py index 9209d5e0..83bcc94f 100644 --- a/simpa/core/device_digital_twins/detection_geometries/linear_array.py +++ b/simpa/core/device_digital_twins/detection_geometries/linear_array.py @@ -60,7 +60,10 @@ def check_settings_prerequisites(self, global_settings: Settings) -> bool: return False return True - def update_settings_for_use_of_model_based_volume_creator(self, global_settings): + def update_settings_for_use_of_model_based_volume_creator(self, global_settings) -> Settings: + pass + + def update_settings_for_use_of_segmentation_based_volume_creator(self, global_settings) -> Settings: pass def get_detector_element_positions_base_mm(self) -> np.ndarray: diff --git a/simpa/core/device_digital_twins/detection_geometries/planar_array.py b/simpa/core/device_digital_twins/detection_geometries/planar_array.py index 1fe2fb93..9cba4d55 100644 --- a/simpa/core/device_digital_twins/detection_geometries/planar_array.py +++ b/simpa/core/device_digital_twins/detection_geometries/planar_array.py @@ -78,7 +78,10 @@ def check_settings_prerequisites(self, global_settings: Settings) -> bool: return False return True - def update_settings_for_use_of_model_based_volume_creator(self, global_settings): + def update_settings_for_use_of_model_based_volume_creator(self, global_settings) -> Settings: + pass + + def update_settings_for_use_of_segmentation_based_volume_creator(self, global_settings) -> Settings: pass def get_detector_element_positions_base_mm(self) -> np.ndarray: diff --git a/simpa/core/device_digital_twins/digital_device_twin_base.py b/simpa/core/device_digital_twins/digital_device_twin_base.py index 478361a0..2095b0d7 100644 --- a/simpa/core/device_digital_twins/digital_device_twin_base.py +++ b/simpa/core/device_digital_twins/digital_device_twin_base.py @@ -9,6 +9,7 @@ import uuid from simpa.utils.serializer import SerializableSIMPAClass from simpa.utils.calculate import are_equal +from simpa.utils import Settings class DigitalDeviceTwinBase(SerializableSIMPAClass): @@ -72,7 +73,7 @@ def check_settings_prerequisites(self, global_settings) -> bool: pass @abstractmethod - def update_settings_for_use_of_model_based_volume_creator(self, global_settings): + def update_settings_for_use_of_model_based_volume_creator(self, global_settings) -> Settings: """ This method can be overwritten by a PA device if the device poses special constraints to the volume that should be considered by the model-based volume creator. @@ -82,6 +83,16 @@ def update_settings_for_use_of_model_based_volume_creator(self, global_settings) """ pass + @abstractmethod + def update_settings_for_use_of_segmentation_based_volume_creator(self, global_settings) -> Settings: + """ + This method can be overwritten by a PA device if the device poses special constraints to the + volume that should be considered by the segmentation-based volume creator. + :param global_settings: Settings for the entire simulation pipeline. + :type global_settings: Settings + """ + pass + def get_field_of_view_mm(self) -> np.ndarray: """ Returns the absolute field of view in mm where the probe position is already diff --git a/simpa/core/device_digital_twins/illumination_geometries/illumination_geometry_base.py b/simpa/core/device_digital_twins/illumination_geometries/illumination_geometry_base.py index 6c5780bf..92abe70d 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/illumination_geometry_base.py +++ b/simpa/core/device_digital_twins/illumination_geometries/illumination_geometry_base.py @@ -57,6 +57,9 @@ def check_settings_prerequisites(self, global_settings) -> bool: def update_settings_for_use_of_model_based_volume_creator(self, global_settings) -> Settings: return global_settings + def update_settings_for_use_of_segmentation_based_volume_creator(self, global_settings) -> Settings: + pass + def serialize(self) -> dict: serialized_device = self.__dict__ device_dict = {"IlluminationGeometryBase": serialized_device} diff --git a/simpa/core/device_digital_twins/pa_devices/ithera_msot_acuity.py b/simpa/core/device_digital_twins/pa_devices/ithera_msot_acuity.py index dd61dbf9..ceb282e8 100644 --- a/simpa/core/device_digital_twins/pa_devices/ithera_msot_acuity.py +++ b/simpa/core/device_digital_twins/pa_devices/ithera_msot_acuity.py @@ -14,7 +14,6 @@ from simpa.utils.libraries.tissue_library import TissueLibrary import numpy as np - class MSOTAcuityEcho(PhotoacousticDevice): """ This class represents a digital twin of the MSOT Acuity Echo, manufactured by iThera Medical, Munich, Germany @@ -218,6 +217,167 @@ def update_settings_for_use_of_model_based_volume_creator(self, global_settings: }) volume_creator_settings[Tags.STRUCTURES][Tags.BACKGROUND] = background_settings + + def update_settings_for_use_of_segmentation_based_volume_creator( + self, + global_settings: Settings, + ): + """ + Updates the volume creation settings of the segmentation based volume creator according to the size of the + device. Adds a heavy-water layer, a mediprene membrane layer, and optionally a US-gel layer (controlled via Tags.US_GEL in the volume creation settings). + Only necessary if the segmentation map does not already include these layers. + + Supports both input formats accepted by ``SegmentationBasedAdapter``: + + * 3D integer label map``[X, Y, Z]``: new z-voxels are filled with the new integer label. + * 4D float fraction map ``[C, X, Y, Z]``: existing channels are zero-padded in z and a new + channel (index == new label) is appended with fraction 1.0 in the new z-slice. + This requires labels to be dense (0...N-1) so that the channel index equals the label. + + :param global_settings: Settings for the entire simulation pipeline. + """ + try: + volume_creator_settings = Settings( + global_settings.get_volume_creation_settings() + ) + except KeyError: + self.logger.warning( + "You called the update_settings_for_use_of_segmentation_based_volume_creator method " + "even though there are no volume creation settings defined in the " + "settings dictionary." + ) + return + + segmentation_map = volume_creator_settings[Tags.INPUT_SEGMENTATION_VOLUME] + segmentation_class_mapping = volume_creator_settings[Tags.SEGMENTATION_CLASS_MAPPING] + spacing_mm = global_settings[Tags.SPACING_MM] + z_dim_position_shift_mm = 0 + + def _prepend_layer(seg, n_px, label): + """Prepend n_px voxels along z for a given label. + + 3D [X,Y,Z]: fills the new voxels with the integer label. + 4D [C,X,Y,Z]: pads existing channels with 0 in z, then appends a new channel + (index == label) with 1.0 in the new z-slice. Labels must be dense so that + label == seg.shape[0] (i.e. the new channel index matches the label value). + """ + if seg.ndim == 3: + return np.pad(seg, ((0, 0), (0, 0), (n_px, 0)), mode="constant", constant_values=label) + else: # 4D [C, X, Y, Z] + # Zero-pad all existing class channels in z (new region belongs to a new class) + padded = np.pad(seg, ((0, 0), (0, 0), (0, 0), (n_px, 0)), mode="constant", constant_values=0) + # Append new channel: fraction 1.0 in the prepended slice, 0.0 everywhere else + new_ch = np.zeros((1, seg.shape[1], seg.shape[2], padded.shape[3]), dtype=padded.dtype) + new_ch[0, :, :, :n_px] = 1.0 + return np.concatenate([padded, new_ch], axis=0) + + # Allocate new integer labels for the device layers by incrementing past the highest + # existing label, so they never collide with user-defined segmentation classes. + used_labels = set(int(l) for l in segmentation_class_mapping.keys()) + + def _next_unused_label(): + label = max(used_labels) + 1 + used_labels.add(label) + return label + + # --- US gel --- + # If specified, add layer of US gel + # Thickness is drawn from a normal distribution (like in model-based volume creation case). + if volume_creator_settings.get(Tags.US_GEL, False): + us_gel_thickness_mm = np.random.normal(0.4, 0.1) + us_gel_thickness_px = int(round(us_gel_thickness_mm / spacing_mm)) + us_gel_label = _next_unused_label() + segmentation_map = _prepend_layer(segmentation_map, us_gel_thickness_px, us_gel_label) + segmentation_class_mapping[us_gel_label] = TissueLibrary.ultrasound_gel() + z_dim_position_shift_mm += us_gel_thickness_px * spacing_mm + self.logger.debug("Added an ultrasound gel layer to the segmentation map.") + + # --- Mediprene membrane --- + # Fixed physical thickness of the membrane (self.mediprene_membrane_height_mm). + mediprene_thickness_px = int(round(self.mediprene_membrane_height_mm / spacing_mm)) + mediprene_label = _next_unused_label() + segmentation_map = _prepend_layer(segmentation_map, mediprene_thickness_px, mediprene_label) + segmentation_class_mapping[mediprene_label] = TissueLibrary.mediprene() + z_dim_position_shift_mm += mediprene_thickness_px * spacing_mm + self.logger.debug("Added a mediprene layer to the segmentation map.") + + # --- Heavy water --- + # Fills the rest of the probe with heavy water (probe height minus the membrane thickness). + heavy_water_label = _next_unused_label() + segmentation_class_mapping[heavy_water_label] = TissueLibrary.heavy_water() + heavy_water_layer_height_mm = self.probe_height_mm - self.mediprene_membrane_height_mm + heavy_water_layer_height_px = int(round(heavy_water_layer_height_mm / spacing_mm)) + segmentation_map = _prepend_layer(segmentation_map, heavy_water_layer_height_px, heavy_water_label) + z_dim_position_shift_mm += heavy_water_layer_height_px * spacing_mm + self.logger.debug( + f"Added a {heavy_water_layer_height_px * spacing_mm:.2f} mm heavy water layer to the segmentation map." + ) + + # Update global z-dimension + global_settings[Tags.DIM_VOLUME_Z_MM] += z_dim_position_shift_mm + self.logger.debug(f"Changed Tags.DIM_VOLUME_Z_MM to {global_settings[Tags.DIM_VOLUME_Z_MM]}") + + # Widen the volume along x if narrower than the probe footprint + # One extra voxel (0.5 on each side) guards against rounding errors + if global_settings[Tags.DIM_VOLUME_X_MM] < round(self.detection_geometry.probe_width_mm) + spacing_mm: + width_shift_for_structures_mm = ( + round(self.detection_geometry.probe_width_mm) + + spacing_mm + - global_settings[Tags.DIM_VOLUME_X_MM] + ) / 2 + # Split the total pixel shift symmetrically left/right. + total_shift_pixels = int(round(2 * width_shift_for_structures_mm / spacing_mm)) + left_shift_pixels = total_shift_pixels // 2 + right_shift_pixels = total_shift_pixels - left_shift_pixels + # For the segmentation map itself, x is axis 0 for 3D but axis 1 for 4D [C,X,Y,Z]. + seg_padding_width = ( + ((0, 0), (left_shift_pixels, right_shift_pixels), (0, 0), (0, 0)) + if np.ndim(segmentation_map) == 4 + else ((left_shift_pixels, right_shift_pixels), (0, 0), (0, 0)) + ) + segmentation_map = np.pad(segmentation_map, seg_padding_width, mode="edge") + global_settings[Tags.DIM_VOLUME_X_MM] = int(round(self.detection_geometry.probe_width_mm)) + spacing_mm + self.logger.debug( + f"Changed Tags.DIM_VOLUME_X_MM to {global_settings[Tags.DIM_VOLUME_X_MM]}, and expanded " + f"the segmentation map accordingly using edge padding" + ) + else: + width_shift_for_structures_mm = 0 + + # Write the modified segmentation map back into the settings. + global_settings[Tags.VOLUME_CREATION_MODEL_SETTINGS][Tags.INPUT_SEGMENTATION_VOLUME] = segmentation_map + self.logger.debug("The segmentation volume has been adjusted to fit the MSOT device") + + # Shift device, detector, and illuminator positions + # The probe layers are prepended at z=0, pushing the tissue surface downward by probe_height_mm. + # The x-shift centres the device over the (potentially widened) volume. + self.device_position_mm = np.add( + self.device_position_mm, + np.array([width_shift_for_structures_mm, 0, self.probe_height_mm]), + ) + self.detection_geometry_position_vector = np.add( + self.device_position_mm, np.array([0, 0, self.focus_in_field_of_view_mm]) + ) + detection_geometry = CurvedArrayDetectionGeometry( + pitch_mm=0.34, + radius_mm=40, + number_detector_elements=256, + detector_element_width_mm=0.24, + detector_element_length_mm=13, + center_frequency_hz=3.96e6, + bandwidth_percent=153, + sampling_frequency_mhz=40, + angular_origin_offset=np.pi, + device_position_mm=self.detection_geometry_position_vector, + field_of_view_extent_mm=self.field_of_view_extent_mm, + ) + self.set_detection_geometry(detection_geometry) + for illumination_geom in self.illumination_geometries: + illumination_geom.device_position_mm = np.add( + illumination_geom.device_position_mm, + np.array([width_shift_for_structures_mm, 0, self.probe_height_mm]), + ) + def serialize(self) -> dict: serialized_device = self.__dict__ device_dict = {"MSOTAcuityEcho": serialized_device} diff --git a/simpa/core/device_digital_twins/pa_devices/photoacoustic_device.py b/simpa/core/device_digital_twins/pa_devices/photoacoustic_device.py index 38b0d202..a2cf3add 100644 --- a/simpa/core/device_digital_twins/pa_devices/photoacoustic_device.py +++ b/simpa/core/device_digital_twins/pa_devices/photoacoustic_device.py @@ -5,6 +5,7 @@ import numpy as np from abc import ABC from simpa.core.device_digital_twins import DigitalDeviceTwinBase +from simpa.utils import Settings class PhotoacousticDevice(DigitalDeviceTwinBase, ABC): @@ -138,7 +139,10 @@ def check_settings_prerequisites(self, global_settings) -> bool: _result = False return _result - def update_settings_for_use_of_model_based_volume_creator(self, global_settings): + def update_settings_for_use_of_model_based_volume_creator(self, global_settings) -> Settings: + pass + + def update_settings_for_use_of_segmentation_based_volume_creator(self, global_settings) -> Settings: pass def serialize(self) -> dict: diff --git a/simpa/core/simulation_modules/volume_creation_module/segmentation_based_adapter.py b/simpa/core/simulation_modules/volume_creation_module/segmentation_based_adapter.py index bf2358d8..055f0553 100644 --- a/simpa/core/simulation_modules/volume_creation_module/segmentation_based_adapter.py +++ b/simpa/core/simulation_modules/volume_creation_module/segmentation_based_adapter.py @@ -2,61 +2,145 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT -from simpa.core.simulation_modules.volume_creation_module import VolumeCreationAdapterBase +from simpa.core.simulation_modules.volume_creation_module import ( + VolumeCreationAdapterBase, +) from simpa.utils import Tags import numpy as np import torch +from scipy.ndimage import uniform_filter class SegmentationBasedAdapter(VolumeCreationAdapterBase): """ - This volume creator expects a np.ndarray to be in the settings - under the Tags.INPUT_SEGMENTATION_VOLUME tag and uses this array - together with a SegmentationClass mapping which is a dict defined in - the settings under Tags.SEGMENTATION_CLASS_MAPPING. + Creates a volume based on a segmentation map. Inputs can be either a 3D integer label map (each voxel has a single class label) + or a 4D fraction map (each voxel is 0-1 fraction of each class). + 3D label maps will be one-hot encoded into 4D tensor with dimensions [C, X, Y, Z] + This creates some memory overhead but allows for easy handling of (pre-computed) fraction maps and partial-volume effects. - With this, an even greater utility is warranted. + If ``Tags.CONSIDER_PARTIAL_VOLUME`` is ``True``, boundary fractions are smoothed via scipy.ndimage.uniform_filter + Smoothing can be controlled via ``Tags.PARTIAL_VOLUME_KERNEL_SIZE`` (default 3 -> 1-voxel transition). + + Fractions are normalised to sum to 1 at each voxel, and blended according to the class properties. + For the resulting segmentation map, the class with the highest fraction at each voxel determines the label. + Oxygenation is blended only where blood is present, and set to NaN where blood volume fraction is zero. """ def create_simulation_volume(self) -> dict: volumes, x_dim_px, y_dim_px, z_dim_px = self.create_empty_volumes() wavelength = self.global_settings[Tags.WAVELENGTH] - for key in volumes.keys(): - volumes[key] = volumes[key].to('cpu') - - segmentation_volume = self.component_settings[Tags.INPUT_SEGMENTATION_VOLUME] - segmentation_classes = np.unique(segmentation_volume, return_counts=False) - x_dim_seg_px, y_dim_seg_px, z_dim_seg_px = np.shape(segmentation_volume) - - if x_dim_px != x_dim_seg_px: - raise ValueError("x_dim of volumes and segmentation must perfectly match but was {} and {}" - .format(x_dim_px, x_dim_seg_px)) - if y_dim_px != y_dim_seg_px: - raise ValueError("y_dim of volumes and segmentation must perfectly match but was {} and {}" - .format(y_dim_px, y_dim_seg_px)) - if z_dim_px != z_dim_seg_px: - raise ValueError("z_dim of volumes and segmentation must perfectly match but was {} and {}" - .format(z_dim_px, z_dim_seg_px)) + raw = self.component_settings[Tags.INPUT_SEGMENTATION_VOLUME] class_mapping = self.component_settings[Tags.SEGMENTATION_CLASS_MAPPING] + # Sort so that blending order is deterministic; lower label = lower priority. + segmentation_classes = sorted(class_mapping.keys()) + # Channel index == class label, so we need max_label + 1 channels. + num_channels = max(segmentation_classes) + 1 + + # normalise input to 4D float [C, X, Y, Z] + if np.issubdtype(np.asarray(raw).dtype, np.floating) and np.ndim(raw) == 4: + # Already 4D float [C, X, Y, Z]: use fractions directly. + # Channel index must equal label value, so the array needs at least max_label + 1 channels + if raw.shape[0] < num_channels: + raise ValueError( + f"4D fraction map has {raw.shape[0]} channels but the class mapping " + f"requires at least {num_channels} (max label = {max(segmentation_classes)})." + ) + fractions = torch.tensor(raw, dtype=torch.float, device=self.torch_device) + elif np.issubdtype(np.asarray(raw).dtype, np.integer) and np.ndim(raw) == 3: + # 3D integer label map → one-hot encode to 0/1 fractions per channel. + seg = torch.tensor(np.asarray(raw), device=self.torch_device) + fractions = torch.zeros( + num_channels, + x_dim_px, + y_dim_px, + z_dim_px, + dtype=torch.float, + device=self.torch_device, + ) + for c in segmentation_classes: + fractions[c] = (seg == c).float() + + if self.component_settings.get(Tags.CONSIDER_PARTIAL_VOLUME, False): + # apply partial volume smoothing + # (kernel_size - 1) / 2 voxels on each side of a boundary, + kernel_size = self.component_settings.get( + Tags.PARTIAL_VOLUME_KERNEL_SIZE, 3 + ) + for c in segmentation_classes: + smoothed = uniform_filter( + fractions[c].cpu().numpy(), size=kernel_size + ) + fractions[c] = torch.tensor( + smoothed, dtype=torch.float, device=self.torch_device + ) + else: + raise ValueError( + f"Unsupported input segmentation volume with dtype {raw.dtype} and shape {raw.shape}. " + "Must be either 3D integer label map or 4D float fraction map." + ) + + # normalise fractions so they sum to 1 at each voxel (Should always hold for 3D label maps) + # But might not be for 4D fraction maps, and can also be slightly off after smoothing due to boundary effects. + fraction_sum = torch.zeros( + x_dim_px, y_dim_px, z_dim_px, dtype=torch.float, device=self.torch_device + ) + for c in segmentation_classes: + fraction_sum += fractions[c] + + if (fraction_sum > 0).any(): + safe = fraction_sum.clone() + safe[safe == 0] = 1.0 + for c in segmentation_classes: + fractions[c] /= safe + + # Blend properties together. + # max_class_fractions tracks the dominant class for the segmentation label. + max_class_fractions = torch.zeros( + x_dim_px, y_dim_px, z_dim_px, dtype=torch.float, device=self.torch_device + ) for seg_class in segmentation_classes: - class_properties = class_mapping[seg_class].get_properties_for_wavelength(self.global_settings, wavelength) - for volume_key in volumes.keys(): - if isinstance(class_properties[volume_key], (int, float)) or class_properties[volume_key] == None: # scalar - assigned_prop = class_properties[volume_key] - if assigned_prop is None: - assigned_prop = torch.nan - volumes[volume_key][segmentation_volume == seg_class] = assigned_prop - elif len(torch.Tensor.size(class_properties[volume_key])) == 3: # 3D map - assigned_prop = class_properties[volume_key][torch.tensor(segmentation_volume == seg_class)] - assigned_prop[assigned_prop is None] = torch.nan - volumes[volume_key][torch.tensor(segmentation_volume == seg_class)] = assigned_prop + class_fractions = fractions[seg_class] + class_props = class_mapping[seg_class].get_properties_for_wavelength( + self.global_settings, wavelength + ) + + mask = class_fractions > 0 + if not mask.any(): + continue + + for key in volumes: + prop = class_props[key] + if prop is None: + continue + + if key == Tags.DATA_FIELD_SEGMENTATION: + # Assign the label of the class with the highest fraction at each voxel. + better = class_fractions > max_class_fractions + if better.any(): + volumes[key][better] = prop + max_class_fractions[better] = class_fractions[better] + else: - raise AssertionError("Properties need to either be a scalar or a 3D map.") + # Linear blending, property is weighted by the voxel fraction. + if isinstance(prop, (int, float)): + volumes[key][mask] += class_fractions[mask] * prop + elif isinstance(prop, torch.Tensor) and prop.dim() == 3: + volumes[key][mask] += ( + class_fractions[mask] * prop.to(self.torch_device)[mask] + ) + else: + raise ValueError( + f"Unsupported property type for '{key}': {type(prop)}" + ) + + # if no blood is present, oxygenation is undefined + bvf = volumes.get(Tags.DATA_FIELD_BLOOD_VOLUME_FRACTION) + if bvf is not None: + volumes[Tags.DATA_FIELD_OXYGENATION][bvf == 0] = torch.nan - # convert volumes back to CPU - for key in volumes.keys(): - volumes[key] = volumes[key].numpy().astype(np.float64, copy=False) + for key in volumes: + volumes[key] = volumes[key].cpu().numpy().astype(np.float64, copy=False) return volumes diff --git a/simpa/utils/tags.py b/simpa/utils/tags.py index 921132bd..e55af954 100644 --- a/simpa/utils/tags.py +++ b/simpa/utils/tags.py @@ -208,6 +208,12 @@ class Tags: only occupying a partial volume of the voxel. \n Usage: adapter versatile_volume_creation """ + PARTIAL_VOLUME_KERNEL_SIZE = ("partial_volume_kernel_size", int) + """ + Integer kernel size for the uniform_filter used to smooth partial-volume boundary fractions + in SegmentationBasedAdapter. Default: 3 (1-voxel-wide on each side of a boundary, must be odd number).\n + Usage: adapter segmentation_based_adapter + """ STRUCTURE_START_MM = ("structure_start", (list, tuple, np.ndarray)) """ diff --git a/simpa_examples/segmentation_loader.py b/simpa_examples/segmentation_loader.py index 9fa4174b..d8e4fea4 100644 --- a/simpa_examples/segmentation_loader.py +++ b/simpa_examples/segmentation_loader.py @@ -7,41 +7,32 @@ import numpy as np from skimage.data import shepp_logan_phantom from scipy.ndimage import zoom -from skimage.transform import resize -# FIXME temporary workaround for newest Intel architectures import os from argparse import ArgumentParser from simpa.utils.profiling import profile -os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE" + # TODO: Please make sure that a valid path_config.env file is located in your home directory, or that you # point to the correct file in the PathManager(). -@profile -def run_segmentation_loader(spacing: float | int = 1.0, input_spacing: float | int = 0.2, path_manager=None, - visualise: bool = True): - """ +DEVICE_CHOICES = ["RSOMExplorerP50", "MSOTAcuityEcho"] - :param spacing: The simulation spacing between voxels in mm - :param input_spacing: The input spacing between voxels in mm - :param path_manager: the path manager to be used, typically sp.PathManager - :param visualise: If VISUALIZE is set to True, the reconstruction result will be plotted - :return: a run through of the example - """ - if path_manager is None: - path_manager = sp.PathManager() +def create_rsomexplorer_sample_segmentation(spacing, input_spacing=0.2): + """ + Creates a segmentation volume from the Shepp-Logan phantom. + """ label_mask = shepp_logan_phantom() - label_mask = np.digitize(label_mask, bins=np.linspace(0.0, 1.0, 11), right=True) label_mask = label_mask[100:300, 100:300] label_mask = np.reshape(label_mask, (label_mask.shape[0], 1, label_mask.shape[1])) segmentation_volume_tiled = np.tile(label_mask, (1, 128, 1)) - segmentation_volume_mask = sp.round_x5_away_from_zero(zoom(segmentation_volume_tiled, input_spacing/spacing, - order=0)).astype(int) + segmentation_volume_mask = sp.round_x5_away_from_zero( + zoom(segmentation_volume_tiled, input_spacing / spacing, order=0) + ).astype(int) def segmentation_class_mapping(): ret_dict = dict() @@ -52,45 +43,148 @@ def segmentation_class_mapping(): ret_dict[4] = sp.TissueLibrary.mediprene() ret_dict[5] = sp.TissueLibrary.ultrasound_gel() ret_dict[6] = sp.TissueLibrary.heavy_water() - ret_dict[7] = (sp.MolecularCompositionGenerator() - .append(sp.MoleculeLibrary.oxyhemoglobin(0.01)) - .append(sp.MoleculeLibrary.deoxyhemoglobin(0.01)) - .append(sp.MoleculeLibrary.water(0.98)) - .get_molecular_composition(sp.SegmentationClasses.COUPLING_ARTIFACT)) + ret_dict[7] = ( + sp.MolecularCompositionGenerator() + .append(sp.MoleculeLibrary.oxyhemoglobin(0.01)) + .append(sp.MoleculeLibrary.deoxyhemoglobin(0.01)) + .append(sp.MoleculeLibrary.water(0.98)) + .get_molecular_composition(sp.SegmentationClasses.COUPLING_ARTIFACT) + ) ret_dict[8] = sp.TissueLibrary.heavy_water() ret_dict[9] = sp.TissueLibrary.heavy_water() ret_dict[10] = sp.TissueLibrary.heavy_water() return ret_dict + return segmentation_volume_mask, segmentation_class_mapping() + + +def create_msotacuity_sample_segmentation( + volume_x_mm, volume_y_mm, volume_z_mm, spacing_mm +): + """ + Use the sample segmentation mask from manual tests, with a muscular background, an epidermis layer, + and two blood vessels at different depths. + """ + + from simpa_tests.manual_tests.volume_creation.SegmentationLoader import ( + create_segmentation_mask, + ) + + seg = create_segmentation_mask(volume_x_mm, volume_y_mm, volume_z_mm, spacing_mm) + + def segmentation_class_mapping(): + ret_dict = dict() + ret_dict[0] = sp.TissueLibrary.muscle() + ret_dict[1] = sp.TissueLibrary.epidermis() + ret_dict[2] = sp.TissueLibrary.blood() + ret_dict[3] = sp.TissueLibrary.blood() + return ret_dict + + return seg, segmentation_class_mapping() + + +def setup_rsom_device(settings, spacing, partial_volume=False): + """Set up RSOMExplorerP50 device.""" + seg, mapping = create_rsomexplorer_sample_segmentation(spacing) + + settings[Tags.VOLUME_NAME] = "SegmentationRSOMExplorer" + settings[Tags.DIM_VOLUME_X_MM] = seg.shape[0] * spacing + settings[Tags.DIM_VOLUME_Y_MM] = seg.shape[1] * spacing + settings[Tags.DIM_VOLUME_Z_MM] = seg.shape[2] * spacing + + volume_creation_settings = { + Tags.INPUT_SEGMENTATION_VOLUME: seg, + Tags.SEGMENTATION_CLASS_MAPPING: mapping, + Tags.CONSIDER_PARTIAL_VOLUME: partial_volume, + } + + settings.set_volume_creation_settings(volume_creation_settings) + + device = sp.RSOMExplorerP50(element_spacing_mm=1.0) + return device + + +def setup_msot_device(settings, spacing, partial_volume=False): + """Set up MSOTAcuityEcho device.""" + volume_x_mm = 40.0 + volume_y_mm = 20.0 + volume_z_mm = 25.0 + + seg, mapping = create_msotacuity_sample_segmentation( + volume_x_mm, volume_y_mm, volume_z_mm, spacing + ) + + settings[Tags.VOLUME_NAME] = "SegmentationMSOTAcuityEcho" + settings[Tags.DIM_VOLUME_X_MM] = volume_x_mm + settings[Tags.DIM_VOLUME_Y_MM] = volume_y_mm + settings[Tags.DIM_VOLUME_Z_MM] = volume_z_mm + + volume_creation_settings = { + Tags.INPUT_SEGMENTATION_VOLUME: seg, + Tags.SEGMENTATION_CLASS_MAPPING: mapping, + Tags.US_GEL: True, + Tags.CONSIDER_PARTIAL_VOLUME: partial_volume, + } + + settings.set_volume_creation_settings(volume_creation_settings) + + device = sp.MSOTAcuityEcho( + device_position_mm=np.array([volume_x_mm / 2, volume_y_mm / 2, 0]) + ) + device.update_settings_for_use_of_segmentation_based_volume_creator(settings) + return device + + +@profile +def run_segmentation_loader( + device_name: str = "RSOMExplorerP50", + spacing: float | int = 0.5, + partial_volume: bool = False, + path_manager=None, + visualise: bool = True, +): + """ + Runs a segmentation-based optical simulation. + + :param device_name: Device to use, one of "RSOMExplorerP50" or "MSOTAcuityEcho" + :param spacing: The simulation spacing between voxels in mm + :param partial_volume: If True, enables partial-volume blending at tissue boundaries + :param path_manager: the path manager to be used, typically sp.PathManager + :param visualise: If True, the simulation result will be plotted + """ + if path_manager is None: + path_manager = sp.PathManager() + + np.random.seed(4711) + settings = sp.Settings() settings[Tags.SIMULATION_PATH] = path_manager.get_hdf5_file_save_path() - settings[Tags.VOLUME_NAME] = "SegmentationTest" - settings[Tags.RANDOM_SEED] = 1234 + settings[Tags.RANDOM_SEED] = 4711 settings[Tags.WAVELENGTHS] = [700, 800] settings[Tags.SPACING_MM] = spacing - settings[Tags.DIM_VOLUME_X_MM] = segmentation_volume_mask.shape[0] * spacing - settings[Tags.DIM_VOLUME_Y_MM] = segmentation_volume_mask.shape[1] * spacing - settings[Tags.DIM_VOLUME_Z_MM] = segmentation_volume_mask.shape[2] * spacing - - settings.set_volume_creation_settings({ - Tags.INPUT_SEGMENTATION_VOLUME: segmentation_volume_mask, - Tags.SEGMENTATION_CLASS_MAPPING: segmentation_class_mapping(), - }) + if device_name == "MSOTAcuityEcho": + device = setup_msot_device(settings, spacing, partial_volume) + elif device_name == "RSOMExplorerP50": + device = setup_rsom_device(settings, spacing, partial_volume) + else: + raise ValueError(f"Unknown device: {device_name}. Choose from {DEVICE_CHOICES}") - settings.set_optical_settings({ - Tags.OPTICAL_MODEL_NUMBER_PHOTONS: 1e8, - Tags.OPTICAL_MODEL_BINARY_PATH: path_manager.get_mcx_binary_path(), - Tags.ILLUMINATION_TYPE: Tags.ILLUMINATION_TYPE_MSOT_ACUITY_ECHO, - Tags.LASER_PULSE_ENERGY_IN_MILLIJOULE: 50, - }) + settings.set_optical_settings( + { + Tags.OPTICAL_MODEL_NUMBER_PHOTONS: 1e7, + Tags.OPTICAL_MODEL_BINARY_PATH: path_manager.get_mcx_binary_path(), + Tags.ILLUMINATION_TYPE: Tags.ILLUMINATION_TYPE_MSOT_ACUITY_ECHO, + Tags.LASER_PULSE_ENERGY_IN_MILLIJOULE: 50, + } + ) pipeline = [ sp.SegmentationBasedAdapter(settings), - sp.MCXAdapter(settings) + sp.MCXAdapter(settings), ] - sp.simulate(pipeline, settings, sp.RSOMExplorerP50(element_spacing_mm=1.0)) + sp.simulate(pipeline, settings, device) if Tags.WAVELENGTH in settings: WAVELENGTH = settings[Tags.WAVELENGTH] @@ -98,19 +192,53 @@ def segmentation_class_mapping(): WAVELENGTH = 700 if visualise: - sp.visualise_data(path_to_hdf5_file=path_manager.get_hdf5_file_save_path() + "/" + "SegmentationTest" + ".hdf5", - wavelength=WAVELENGTH, - show_initial_pressure=True, - show_segmentation_map=True) + sp.visualise_data( + path_to_hdf5_file=settings[Tags.SIMPA_OUTPUT_FILE_PATH], + wavelength=WAVELENGTH, + show_initial_pressure=True, + show_segmentation_map=True, + show_absorption=True, + show_fluence=True, + save_path=path_manager.get_hdf5_file_save_path() + + "/" + + settings[Tags.VOLUME_NAME] + + ".png", + ) if __name__ == "__main__": - parser = ArgumentParser(description='Run the segmentation loader example') - parser.add_argument("--spacing", default=1, type=float, help='the voxel spacing in mm') - parser.add_argument("--input_spacing", default=0.2, type=float, help='the input spacing in mm') - parser.add_argument("--path_manager", default=None, help='the path manager, None uses sp.PathManager') - parser.add_argument("--visualise", default=True, type=bool, help='whether to visualise the result') + parser = ArgumentParser( + description="Run segmentation-based simulation with a configurable device" + ) + parser.add_argument( + "--device", + default="MSOTAcuityEcho", + choices=DEVICE_CHOICES, + help="the device to use", + ) + parser.add_argument( + "--spacing", default=0.5, type=float, help="the voxel spacing in mm" + ) + parser.add_argument( + "--partial_volume", + action="store_true", + default=True, + help="enable partial-volume blending at tissue boundaries", + ) + parser.add_argument( + "--path_manager", + default=None, + help="the path manager, None uses sp.PathManager", + ) + parser.add_argument( + "--visualise", default=True, type=bool, help="whether to visualise the result" + ) config = parser.parse_args() - run_segmentation_loader(spacing=config.spacing, input_spacing=config.input_spacing, - path_manager=config.path_manager, visualise=config.visualise) + run_segmentation_loader( + device_name=config.device, + spacing=config.spacing, + partial_volume=config.partial_volume, + path_manager=config.path_manager, + visualise=config.visualise, + ) diff --git a/simpa_tests/manual_tests/volume_creation/SegmentationLoader.py b/simpa_tests/manual_tests/volume_creation/SegmentationLoader.py index e4af5965..2fb3e0fb 100644 --- a/simpa_tests/manual_tests/volume_creation/SegmentationLoader.py +++ b/simpa_tests/manual_tests/volume_creation/SegmentationLoader.py @@ -6,26 +6,29 @@ import simpa as sp import numpy as np from skimage.data import shepp_logan_phantom -from scipy.ndimage import zoom +from scipy.ndimage import zoom, uniform_filter from simpa_tests.manual_tests import ManualIntegrationTestClass # FIXME temporary workaround for newest Intel architectures import os + os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE" -class SegmentationLoaderTest(ManualIntegrationTestClass): +class SegmentationLoaderRSOMExplorerTest(ManualIntegrationTestClass): def setup(self): self.path_manager = sp.PathManager() + self.volume_name = "SegmentationRSOMExplorerTest" target_spacing = 1.0 label_mask = shepp_logan_phantom() label_mask = np.digitize(label_mask, bins=np.linspace(0.0, 1.0, 11), right=True) label_mask = np.reshape(label_mask, (400, 1, 400)) input_spacing = 0.2 segmentation_volume_tiled = np.tile(label_mask, (1, 128, 1)) - segmentation_volume_mask = sp.round_x5_away_from_zero(zoom(segmentation_volume_tiled, input_spacing/target_spacing, - order=0)).astype(int) + segmentation_volume_mask = sp.round_x5_away_from_zero( + zoom(segmentation_volume_tiled, input_spacing / target_spacing, order=0) + ).astype(int) def segmentation_class_mapping(): ret_dict = dict() @@ -36,11 +39,13 @@ def segmentation_class_mapping(): ret_dict[4] = sp.TissueLibrary.mediprene() ret_dict[5] = sp.TissueLibrary.ultrasound_gel() ret_dict[6] = sp.TissueLibrary.heavy_water() - ret_dict[7] = (sp.MolecularCompositionGenerator() - .append(sp.MoleculeLibrary.oxyhemoglobin(0.01)) - .append(sp.MoleculeLibrary.deoxyhemoglobin(0.01)) - .append(sp.MoleculeLibrary.water(0.98)) - .get_molecular_composition(sp.SegmentationClasses.COUPLING_ARTIFACT)) + ret_dict[7] = ( + sp.MolecularCompositionGenerator() + .append(sp.MoleculeLibrary.oxyhemoglobin(0.01)) + .append(sp.MoleculeLibrary.deoxyhemoglobin(0.01)) + .append(sp.MoleculeLibrary.water(0.98)) + .get_molecular_composition(sp.SegmentationClasses.COUPLING_ARTIFACT) + ) ret_dict[8] = sp.TissueLibrary.heavy_water() ret_dict[9] = sp.TissueLibrary.heavy_water() ret_dict[10] = sp.TissueLibrary.heavy_water() @@ -48,8 +53,10 @@ def segmentation_class_mapping(): return ret_dict self.settings = sp.Settings() - self.settings[Tags.SIMULATION_PATH] = self.path_manager.get_hdf5_file_save_path() - self.settings[Tags.VOLUME_NAME] = "SegmentationTest" + self.settings[Tags.SIMULATION_PATH] = ( + self.path_manager.get_hdf5_file_save_path() + ) + self.settings[Tags.VOLUME_NAME] = self.volume_name self.settings[Tags.RANDOM_SEED] = 1234 self.settings[Tags.WAVELENGTHS] = [700] self.settings[Tags.SPACING_MM] = target_spacing @@ -58,29 +65,38 @@ def segmentation_class_mapping(): self.settings[Tags.DIM_VOLUME_Z_MM] = 400 / (target_spacing / input_spacing) # self.settings[Tags.IGNORE_QA_ASSERTIONS] = True - self.settings.set_volume_creation_settings({ - Tags.INPUT_SEGMENTATION_VOLUME: segmentation_volume_mask, - Tags.SEGMENTATION_CLASS_MAPPING: segmentation_class_mapping(), - - }) + self.settings.set_volume_creation_settings( + { + Tags.INPUT_SEGMENTATION_VOLUME: segmentation_volume_mask, + Tags.SEGMENTATION_CLASS_MAPPING: segmentation_class_mapping(), + } + ) - self.settings.set_optical_settings({ - Tags.OPTICAL_MODEL_NUMBER_PHOTONS: 1e7, - Tags.OPTICAL_MODEL_BINARY_PATH: self.path_manager.get_mcx_binary_path(), - Tags.ILLUMINATION_TYPE: Tags.ILLUMINATION_TYPE_MSOT_ACUITY_ECHO, - Tags.LASER_PULSE_ENERGY_IN_MILLIJOULE: 50, - }) + self.settings.set_optical_settings( + { + Tags.OPTICAL_MODEL_NUMBER_PHOTONS: 1e7, + Tags.OPTICAL_MODEL_BINARY_PATH: self.path_manager.get_mcx_binary_path(), + Tags.ILLUMINATION_TYPE: Tags.ILLUMINATION_TYPE_MSOT_ACUITY_ECHO, + Tags.LASER_PULSE_ENERGY_IN_MILLIJOULE: 50, + } + ) self.pipeline = [ sp.SegmentationBasedAdapter(self.settings), - sp.MCXAdapter(self.settings) + sp.MCXAdapter(self.settings), ] def perform_test(self): - sp.simulate(self.pipeline, self.settings, sp.RSOMExplorerP50(element_spacing_mm=2.0, - number_elements_y=10, - number_elements_x=20, - device_position_mm=np.asarray([20, 10, 0]))) + sp.simulate( + self.pipeline, + self.settings, + sp.RSOMExplorerP50( + element_spacing_mm=2.0, + number_elements_y=10, + number_elements_x=20, + device_position_mm=np.asarray([20, 10, 0]), + ), + ) def tear_down(self): os.remove(self.settings[Tags.SIMPA_OUTPUT_FILE_PATH]) @@ -90,22 +106,219 @@ def visualise_result(self, show_figure_on_screen=True, save_path=None): if show_figure_on_screen: save_path = None else: - save_path = save_path + "SegmentationLoaderExample.png" - - sp.visualise_data(path_to_hdf5_file=self.path_manager.get_hdf5_file_save_path() + "/" + "SegmentationTest" + ".hdf5", - wavelength=700, - show_initial_pressure=True, - show_segmentation_map=True, - show_absorption=True, - show_fluence=True, - show_tissue_density=True, - show_speed_of_sound=True, - show_anisotropy=True, - show_scattering=True, - save_path=save_path, - log_scale=False) + save_path = save_path + f"{self.volume_name}.png" + + sp.visualise_data( + path_to_hdf5_file=self.path_manager.get_hdf5_file_save_path() + + "/" + + self.volume_name + + ".hdf5", + wavelength=700, + show_initial_pressure=True, + show_segmentation_map=True, + show_absorption=True, + show_fluence=True, + show_tissue_density=True, + show_speed_of_sound=True, + show_anisotropy=True, + show_scattering=True, + save_path=save_path, + log_scale=False, + ) + + +def create_segmentation_mask(volume_x_mm, volume_y_mm, volume_z_mm, spacing_mm): + """ + Creates a simple segmentation mask with a muscular background, an epidermis layer, + and two blood vessels at different depths. + """ + nx = int(round(volume_x_mm / spacing_mm)) + ny = int(round(volume_y_mm / spacing_mm)) + nz = int(round(volume_z_mm / spacing_mm)) + + # background (muscle): class 0 + seg = np.zeros((nx, ny, nz), dtype=int) + + # epidermis: class 1, 1mm at the top of the volume + epidermis_end_pix = max(1, int(round(1.0 / spacing_mm))) + seg[:, :, :epidermis_end_pix] = 1 + + # shallow blood vessel: class 2, cylinder at z=5mm, radius 2mm + vessel1_x_center = volume_x_mm / 2 - 8 + vessel1_z_center = 5.0 + vessel1_radius = 2.0 + + # deep blood vessel: class 3, cylinder at z=12mm, radius 3mm + vessel2_x_center = volume_x_mm / 2 + 5 + vessel2_z_center = 12.0 + vessel2_radius = 3.0 + + x_coords = (np.arange(nx) + 0.5) * spacing_mm + z_coords = (np.arange(nz) + 0.5) * spacing_mm + xx, zz = np.meshgrid(x_coords, z_coords, indexing="ij") + + dist1 = np.sqrt((xx - vessel1_x_center) ** 2 + (zz - vessel1_z_center) ** 2) + dist2 = np.sqrt((xx - vessel2_x_center) ** 2 + (zz - vessel2_z_center) ** 2) + + vessel1_mask_2d = dist1 <= vessel1_radius + vessel2_mask_2d = dist2 <= vessel2_radius + + # repeat accross the y-dimension + for iy in range(ny): + seg[:, iy, :][vessel1_mask_2d] = 2 + seg[:, iy, :][vessel2_mask_2d] = 3 + + return seg + + +class SegmentationLoaderMSOTAcuityEchoTest(ManualIntegrationTestClass): + """ + Tests the MSOTAcuityEcho device with the segmentation-based volume creator, + including the update_settings_for_use_of_segmentation_based_volume_creator method + which adds mediprene, heavy water, and optionally US gel layers. + PV_EFFECTS enables partial-volume effects in the SegmentationBasedAdapter + FRACTION_MAP can be used to simulate a pre-computed 4D fraction map, where partial volume effects can already be included + """ + + PV_EFFECTS = False + FRACTION_MAP = False + + def setup(self): + self.path_manager = sp.PathManager() + + self.volume_name = f"SegmentationMSOTAcuityEchoTest{'_PV' if self.PV_EFFECTS else ''}{'_4D' if self.FRACTION_MAP else ''}" + + target_spacing = 0.5 + volume_x_mm = 40.0 + volume_y_mm = 20.0 + volume_z_mm = 25.0 + + np.random.seed(4711) + + segmentation_volume_mask = create_segmentation_mask( + volume_x_mm, volume_y_mm, volume_z_mm, target_spacing + ) + + if self.FRACTION_MAP: + # transform into 4d fraction map and simulate existing partial volume effects via smoothing + # One-hot encode to [C, X, Y, Z] + # Labels must be dense (0...N-1) so that channel index == label value, which is + # required by update_settings_for_use_of_segmentation_based_volume_creator. + num_classes = int(segmentation_volume_mask.max()) + 1 + seg_4d = np.zeros( + (num_classes, *segmentation_volume_mask.shape), dtype=np.float32 + ) + for c in range(num_classes): + seg_4d[c] = (segmentation_volume_mask == c).astype(np.float32) + if self.PV_EFFECTS: + for c in range(num_classes): + seg_4d[c] = uniform_filter(seg_4d[c], size=3).astype(np.float32) + segmentation_volume_mask = seg_4d + + def segmentation_class_mapping(): + ret_dict = dict() + ret_dict[0] = sp.TissueLibrary.muscle() + ret_dict[1] = sp.TissueLibrary.epidermis() + ret_dict[2] = sp.TissueLibrary.blood() + ret_dict[3] = sp.TissueLibrary.blood() + return ret_dict + + self.settings = sp.Settings() + self.settings[Tags.SIMULATION_PATH] = ( + self.path_manager.get_hdf5_file_save_path() + ) + self.settings[Tags.VOLUME_NAME] = self.volume_name + self.settings[Tags.RANDOM_SEED] = 4711 + self.settings[Tags.WAVELENGTHS] = [700, 800] + self.settings[Tags.SPACING_MM] = target_spacing + self.settings[Tags.DIM_VOLUME_X_MM] = volume_x_mm + self.settings[Tags.DIM_VOLUME_Y_MM] = volume_y_mm + self.settings[Tags.DIM_VOLUME_Z_MM] = volume_z_mm + + self.settings.set_volume_creation_settings( + { + Tags.INPUT_SEGMENTATION_VOLUME: segmentation_volume_mask, + Tags.SEGMENTATION_CLASS_MAPPING: segmentation_class_mapping(), + Tags.US_GEL: True, + Tags.CONSIDER_PARTIAL_VOLUME: self.PV_EFFECTS + and not self.FRACTION_MAP, # only for 3D map. 4D map already smoothed if PV_EFFECTS is True. + } + ) + + self.settings.set_optical_settings( + { + Tags.OPTICAL_MODEL_NUMBER_PHOTONS: 1e7, + Tags.OPTICAL_MODEL_BINARY_PATH: self.path_manager.get_mcx_binary_path(), + Tags.ILLUMINATION_TYPE: Tags.ILLUMINATION_TYPE_MSOT_ACUITY_ECHO, + Tags.LASER_PULSE_ENERGY_IN_MILLIJOULE: 50, + } + ) + + # Set up the MSOTAcuityEcho device + self.device = sp.MSOTAcuityEcho( + device_position_mm=np.array([volume_x_mm / 2, volume_y_mm / 2, 0]) + ) + + # Use update_settings_for_use_of_segmentation_based_volume_creator to add + # mediprene, heavy water, and us gel (US gel controlled via Tags.US_GEL). + self.device.update_settings_for_use_of_segmentation_based_volume_creator( + self.settings, + ) + + self.pipeline = [ + sp.SegmentationBasedAdapter(self.settings), + sp.MCXAdapter(self.settings), + ] + + def perform_test(self): + sp.simulate(self.pipeline, self.settings, self.device) + + def tear_down(self): + os.remove(self.settings[Tags.SIMPA_OUTPUT_FILE_PATH]) + + def visualise_result(self, show_figure_on_screen=True, save_path=None): + if show_figure_on_screen: + save_path = None + else: + save_path = save_path + f"{self.volume_name}.png" + + sp.visualise_data( + path_to_hdf5_file=self.path_manager.get_hdf5_file_save_path() + + "/" + + self.volume_name + + ".hdf5", + wavelength=700, + show_initial_pressure=True, + show_segmentation_map=True, + show_absorption=True, + show_fluence=True, + show_tissue_density=True, + show_speed_of_sound=True, + show_anisotropy=True, + show_scattering=True, + save_path=save_path, + log_scale=False, + ) if __name__ == "__main__": - test = SegmentationLoaderTest() - test.run_test(show_figure_on_screen=False) + print("Running SegmentationLoaderRSOMExplorerTest...") + test_rsom = SegmentationLoaderRSOMExplorerTest() + test_rsom.run_test(show_figure_on_screen=False) + + print("Running SegmentationLoaderMSOTAcuityEchoTest (3D label map)...") + test_msot = SegmentationLoaderMSOTAcuityEchoTest() + test_msot.run_test(show_figure_on_screen=False) + + print("Running SegmentationLoaderMSOTAcuityEchoTest (3D label map, PV effects)...") + test_msot_pv = SegmentationLoaderMSOTAcuityEchoTest() + test_msot_pv.PV_EFFECTS = True + test_msot_pv.run_test(show_figure_on_screen=False) + + print( + "Running SegmentationLoaderMSOTAcuityEchoTest (4D fraction map including PV effects)..." + ) + test_msot_4d = SegmentationLoaderMSOTAcuityEchoTest() + test_msot_4d.PV_EFFECTS = True + test_msot_4d.FRACTION_MAP = True + test_msot_4d.run_test(show_figure_on_screen=False)