Merge remote-tracking branch 'refs/remotes/origin/develop' into T323_…

…clean_logging

# Conflicts:
#	docs/source/simpa.core.device_digital_twins.illumination_geometries.rst
#	docs/source/simpa_examples.rst

docs/source/clean_up_rst_files.py

@@ -42,7 +42,7 @@
 simpa_examples_rst_file = open(os.path.join(current_dir, "simpa_examples.rst"), "w")
 simpa_examples_rst_file.write(
     "simpa\_examples\n=========================================\n\n.. toctree::\n   :maxdepth: 2\n\n")
-examples = glob.glob(os.path.join(current_dir, "../" + folder_level + "simpa_examples/*.py"))
+examples = sorted(glob.glob(os.path.join(current_dir, "../" + folder_level + "simpa_examples/*.py")))
 for example in examples:
     example_file_name = example.split("/")[-1]
     if example_file_name == "__init__.py":

docs/source/simpa.core.device_digital_twins.illumination_geometries.rst

@@ -48,6 +48,12 @@ illumination\_geometries
    :show-inheritance:
 
 
+.. automodule:: simpa.core.device_digital_twins.illumination_geometries.ring_illumination
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+
 .. automodule:: simpa.core.device_digital_twins.illumination_geometries.slit_illumination
    :members:
    :undoc-members:

docs/source/simpa_examples.rst

@@ -9,8 +9,8 @@ simpa\_examples
    linear_unmixing
    minimal_optical_simulation
    minimal_optical_simulation_uniform_cube
-   perform_iterative_qPAI_reconstruction
-   optical_and_acoustic_simulation
-   segmentation_loader
    msot_invision_simulation
+   optical_and_acoustic_simulation
    perform_image_reconstruction
+   perform_iterative_qPAI_reconstruction
+   segmentation_loader

simpa/core/device_digital_twins/__init__.py

@@ -131,6 +131,7 @@ def deserialize(dictionary_to_deserialize):
             field_of_view_extent_mm=dictionary_to_deserialize["field_of_view_extent_mm"])
         return deserialized_device
 
+
 """
 It is important to have these relative imports after the definition of the DigitalDeviceTwinBase class to avoid circular imports triggered by imported child classes
 """
@@ -146,6 +147,7 @@ def deserialize(dictionary_to_deserialize):
 from .illumination_geometries.pencil_beam_illumination import PencilBeamIlluminationGeometry  # nopep8
 from .illumination_geometries.disk_illumination import DiskIlluminationGeometry  # nopep8
 from .illumination_geometries.rectangle_illumination import RectangleIlluminationGeometry  # nopep8
+from .illumination_geometries.ring_illumination import RingIlluminationGeometry  # nopep8
 from .illumination_geometries.ithera_msot_acuity_illumination import MSOTAcuityIlluminationGeometry  # nopep8
 from .illumination_geometries.ithera_msot_invision_illumination import MSOTInVisionIlluminationGeometry  # nopep8
 from .pa_devices.ithera_msot_invision import InVision256TF  # nopep8

simpa/core/device_digital_twins/illumination_geometries/ring_illumination.py

