Merge pull request #290 from IMSY-DKFZ/T154_heterogeneous_tissue

T154 heterogeneous tissue

.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/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 = "[email protected]"},
-    {name = "Janek Groehl"}
+    {name = "Janek Groehl <[email protected]>"}
 ]
 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/reconstruction_module/reconstruction_module_time_reversal_adapter.py

@@ -170,7 +170,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)

simpa/core/simulation_modules/volume_creation_module/volume_creation_module_segmentation_based_adapter.py

@@ -12,7 +12,7 @@
 
 class SegmentationBasedVolumeCreationAdapter(VolumeCreatorModuleBase):
     """
-    This volume creator expects a np.ndarray to be in the settigs
+    This volume creator expects a np.ndarray to be in the settings
     under the Tags.INPUT_SEGMENTATION_VOLUME tag and uses this array
     together with a SegmentationClass mapping which is a dict defined in
     the settings under Tags.SEGMENTATION_CLASS_MAPPING.
@@ -23,6 +23,8 @@ class SegmentationBasedVolumeCreationAdapter(VolumeCreatorModuleBase):
     def create_simulation_volume(self) -> dict:
         volumes, x_dim_px, y_dim_px, z_dim_px = self.create_empty_volumes()
         wavelength = self.global_settings[Tags.WAVELENGTH]
+        for key in volumes.keys():
+            volumes[key] = volumes[key].to('cpu')
 
         segmentation_volume = self.component_settings[Tags.INPUT_SEGMENTATION_VOLUME]
         segmentation_classes = np.unique(segmentation_volume, return_counts=False)
@@ -41,17 +43,24 @@ def create_simulation_volume(self) -> dict:
         class_mapping = self.component_settings[Tags.SEGMENTATION_CLASS_MAPPING]
 
         for seg_class in segmentation_classes:
-            class_properties = class_mapping[seg_class].get_properties_for_wavelength(wavelength)
+            class_properties = class_mapping[seg_class].get_properties_for_wavelength(self.global_settings, wavelength)
             for prop_tag in property_tags:
-                assigned_prop = class_properties[prop_tag]
-                if assigned_prop is None:
-                    assigned_prop = torch.nan
-                volumes[prop_tag][segmentation_volume == seg_class] = assigned_prop
+                if isinstance(class_properties[prop_tag], (int, float)) or class_properties[prop_tag] == None:  # scalar
+                    assigned_prop = class_properties[prop_tag]
+                    if assigned_prop is None:
+                        assigned_prop = torch.nan
+                    volumes[prop_tag][segmentation_volume == seg_class] = assigned_prop
+                elif len(torch.Tensor.size(class_properties[prop_tag])) == 3:  # 3D map
+                    assigned_prop = class_properties[prop_tag][torch.tensor(segmentation_volume == seg_class)]
+                    assigned_prop[assigned_prop is None] = torch.nan
+                    volumes[prop_tag][torch.tensor(segmentation_volume == seg_class)] = assigned_prop
+                else:
+                    raise AssertionError("Properties need to either be a scalar or a 3D map.")
 
         save_hdf5(self.global_settings, self.global_settings[Tags.SIMPA_OUTPUT_PATH], "/settings/")
 
         # convert volumes back to CPU
         for key in volumes.keys():
-            volumes[key] = volumes[key].cpu().numpy().astype(np.float64, copy=False)
+            volumes[key] = volumes[key].numpy().astype(np.float64, copy=False)
 
         return volumes

simpa/utils/__init__.py

@@ -4,11 +4,19 @@
 
 # First load everything without internal dependencies
 from .tags import Tags
+from .settings import Settings
+
 from .libraries.literature_values import MorphologicalTissueProperties
 from .libraries.literature_values import StandardProperties
 from .libraries.literature_values import OpticalTissueProperties
 from .constants import SegmentationClasses
 
+# Heterogeneity
+
+from .libraries.heterogeneity_generator import RandomHeterogeneity
+from .libraries.heterogeneity_generator import BlobHeterogeneity
+from .libraries.heterogeneity_generator import ImageHeterogeneity
+
 # Then load classes and methods with an <b>increasing</b> amount of internal dependencies.
 # If there are import errors in the tests, it is probably due to an incorrect
 # initialization order
@@ -33,8 +41,6 @@
 from .deformation_manager import create_deformation_settings
 from .deformation_manager import get_functional_from_deformation_settings
 
-from .settings import Settings
-
 from .dict_path_manager import generate_dict_path
 from .dict_path_manager import get_data_field_from_simpa_output
 

simpa/utils/calculate.py

@@ -12,7 +12,7 @@
 def extract_hemoglobin_fractions(molecule_list: List) -> Dict[str, float]:
     """
     Extract hemoglobin volume fractions from a list of molecules.
-    
+
     :param molecule_list: List of molecules with their spectrum information and volume fractions.
     :return: A dictionary with hemoglobin types as keys and their volume fractions as values.
     """
@@ -34,7 +34,7 @@ def extract_hemoglobin_fractions(molecule_list: List) -> Dict[str, float]:
 def calculate_oxygenation(molecule_list: List) -> Optional[float]:
     """
     Calculate the oxygenation level based on the volume fractions of deoxyhemoglobin and oxyhemoglobin.
-    
+
     :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.
     """
@@ -44,16 +44,20 @@ def calculate_oxygenation(molecule_list: List) -> Optional[float]:
     total = hb + hbO2
 
     # Avoid division by zero. If none of the hemoglobin types are present, the oxygenation level is not computable.
-    if total < 1e-10:
-        return None
+    if isinstance(hb, torch.Tensor) or isinstance(hbO2, torch.Tensor):
+        return torch.where(total < 1e-10, 0, hbO2 / total)
 
-    return hbO2 / total
+    else:
+        if total < 1e-10:
+            return None
+        else:
+            return hbO2 / total
 
 
 def calculate_bvf(molecule_list: List) -> Union[float, int]:
     """
     Calculate the blood volume fraction based on the volume fractions of deoxyhemoglobin and oxyhemoglobin.
-    
+
     :param molecule_list: List of molecules with their spectrum information and volume fractions.
     :return: The blood volume fraction value between 0 and 1, or 0, if oxy and deoxy not present.
     """