Merge branch 'develop' into T364_update_segmentation_loader

.gitignore

@@ -15,6 +15,8 @@ docs/objects.inv
 
 simpa_tests/figures
 simpa_tests/figures/*
+simpa_tests/**/figures
+simpa_tests/**/figures/*
 simpa_tests/manual_tests/figures
 simpa_tests/manual_tests/reference_figures
 simpa_tests/manual_tests/manual_tests_*
@@ -154,4 +156,4 @@ dmypy.json
 
 # numpy files
 *npy
-/.mailmap
+/.mailmap

docs/source/benchmarking.md

@@ -7,15 +7,20 @@ It allows customization of initial spacing, final spacing, step size, output fil
 To check and amend which simulations are run, see the [performance check script](../../simpa_examples/benchmarking/performance_check.py).
 
 ## Usage
-To use this script, run it from the command line with the desired options. Please ensure you check two things before
-running the script: First, ensure that the device will not be in use for the duration - ideally restart before 
-benchmarking - of the benchmarking process, 
+In order to be able to use this script, please first ensure that you have the dependencies required for the benchmarking
+scripts. To do this, please navigate to the simpa directory and execute `pip install .[profile]`.
+
+Now, you can run [performance_check.py](../../simpa_examples/benchmarking/performance_check.py) and [run_benchmarking.sh](../../simpa_examples/benchmarking/run_benchmarking.sh) from the command line with the desired
+options. Please ensure you check two things before running the script: First, ensure that the device will not be in use
+for the duration - ideally restart before benchmarking - of the benchmarking process, 
 as this will create large uncertainties within the outcomes. Second, ensure that you don't accidentally write over any 
 existing file by saving the files created by this script after runtime to different location.
 
-The script will create multiple text files (eg. benchmarking_data_TIME_0.2.txt), showing the line by line profiling of
-the most recent runs, as well as two csv's with the data from all the runs (benchmarking_data_frame.csv) and the means
-and standard deviations of all the runs (benchmarking_data_frame_mean.csv).
+The both scripts create text files (eg. benchmarking_data_TIME_0.2.txt), showing the line by line profiling of
+the most recent runs. The [run_benchmarking.sh](../../simpa_examples/benchmarking/run_benchmarking.sh) also creates two csv's: one with the data from all the runs
+(benchmarking_data_frame.csv) and one with the means and standard deviations of all the runs
+(benchmarking_data_frame_mean.csv). With both scripts, unless the user intentionally changes the save folder name,
+the output files will be overwritten.
 
 Below is a description of the available options and how to use them.
 
@@ -28,8 +33,8 @@ Please put this in the conversation of the pull request, and not add it to the f
 
 ## Options
 - **`-i, --init`**: First spacing to benchmark (default = 0.2mm).
-- **`-c, --cease`**: Final spacing to benchmark (default = 0.4mm).
-- **`-s, --step`**: Step between spacings (default = 0.1mm).
+- **`-c, --cease`**: Final spacing to benchmark (default = 0.25mm).
+- **`-s, --step`**: Step between spacings (default = 0.05mm).
 - **`-f, --file`**: Where to store the output files (default = save in current directory; 'print' prints it in console).
 - **`-t, --time`**: Profile times taken (if no profile is specified, all are set).
 - **`-g, --gpu`**: Profile GPU usage (if no profile is specified, all are set).
@@ -40,9 +45,9 @@ Please put this in the conversation of the pull request, and not add it to the f
 ## Default Values
 If no options are provided for initial spacing, final spacing, or step size, the script uses the following default
 values:
-- **Initial Spacing**: 0.2mm
-- **Final Spacing**: 0.4mm
-- **Step Size**: 0.1mm
+- **Initial Spacing**: 0.15mm
+- **Final Spacing**: 0.25mm
+- **Step Size**: 0.05mm
 
 If no profiling options are specified, all three profilers (time, GPU memory, and memory) are used by default.
 

docs/source/minimal_optical_simulation_heterogeneous_tissue.rst

@@ -0,0 +1,7 @@
+minimal_optical_simulation_heterogeneous_tissue
+=========================================
+
+.. literalinclude:: ../../simpa_examples/minimal_optical_simulation_heterogeneous_tissue.py
+   :language: python
+   :lines: 1-
+

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.heterogeneity_generator
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+
 .. automodule:: simpa.utils.libraries.literature_values
    :members:
    :undoc-members:

docs/source/simpa_examples.rst

@@ -8,6 +8,7 @@ simpa\_examples
    create_custom_tissues
    linear_unmixing
    minimal_optical_simulation
+   minimal_optical_simulation_heterogeneous_tissue
    minimal_optical_simulation_uniform_cube
    msot_invision_simulation
    optical_and_acoustic_simulation

pyproject.toml

@@ -7,7 +7,7 @@ name = "simpa"
 dynamic = ["version"]
 authors = [
     {name = "Division of Intelligent Medical Systems (IMSY), DKFZ", email = "k.dreher@dkfz-heidelberg.de"},
-    {name = "Janek Groehl"}
+    {name = "Janek Groehl <janekgroehl@live.de>"}
 ]
 description = "Simulation and Image Processing for Photonics and Acoustics"
 license = {text = "MIT"}
@@ -32,7 +32,8 @@ dependencies = [
     "wget>=3.2",               # Is Public Domain (MIT compatible)
     "jdata>=0.5.2",            # Uses Apache 2.0-License (MIT compatible)
     "pre-commit>=3.2.2",       # Uses MIT-License (MIT compatible)
-    "PyWavelets"               # Uses MIT-License (MIT compatible)
+    "PyWavelets",              # Uses MIT-License (MIT compatible)
+    "scikit-learn>=1.1.0",     # Uses BSD-License (MIT compatible)
 ]
 
 [project.optional-dependencies]

simpa/__init__.py

@@ -10,7 +10,7 @@
     __version__ = version("simpa")
 except PackageNotFoundError:
     __version__ = "unknown version"
-    
+
 from .core.simulation_modules.volume_creation_module.volume_creation_module_model_based_adapter import \
     ModelBasedVolumeCreationAdapter
 from .core.simulation_modules.volume_creation_module.volume_creation_module_segmentation_based_adapter import \

simpa/core/device_digital_twins/pa_devices/ithera_msot_acuity.py

@@ -1,6 +1,8 @@
 # SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ
 # SPDX-FileCopyrightText: 2021 Janek Groehl
 # SPDX-License-Identifier: MIT
+import torch
+import torch.nn.functional as F
 
 from simpa.core.device_digital_twins import PhotoacousticDevice, \
     CurvedArrayDetectionGeometry, MSOTAcuityIlluminationGeometry
@@ -88,7 +90,7 @@ def __init__(self, device_position_mm: np.ndarray = None,
                                                                                             -y_pos_relative_to_membrane,
                                                                                             -43.2]))
 
-    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):
         """
         Updates the volume creation settings of the model based volume creator according to the size of the device.
         :param global_settings: Settings for the entire simulation pipeline.
@@ -105,6 +107,8 @@ def update_settings_for_use_of_model_based_volume_creator(self, global_settings)
         probe_size_mm = self.probe_height_mm
         mediprene_layer_height_mm = self.mediprene_membrane_height_mm
         heavy_water_layer_height_mm = probe_size_mm - mediprene_layer_height_mm
+        spacing_mm = global_settings[Tags.SPACING_MM]
+        old_volume_height_pixels = round(global_settings[Tags.DIM_VOLUME_Z_MM] / spacing_mm)
 
         if Tags.US_GEL in volume_creator_settings and volume_creator_settings[Tags.US_GEL]:
             us_gel_thickness = np.random.normal(0.4, 0.1)
@@ -122,13 +126,10 @@ def update_settings_for_use_of_model_based_volume_creator(self, global_settings)
         # 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) + \
-                global_settings[Tags.SPACING_MM]:
-            width_shift_for_structures_mm = (round(self.detection_geometry.probe_width_mm) +
-                                             global_settings[Tags.SPACING_MM] -
+        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
-            global_settings[Tags.DIM_VOLUME_X_MM] = round(self.detection_geometry.probe_width_mm) + \
-                global_settings[Tags.SPACING_MM]
+            global_settings[Tags.DIM_VOLUME_X_MM] = 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
@@ -139,6 +140,17 @@ def update_settings_for_use_of_model_based_volume_creator(self, global_settings)
             self.logger.debug("Adjusting " + str(structure_key))
             structure_dict = volume_creator_settings[Tags.STRUCTURES][structure_key]
             if Tags.STRUCTURE_START_MM in structure_dict:
+                for molecule in structure_dict[Tags.MOLECULE_COMPOSITION]:
+                    old_volume_fraction = getattr(molecule, "volume_fraction")
+                    if isinstance(old_volume_fraction, torch.Tensor):
+                        if old_volume_fraction.shape[2] == old_volume_height_pixels:
+                            width_shift_pixels = round(width_shift_for_structures_mm / spacing_mm)
+                            z_shift_pixels = round(z_dim_position_shift_mm / spacing_mm)
+                            padding_height = (z_shift_pixels, 0, 0, 0, 0, 0)
+                            padding_width = ((width_shift_pixels, width_shift_pixels), (0, 0), (0, 0))
+                            padded_up = F.pad(old_volume_fraction, padding_height, mode='constant', value=0)
+                            padded_vol = np.pad(padded_up.numpy(), padding_width, mode='edge')
+                            setattr(molecule, "volume_fraction", torch.from_numpy(padded_vol))
                 structure_dict[Tags.STRUCTURE_START_MM][0] = structure_dict[Tags.STRUCTURE_START_MM][
                     0] + width_shift_for_structures_mm
                 structure_dict[Tags.STRUCTURE_START_MM][2] = structure_dict[Tags.STRUCTURE_START_MM][

simpa/core/processing_components/monospectral/iterative_qPAI_algorithm.py

@@ -369,7 +369,7 @@ def standard_optical_properties(self, image_data: np.ndarray) -> dict:
             scattering = float(self.global_settings[Tags.DATA_FIELD_SCATTERING_PER_CM]) * np.ones(shape)
         else:
             background_dict = TISSUE_LIBRARY.muscle()
-            scattering = float(MolecularComposition.get_properties_for_wavelength(background_dict,
+            scattering = float(MolecularComposition.get_properties_for_wavelength(background_dict, self.global_settings,
                                                                                   wavelength=800)["mus"])
             scattering = scattering * np.ones(shape)
 

simpa/core/simulation.py

@@ -56,7 +56,7 @@ def simulate(simulation_pipeline: list, settings: Settings, digital_device_twin:
         simpa_output_path = path + settings[Tags.VOLUME_NAME]
 
     settings[Tags.SIMPA_OUTPUT_PATH] = simpa_output_path + ".hdf5"
-    
+
     simpa_output[Tags.SIMPA_VERSION] = __version__
     simpa_output[Tags.SETTINGS] = settings
     simpa_output[Tags.DIGITAL_DEVICE] = digital_device_twin

simpa/core/simulation_modules/__init__.py

@@ -41,4 +41,4 @@ def get_additional_flags(self) -> List[str]:
         if Tags.ADDITIONAL_FLAGS in self.component_settings:
             for flag in self.component_settings[Tags.ADDITIONAL_FLAGS]:
                 cmd.append(str(flag))
-        return cmd
+        return cmd

simpa/core/simulation_modules/acoustic_forward_module/acoustic_forward_module_k_wave_adapter.py

@@ -5,7 +5,7 @@
 import gc
 import os
 import subprocess
-from typing import List 
+from typing import List
 
 import numpy as np
 import scipy.io as sio
@@ -99,16 +99,16 @@ def forward_model(self, detection_geometry: DetectionGeometryBase) -> np.ndarray
 
         detectors_are_aligned_along_x_axis = field_of_view_extent[2] == 0 and field_of_view_extent[3] == 0
         detectors_are_aligned_along_y_axis = field_of_view_extent[0] == 0 and field_of_view_extent[1] == 0
-        if detectors_are_aligned_along_x_axis or detectors_are_aligned_along_y_axis:
-            axes = (0, 1)
+        if not (Tags.ACOUSTIC_SIMULATION_3D in self.component_settings
+                and self.component_settings[Tags.ACOUSTIC_SIMULATION_3D]) and \
+                (detectors_are_aligned_along_x_axis or detectors_are_aligned_along_y_axis):
             if detectors_are_aligned_along_y_axis:
                 transducer_plane = int(round((detector_positions_mm[0, 0] / self.global_settings[Tags.SPACING_MM]))) - 1
                 image_slice = np.s_[transducer_plane, :, :]
             else:
                 transducer_plane = int(round((detector_positions_mm[0, 1] / self.global_settings[Tags.SPACING_MM]))) - 1
                 image_slice = np.s_[:, transducer_plane, :]
         else:
-            axes = (0, 2)
             image_slice = np.s_[:]
 
         data_dict[Tags.DATA_FIELD_SPEED_OF_SOUND] = data_dict[Tags.DATA_FIELD_SPEED_OF_SOUND][image_slice].T
@@ -156,7 +156,8 @@ def k_wave_acoustic_forward_model(self, detection_geometry: DetectionGeometryBas
         field_of_view = pa_device.get_field_of_view_mm()
         detector_positions_mm = pa_device.get_detector_element_positions_accounting_for_device_position_mm()
 
-        if not self.component_settings.get(Tags.ACOUSTIC_SIMULATION_3D):
+        if not (Tags.ACOUSTIC_SIMULATION_3D in self.component_settings
+                and self.component_settings[Tags.ACOUSTIC_SIMULATION_3D]):
             detectors_are_aligned_along_x_axis = np.abs(field_of_view[2] - field_of_view[3]) < 1e-5
             detectors_are_aligned_along_y_axis = np.abs(field_of_view[0] - field_of_view[1]) < 1e-5
             if detectors_are_aligned_along_x_axis or detectors_are_aligned_along_y_axis:

simpa/core/simulation_modules/reconstruction_module/reconstruction_module_time_reversal_adapter.py

@@ -55,7 +55,8 @@ def get_acoustic_properties(self, input_data: dict, detection_geometry):
             raise AttributeError("Please specify a value for SPACING_MM")
 
         detector_positions = detection_geometry.get_detector_element_positions_accounting_for_device_position_mm()
-        detector_positions_voxels = np.round(detector_positions / spacing_in_mm).astype(int)
+        # we add eps of 1e-10 because numpy rounds 0.5 to the next even number
+        detector_positions_voxels = np.round(detector_positions / spacing_in_mm + 1e-10).astype(int)
 
         # plus 2 because of off-
         volume_x_dim = int(np.ceil(self.global_settings[Tags.DIM_VOLUME_X_MM] / spacing_in_mm) + 1)
@@ -170,7 +171,7 @@ def reconstruction_algorithm(self, time_series_sensor_data, detection_geometry):
 
         matlab_binary_path = self.component_settings[Tags.ACOUSTIC_MODEL_BINARY_PATH]
         cmd = generate_matlab_cmd(matlab_binary_path, time_reversal_script, acoustic_path, self.get_additional_flags())
-        
+
         cur_dir = os.getcwd()
         os.chdir(self.global_settings[Tags.SIMULATION_PATH])
         self.logger.info(cmd)

simpa/core/simulation_modules/volume_creation_module/__init__.py

@@ -30,10 +30,7 @@ def load_component_settings(self) -> Settings:
 
     def create_empty_volumes(self):
         volumes = dict()
-        voxel_spacing = self.global_settings[Tags.SPACING_MM]
-        volume_x_dim = int(round(self.global_settings[Tags.DIM_VOLUME_X_MM] / voxel_spacing))
-        volume_y_dim = int(round(self.global_settings[Tags.DIM_VOLUME_Y_MM] / voxel_spacing))
-        volume_z_dim = int(round(self.global_settings[Tags.DIM_VOLUME_Z_MM] / voxel_spacing))
+        volume_x_dim, volume_y_dim, volume_z_dim = self.global_settings.get_volume_dimensions_voxels()
         sizes = (volume_x_dim, volume_y_dim, volume_z_dim)
 
         wavelength = self.global_settings[Tags.WAVELENGTH]

simpa/core/simulation_modules/volume_creation_module/volume_creation_module_model_based_adapter.py

@@ -68,7 +68,7 @@ def create_simulation_volume(self) -> dict:
         for structure in priority_sorted_structures(self.global_settings, self.component_settings):
             self.logger.debug(type(structure))
 
-            structure_properties = structure.properties_for_wavelength(wavelength)
+            structure_properties = structure.properties_for_wavelength(self.global_settings, wavelength)
 
             structure_volume_fractions = torch.as_tensor(
                 structure.geometrical_volume, dtype=torch.float, device=self.torch_device)
@@ -95,10 +95,21 @@ def create_simulation_volume(self) -> dict:
                     max_added_fractions[added_fraction_greater_than_any_added_fraction & mask] = \
                         added_volume_fraction[added_fraction_greater_than_any_added_fraction & mask]
                 else:
-                    volumes[key][mask] += added_volume_fraction[mask] * structure_properties[key]
+                    if isinstance(structure_properties[key], torch.Tensor):
+                        volumes[key][mask] += added_volume_fraction[mask] * \
+                            structure_properties[key].to(self.torch_device)[mask]
+                    elif isinstance(structure_properties[key], (float, np.float64, int, np.int64)):
+                        volumes[key][mask] += added_volume_fraction[mask] * structure_properties[key]
+                    else:
+                        raise ValueError(f"Unsupported type of structure property. "
+                                         f"Was {type(structure_properties[key])}.")
 
             global_volume_fractions[mask] += added_volume_fraction[mask]
 
+        if (torch.abs(global_volume_fractions[global_volume_fractions > 1]) < 1e-5).any():
+            raise AssertionError("Invalid Molecular composition! The volume fractions of all molecules must be"
+                                 "exactly 100%!")
+
         # convert volumes back to CPU
         for key in volumes.keys():
             volumes[key] = volumes[key].cpu().numpy().astype(np.float64, copy=False)