@@ -0,0 +1,124 @@
+# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ
+# SPDX-FileCopyrightText: 2021 Janek Groehl
+# SPDX-License-Identifier: MIT
+
+import typing
+
+from simpa.core.device_digital_twins import IlluminationGeometryBase
+from simpa.utils import Settings, Tags
+import numpy as np
+from simpa.utils.serializer import SerializableSIMPAClass
+
+
+class RingIlluminationGeometry(IlluminationGeometryBase):
+    """
+    Defines a ring illumination geometry.
+    The device position is defined as the center of the ring.
+
+    Note: To create a ring light which illuminates a square tissue with the same center and any inner radius r,
+    create the following geometry using the tissue width:
+
+    >>> ring_light = RingIlluminationGeometry(inner_radius_in_mm=r,
+                                              outer_radius_in_mm=tissue_width / 2.,
+                                              device_position_mm=np.array([tissue_width / 2., tissue_width / 2., 0]))
+    """
+
+    def __init__(self,
+                 outer_radius_in_mm: float = 1,
+                 inner_radius_in_mm: float = 0,
+                 lower_angular_bound: float = 0,
+                 upper_angular_bound: float = 0,
+                 device_position_mm: typing.Optional[np.ndarray] = None,
+                 source_direction_vector: typing.Optional[np.ndarray] = None,
+                 field_of_view_extent_mm: typing.Optional[np.ndarray] = None):
+        """
+        :param outer_radius_in_mm: The outer radius of the ring in mm.
+        :param inner_radius_in_mm: The inner radius of the ring in mm. If 0, should match the disk illumination.
+        :param lower_angular_bound: The lower angular bound in radians. If both bounds are 0, than no bound is applied.
+        Note that the bound of 0 starts from the x-axis on the right side and is applied clockwise.
+        :param upper_angular_bound: The upper angular bound in radians. If both bounds are 0, than no bound is applied.
+        Note that the bound is applied clockwise in relation to the lower bound.
+        :param device_position_mm: The device position in mm, the center of the ring.
+        If None, the position is defined as [0, 0, 0].
+        :param source_direction_vector: Direction of the illumination source.
+        If None, the direction is defined as [0, 0, 1].
+        :param field_of_view_extent_mm: Field of view which is defined as a numpy array of the shape \
+        [xs, xe, ys, ye, zs, ze], where x, y, and z denote the coordinate axes and s and e denote the start and end \
+        positions.
+        """
+        if device_position_mm is None:
+            device_position_mm = np.zeros(3)
+
+        if source_direction_vector is None:
+            source_direction_vector = np.array([0, 0, 1])
+
+        super(RingIlluminationGeometry, self).__init__(device_position_mm=device_position_mm,
+                                                       source_direction_vector=source_direction_vector,
+                                                       field_of_view_extent_mm=field_of_view_extent_mm)
+
+        assert inner_radius_in_mm >= 0, f"The inner radius has to be 0 or positive, not {inner_radius_in_mm}!"
+        assert outer_radius_in_mm >= inner_radius_in_mm, \
+            f"The outer radius ({outer_radius_in_mm}) has to be at least as large " \
+            f"as the inner radius ({inner_radius_in_mm})!"
+
+        assert lower_angular_bound >= 0, f"The lower angular bound has to be 0 or positive, not {lower_angular_bound}!"
+        assert upper_angular_bound >= lower_angular_bound, \
+            f"The outer radius ({upper_angular_bound}) has to be at least as large " \
+            f"as the inner radius ({lower_angular_bound})!"
+
+        self.outer_radius_in_mm = outer_radius_in_mm
+        self.inner_radius_in_mm = inner_radius_in_mm
+        self.lower_angular_bound = lower_angular_bound
+        self.upper_angular_bound = upper_angular_bound
+
+    def get_mcx_illuminator_definition(self, global_settings: Settings) -> dict:
+        """
+        Returns the illumination parameters for MCX simulations.
+        :param global_settings: The global settings.
+        :return: The illumination parameters as a dictionary.
+        """
+        assert isinstance(global_settings, Settings), type(global_settings)
+
+        source_type = Tags.ILLUMINATION_TYPE_RING
+
+        spacing = global_settings[Tags.SPACING_MM]
+
+        device_position = list(self.device_position_mm / spacing + 1)  # No need to round
+
+        source_direction = list(self.normalized_source_direction_vector)
+
+        source_param1 = [self.outer_radius_in_mm / spacing,
+                         self.inner_radius_in_mm / spacing,
+                         self.lower_angular_bound,
+                         self.upper_angular_bound]
+
+        return {
+            "Type": source_type,
+            "Pos": device_position,
+            "Dir": source_direction,
+            "Param1": source_param1
+        }
+
+    def serialize(self) -> dict:
+        """
+        Serializes the object into a dictionary.
+        :return: The dictionary representing the serialized object.
+        """
+        serialized_device = self.__dict__
+        device_dict = {RingIlluminationGeometry.__name__: serialized_device}
+        return device_dict
+
+    @staticmethod
+    def deserialize(dictionary_to_deserialize: dict) -> SerializableSIMPAClass:
+        """
+        Deserializes the provided dict into an object of this type.
+        :param dictionary_to_deserialize: The dictionary to deserialize.
+        :return: The deserialized object from the dictionary.
+        """
+        assert isinstance(dictionary_to_deserialize, dict), type(dictionary_to_deserialize)
+
+        deserialized_device = RingIlluminationGeometry()
+        for key, value in dictionary_to_deserialize.items():
+            deserialized_device.__dict__[key] = value
+
+        return deserialized_device

simpa/utils/calculate.py

@@ -9,12 +9,18 @@
 from scipy.interpolate import interp1d
 
 
-def calculate_oxygenation(molecule_list):
+def calculate_oxygenation(molecule_list: list) -> Union[float, int, torch.Tensor]:
     """
-    :return: an oxygenation value between 0 and 1 if possible, or None, if not computable.
+    Calculate the oxygenation level based on the volume fractions of deoxyhaemoglobin and oxyhaemoglobin.
+
+    This function takes a list of molecules and returns an oxygenation value between 0 and 1 if computable,
+    otherwise returns None.
+
+    :param molecule_list: List of molecules with their spectrum information and volume fractions.
+    :return: An oxygenation value between 0 and 1 if possible, or None if not computable.
     """
-    hb = None
-    hbO2 = None
+    hb = None    # Volume fraction of deoxyhaemoglobin
+    hbO2 = None  # Volume fraction of oxyhaemoglobin
 
     for molecule in molecule_list:
         if molecule.spectrum.spectrum_name == "Deoxyhemoglobin":
@@ -36,7 +42,8 @@ def calculate_oxygenation(molecule_list):
     return hbO2 / (hb + hbO2)
 
 
