From 9d21a70541549702d15b3a3063d2ffa3ad295225 Mon Sep 17 00:00:00 2001 From: Friso Grace Date: Fri, 12 Jul 2024 14:33:42 +0200 Subject: [PATCH 01/12] Segmentation fully works with the MSOT device and updates the settings properly --- .../pa_devices/ithera_msot_acuity.py | 131 ++++++++++++++++++ simpa/utils/tags.py | 19 +++ 2 files changed, 150 insertions(+) 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 091c9396..c0fd3052 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 @@ -1,6 +1,7 @@ # SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT +import torch from simpa.core.device_digital_twins import PhotoacousticDevice, \ CurvedArrayDetectionGeometry, MSOTAcuityIlluminationGeometry @@ -207,6 +208,136 @@ 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, + add_layers: list = [Tags.ADD_US_GEL, + Tags.ADD_MEDIPRENE, + Tags.ADD_HEAVY_WATER], + current_heavy_water_depth: (float, int) = 0, + heavy_water_tag: int = None): + """ + Updates the volume creation settings of the segmentation based volume creator according to the size of the + device. On the occasion that your segmentation already includes the mediprene, ultrasound gel and some of the + heavy water, you may specify the existing depth of the heavy water so that it can be adapted to the depth of + device. + :param add_layers: The layers to add to the volume, all configured to the typical thicknesses for MSOT acuity + echo. + :param current_heavy_water_depth: the current heavy water depth (mm). + :param heavy_water_tag: the existing heavy water tag in the segmentation map. + :param global_settings: Settings for the entire simulation pipeline. + """ + try: + volume_creator_settings = Settings(global_settings.get_volume_creation_settings()) + except KeyError as e: + 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 + mediprene_layer_height_mm = 0 + heavy_water_layer_height_mm = 0 + + if Tags.ADD_US_GEL in add_layers: + us_gel_thickness_mm = np.random.normal(0.4, 0.1) + us_gel_thickness_pix = int(round(us_gel_thickness_mm/spacing_mm)) + padding_dims = ((0, 0), (0, 0), (us_gel_thickness_pix, 0)) + segmentation_map = np.pad(segmentation_map, padding_dims, mode='constant', constant_values=64) + segmentation_class_mapping[64] = TISSUE_LIBRARY.ultrasound_gel() + z_dim_position_shift_mm += us_gel_thickness_pix * spacing_mm + + if Tags.ADD_MEDIPRENE in add_layers: + mediprene_layer_height_mm = self.mediprene_membrane_height_mm + mediprene_layer_height_pix = int(round(mediprene_layer_height_mm/spacing_mm)) + padding_dims = ((0, 0), (0, 0), (mediprene_layer_height_pix, 0)) + segmentation_map = np.pad(segmentation_map, padding_dims, mode='constant', constant_values=128) + segmentation_class_mapping[128] = TISSUE_LIBRARY.mediprene() + z_dim_position_shift_mm += mediprene_layer_height_pix * spacing_mm + + if Tags.ADD_HEAVY_WATER in add_layers: + if heavy_water_tag is None: + heavy_water_tag = 256 + segmentation_class_mapping[256] = TISSUE_LIBRARY.heavy_water() + probe_size_mm = self.probe_height_mm + mediprene_layer_height_mm = self.mediprene_membrane_height_mm + heavy_water_layer_height_mm = probe_size_mm - current_heavy_water_depth - mediprene_layer_height_mm + heavy_water_layer_height_pix = int(round(heavy_water_layer_height_mm / spacing_mm)) + padding_dims = ((0, 0), (0, 0), (heavy_water_layer_height_pix, 0)) + segmentation_map = np.pad(segmentation_map, padding_dims, mode='constant', constant_values=heavy_water_tag) + segmentation_class_mapping[heavy_water_tag] = TISSUE_LIBRARY.heavy_water() + z_dim_position_shift_mm += heavy_water_layer_height_pix * spacing_mm + + new_volume_height_mm = global_settings[Tags.DIM_VOLUME_Z_MM] + z_dim_position_shift_mm + + # adjust the z-dim to msot probe height + global_settings[Tags.DIM_VOLUME_Z_MM] = new_volume_height_mm + self.logger.debug(f"Changed Tags.DIM_VOLUME_Z_MM to {global_settings[Tags.DIM_VOLUME_Z_MM]}") + + # adjust the x-dim to msot probe width + # 1 voxel is added (0.5 on both sides) to make sure no rounding errors lead to a detector element being outside + # of the simulated volume. + + 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 + # specific left and right to avoid rounding errors + left_shift_pixels = int(round(width_shift_for_structures_mm / spacing_mm)) + right_shift_pixels = int(round((round(self.detection_geometry.probe_width_mm) + spacing_mm - + global_settings[Tags.DIM_VOLUME_X_MM])/spacing_mm)) - left_shift_pixels + padding_width = ((left_shift_pixels, right_shift_pixels), (0, 0), (0, 0)) + segmentation_map = np.pad(segmentation_map, 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]}") + + else: + width_shift_for_structures_mm = 0 + padding_width = ((0, 0), (0, 0), (0, 0)) + + 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") + + for structure_key in volume_creator_settings[Tags.SEGMENTATION_CLASS_MAPPING]: + self.logger.debug("Adjusting " + str(structure_key)) + structure_dict = volume_creator_settings[Tags.SEGMENTATION_CLASS_MAPPING][structure_key] + for molecule in structure_dict: + try: + old_volume_fraction = getattr(molecule, Tags.VOLUME_FRACTION) + except AttributeError: + continue + if isinstance(old_volume_fraction, torch.Tensor): + if old_volume_fraction.shape != segmentation_map.shape: + z_shift_pixels = int(round(z_dim_position_shift_mm / spacing_mm)) + padding_height = ((0, 0), (0, 0), (z_shift_pixels, 0)) + padded_up = np.pad(old_volume_fraction.numpy(), padding_height, mode='edge') + padded_vol = np.pad(padded_up, padding_width, mode='edge') + setattr(molecule, Tags.VOLUME_FRACTION, torch.tensor(padded_vol, dtype=torch.float32)) + + device_change_in_height = mediprene_layer_height_mm + heavy_water_layer_height_mm + self.device_position_mm = np.add(self.device_position_mm, np.array([width_shift_for_structures_mm, 0, + device_change_in_height])) + 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=55, + 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, + device_change_in_height])) + def serialize(self) -> dict: serialized_device = self.__dict__ device_dict = {"MSOTAcuityEcho": serialized_device} diff --git a/simpa/utils/tags.py b/simpa/utils/tags.py index 39f3c4f4..8e932c4d 100644 --- a/simpa/utils/tags.py +++ b/simpa/utils/tags.py @@ -1496,3 +1496,22 @@ class Tags: """ Identifier for the environment varibale that defines the path the the matlab executable. """ + + ADD_US_GEL = "add_us_gel" + """ + Identifier to add ultrasound gel to the segmentation + Usage:simpa.core.digital_device_twins.pa_devices.ithera_msot_acuity + """ + + ADD_MEDIPRENE = "add_mediprene" + """ + Identifier to add mediprene to the segmentation + Usage:simpa.core.digital_device_twins.pa_devices.ithera_msot_acuity + """ + + ADD_HEAVY_WATER = "add_heavy_water" + """ + Identifier to add heavy water to the segmentation + Usage:simpa.core.digital_device_twins.pa_devices.ithera_msot_acuity + + """ From fcf484934d8a0e76af28cf565fbf56820c160a7a Mon Sep 17 00:00:00 2001 From: Friso Grace Date: Fri, 12 Jul 2024 15:44:57 +0200 Subject: [PATCH 02/12] Lots of extra type hints :) --- simpa/core/device_digital_twins/__init__.py | 20 +++++++++++++++---- .../detection_geometries/__init__.py | 9 +++++---- .../detection_geometries/curved_array.py | 14 ++++++------- .../detection_geometries/linear_array.py | 7 ++++--- .../detection_geometries/planar_array.py | 3 ++- .../illumination_geometries/__init__.py | 5 +++-- .../disk_illumination.py | 3 ++- .../gaussian_beam_illumination.py | 5 +++-- .../ithera_msot_invision_illumination.py | 7 +------ .../pencil_array_illumination.py | 4 ++-- .../slit_illumination.py | 2 +- .../pa_devices/ithera_rsom.py | 6 +++--- 12 files changed, 49 insertions(+), 36 deletions(-) diff --git a/simpa/core/device_digital_twins/__init__.py b/simpa/core/device_digital_twins/__init__.py index 13bd5bcf..c34d1760 100644 --- a/simpa/core/device_digital_twins/__init__.py +++ b/simpa/core/device_digital_twins/__init__.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): @@ -17,7 +18,7 @@ class DigitalDeviceTwinBase(SerializableSIMPAClass): which has representations of both. """ - def __init__(self, device_position_mm=None, field_of_view_extent_mm=None): + def __init__(self, device_position_mm: np.ndarray = None, field_of_view_extent_mm: np.ndarray = None): """ :param device_position_mm: Each device has an internal position which serves as origin for internal \ representations of e.g. detector element positions or illuminator positions. @@ -54,7 +55,7 @@ def __eq__(self, other): return False @abstractmethod - def check_settings_prerequisites(self, global_settings) -> bool: + def check_settings_prerequisites(self, global_settings: Settings) -> bool: """ It might be that certain device geometries need a certain dimensionality of the simulated PAI volume, or that it requires the existence of certain Tags in the global global_settings. @@ -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 + """ + def get_field_of_view_mm(self) -> np.ndarray: """ Returns the absolute field of view in mm where the probe position is already @@ -133,7 +144,8 @@ def deserialize(dictionary_to_deserialize): """ -It is important to have these relative imports after the definition of the DigitalDeviceTwinBase class to avoid circular imports triggered by imported child classes +It is important to have these relative imports after the definition of the DigitalDeviceTwinBase class to avoid circular +imports triggered by imported child classes """ from .pa_devices import PhotoacousticDevice # nopep8 from simpa.core.device_digital_twins.detection_geometries import DetectionGeometryBase # nopep8 diff --git a/simpa/core/device_digital_twins/detection_geometries/__init__.py b/simpa/core/device_digital_twins/detection_geometries/__init__.py index eba6bdd5..8afa7188 100644 --- a/simpa/core/device_digital_twins/detection_geometries/__init__.py +++ b/simpa/core/device_digital_twins/detection_geometries/__init__.py @@ -5,6 +5,7 @@ from abc import abstractmethod from simpa.core.device_digital_twins import DigitalDeviceTwinBase import numpy as np +from typing import Union class DetectionGeometryBase(DigitalDeviceTwinBase): @@ -12,10 +13,10 @@ class DetectionGeometryBase(DigitalDeviceTwinBase): This class is the base class for representing all detector geometries. """ - def __init__(self, number_detector_elements, detector_element_width_mm, - detector_element_length_mm, center_frequency_hz, bandwidth_percent, - sampling_frequency_mhz, device_position_mm: np.ndarray = None, - field_of_view_extent_mm: np.ndarray = None): + def __init__(self, number_detector_elements: int, detector_element_width_mm: Union[float, int], + detector_element_length_mm: Union[float, int], center_frequency_hz: Union[float, int], + bandwidth_percent: Union[float, int], sampling_frequency_mhz: Union[float, int], + device_position_mm: np.ndarray = None, field_of_view_extent_mm: np.ndarray = None): """ :param number_detector_elements: Total number of detector elements. 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..2ba7132a 100644 --- a/simpa/core/device_digital_twins/detection_geometries/curved_array.py +++ b/simpa/core/device_digital_twins/detection_geometries/curved_array.py @@ -2,9 +2,9 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT import numpy as np +from typing import Union from simpa.core.device_digital_twins import DetectionGeometryBase -from simpa.utils import Tags class CurvedArrayDetectionGeometry(DetectionGeometryBase): @@ -13,17 +13,17 @@ class CurvedArrayDetectionGeometry(DetectionGeometryBase): with a curved detection geometry. The origin for this device is the center (focus) of the curved array. """ - def __init__(self, pitch_mm=0.5, - radius_mm=40, - number_detector_elements=256, + def __init__(self, pitch_mm: Union[float, int] = 0.5, + radius_mm: Union[float, int] = 40, + number_detector_elements: int = 256, detector_element_width_mm=0.24, detector_element_length_mm=13, center_frequency_hz=3.96e6, bandwidth_percent=55, sampling_frequency_mhz=40, - angular_origin_offset=np.pi, - device_position_mm=None, - field_of_view_extent_mm=None): + angular_origin_offset: Union[float, int] = np.pi, + device_position_mm: np.ndarray = None, + field_of_view_extent_mm: np.ndarray = None): """ :param pitch_mm: In-plane distance between the beginning of one detector element to the next detector element. 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..18bb12e7 100644 --- a/simpa/core/device_digital_twins/detection_geometries/linear_array.py +++ b/simpa/core/device_digital_twins/detection_geometries/linear_array.py @@ -2,9 +2,10 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT import numpy as np +from typing import Union from simpa.core.device_digital_twins import DetectionGeometryBase -from simpa.utils import Settings, Tags +from simpa.utils import Tags class LinearArrayDetectionGeometry(DetectionGeometryBase): @@ -15,7 +16,7 @@ class LinearArrayDetectionGeometry(DetectionGeometryBase): """ - def __init__(self, pitch_mm=0.5, + def __init__(self, pitch_mm: Union[float, int] = 0.5, number_detector_elements=100, detector_element_width_mm=0.24, detector_element_length_mm=0.5, @@ -51,7 +52,7 @@ def __init__(self, pitch_mm=0.5, self.pitch_mm = pitch_mm self.probe_width_mm = (number_detector_elements - 1) * self.pitch_mm - def check_settings_prerequisites(self, global_settings: Settings) -> bool: + def check_settings_prerequisites(self, global_settings) -> bool: if global_settings[Tags.DIM_VOLUME_X_MM] < self.probe_width_mm + global_settings[Tags.SPACING_MM]: self.logger.error("Volume x dimension is too small to encompass MSOT device in simulation!" "Must be at least {} mm but was {} mm" 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..886bb4d7 100644 --- a/simpa/core/device_digital_twins/detection_geometries/planar_array.py +++ b/simpa/core/device_digital_twins/detection_geometries/planar_array.py @@ -5,6 +5,7 @@ from simpa.core.device_digital_twins import DetectionGeometryBase from simpa.utils import Settings, Tags +from typing import Union class PlanarArrayDetectionGeometry(DetectionGeometryBase): @@ -14,7 +15,7 @@ class PlanarArrayDetectionGeometry(DetectionGeometryBase): """ - def __init__(self, pitch_mm=0.5, + def __init__(self, pitch_mm: Union[float, int] = 0.5, number_detector_elements_x=100, number_detector_elements_y=100, detector_element_width_mm=0.24, diff --git a/simpa/core/device_digital_twins/illumination_geometries/__init__.py b/simpa/core/device_digital_twins/illumination_geometries/__init__.py index 6c5780bf..a19fcfa9 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/__init__.py +++ b/simpa/core/device_digital_twins/illumination_geometries/__init__.py @@ -13,7 +13,8 @@ class IlluminationGeometryBase(DigitalDeviceTwinBase): This class is the base class for representing all illumination geometries. """ - def __init__(self, device_position_mm=None, source_direction_vector=None, field_of_view_extent_mm=None): + def __init__(self, device_position_mm=None, source_direction_vector: np.ndarray = None, + field_of_view_extent_mm=None): """ :param device_position_mm: Each device has an internal position which serves as origin for internal \ representations of illuminator positions. @@ -38,7 +39,7 @@ def __init__(self, device_position_mm=None, source_direction_vector=None, field_ self.source_direction_vector) @abstractmethod - def get_mcx_illuminator_definition(self, global_settings) -> dict: + def get_mcx_illuminator_definition(self, global_settings: Settings) -> dict: """ IMPORTANT: This method creates a dictionary that contains tags as they are expected for the mcx simulation tool to represent the illumination geometry of this device. diff --git a/simpa/core/device_digital_twins/illumination_geometries/disk_illumination.py b/simpa/core/device_digital_twins/illumination_geometries/disk_illumination.py index 42b12909..51e2eddc 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/disk_illumination.py +++ b/simpa/core/device_digital_twins/illumination_geometries/disk_illumination.py @@ -4,6 +4,7 @@ from simpa.core.device_digital_twins import IlluminationGeometryBase from simpa.utils import Tags +from typing import Union class DiskIlluminationGeometry(IlluminationGeometryBase): @@ -12,7 +13,7 @@ class DiskIlluminationGeometry(IlluminationGeometryBase): The device position is defined as the middle of the disk. """ - def __init__(self, beam_radius_mm=None, device_position_mm=None, field_of_view_extent_mm=None): + def __init__(self, beam_radius_mm: Union[int, float] = None, device_position_mm=None, field_of_view_extent_mm=None): """ :param beam_radius_mm: Radius of the disk in mm. :type beam_radius_mm: int, float diff --git a/simpa/core/device_digital_twins/illumination_geometries/gaussian_beam_illumination.py b/simpa/core/device_digital_twins/illumination_geometries/gaussian_beam_illumination.py index b3a9b8a6..1f0a14dd 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/gaussian_beam_illumination.py +++ b/simpa/core/device_digital_twins/illumination_geometries/gaussian_beam_illumination.py @@ -6,6 +6,7 @@ from collections.abc import Sized from simpa.core.device_digital_twins import IlluminationGeometryBase from simpa.utils import Tags +from typing import Union class GaussianBeamIlluminationGeometry(IlluminationGeometryBase): @@ -14,8 +15,8 @@ class GaussianBeamIlluminationGeometry(IlluminationGeometryBase): The position is defined as the middle of the beam. """ - def __init__(self, beam_radius_mm=None, focal_length_mm=None, device_position_mm=None, - field_of_view_extent_mm=None): + def __init__(self, beam_radius_mm: Union[int, float] = None, focal_length_mm: Union[int, float] = None, + device_position_mm=None, field_of_view_extent_mm=None): """ :param beam_radius_mm: Initial radius of the gaussian beam at half maximum (full width at half maximum (FWHM)) in mm. diff --git a/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_invision_illumination.py b/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_invision_illumination.py index 0c22a346..a31f4c6a 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_invision_illumination.py +++ b/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_invision_illumination.py @@ -13,14 +13,9 @@ class MSOTInVisionIlluminationGeometry(IlluminationGeometryBase): This class represents the illumination geometry of the MSOT InVision photoacoustic device. """ - def __init__(self, invision_position=None): + def __init__(self, invision_position: list = [0, 0, 0]): super().__init__() - if invision_position is None: - self.invision_position = [0, 0, 0] - else: - self.invision_position = invision_position - det_sep_half = 24.74 / 2 detector_iso_distance = 74.05 / 2 detector_width = 2 * 6.12 diff --git a/simpa/core/device_digital_twins/illumination_geometries/pencil_array_illumination.py b/simpa/core/device_digital_twins/illumination_geometries/pencil_array_illumination.py index e5a4deab..0712bc86 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/pencil_array_illumination.py +++ b/simpa/core/device_digital_twins/illumination_geometries/pencil_array_illumination.py @@ -14,8 +14,8 @@ class PencilArrayIlluminationGeometry(IlluminationGeometryBase): The device position is defined as the middle of the array. """ - def __init__(self, pitch_mm=0.5, number_illuminators_x=100, number_illuminators_y=100, device_position_mm=None, - field_of_view_extent_mm=None): + def __init__(self, pitch_mm: float = 0.5, number_illuminators_x: int = 100, number_illuminators_y: int = 100, + device_position_mm=None, field_of_view_extent_mm=None): """ :param pitch_mm: Defines the x and y distance between the illumination positions :type pitch_mm: float diff --git a/simpa/core/device_digital_twins/illumination_geometries/slit_illumination.py b/simpa/core/device_digital_twins/illumination_geometries/slit_illumination.py index 4dc58eba..87b05bf0 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/slit_illumination.py +++ b/simpa/core/device_digital_twins/illumination_geometries/slit_illumination.py @@ -14,7 +14,7 @@ class SlitIlluminationGeometry(IlluminationGeometryBase): The device position is defined as the middle of the slit. """ - def __init__(self, slit_vector_mm=None, direction_vector_mm=None, device_position_mm=None, + def __init__(self, slit_vector_mm: list = None, direction_vector_mm: list = None, device_position_mm=None, field_of_view_extent_mm=None): """ :param slit_vector_mm: Defines the slit in vector form. For example a slit along the x-axis with length 5mm diff --git a/simpa/core/device_digital_twins/pa_devices/ithera_rsom.py b/simpa/core/device_digital_twins/pa_devices/ithera_rsom.py index eba16ba3..152eb27f 100644 --- a/simpa/core/device_digital_twins/pa_devices/ithera_rsom.py +++ b/simpa/core/device_digital_twins/pa_devices/ithera_rsom.py @@ -35,9 +35,9 @@ class RSOMExplorerP50(PhotoacousticDevice): """ - def __init__(self, element_spacing_mm=0.02, - number_elements_x=10, - number_elements_y=10, + def __init__(self, element_spacing_mm: float = 0.02, + number_elements_x: int = 10, + number_elements_y: int = 10, device_position_mm: np.ndarray = None, field_of_view_extent_mm: np.ndarray = None): """ From 34fdf5ffd14e65be1194e060715a9faafc3efebc Mon Sep 17 00:00:00 2001 From: Friso Grace Date: Fri, 12 Jul 2024 15:58:44 +0200 Subject: [PATCH 03/12] Lots of extra type hints :) --- .../detection_geometries/curved_array.py | 4 ++++ .../detection_geometries/linear_array.py | 4 ++++ .../detection_geometries/planar_array.py | 3 +++ .../device_digital_twins/illumination_geometries/__init__.py | 5 ++++- .../ithera_msot_invision_illumination.py | 2 ++ simpa/core/device_digital_twins/pa_devices/__init__.py | 5 +++++ 6 files changed, 22 insertions(+), 1 deletion(-) 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 2ba7132a..f281516e 100644 --- a/simpa/core/device_digital_twins/detection_geometries/curved_array.py +++ b/simpa/core/device_digital_twins/detection_geometries/curved_array.py @@ -5,6 +5,7 @@ from typing import Union from simpa.core.device_digital_twins import DetectionGeometryBase +from simpa import Settings, Tags class CurvedArrayDetectionGeometry(DetectionGeometryBase): @@ -85,6 +86,9 @@ def check_settings_prerequisites(self, global_settings) -> bool: def update_settings_for_use_of_model_based_volume_creator(self, global_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: pitch_angle = self.pitch_mm / self.radius_mm 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 18bb12e7..b3a1821c 100644 --- a/simpa/core/device_digital_twins/detection_geometries/linear_array.py +++ b/simpa/core/device_digital_twins/detection_geometries/linear_array.py @@ -4,6 +4,7 @@ import numpy as np from typing import Union +from simpa import Settings from simpa.core.device_digital_twins import DetectionGeometryBase from simpa.utils import Tags @@ -64,6 +65,9 @@ def check_settings_prerequisites(self, global_settings) -> bool: def update_settings_for_use_of_model_based_volume_creator(self, global_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: detector_positions = np.zeros((self.number_detector_elements, 3)) 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 886bb4d7..7dc7fbec 100644 --- a/simpa/core/device_digital_twins/detection_geometries/planar_array.py +++ b/simpa/core/device_digital_twins/detection_geometries/planar_array.py @@ -82,6 +82,9 @@ def check_settings_prerequisites(self, global_settings: Settings) -> bool: def update_settings_for_use_of_model_based_volume_creator(self, global_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: detector_element_positions_mm = np.zeros((self.number_detector_elements, 3)) for x in range(self.number_detector_elements_x): diff --git a/simpa/core/device_digital_twins/illumination_geometries/__init__.py b/simpa/core/device_digital_twins/illumination_geometries/__init__.py index a19fcfa9..7ccf31fa 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/__init__.py +++ b/simpa/core/device_digital_twins/illumination_geometries/__init__.py @@ -56,7 +56,10 @@ def check_settings_prerequisites(self, global_settings) -> bool: return True def update_settings_for_use_of_model_based_volume_creator(self, global_settings) -> Settings: - return global_settings + pass + + def update_settings_for_use_of_segmentation_based_volume_creator(self, global_settings: Settings): + pass def serialize(self) -> dict: serialized_device = self.__dict__ diff --git a/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_invision_illumination.py b/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_invision_illumination.py index a31f4c6a..8cd54a1a 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_invision_illumination.py +++ b/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_invision_illumination.py @@ -16,6 +16,8 @@ class MSOTInVisionIlluminationGeometry(IlluminationGeometryBase): def __init__(self, invision_position: list = [0, 0, 0]): super().__init__() + self.invision_position = invision_position + det_sep_half = 24.74 / 2 detector_iso_distance = 74.05 / 2 detector_width = 2 * 6.12 diff --git a/simpa/core/device_digital_twins/pa_devices/__init__.py b/simpa/core/device_digital_twins/pa_devices/__init__.py index 38b0d202..b35366ef 100644 --- a/simpa/core/device_digital_twins/pa_devices/__init__.py +++ b/simpa/core/device_digital_twins/pa_devices/__init__.py @@ -4,6 +4,8 @@ import numpy as np from abc import ABC + +from simpa import Settings from simpa.core.device_digital_twins import DigitalDeviceTwinBase @@ -141,6 +143,9 @@ def check_settings_prerequisites(self, global_settings) -> bool: def update_settings_for_use_of_model_based_volume_creator(self, global_settings): pass + 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 = {"PhotoacousticDevice": serialized_device} From c3564f02a8228131bc4b3381d382f660a1ddb22a Mon Sep 17 00:00:00 2001 From: Friso Grace <155451092+frisograce@users.noreply.github.com> Date: Fri, 12 Jul 2024 16:07:22 +0200 Subject: [PATCH 04/12] Updateremove accidental update to tags --- simpa/utils/tags.py | 19 ------------------- 1 file changed, 19 deletions(-) diff --git a/simpa/utils/tags.py b/simpa/utils/tags.py index 8e932c4d..39f3c4f4 100644 --- a/simpa/utils/tags.py +++ b/simpa/utils/tags.py @@ -1496,22 +1496,3 @@ class Tags: """ Identifier for the environment varibale that defines the path the the matlab executable. """ - - ADD_US_GEL = "add_us_gel" - """ - Identifier to add ultrasound gel to the segmentation - Usage:simpa.core.digital_device_twins.pa_devices.ithera_msot_acuity - """ - - ADD_MEDIPRENE = "add_mediprene" - """ - Identifier to add mediprene to the segmentation - Usage:simpa.core.digital_device_twins.pa_devices.ithera_msot_acuity - """ - - ADD_HEAVY_WATER = "add_heavy_water" - """ - Identifier to add heavy water to the segmentation - Usage:simpa.core.digital_device_twins.pa_devices.ithera_msot_acuity - - """ From 72fada557d96141df89e3a6fcbef3767c4e4e314 Mon Sep 17 00:00:00 2001 From: Friso Grace Date: Fri, 12 Jul 2024 16:18:13 +0200 Subject: [PATCH 05/12] oops shouldnt have removed the tags... --- simpa/utils/tags.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/simpa/utils/tags.py b/simpa/utils/tags.py index 39f3c4f4..308559e3 100644 --- a/simpa/utils/tags.py +++ b/simpa/utils/tags.py @@ -1496,3 +1496,21 @@ class Tags: """ Identifier for the environment varibale that defines the path the the matlab executable. """ + + ADD_US_GEL = "add_us_gel" + """ + Identifier to add ultrasound gel to the segmentation + Usage:simpa.core.digital_device_twins.pa_devices.ithera_msot_acuity + """ + + ADD_MEDIPRENE = "add_mediprene" + """ + Identifier to add mediprene to the segmentation + Usage:simpa.core.digital_device_twins.pa_devices.ithera_msot_acuity + """ + + ADD_HEAVY_WATER = "add_heavy_water" + """ + Identifier to add heavy water to the segmentation + Usage:simpa.core.digital_device_twins.pa_devices.ithera_msot_acuity + """ From 271d029b299879a39daecb55bed9ec58fc5bdef8 Mon Sep 17 00:00:00 2001 From: Friso Grace Date: Mon, 15 Jul 2024 10:06:57 +0200 Subject: [PATCH 06/12] Extra logging fo the new feature --- .../detection_geometries/linear_array.py | 3 +-- .../pa_devices/ithera_msot_acuity.py | 11 ++++++++--- simpa/utils/tags.py | 6 ++++++ 3 files changed, 15 insertions(+), 5 deletions(-) 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 b3a1821c..c1f91902 100644 --- a/simpa/core/device_digital_twins/detection_geometries/linear_array.py +++ b/simpa/core/device_digital_twins/detection_geometries/linear_array.py @@ -4,9 +4,8 @@ import numpy as np from typing import Union -from simpa import Settings from simpa.core.device_digital_twins import DetectionGeometryBase -from simpa.utils import Tags +from simpa.utils import Settings, Tags class LinearArrayDetectionGeometry(DetectionGeometryBase): 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 c0fd3052..4592116c 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 @@ -6,8 +6,6 @@ from simpa.core.device_digital_twins import PhotoacousticDevice, \ CurvedArrayDetectionGeometry, MSOTAcuityIlluminationGeometry -from simpa.core.device_digital_twins.pa_devices import PhotoacousticDevice -from simpa.core.device_digital_twins.detection_geometries.curved_array import CurvedArrayDetectionGeometry from simpa.utils.settings import Settings from simpa.utils import Tags from simpa.utils.libraries.tissue_library import TISSUE_LIBRARY @@ -246,7 +244,10 @@ def update_settings_for_use_of_segmentation_based_volume_creator(self, global_se padding_dims = ((0, 0), (0, 0), (us_gel_thickness_pix, 0)) segmentation_map = np.pad(segmentation_map, padding_dims, mode='constant', constant_values=64) segmentation_class_mapping[64] = TISSUE_LIBRARY.ultrasound_gel() + # Whilst it may seem strange, to avoid rounding differences through the pipline it is best to let the + # z dimension shift us the same rounding as the pixel thickness. (This is the same for all three layers) z_dim_position_shift_mm += us_gel_thickness_pix * spacing_mm + self.logger.debug("Added a ultrasound gel layer to the segmentation map.") if Tags.ADD_MEDIPRENE in add_layers: mediprene_layer_height_mm = self.mediprene_membrane_height_mm @@ -255,6 +256,7 @@ def update_settings_for_use_of_segmentation_based_volume_creator(self, global_se segmentation_map = np.pad(segmentation_map, padding_dims, mode='constant', constant_values=128) segmentation_class_mapping[128] = TISSUE_LIBRARY.mediprene() z_dim_position_shift_mm += mediprene_layer_height_pix * spacing_mm + self.logger.debug("Added a mediprene layer to the segmentation map.") if Tags.ADD_HEAVY_WATER in add_layers: if heavy_water_tag is None: @@ -268,6 +270,8 @@ def update_settings_for_use_of_segmentation_based_volume_creator(self, global_se segmentation_map = np.pad(segmentation_map, padding_dims, mode='constant', constant_values=heavy_water_tag) segmentation_class_mapping[heavy_water_tag] = TISSUE_LIBRARY.heavy_water() z_dim_position_shift_mm += heavy_water_layer_height_pix * spacing_mm + self.logger.debug(f"Added a {heavy_water_layer_height_pix * spacing_mm}mm heavy water layer to the" + f"segmentation map.") new_volume_height_mm = global_settings[Tags.DIM_VOLUME_Z_MM] + z_dim_position_shift_mm @@ -289,7 +293,8 @@ def update_settings_for_use_of_segmentation_based_volume_creator(self, global_se padding_width = ((left_shift_pixels, right_shift_pixels), (0, 0), (0, 0)) segmentation_map = np.pad(segmentation_map, 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]}") + 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 diff --git a/simpa/utils/tags.py b/simpa/utils/tags.py index 308559e3..030b7aa1 100644 --- a/simpa/utils/tags.py +++ b/simpa/utils/tags.py @@ -1497,6 +1497,12 @@ class Tags: Identifier for the environment varibale that defines the path the the matlab executable. """ + VOLUME_FRACTION = "volume_fraction" + """ + Identifier for the volume fraction for the simulation + Usage: simpa.core.digital_device_twins.pa_devices.ithera_msot_acuity + """ + ADD_US_GEL = "add_us_gel" """ Identifier to add ultrasound gel to the segmentation From 766e36968ef8eb8e1e3e03e015eaff9479795e61 Mon Sep 17 00:00:00 2001 From: Friso Grace Date: Tue, 16 Jul 2024 10:18:40 +0200 Subject: [PATCH 07/12] Pep8 updates --- simpa/core/device_digital_twins/__init__.py | 10 +++++----- simpa/core/device_digital_twins/pa_devices/__init__.py | 2 +- .../pa_devices/ithera_msot_acuity.py | 2 +- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/simpa/core/device_digital_twins/__init__.py b/simpa/core/device_digital_twins/__init__.py index c34d1760..cb5c5f3b 100644 --- a/simpa/core/device_digital_twins/__init__.py +++ b/simpa/core/device_digital_twins/__init__.py @@ -86,12 +86,12 @@ def update_settings_for_use_of_model_based_volume_creator(self, global_settings: @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. + 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 - """ + :param global_settings: Settings for the entire simulation pipeline. + :type global_settings: Settings + """ def get_field_of_view_mm(self) -> np.ndarray: """ diff --git a/simpa/core/device_digital_twins/pa_devices/__init__.py b/simpa/core/device_digital_twins/pa_devices/__init__.py index b35366ef..2755ea88 100644 --- a/simpa/core/device_digital_twins/pa_devices/__init__.py +++ b/simpa/core/device_digital_twins/pa_devices/__init__.py @@ -143,7 +143,7 @@ def check_settings_prerequisites(self, global_settings) -> bool: def update_settings_for_use_of_model_based_volume_creator(self, global_settings): pass - def update_settings_for_use_of_segmentation_based_volume_creator(self, global_settings: Settings): + def update_settings_for_use_of_segmentation_based_volume_creator(self, global_settings): pass def serialize(self) -> dict: 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 4592116c..4a6bf23d 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 @@ -206,7 +206,7 @@ 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, + def update_settings_for_use_of_segmentation_based_volume_creator(self, global_settings, add_layers: list = [Tags.ADD_US_GEL, Tags.ADD_MEDIPRENE, Tags.ADD_HEAVY_WATER], From 20714d625d008e89a971ef6cc6e4e0b4065a12f0 Mon Sep 17 00:00:00 2001 From: Friso Grace Date: Wed, 17 Jul 2024 15:41:09 +0200 Subject: [PATCH 08/12] Add diffusive model for oxy --- docs/source/simpa.utils.libraries.rst | 6 +++ .../libraries/heterogeneity_generator.py | 39 +++++++++++++++++++ 2 files changed, 45 insertions(+) diff --git a/docs/source/simpa.utils.libraries.rst b/docs/source/simpa.utils.libraries.rst index fdff5eef..5ec913e7 100644 --- a/docs/source/simpa.utils.libraries.rst +++ b/docs/source/simpa.utils.libraries.rst @@ -14,6 +14,12 @@ libraries simpa.utils.libraries.scattering_spectra_data simpa.utils.libraries.structure_library +.. automodule:: simpa.utils.libraries.diffusion model + :members: + :undoc-members: + :show-inheritance: + + .. automodule:: simpa.utils.libraries.heterogeneity_generator :members: :undoc-members: diff --git a/simpa/utils/libraries/heterogeneity_generator.py b/simpa/utils/libraries/heterogeneity_generator.py index 8d01b3f1..5537f37b 100644 --- a/simpa/utils/libraries/heterogeneity_generator.py +++ b/simpa/utils/libraries/heterogeneity_generator.py @@ -3,12 +3,14 @@ # SPDX-License-Identifier: MIT import numpy as np +from rich import measure from sklearn.datasets import make_blobs from scipy.ndimage.filters import gaussian_filter from skimage import transform from simpa.utils import Tags from typing import Union, Optional from simpa.log import Logger +from simpa.utils import Settings class HeterogeneityGeneratorBase(object): @@ -151,6 +153,43 @@ def __init__(self, xdim, ydim, zdim, spacing_mm, num_centers=None, cluster_std=N self.map = gaussian_filter(self.map, 5) +class SpecificDiffusiveBlobs(object): + """ + Here the aim is to create diffusion blobs based on the positions of the sources/sinks. + """ + + def __init__(self, xdim: int, ydim: int, zdim: int, spacing_mm: Union[int, float], segmentation_area: np.ndarray, + segmentation_mapping: dict, segments_to_include: list, diffusion_const: Union[int, float], + n_steps: int = 10000, ): + ds2 = spacing_mm**2 + dt = ds2 ** 2 / (4 * diffusion_const * ds2) + + initial_conditions = np.zeros(segmentation_area.shape) + for segment in segments_to_include: + initial_conditions[np.where(segmentation_area == segment)] = segmentation_mapping[segment] + + u0 = initial_conditions.copy() + u0[np.where(initial_conditions == 0)] = 0.7 + u = u0.copy() + + for m in range(n_steps): + u0, u = self.do_timestep(u0, u) + + self.map = u0 + + def do_timestep(self, u0, u, initial_conditions, D, dt, ds2): + # Propagate with forward-difference in time, central-difference in space + u[1:-1, 1:-1] = u0[1:-1, 1:-1] + D * dt * ( + (u0[2:, 1:-1] - 2 * u0[1:-1, 1:-1] + u0[:-2, 1:-1]) / ds2 + + (u0[1:-1, 2:] - 2 * u0[1:-1, 1:-1] + u0[1:-1, :-2]) / ds2) + u[np.where(initial_conditions != 0)] = initial_conditions[np.where(initial_conditions != 0)] + u0 = u.copy() + return u0, u + + def get_map(self): + return self.map.astype(float) + + class ImageHeterogeneity(HeterogeneityGeneratorBase): """ This heterogeneity generator takes a pre-specified 2D image, currently only supporting numpy arrays, and uses them From e0cd50b4c545ee22144579503d45dbfcb0866c02 Mon Sep 17 00:00:00 2001 From: Friso Grace Date: Thu, 18 Jul 2024 14:22:10 +0200 Subject: [PATCH 09/12] Add diffusive model for oxy, using 2d and 3d models --- simpa/utils/__init__.py | 2 + .../libraries/heterogeneity_generator.py | 227 +++++++++++++++--- 2 files changed, 190 insertions(+), 39 deletions(-) diff --git a/simpa/utils/__init__.py b/simpa/utils/__init__.py index 2393f806..0ecfa5d5 100644 --- a/simpa/utils/__init__.py +++ b/simpa/utils/__init__.py @@ -16,6 +16,8 @@ from .libraries.heterogeneity_generator import RandomHeterogeneity from .libraries.heterogeneity_generator import BlobHeterogeneity from .libraries.heterogeneity_generator import ImageHeterogeneity +from .libraries.heterogeneity_generator import OxygenationDiffusionMap3D +from .libraries.heterogeneity_generator import OxygenationDiffusionMap2D # Then load classes and methods with an increasing amount of internal dependencies. # If there are import errors in the tests, it is probably due to an incorrect diff --git a/simpa/utils/libraries/heterogeneity_generator.py b/simpa/utils/libraries/heterogeneity_generator.py index 5537f37b..13e1fb57 100644 --- a/simpa/utils/libraries/heterogeneity_generator.py +++ b/simpa/utils/libraries/heterogeneity_generator.py @@ -2,15 +2,15 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT +from abc import abstractmethod, ABC import numpy as np -from rich import measure from sklearn.datasets import make_blobs from scipy.ndimage.filters import gaussian_filter from skimage import transform from simpa.utils import Tags from typing import Union, Optional from simpa.log import Logger -from simpa.utils import Settings +from simpa.utils.calculate import calculate_oxygenation class HeterogeneityGeneratorBase(object): @@ -153,43 +153,6 @@ def __init__(self, xdim, ydim, zdim, spacing_mm, num_centers=None, cluster_std=N self.map = gaussian_filter(self.map, 5) -class SpecificDiffusiveBlobs(object): - """ - Here the aim is to create diffusion blobs based on the positions of the sources/sinks. - """ - - def __init__(self, xdim: int, ydim: int, zdim: int, spacing_mm: Union[int, float], segmentation_area: np.ndarray, - segmentation_mapping: dict, segments_to_include: list, diffusion_const: Union[int, float], - n_steps: int = 10000, ): - ds2 = spacing_mm**2 - dt = ds2 ** 2 / (4 * diffusion_const * ds2) - - initial_conditions = np.zeros(segmentation_area.shape) - for segment in segments_to_include: - initial_conditions[np.where(segmentation_area == segment)] = segmentation_mapping[segment] - - u0 = initial_conditions.copy() - u0[np.where(initial_conditions == 0)] = 0.7 - u = u0.copy() - - for m in range(n_steps): - u0, u = self.do_timestep(u0, u) - - self.map = u0 - - def do_timestep(self, u0, u, initial_conditions, D, dt, ds2): - # Propagate with forward-difference in time, central-difference in space - u[1:-1, 1:-1] = u0[1:-1, 1:-1] + D * dt * ( - (u0[2:, 1:-1] - 2 * u0[1:-1, 1:-1] + u0[:-2, 1:-1]) / ds2 - + (u0[1:-1, 2:] - 2 * u0[1:-1, 1:-1] + u0[1:-1, :-2]) / ds2) - u[np.where(initial_conditions != 0)] = initial_conditions[np.where(initial_conditions != 0)] - u0 = u.copy() - return u0, u - - def get_map(self): - return self.map.astype(float) - - class ImageHeterogeneity(HeterogeneityGeneratorBase): """ This heterogeneity generator takes a pre-specified 2D image, currently only supporting numpy arrays, and uses them @@ -319,6 +282,9 @@ def crop_image(self, xdim: int, zdim: int, image: np.ndarray, crop_horizontal = round((image_width_pixels - crop_width) / 2) elif isinstance(crop_placement[0], int): crop_horizontal = crop_placement[0] + else: + raise ValueError(f'Unrecognised crop placement {crop_placement[0]}, please see the available options in' + f'the documentation') if crop_placement[1] == Tags.CROP_POSITION_TOP: crop_vertical = 0 @@ -328,6 +294,9 @@ def crop_image(self, xdim: int, zdim: int, image: np.ndarray, crop_vertical = round((image_height_pixels - crop_height) / 2) elif isinstance(crop_placement[1], int): crop_vertical = crop_placement[1] + else: + raise ValueError(f'Unrecognised crop placement {crop_placement[1]}, please see the available options in' + f'the documentation') elif isinstance(crop_placement, str): if crop_placement == Tags.CROP_POSITION_CENTRE: @@ -340,6 +309,12 @@ def crop_image(self, xdim: int, zdim: int, image: np.ndarray, crop_vertical = image_height_pixels - crop_height if crop_vertical != 0: crop_vertical = np.random.randint(0, crop_vertical) + else: + raise ValueError(f'Unrecognised crop placement {crop_placement}, please see the available options in' + f'the documentation') + else: + raise ValueError(f'Unrecognised crop placement {crop_placement}, please see the available options in' + f'the documentation') cropped_image = image[crop_horizontal: crop_horizontal + crop_width, crop_vertical: crop_vertical + crop_height] @@ -364,3 +339,177 @@ def change_resolution(self, image: np.ndarray, spacing_mm: Union[int, float], self.logger.warning( "The input image has changed pixel spacing to {} to match the simulation volume".format(spacing_mm)) return transform.resize(image, (new_image_pixel_width, new_image_pixel_height)) + + +class OxygenationDiffusionMap(object): + """ + Here the aim is to create diffusion maps for oxygenation based on the positions and values of pre-set of the + sources/sinks. This predefined oxygenation map must be 2D. The user should understand that this technique requires + knowing when the diffusion equations have approached a stable point. This works using a finite difference method on + the diffusion equations. + + Some suggestions for how to optimise the time of convergence: + - We recommend that the user try a number of time steps to begin with to gauge reasonable values for your + simulation. + - We recommend that the user amends the diffusion constant. At its essence, this value determines the size of + the steps in this numerical method. Large diffusion constant will lead to faster convergence, but will reduce + the accuracy + - We recommend caution around the value decide for the background_oxygenation. This value should represent what + you expect the long term oxygenation to look like, as this will reduce the amount each value needs to change in + order to reach ots converged value. In addition, this value will be used as the boundary. It represents + condition which must be fulfilled throughout the simulation and will therefore have a large baring on the final + results + + Attributes: + map: The final oxygenation map + """ + + def __init__(self, segmentation_array: np.ndarray, segmentation_mapping: dict, + segments_to_include: list, diffusion_const: Union[int, float], dt: Union[int, float], + ds2: Union[int, float], n_steps: int, background_oxygenation: Union[int, float]): + """ + Upon initialisation of this class, the diffusion map will be created. + :param segmentation_array: An array segmenting the medium + :param segmentation_mapping: A mapping from the segmentations to the oxygenations of the blood vessels + :param segments_to_include: Which segments to include in the initial oxygenation map + :param diffusion_const: The diffusion constant of the oxygenation + :param n_steps: The number of steps the simulation should use + :param background_oxygenation: The oxygenation that will be used as the starting value everywhere where the + oxygenation is not defined, and as the edge value + """ + initial_conditions = np.zeros(segmentation_array.shape) + for segment in segments_to_include: + initial_conditions[np.where(segmentation_array == segment) + ] = calculate_oxygenation(segmentation_mapping[segment]) + + oxygenation = initial_conditions.copy() + oxygenation[np.where(initial_conditions == 0)] = background_oxygenation + + for m in range(n_steps): + oxygenation = self.diffusion_timestep(oxygenation, initial_conditions, diffusion_const, dt, ds2) + + self.map = oxygenation + + @abstractmethod + def diffusion_timestep(self, u, boundary_conditions, diffusion_const, dt, ds2): + """ + Here, a diffusion equation can be applied to determine how the oxygenation should evolve. + :param u: The initial map to perform a time step of the diffusion model + :param boundary_conditions: The conditions that you wish to remain constant throughout the simulation + :param diffusion_const: The diffusion constant of the oxygenation + :param dt: The time step + :param ds2: The second order spatial step + :return: The map after a time step of the diffusion model + """ + pass + + def get_map(self): + return self.map + + +class OxygenationDiffusionMap2D(OxygenationDiffusionMap): + """ + This uses the OxygenationDiffusionMap class and applies it to a volume with a constant cross-sectoin. The predefined + oxygenation map must be 2D. + + Attributes: + map: The final oxygenation map + """ + + def __init__(self, spacing_mm, segmentation_area: np.ndarray, segmentation_mapping, segments_to_include, + diffusion_const, ydim: int, n_steps=10000, background_oxygenation=0.8): + """ + Upon initialisation of this class, the diffusion map will be created. + :param spacing_mm: the spacing of the volume in mm + :param segmentation_area: An array segmenting the cross-sectional area of the medium + :param segmentation_mapping: A mapping from the segmentations to the oxygenations of the blood vessels + :param segments_to_include: Which segments to include in the initial oxygenation map + :param diffusion_const: The diffusion constant of the oxygenation + :param ydim: the dimension to scale the image to make it into a volume + :param n_steps: The number of steps the simulation should use + :param background_oxygenation: The oxygenation that will be used as the starting value everywhere where the + oxygenation is not defined, and as the edge value + """ + ds2 = spacing_mm**2 + dt = ds2 ** 2 / (4 * diffusion_const * ds2) + super().__init__(segmentation_area, segmentation_mapping, segments_to_include, diffusion_const, dt, ds2, + n_steps, background_oxygenation) + + self.map = np.repeat(self.map[:, np.newaxis, :], ydim, axis=1) + + def diffusion_timestep(self, u, boundary_conditions, diffusion_const, dt, ds2): + """ + :param u: The initial map to perform a time step of the diffusion model + :param boundary_conditions: The conditions that you wish to remain constant throughout the simulation + :param diffusion_const: The diffusion constant of the oxygenation + :param dt: The time step + :param ds2: The second order spatial step + :return: The map after a time step of the diffusion model + """ + # Propagate with forward-difference in time, central-difference in space + u[1:-1, 1:-1] = u[1:-1, 1:-1] + diffusion_const * dt * ( + (u[2:, 1:-1] - 2 * u[1:-1, 1:-1] + u[:-2, 1:-1]) / ds2 + + (u[1:-1, 2:] - 2 * u[1:-1, 1:-1] + u[1:-1, :-2]) / ds2) + u[np.where(boundary_conditions != 0)] = boundary_conditions[np.where(boundary_conditions != 0)] + return u + + +class OxygenationDiffusionMap3D(OxygenationDiffusionMap): + """ + Here the aim is to create diffusion maps for oxygenation based on the positions and values of pre-set of the + sources/sinks. This predefined oxygenation map must be 3D. The user should understand that this technique requires + knowing when the diffusion equations have approached a stable point. This works using a finite difference method on + the diffusion equations. + + Some suggestions for how to optimise the time of convergence: + - We recommend that the user try a number of time steps to begin with to gauge reasonable values for your + simulation. + - We recommend that the user amends the diffusion constant. At its essence, this value determines the size of + the steps in this numerical method. Large diffusion constant will lead to faster convergence, but will reduce + the accuracy + - We recommend caution around the value decide for the background_oxygenation. This value should represent what + you expect the long term oxygenation to look like, as this will reduce the amount each value needs to change in + order to reach ots converged value. In addition, this value will be used as the boundary. It represents + condition which must be fulfilled throughout the simulation and will therefore have a large baring on the final + results + + Attributes: + map: The final oxygenation map + """ + + def __init__(self, spacing_mm: Union[int, float], segmentation_volume: np.ndarray, segmentation_mapping: dict, + segments_to_include: list, diffusion_const: Union[int, float], + n_steps: int = 10000, background_oxygenation: Union[int, float] = 0.8): + """ + Upon initialisation of this class, the diffusion map will be created. + :param spacing_mm: the spacing of the volume in mm + :param segmentation_volume: An array segmenting the volume of the medium + :param segmentation_mapping: A mapping from the segmentations to the oxygenations of the blood vessels + :param segments_to_include: Which segments to include in the initial oxygenation map + :param diffusion_const: The diffusion constant of the oxygenation + :param n_steps: The number of steps the simulation should use + :param background_oxygenation: The oxygenation that will be used: as the starting value everywhere where the + oxygenation is not defined, and as the edge value + """ + ds2 = spacing_mm**2 + dt = ds2 ** 2 / (6 * diffusion_const * ds2) + + super().__init__(segmentation_volume, segmentation_mapping, segments_to_include, diffusion_const, dt, + ds2, n_steps, background_oxygenation) + + def diffusion_timestep(self, u, boundary_conditions, diffusion_const, dt, ds2): + """ + :param u: The initial map to perform a time step of the diffusion model + :param boundary_conditions: The conditions that you wish to remain constant throughout the simulation + :param diffusion_const: The diffusion constant of the oxygenation + :param dt: The time step + :param ds2: The second order spatial step + :return: The map after a time step of the diffusion model + """ + # Propagate with forward-difference in time, central-difference in space + u[1:-1, 1:-1, 1:-1] = u[1:-1, 1:-1, 1:-1] + diffusion_const * dt * ( + (u[2:, 1:-1, 1:-1] - 2 * u[1:-1, 1:-1, 1:-1] + u[:-2, 1:-1, 1:-1]) / ds2 + + (u[1:-1, 2:, 1:-1] - 2 * u[1:-1, 1:-1, 1:-1] + u[1:-1, :-2, 1:-1]) / ds2 + + (u[1:-1, 1:-1, 2:] - 2 * u[1:-1, 1:-1, 1:-1] + u[1:-1, 1:-1, :-2]) / ds2) + u[np.where(boundary_conditions != 0)] = boundary_conditions[np.where(boundary_conditions != 0)] + return u From c92360124ab892148d640c8c1930d3ecc212d2c1 Mon Sep 17 00:00:00 2001 From: Friso Grace Date: Thu, 18 Jul 2024 14:28:17 +0200 Subject: [PATCH 10/12] Remove T348 updates --- simpa/core/device_digital_twins/__init__.py | 20 +-- .../detection_geometries/__init__.py | 9 +- .../detection_geometries/curved_array.py | 18 +-- .../detection_geometries/linear_array.py | 8 +- .../detection_geometries/planar_array.py | 6 +- .../illumination_geometries/__init__.py | 10 +- .../disk_illumination.py | 3 +- .../gaussian_beam_illumination.py | 5 +- .../ithera_msot_invision_illumination.py | 7 +- .../pencil_array_illumination.py | 4 +- .../slit_illumination.py | 2 +- .../pa_devices/__init__.py | 5 - .../pa_devices/ithera_msot_acuity.py | 139 +----------------- .../pa_devices/ithera_rsom.py | 6 +- 14 files changed, 37 insertions(+), 205 deletions(-) diff --git a/simpa/core/device_digital_twins/__init__.py b/simpa/core/device_digital_twins/__init__.py index cb5c5f3b..13bd5bcf 100644 --- a/simpa/core/device_digital_twins/__init__.py +++ b/simpa/core/device_digital_twins/__init__.py @@ -9,7 +9,6 @@ import uuid from simpa.utils.serializer import SerializableSIMPAClass from simpa.utils.calculate import are_equal -from simpa.utils import Settings class DigitalDeviceTwinBase(SerializableSIMPAClass): @@ -18,7 +17,7 @@ class DigitalDeviceTwinBase(SerializableSIMPAClass): which has representations of both. """ - def __init__(self, device_position_mm: np.ndarray = None, field_of_view_extent_mm: np.ndarray = None): + def __init__(self, device_position_mm=None, field_of_view_extent_mm=None): """ :param device_position_mm: Each device has an internal position which serves as origin for internal \ representations of e.g. detector element positions or illuminator positions. @@ -55,7 +54,7 @@ def __eq__(self, other): return False @abstractmethod - def check_settings_prerequisites(self, global_settings: Settings) -> bool: + def check_settings_prerequisites(self, global_settings) -> bool: """ It might be that certain device geometries need a certain dimensionality of the simulated PAI volume, or that it requires the existence of certain Tags in the global global_settings. @@ -73,7 +72,7 @@ def check_settings_prerequisites(self, global_settings: Settings) -> bool: pass @abstractmethod - def update_settings_for_use_of_model_based_volume_creator(self, global_settings: Settings): + def update_settings_for_use_of_model_based_volume_creator(self, global_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. @@ -83,16 +82,6 @@ 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 - """ - def get_field_of_view_mm(self) -> np.ndarray: """ Returns the absolute field of view in mm where the probe position is already @@ -144,8 +133,7 @@ def deserialize(dictionary_to_deserialize): """ -It is important to have these relative imports after the definition of the DigitalDeviceTwinBase class to avoid circular -imports triggered by imported child classes +It is important to have these relative imports after the definition of the DigitalDeviceTwinBase class to avoid circular imports triggered by imported child classes """ from .pa_devices import PhotoacousticDevice # nopep8 from simpa.core.device_digital_twins.detection_geometries import DetectionGeometryBase # nopep8 diff --git a/simpa/core/device_digital_twins/detection_geometries/__init__.py b/simpa/core/device_digital_twins/detection_geometries/__init__.py index 8afa7188..eba6bdd5 100644 --- a/simpa/core/device_digital_twins/detection_geometries/__init__.py +++ b/simpa/core/device_digital_twins/detection_geometries/__init__.py @@ -5,7 +5,6 @@ from abc import abstractmethod from simpa.core.device_digital_twins import DigitalDeviceTwinBase import numpy as np -from typing import Union class DetectionGeometryBase(DigitalDeviceTwinBase): @@ -13,10 +12,10 @@ class DetectionGeometryBase(DigitalDeviceTwinBase): This class is the base class for representing all detector geometries. """ - def __init__(self, number_detector_elements: int, detector_element_width_mm: Union[float, int], - detector_element_length_mm: Union[float, int], center_frequency_hz: Union[float, int], - bandwidth_percent: Union[float, int], sampling_frequency_mhz: Union[float, int], - device_position_mm: np.ndarray = None, field_of_view_extent_mm: np.ndarray = None): + def __init__(self, number_detector_elements, detector_element_width_mm, + detector_element_length_mm, center_frequency_hz, bandwidth_percent, + sampling_frequency_mhz, device_position_mm: np.ndarray = None, + field_of_view_extent_mm: np.ndarray = None): """ :param number_detector_elements: Total number of detector elements. 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 f281516e..c530e4f7 100644 --- a/simpa/core/device_digital_twins/detection_geometries/curved_array.py +++ b/simpa/core/device_digital_twins/detection_geometries/curved_array.py @@ -2,10 +2,9 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT import numpy as np -from typing import Union from simpa.core.device_digital_twins import DetectionGeometryBase -from simpa import Settings, Tags +from simpa.utils import Tags class CurvedArrayDetectionGeometry(DetectionGeometryBase): @@ -14,17 +13,17 @@ class CurvedArrayDetectionGeometry(DetectionGeometryBase): with a curved detection geometry. The origin for this device is the center (focus) of the curved array. """ - def __init__(self, pitch_mm: Union[float, int] = 0.5, - radius_mm: Union[float, int] = 40, - number_detector_elements: int = 256, + def __init__(self, pitch_mm=0.5, + 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=55, sampling_frequency_mhz=40, - angular_origin_offset: Union[float, int] = np.pi, - device_position_mm: np.ndarray = None, - field_of_view_extent_mm: np.ndarray = None): + angular_origin_offset=np.pi, + device_position_mm=None, + field_of_view_extent_mm=None): """ :param pitch_mm: In-plane distance between the beginning of one detector element to the next detector element. @@ -86,9 +85,6 @@ def check_settings_prerequisites(self, global_settings) -> bool: def update_settings_for_use_of_model_based_volume_creator(self, global_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: pitch_angle = self.pitch_mm / self.radius_mm 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 c1f91902..9209d5e0 100644 --- a/simpa/core/device_digital_twins/detection_geometries/linear_array.py +++ b/simpa/core/device_digital_twins/detection_geometries/linear_array.py @@ -2,7 +2,6 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT import numpy as np -from typing import Union from simpa.core.device_digital_twins import DetectionGeometryBase from simpa.utils import Settings, Tags @@ -16,7 +15,7 @@ class LinearArrayDetectionGeometry(DetectionGeometryBase): """ - def __init__(self, pitch_mm: Union[float, int] = 0.5, + def __init__(self, pitch_mm=0.5, number_detector_elements=100, detector_element_width_mm=0.24, detector_element_length_mm=0.5, @@ -52,7 +51,7 @@ def __init__(self, pitch_mm: Union[float, int] = 0.5, self.pitch_mm = pitch_mm self.probe_width_mm = (number_detector_elements - 1) * self.pitch_mm - def check_settings_prerequisites(self, global_settings) -> bool: + def check_settings_prerequisites(self, global_settings: Settings) -> bool: if global_settings[Tags.DIM_VOLUME_X_MM] < self.probe_width_mm + global_settings[Tags.SPACING_MM]: self.logger.error("Volume x dimension is too small to encompass MSOT device in simulation!" "Must be at least {} mm but was {} mm" @@ -64,9 +63,6 @@ def check_settings_prerequisites(self, global_settings) -> bool: def update_settings_for_use_of_model_based_volume_creator(self, global_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: detector_positions = np.zeros((self.number_detector_elements, 3)) 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 7dc7fbec..1fe2fb93 100644 --- a/simpa/core/device_digital_twins/detection_geometries/planar_array.py +++ b/simpa/core/device_digital_twins/detection_geometries/planar_array.py @@ -5,7 +5,6 @@ from simpa.core.device_digital_twins import DetectionGeometryBase from simpa.utils import Settings, Tags -from typing import Union class PlanarArrayDetectionGeometry(DetectionGeometryBase): @@ -15,7 +14,7 @@ class PlanarArrayDetectionGeometry(DetectionGeometryBase): """ - def __init__(self, pitch_mm: Union[float, int] = 0.5, + def __init__(self, pitch_mm=0.5, number_detector_elements_x=100, number_detector_elements_y=100, detector_element_width_mm=0.24, @@ -82,9 +81,6 @@ def check_settings_prerequisites(self, global_settings: Settings) -> bool: def update_settings_for_use_of_model_based_volume_creator(self, global_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: detector_element_positions_mm = np.zeros((self.number_detector_elements, 3)) for x in range(self.number_detector_elements_x): diff --git a/simpa/core/device_digital_twins/illumination_geometries/__init__.py b/simpa/core/device_digital_twins/illumination_geometries/__init__.py index 7ccf31fa..6c5780bf 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/__init__.py +++ b/simpa/core/device_digital_twins/illumination_geometries/__init__.py @@ -13,8 +13,7 @@ class IlluminationGeometryBase(DigitalDeviceTwinBase): This class is the base class for representing all illumination geometries. """ - def __init__(self, device_position_mm=None, source_direction_vector: np.ndarray = None, - field_of_view_extent_mm=None): + def __init__(self, device_position_mm=None, source_direction_vector=None, field_of_view_extent_mm=None): """ :param device_position_mm: Each device has an internal position which serves as origin for internal \ representations of illuminator positions. @@ -39,7 +38,7 @@ def __init__(self, device_position_mm=None, source_direction_vector: np.ndarray self.source_direction_vector) @abstractmethod - def get_mcx_illuminator_definition(self, global_settings: Settings) -> dict: + def get_mcx_illuminator_definition(self, global_settings) -> dict: """ IMPORTANT: This method creates a dictionary that contains tags as they are expected for the mcx simulation tool to represent the illumination geometry of this device. @@ -56,10 +55,7 @@ def check_settings_prerequisites(self, global_settings) -> bool: return True 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 + return global_settings def serialize(self) -> dict: serialized_device = self.__dict__ diff --git a/simpa/core/device_digital_twins/illumination_geometries/disk_illumination.py b/simpa/core/device_digital_twins/illumination_geometries/disk_illumination.py index 51e2eddc..42b12909 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/disk_illumination.py +++ b/simpa/core/device_digital_twins/illumination_geometries/disk_illumination.py @@ -4,7 +4,6 @@ from simpa.core.device_digital_twins import IlluminationGeometryBase from simpa.utils import Tags -from typing import Union class DiskIlluminationGeometry(IlluminationGeometryBase): @@ -13,7 +12,7 @@ class DiskIlluminationGeometry(IlluminationGeometryBase): The device position is defined as the middle of the disk. """ - def __init__(self, beam_radius_mm: Union[int, float] = None, device_position_mm=None, field_of_view_extent_mm=None): + def __init__(self, beam_radius_mm=None, device_position_mm=None, field_of_view_extent_mm=None): """ :param beam_radius_mm: Radius of the disk in mm. :type beam_radius_mm: int, float diff --git a/simpa/core/device_digital_twins/illumination_geometries/gaussian_beam_illumination.py b/simpa/core/device_digital_twins/illumination_geometries/gaussian_beam_illumination.py index 1f0a14dd..b3a9b8a6 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/gaussian_beam_illumination.py +++ b/simpa/core/device_digital_twins/illumination_geometries/gaussian_beam_illumination.py @@ -6,7 +6,6 @@ from collections.abc import Sized from simpa.core.device_digital_twins import IlluminationGeometryBase from simpa.utils import Tags -from typing import Union class GaussianBeamIlluminationGeometry(IlluminationGeometryBase): @@ -15,8 +14,8 @@ class GaussianBeamIlluminationGeometry(IlluminationGeometryBase): The position is defined as the middle of the beam. """ - def __init__(self, beam_radius_mm: Union[int, float] = None, focal_length_mm: Union[int, float] = None, - device_position_mm=None, field_of_view_extent_mm=None): + def __init__(self, beam_radius_mm=None, focal_length_mm=None, device_position_mm=None, + field_of_view_extent_mm=None): """ :param beam_radius_mm: Initial radius of the gaussian beam at half maximum (full width at half maximum (FWHM)) in mm. diff --git a/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_invision_illumination.py b/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_invision_illumination.py index 8cd54a1a..0c22a346 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_invision_illumination.py +++ b/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_invision_illumination.py @@ -13,10 +13,13 @@ class MSOTInVisionIlluminationGeometry(IlluminationGeometryBase): This class represents the illumination geometry of the MSOT InVision photoacoustic device. """ - def __init__(self, invision_position: list = [0, 0, 0]): + def __init__(self, invision_position=None): super().__init__() - self.invision_position = invision_position + if invision_position is None: + self.invision_position = [0, 0, 0] + else: + self.invision_position = invision_position det_sep_half = 24.74 / 2 detector_iso_distance = 74.05 / 2 diff --git a/simpa/core/device_digital_twins/illumination_geometries/pencil_array_illumination.py b/simpa/core/device_digital_twins/illumination_geometries/pencil_array_illumination.py index 0712bc86..e5a4deab 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/pencil_array_illumination.py +++ b/simpa/core/device_digital_twins/illumination_geometries/pencil_array_illumination.py @@ -14,8 +14,8 @@ class PencilArrayIlluminationGeometry(IlluminationGeometryBase): The device position is defined as the middle of the array. """ - def __init__(self, pitch_mm: float = 0.5, number_illuminators_x: int = 100, number_illuminators_y: int = 100, - device_position_mm=None, field_of_view_extent_mm=None): + def __init__(self, pitch_mm=0.5, number_illuminators_x=100, number_illuminators_y=100, device_position_mm=None, + field_of_view_extent_mm=None): """ :param pitch_mm: Defines the x and y distance between the illumination positions :type pitch_mm: float diff --git a/simpa/core/device_digital_twins/illumination_geometries/slit_illumination.py b/simpa/core/device_digital_twins/illumination_geometries/slit_illumination.py index 87b05bf0..4dc58eba 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/slit_illumination.py +++ b/simpa/core/device_digital_twins/illumination_geometries/slit_illumination.py @@ -14,7 +14,7 @@ class SlitIlluminationGeometry(IlluminationGeometryBase): The device position is defined as the middle of the slit. """ - def __init__(self, slit_vector_mm: list = None, direction_vector_mm: list = None, device_position_mm=None, + def __init__(self, slit_vector_mm=None, direction_vector_mm=None, device_position_mm=None, field_of_view_extent_mm=None): """ :param slit_vector_mm: Defines the slit in vector form. For example a slit along the x-axis with length 5mm diff --git a/simpa/core/device_digital_twins/pa_devices/__init__.py b/simpa/core/device_digital_twins/pa_devices/__init__.py index 2755ea88..38b0d202 100644 --- a/simpa/core/device_digital_twins/pa_devices/__init__.py +++ b/simpa/core/device_digital_twins/pa_devices/__init__.py @@ -4,8 +4,6 @@ import numpy as np from abc import ABC - -from simpa import Settings from simpa.core.device_digital_twins import DigitalDeviceTwinBase @@ -143,9 +141,6 @@ def check_settings_prerequisites(self, global_settings) -> bool: def update_settings_for_use_of_model_based_volume_creator(self, global_settings): pass - def update_settings_for_use_of_segmentation_based_volume_creator(self, global_settings): - pass - def serialize(self) -> dict: serialized_device = self.__dict__ device_dict = {"PhotoacousticDevice": 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 0fa4f466..749e0f19 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 @@ -7,6 +7,8 @@ from simpa.core.device_digital_twins import PhotoacousticDevice, \ CurvedArrayDetectionGeometry, MSOTAcuityIlluminationGeometry +from simpa.core.device_digital_twins.pa_devices import PhotoacousticDevice +from simpa.core.device_digital_twins.detection_geometries.curved_array import CurvedArrayDetectionGeometry from simpa.utils.settings import Settings from simpa.utils import Tags from simpa.utils.libraries.tissue_library import TISSUE_LIBRARY @@ -217,143 +219,6 @@ 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, - add_layers: list = [Tags.ADD_US_GEL, - Tags.ADD_MEDIPRENE, - Tags.ADD_HEAVY_WATER], - current_heavy_water_depth: (float, int) = 0, - heavy_water_tag: int = None): - """ - Updates the volume creation settings of the segmentation based volume creator according to the size of the - device. On the occasion that your segmentation already includes the mediprene, ultrasound gel and some of the - heavy water, you may specify the existing depth of the heavy water so that it can be adapted to the depth of - device. - :param add_layers: The layers to add to the volume, all configured to the typical thicknesses for MSOT acuity - echo. - :param current_heavy_water_depth: the current heavy water depth (mm). - :param heavy_water_tag: the existing heavy water tag in the segmentation map. - :param global_settings: Settings for the entire simulation pipeline. - """ - try: - volume_creator_settings = Settings(global_settings.get_volume_creation_settings()) - except KeyError as e: - 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 - mediprene_layer_height_mm = 0 - heavy_water_layer_height_mm = 0 - - if Tags.ADD_US_GEL in add_layers: - us_gel_thickness_mm = np.random.normal(0.4, 0.1) - us_gel_thickness_pix = int(round(us_gel_thickness_mm/spacing_mm)) - padding_dims = ((0, 0), (0, 0), (us_gel_thickness_pix, 0)) - segmentation_map = np.pad(segmentation_map, padding_dims, mode='constant', constant_values=64) - segmentation_class_mapping[64] = TISSUE_LIBRARY.ultrasound_gel() - # Whilst it may seem strange, to avoid rounding differences through the pipline it is best to let the - # z dimension shift us the same rounding as the pixel thickness. (This is the same for all three layers) - z_dim_position_shift_mm += us_gel_thickness_pix * spacing_mm - self.logger.debug("Added a ultrasound gel layer to the segmentation map.") - - if Tags.ADD_MEDIPRENE in add_layers: - mediprene_layer_height_mm = self.mediprene_membrane_height_mm - mediprene_layer_height_pix = int(round(mediprene_layer_height_mm/spacing_mm)) - padding_dims = ((0, 0), (0, 0), (mediprene_layer_height_pix, 0)) - segmentation_map = np.pad(segmentation_map, padding_dims, mode='constant', constant_values=128) - segmentation_class_mapping[128] = TISSUE_LIBRARY.mediprene() - z_dim_position_shift_mm += mediprene_layer_height_pix * spacing_mm - self.logger.debug("Added a mediprene layer to the segmentation map.") - - if Tags.ADD_HEAVY_WATER in add_layers: - if heavy_water_tag is None: - heavy_water_tag = 256 - segmentation_class_mapping[256] = TISSUE_LIBRARY.heavy_water() - probe_size_mm = self.probe_height_mm - mediprene_layer_height_mm = self.mediprene_membrane_height_mm - heavy_water_layer_height_mm = probe_size_mm - current_heavy_water_depth - mediprene_layer_height_mm - heavy_water_layer_height_pix = int(round(heavy_water_layer_height_mm / spacing_mm)) - padding_dims = ((0, 0), (0, 0), (heavy_water_layer_height_pix, 0)) - segmentation_map = np.pad(segmentation_map, padding_dims, mode='constant', constant_values=heavy_water_tag) - segmentation_class_mapping[heavy_water_tag] = TISSUE_LIBRARY.heavy_water() - z_dim_position_shift_mm += heavy_water_layer_height_pix * spacing_mm - self.logger.debug(f"Added a {heavy_water_layer_height_pix * spacing_mm}mm heavy water layer to the" - f"segmentation map.") - - new_volume_height_mm = global_settings[Tags.DIM_VOLUME_Z_MM] + z_dim_position_shift_mm - - # adjust the z-dim to msot probe height - global_settings[Tags.DIM_VOLUME_Z_MM] = new_volume_height_mm - self.logger.debug(f"Changed Tags.DIM_VOLUME_Z_MM to {global_settings[Tags.DIM_VOLUME_Z_MM]}") - - # adjust the x-dim to msot probe width - # 1 voxel is added (0.5 on both sides) to make sure no rounding errors lead to a detector element being outside - # of the simulated volume. - - 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 - # specific left and right to avoid rounding errors - left_shift_pixels = int(round(width_shift_for_structures_mm / spacing_mm)) - right_shift_pixels = int(round((round(self.detection_geometry.probe_width_mm) + spacing_mm - - global_settings[Tags.DIM_VOLUME_X_MM])/spacing_mm)) - left_shift_pixels - padding_width = ((left_shift_pixels, right_shift_pixels), (0, 0), (0, 0)) - segmentation_map = np.pad(segmentation_map, 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 - padding_width = ((0, 0), (0, 0), (0, 0)) - - 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") - - for structure_key in volume_creator_settings[Tags.SEGMENTATION_CLASS_MAPPING]: - self.logger.debug("Adjusting " + str(structure_key)) - structure_dict = volume_creator_settings[Tags.SEGMENTATION_CLASS_MAPPING][structure_key] - for molecule in structure_dict: - try: - old_volume_fraction = getattr(molecule, Tags.VOLUME_FRACTION) - except AttributeError: - continue - if isinstance(old_volume_fraction, torch.Tensor): - if old_volume_fraction.shape != segmentation_map.shape: - z_shift_pixels = int(round(z_dim_position_shift_mm / spacing_mm)) - padding_height = ((0, 0), (0, 0), (z_shift_pixels, 0)) - padded_up = np.pad(old_volume_fraction.numpy(), padding_height, mode='edge') - padded_vol = np.pad(padded_up, padding_width, mode='edge') - setattr(molecule, Tags.VOLUME_FRACTION, torch.tensor(padded_vol, dtype=torch.float32)) - - device_change_in_height = mediprene_layer_height_mm + heavy_water_layer_height_mm - self.device_position_mm = np.add(self.device_position_mm, np.array([width_shift_for_structures_mm, 0, - device_change_in_height])) - 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=55, - 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, - device_change_in_height])) - def serialize(self) -> dict: serialized_device = self.__dict__ device_dict = {"MSOTAcuityEcho": serialized_device} diff --git a/simpa/core/device_digital_twins/pa_devices/ithera_rsom.py b/simpa/core/device_digital_twins/pa_devices/ithera_rsom.py index 152eb27f..eba16ba3 100644 --- a/simpa/core/device_digital_twins/pa_devices/ithera_rsom.py +++ b/simpa/core/device_digital_twins/pa_devices/ithera_rsom.py @@ -35,9 +35,9 @@ class RSOMExplorerP50(PhotoacousticDevice): """ - def __init__(self, element_spacing_mm: float = 0.02, - number_elements_x: int = 10, - number_elements_y: int = 10, + def __init__(self, element_spacing_mm=0.02, + number_elements_x=10, + number_elements_y=10, device_position_mm: np.ndarray = None, field_of_view_extent_mm: np.ndarray = None): """ From 913259a21db8e45de4f2bf193dae60dc12c34283 Mon Sep 17 00:00:00 2001 From: Friso Grace Date: Thu, 18 Jul 2024 14:31:22 +0200 Subject: [PATCH 11/12] Remove old diffusion model --- docs/source/simpa.utils.libraries.rst | 6 ------ 1 file changed, 6 deletions(-) diff --git a/docs/source/simpa.utils.libraries.rst b/docs/source/simpa.utils.libraries.rst index 5ec913e7..fdff5eef 100644 --- a/docs/source/simpa.utils.libraries.rst +++ b/docs/source/simpa.utils.libraries.rst @@ -14,12 +14,6 @@ libraries simpa.utils.libraries.scattering_spectra_data simpa.utils.libraries.structure_library -.. automodule:: simpa.utils.libraries.diffusion model - :members: - :undoc-members: - :show-inheritance: - - .. automodule:: simpa.utils.libraries.heterogeneity_generator :members: :undoc-members: From 2c33d23a5b733a83008efa3b1765847017d76c23 Mon Sep 17 00:00:00 2001 From: Friso Grace Date: Thu, 18 Jul 2024 14:33:16 +0200 Subject: [PATCH 12/12] Remove unused Tags --- simpa/utils/tags.py | 19 ------------------- 1 file changed, 19 deletions(-) diff --git a/simpa/utils/tags.py b/simpa/utils/tags.py index b2607a21..dee642b6 100644 --- a/simpa/utils/tags.py +++ b/simpa/utils/tags.py @@ -1569,23 +1569,4 @@ class Tags: VOLUME_FRACTION = "volume_fraction" """ Identifier for the volume fraction for the simulation - Usage: simpa.core.digital_device_twins.pa_devices.ithera_msot_acuity - """ - - ADD_US_GEL = "add_us_gel" - """ - Identifier to add ultrasound gel to the segmentation - Usage:simpa.core.digital_device_twins.pa_devices.ithera_msot_acuity - """ - - ADD_MEDIPRENE = "add_mediprene" - """ - Identifier to add mediprene to the segmentation - Usage:simpa.core.digital_device_twins.pa_devices.ithera_msot_acuity - """ - - ADD_HEAVY_WATER = "add_heavy_water" - """ - Identifier to add heavy water to the segmentation - Usage:simpa.core.digital_device_twins.pa_devices.ithera_msot_acuity """