Merge pull request #340 from IMSY-DKFZ/T339_specific_hetero

T339 specific hetero

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
@@ -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/utils/libraries/molecule_library.py

@@ -38,7 +38,7 @@ def __init__(self, segmentation_type: Optional[str] = None, molecular_compositio
         """
         super().__init__()
         self.segmentation_type = segmentation_type
-        self.internal_properties = None
+        self.internal_properties: TissueProperties = None
 
         if molecular_composition_settings is None:
             return
@@ -70,7 +70,7 @@ def update_internal_properties(self, settings):
             self.internal_properties[Tags.DATA_FIELD_ALPHA_COEFF] += molecule.volume_fraction * \
                 molecule.alpha_coefficient
 
-        if (np.abs(self.internal_properties.volume_fraction - 1.0) > 1e-5).any():
+        if not (torch.abs(self.internal_properties.volume_fraction - 1.0) < 1e-5).any():
             raise AssertionError("Invalid Molecular composition! The volume fractions of all molecules must be"
                                  "exactly 100%!")
 

simpa/utils/libraries/structure_library/HorizontalLayerStructure.py

@@ -7,6 +7,9 @@
 from simpa.utils import Tags
 from simpa.utils.libraries.molecule_library import MolecularComposition
 from simpa.utils.libraries.structure_library.StructureBase import GeometricalStructure
+from simpa.log import Logger
+
+logger = Logger()
 
 
 class HorizontalLayerStructure(GeometricalStructure):
@@ -100,6 +103,36 @@ def get_enclosed_indices(self):
 
         return bools_all_layers.cpu().numpy(), volume_fractions[bools_all_layers].cpu().numpy()
 
+    def update_molecule_volume_fractions(self, single_structure_settings):
+        for molecule in self.molecule_composition:
+            old_vol = getattr(molecule, "volume_fraction")
+            if isinstance(old_vol, torch.Tensor):
+                structure_start_voxels = (torch.tensor(single_structure_settings[Tags.STRUCTURE_START_MM]) /
+                                          self.voxel_spacing)
+                structure_end_voxels = (torch.tensor(single_structure_settings[Tags.STRUCTURE_END_MM]) /
+                                        self.voxel_spacing)
+                structure_size = structure_end_voxels - structure_start_voxels
+
+                if self.volume_dimensions_voxels[2] != old_vol.shape[2]:
+                    if self.volume_dimensions_voxels[2] == structure_size[2]:
+                        pad_start = structure_start_voxels.flip(dims=[0])
+                        pad_end = ((torch.from_numpy(self.volume_dimensions_voxels) - structure_end_voxels)
+                                   .flip(dims=[0]))
+                        for count, structure_end in enumerate(structure_end_voxels):
+                            if structure_end == 0:
+                                pad_end[2 - count] = 0
+
+                        if (pad_start > 1e-6).any() or (pad_end > 1e-6).any():
+                            padding_list = torch.flatten(torch.stack((pad_start, pad_end), 1)).tolist()
+                            padding = tuple(map(int, padding_list))
+                            padded_vol = torch.nn.functional.pad(old_vol, padding, mode='constant', value=0)
+                            setattr(molecule, "volume_fraction", padded_vol)
+                    else:
+                        logger.critical("Tensor does not have the same dimensionality as the area it should fill" +
+                                        "{} or the dimensions of the entire ".format(old_vol.shape) +
+                                        "simulation volume{}! Please change this.".format(
+                                            self.volume_dimensions_voxels.shape))
+
 
 def define_horizontal_layer_structure_settings(molecular_composition: MolecularComposition,
                                                z_start_mm: float = 0, thickness_mm: float = 0, priority: int = 10,

simpa/utils/libraries/structure_library/StructureBase.py

@@ -66,13 +66,18 @@ def __init__(self, global_settings: Settings,
         else:
             self.partial_volume = False
 
-        self.molecule_composition = single_structure_settings[Tags.MOLECULE_COMPOSITION]
+        self.molecule_composition: MolecularComposition = single_structure_settings[Tags.MOLECULE_COMPOSITION]
+        self.update_molecule_volume_fractions(single_structure_settings)
         self.molecule_composition.update_internal_properties(global_settings)
 
         self.geometrical_volume = np.zeros(self.volume_dimensions_voxels, dtype=np.float32)
         self.params = self.get_params_from_settings(single_structure_settings)
         self.fill_internal_volume()
 
+        assert ((self.molecule_composition.internal_properties.volume_fraction[self.geometrical_volume != 0] - 1 < 1e-5)
+                .all()), ("Invalid Molecular composition! The volume fractions of all molecules in the structure must"
+                          "be exactly 100%!")
+
     def fill_internal_volume(self):
         """
         Fills self.geometrical_volume of the GeometricalStructure.
@@ -95,7 +100,7 @@ def get_enclosed_indices(self):
         pass
 
     @abstractmethod
-    def get_params_from_settings(self, single_structure_settings):
+    def get_params_from_settings(self, single_structure_settings: Settings):
         """
         Gets all the parameters required for the specific GeometricalStructure.
         :param single_structure_settings: Settings which describe the specific GeometricalStructure.
@@ -112,6 +117,17 @@ def properties_for_wavelength(self, settings, wavelength) -> TissueProperties:
         """
         return self.molecule_composition.get_properties_for_wavelength(settings, wavelength)
 
+    def update_molecule_volume_fractions(self, single_structure_settings: Settings):
+        """
+        In particular cases, only molecule volume fractions are determined by tensors that expand the structure. This
+        method allows the tensor to only have the size of the structure, with the rest of the volume filled with volume
+        fractions of 0.
+        In the case where the tensors defined are such that they fill the volume, they will be left. Later, when
+        using priority_sorted_structures the volume used will be that within the boundaries in the shape.
+        :param single_structure_settings: Settings which describe the specific GeometricalStructure.
+        """
+        pass
+
     @abstractmethod
     def to_settings(self) -> Settings:
         """

simpa/utils/libraries/structure_library/__init__.py

@@ -36,6 +36,7 @@ def priority_sorted_structures(settings: Settings, volume_creator_settings: dict
     sorted_structure_settings = sorted(
         [structure_setting for structure_setting in volume_creator_settings[Tags.STRUCTURES].values()],
         key=lambda s: s[Tags.PRIORITY] if Tags.PRIORITY in s else 0, reverse=True)
+
     for structure_setting in sorted_structure_settings:
         try:
             structure_class = globals()[structure_setting[Tags.STRUCTURE_TYPE]]

simpa_examples/minimal_optical_simulation_heterogeneous_tissue.py

@@ -50,8 +50,8 @@ def create_example_tissue(settings):
     muscle_dictionary[Tags.STRUCTURE_START_MM] = [0, 0, 10]
     muscle_dictionary[Tags.STRUCTURE_END_MM] = [0, 0, 100]
     muscle_dictionary[Tags.MOLECULE_COMPOSITION] = tissue_library.muscle(
-        background_oxy=sp.BlobHeterogeneity(dim_x, dim_y, dim_z, SPACING,
-                                            target_min=0.5, target_max=0.8).get_map(),
+        oxygenation=sp.BlobHeterogeneity(dim_x, dim_y, dim_z, SPACING,
+                                         target_min=0.5, target_max=0.8).get_map(),
         blood_volume_fraction=sp.BlobHeterogeneity(dim_x, dim_y, dim_z, SPACING,
                                                    target_min=0.01, target_max=0.1).get_map())
     muscle_dictionary[Tags.CONSIDER_PARTIAL_VOLUME] = True
@@ -78,8 +78,8 @@ def create_example_tissue(settings):
     epidermis_dictionary[Tags.STRUCTURE_START_MM] = [0, 0, 9]
     epidermis_dictionary[Tags.STRUCTURE_END_MM] = [0, 0, 10]
     epidermis_dictionary[Tags.MOLECULE_COMPOSITION] = tissue_library.epidermis(
-        melanosom_volume_fraction=sp.RandomHeterogeneity(dim_x, dim_y, dim_z, SPACING,
-                                                         target_min=0.1, target_max=0.2).get_map())
+        melanin_volume_fraction=sp.RandomHeterogeneity(dim_x, dim_y, dim_z, SPACING,
+                                                       target_min=0.1, target_max=0.2).get_map())
     epidermis_dictionary[Tags.CONSIDER_PARTIAL_VOLUME] = True
     epidermis_dictionary[Tags.ADHERE_TO_DEFORMATION] = True
     epidermis_dictionary[Tags.STRUCTURE_TYPE] = Tags.HORIZONTAL_LAYER_STRUCTURE