-def create_spline_for_range(xmin_mm=0, xmax_mm=10, maximum_y_elevation_mm=1, spacing=0.1):
+def create_spline_for_range(xmin_mm: Union[float, int] = 0, xmax_mm: Union[float, int] = 10,
+                            maximum_y_elevation_mm: Union[float, int] = 1, spacing: Union[float, int] = 0.1) -> tuple:
     """
     Creates a functional that simulates distortion along the y position
     between the minimum and maximum x positions. The elevation can never be
@@ -45,6 +52,7 @@ def create_spline_for_range(xmin_mm=0, xmax_mm=10, maximum_y_elevation_mm=1, spa
     :param xmin_mm: the minimum x axis value the return functional is defined in
     :param xmax_mm: the maximum x axis value the return functional is defined in
     :param maximum_y_elevation_mm: the maximum y axis value the return functional will yield
+    :param spacing: the voxel spacing in the simulation
     :return: a functional that describes a distortion field along the y axis
 
     """
@@ -85,7 +93,22 @@ def create_spline_for_range(xmin_mm=0, xmax_mm=10, maximum_y_elevation_mm=1, spa
     return spline, max_el
 
 
-def spline_evaluator2d_voxel(x, y, spline, offset_voxel, thickness_voxel):
+def spline_evaluator2d_voxel(x: int, y: int, spline: Union[list, np.ndarray], offset_voxel: Union[float, int],
+                             thickness_voxel: int) -> bool:
+    """
+    Evaluate whether a given point (x, y) lies within the thickness bounds around a spline curve.
+
+    This function checks if the y-coordinate of a point lies within a vertical range defined
+    around a spline curve at a specific x-coordinate. The range is determined by the spline elevation,
+    an offset, and a thickness.
+
+    :param x: The x-coordinate of the point to evaluate.
+    :param y: The y-coordinate of the point to evaluate.
+    :param spline: A 1D array or list representing the spline curve elevations at each x-coordinate.
+    :param offset_voxel: The offset to be added to the spline elevation to define the starting y-coordinate of the range.
+    :param thickness_voxel: The vertical thickness of the range around the spline.
+    :return: True if the point (x, y) lies within the range around the spline, False otherwise.
+    """
     elevation = spline[x]
     y_value = np.round(elevation + offset_voxel)
     if y_value <= y < thickness_voxel + y_value:
@@ -94,7 +117,7 @@ def spline_evaluator2d_voxel(x, y, spline, offset_voxel, thickness_voxel):
         return False
 
 
-def calculate_gruneisen_parameter_from_temperature(temperature_in_celcius):
+def calculate_gruneisen_parameter_from_temperature(temperature_in_celcius: Union[float, int]) -> Union[float, int]:
     """
     This function returns the dimensionless gruneisen parameter based on a heuristic formula that
     was determined experimentally::
@@ -114,7 +137,7 @@ def calculate_gruneisen_parameter_from_temperature(temperature_in_celcius):
     return 0.0043 + 0.0053 * temperature_in_celcius
 
 
-def randomize_uniform(min_value: float, max_value: float):
+def randomize_uniform(min_value: float, max_value: float) -> Union[float, int]:
     """
     returns a uniformly drawn random number in [min_value, max_value[
 
@@ -126,7 +149,7 @@ def randomize_uniform(min_value: float, max_value: float):
     return (np.random.random() * (max_value-min_value)) + min_value
 
 
-def rotation_x(theta):
+def rotation_x(theta: Union[float, int]) -> torch.Tensor:
     """
     Rotation matrix around the x-axis with angle theta.
 
@@ -138,7 +161,7 @@ def rotation_x(theta):
                          [0, torch.sin(theta), torch.cos(theta)]])
 
 
-def rotation_y(theta):
+def rotation_y(theta: Union[float, int]) -> torch.Tensor:
     """
     Rotation matrix around the y-axis with angle theta.
 
@@ -150,7 +173,7 @@ def rotation_y(theta):
                          [-torch.sin(theta), 0, torch.cos(theta)]])
 
 
-def rotation_z(theta):
+def rotation_z(theta: Union[float, int]) -> torch.Tensor:
     """
     Rotation matrix around the z-axis with angle theta.
 
@@ -162,7 +185,7 @@ def rotation_z(theta):
                          [0, 0, 1]])
 
 
-def rotation(angles):
+def rotation(angles: Union[list, np.ndarray]) -> torch.Tensor:
     """
     Rotation matrix around the x-, y-, and z-axis with angles [theta_x, theta_y, theta_z].
 
@@ -172,7 +195,7 @@ def rotation(angles):
     return rotation_x(angles[0]) * rotation_y(angles[1]) * rotation_z(angles[2])
 
 
-def rotation_matrix_between_vectors(a, b):
+def rotation_matrix_between_vectors(a: np.ndarray, b: np.ndarray) -> np.ndarray:
     """
     Returns the rotation matrix from a to b