diff --git a/.gitignore b/.gitignore index 7f2ea1d0..79570db7 100644 --- a/.gitignore +++ b/.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 \ No newline at end of file +/.mailmap diff --git a/docs/source/benchmarking.md b/docs/source/benchmarking.md index 12d87d63..dbfa1d8f 100644 --- a/docs/source/benchmarking.md +++ b/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. diff --git a/docs/source/minimal_optical_simulation_heterogeneous_tissue.rst b/docs/source/minimal_optical_simulation_heterogeneous_tissue.rst new file mode 100644 index 00000000..13cb3588 --- /dev/null +++ b/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- + diff --git a/docs/source/simpa.utils.libraries.rst b/docs/source/simpa.utils.libraries.rst index 5c38a7a0..fdff5eef 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.heterogeneity_generator + :members: + :undoc-members: + :show-inheritance: + + .. automodule:: simpa.utils.libraries.literature_values :members: :undoc-members: diff --git a/docs/source/simpa_examples.rst b/docs/source/simpa_examples.rst index e8a54665..9d229be6 100644 --- a/docs/source/simpa_examples.rst +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index fde75ceb..fc64165c 100644 --- a/pyproject.toml +++ b/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 "} ] 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] diff --git a/simpa/__init__.py b/simpa/__init__.py index b336a76f..c25f2ebb 100644 --- a/simpa/__init__.py +++ b/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 \ 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..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 @@ -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][ diff --git a/simpa/core/processing_components/monospectral/iterative_qPAI_algorithm.py b/simpa/core/processing_components/monospectral/iterative_qPAI_algorithm.py index 3a942564..59a6ad73 100644 --- a/simpa/core/processing_components/monospectral/iterative_qPAI_algorithm.py +++ b/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) diff --git a/simpa/core/simulation.py b/simpa/core/simulation.py index 52c90245..157ff678 100644 --- a/simpa/core/simulation.py +++ b/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 diff --git a/simpa/core/simulation_modules/__init__.py b/simpa/core/simulation_modules/__init__.py index 9624110c..a2af462b 100644 --- a/simpa/core/simulation_modules/__init__.py +++ b/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 \ No newline at end of file + return cmd diff --git a/simpa/core/simulation_modules/acoustic_forward_module/acoustic_forward_module_k_wave_adapter.py b/simpa/core/simulation_modules/acoustic_forward_module/acoustic_forward_module_k_wave_adapter.py index 17114651..249c5cbd 100644 --- a/simpa/core/simulation_modules/acoustic_forward_module/acoustic_forward_module_k_wave_adapter.py +++ b/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,8 +99,9 @@ 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, :, :] @@ -108,7 +109,6 @@ def forward_model(self, detection_geometry: DetectionGeometryBase) -> np.ndarray 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: diff --git a/simpa/core/simulation_modules/reconstruction_module/reconstruction_module_time_reversal_adapter.py b/simpa/core/simulation_modules/reconstruction_module/reconstruction_module_time_reversal_adapter.py index a67116e6..e266c7a1 100644 --- a/simpa/core/simulation_modules/reconstruction_module/reconstruction_module_time_reversal_adapter.py +++ b/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) diff --git a/simpa/core/simulation_modules/volume_creation_module/__init__.py b/simpa/core/simulation_modules/volume_creation_module/__init__.py index 6b9f0753..3aeecf14 100644 --- a/simpa/core/simulation_modules/volume_creation_module/__init__.py +++ b/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] diff --git a/simpa/core/simulation_modules/volume_creation_module/volume_creation_module_model_based_adapter.py b/simpa/core/simulation_modules/volume_creation_module/volume_creation_module_model_based_adapter.py index 2379adf3..ef10ef97 100644 --- a/simpa/core/simulation_modules/volume_creation_module/volume_creation_module_model_based_adapter.py +++ b/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) diff --git a/simpa/core/simulation_modules/volume_creation_module/volume_creation_module_segmentation_based_adapter.py b/simpa/core/simulation_modules/volume_creation_module/volume_creation_module_segmentation_based_adapter.py index 490ce2eb..9f6c0cb1 100644 --- a/simpa/core/simulation_modules/volume_creation_module/volume_creation_module_segmentation_based_adapter.py +++ b/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 diff --git a/simpa/utils/__init__.py b/simpa/utils/__init__.py index 73bdd1e9..2393f806 100644 --- a/simpa/utils/__init__.py +++ b/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 increasing 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 diff --git a/simpa/utils/calculate.py b/simpa/utils/calculate.py index bbed1609..31b59299 100644 --- a/simpa/utils/calculate.py +++ b/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. """ diff --git a/simpa/utils/libraries/heterogeneity_generator.py b/simpa/utils/libraries/heterogeneity_generator.py new file mode 100644 index 00000000..8d01b3f1 --- /dev/null +++ b/simpa/utils/libraries/heterogeneity_generator.py @@ -0,0 +1,327 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT + +import numpy as np +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 + + +class HeterogeneityGeneratorBase(object): + """ + This is the base class to define heterogeneous structure maps. + """ + + def __init__(self, xdim, ydim, zdim, spacing_mm, target_mean=None, + target_std=None, target_min=None, target_max=None, + eps=1e-5): + """ + :param xdim: the x dimension of the volume in voxels + :param ydim: the y dimension of the volume in voxels + :param zdim: the z dimension of the volume in voxels + :param spacing_mm: the spacing of the volume in mm + :param target_mean: (optional) the mean of the created heterogeneity map + :param target_std: (optional) the standard deviation of the created heterogeneity map + :param target_min: (optional) the minimum of the created heterogeneity map + :param target_max: (optional) the maximum of the created heterogeneity map + :param eps: (optional) the threshold when a re-normalisation should be triggered (default: 1e-5) + """ + self._xdim = xdim + self._ydim = ydim + self._zdim = zdim + self._spacing_mm = spacing_mm + self._mean = target_mean + self._std = target_std + self._min = target_min + self._max = target_max + self.eps = eps + + self.map = np.ones((self._xdim, self._ydim, self._zdim), dtype=float) + + def get_map(self): + self.normalise_map() + return self.map.astype(float) + + def normalise_map(self): + """ + If mean and std are set, then the data will be normalised to have the desired mean and the + desired standard deviation. + If min and max are set, then the data will be normalised to have the desired minimum and the + desired maximum value. + If all four values are set, then the data will be normalised to have the desired mean and the + desired standard deviation first. afterwards all values smaller than min will be ste to min and + all values larger than max will be set to max. + """ + # Testing mean mean/std normalisation needs to be done + if self._mean is not None and self._std is not None: + if (np.abs(np.mean(self.map) - self._mean) > self.eps or + np.abs(np.std(self.map) - self._std) > self.eps): + mean = np.mean(self.map) + std = np.std(self.map) + self.map = (self.map - mean) / std + self.map = (self.map * self._std) + self._mean + if self._min is not None and self._max is not None: + self.map[self.map < self._min] = self._min + self.map[self.map > self._max] = self._max + + # Testing if min max normalisation needs to be done + if self._min is None or self._max is None: + return + + if (np.abs(np.min(self.map) - self._min) < self.eps and + np.abs(np.max(self.map) - self._max) < self.eps): + return + + _min = np.min(self.map) + _max = np.max(self.map) + self.map = (self.map - _min) / (_max-_min) + self.map = (self.map * (self._max - self._min)) + self._min + + +class RandomHeterogeneity(HeterogeneityGeneratorBase): + """ + This heterogeneity generator represents a uniform random sampling between the given bounds. + Optionally, a Gaussian blur can be specified. Please not that a Gaussian blur will transform the random + distribution to a Gaussian. + """ + + def __init__(self, xdim, ydim, zdim, spacing_mm, gaussian_blur_size_mm=None, target_mean=None, target_std=None, + target_min=None, target_max=None, eps=1e-5): + """ + :param xdim: the x dimension of the volume in voxels + :param ydim: the y dimension of the volume in voxels + :param zdim: the z dimension of the volume in voxels + :param spacing_mm: the spacing of the volume in mm + :param gaussian_blur_size_mm: the size of the standard deviation for the Gaussian blur + :param target_mean: (optional) the mean of the created heterogeneity map + :param target_std: (optional) the standard deviation of the created heterogeneity map + :param target_min: (optional) the minimum of the created heterogeneity map + :param target_max: (optional) the maximum of the created heterogeneity map + :param eps: (optional) the threshold when a re-normalisation should be triggered (default: 1e-5) + """ + super().__init__(xdim, ydim, zdim, spacing_mm, target_mean, target_std, target_min, target_max, eps) + + self.map = np.random.random((xdim, ydim, zdim)) + if gaussian_blur_size_mm is not None: + _gaussian_blur_size_voxels = gaussian_blur_size_mm / spacing_mm + self.map = gaussian_filter(self.map, _gaussian_blur_size_voxels) + + +class BlobHeterogeneity(HeterogeneityGeneratorBase): + """ + This heterogeneity generator representes a blob-like random sampling between the given bounds using the + sklearn.datasets.make_blobs method. Please look into their documentation for optimising the given hyperparameters. + + """ + + def __init__(self, xdim, ydim, zdim, spacing_mm, num_centers=None, cluster_std=None, target_mean=None, + target_std=None, target_min=None, target_max=None, random_state=None): + """ + :param xdim: the x dimension of the volume in voxels + :param ydim: the y dimension of the volume in voxels + :param zdim: the z dimension of the volume in voxels + :param spacing_mm: the spacing of the volume in mm + :param num_centers: the number of blobs + :param cluster_std: the size of the blobs + :param target_mean: (optional) the mean of the created heterogeneity map + :param target_std: (optional) the standard deviation of the created heterogeneity map + :param target_min: (optional) the minimum of the created heterogeneity map + :param target_max: (optional) the maximum of the created heterogeneity map + """ + super().__init__(xdim, ydim, zdim, spacing_mm, target_mean, target_std, target_min, target_max) + + if num_centers is None: + num_centers = int(np.round(np.float_power((xdim * ydim * zdim) * spacing_mm, 1 / 3))) + + if cluster_std is None: + cluster_std = 1 + x, y = make_blobs(n_samples=(xdim * ydim * zdim) * 10, n_features=3, centers=num_centers, + random_state=random_state, cluster_std=cluster_std) + + self.map = np.histogramdd(x, bins=(xdim, ydim, zdim), range=((np.percentile(x[:, 0], 5), + np.percentile(x[:, 0], 95)), + (np.percentile(x[:, 1], 5), + np.percentile(x[:, 1], 95)), + (np.percentile(x[:, 2], 5), + np.percentile(x[:, 2], 95))))[0] + self.map = gaussian_filter(self.map, 5) + + +class ImageHeterogeneity(HeterogeneityGeneratorBase): + """ + This heterogeneity generator takes a pre-specified 2D image, currently only supporting numpy arrays, and uses them + as a map for heterogeneity within the tissue. + + Attributes: + map: the np array of the heterogeneity map transformed and augments to fill the area + """ + + def __init__(self, xdim: int, ydim: int, zdim: int, heterogeneity_image: np.ndarray, + spacing_mm: Union[int, float] = None, image_pixel_spacing_mm: Union[int, float] = None, + scaling_type: [None, str] = None, constant: Union[int, float] = 0, + crop_placement=('centre', 'centre'), target_mean: Union[int, float] = None, + target_std: Union[int, float] = None, target_min: Union[int, float] = None, + target_max: Union[int, float] = None): + """ + :param xdim: the x dimension of the volume in voxels + :param ydim: the y dimension of the volume in voxels + :param zdim: the z dimension of the volume in voxels + :param heterogeneity_image: the 2D prior image of the heterogeneity map + :param spacing_mm: the spacing of the volume in mm + :param image_pixel_spacing_mm: the pixel spacing of the image in mm (pixel spacing) + :param scaling_type: the scaling type of the heterogeneity map, with default being that no scaling occurs + OPTIONS: + TAGS.IMAGE_SCALING_SYMMETRIC: symmetric reflections of the image to span the area + TAGS.IMAGE_SCALING_STRETCH: stretch the image to span the area + TAGS.IMAGE_SCALING_WRAP: multiply the image to span the area + TAGS.IMAGE_SCALING_EDGE: continue the values at the edge of the area to fill the shape + TAGS.IMAGE_SCALING_CONSTANT: span the left-over area with a constant + :param constant: the scaling constant of the heterogeneity map, used only for scaling type 'constant' + WARNING: scaling constant must be in reference to the values in the heterogeneity_image + :param crop_placement: the placement of where the heterogeneity map is cropped + :param target_mean: (optional) the mean of the created heterogeneity map + :param target_std: (optional) the standard deviation of the created heterogeneity map + :param target_min: (optional) the minimum of the created heterogeneity map + :param target_max: (optional) the maximum of the created heterogeneity map + """ + super().__init__(xdim, ydim, zdim, spacing_mm, target_mean, target_std, target_min, target_max) + self.logger = Logger() + + if image_pixel_spacing_mm is None: + image_pixel_spacing_mm = spacing_mm + + (image_width_pixels, image_height_pixels) = heterogeneity_image.shape + [image_width_mm, image_height_mm] = np.array([image_width_pixels, image_height_pixels]) * image_pixel_spacing_mm + (xdim_mm, ydim_mm, zdim_mm) = np.array([xdim, ydim, zdim]) * spacing_mm + + wider = image_width_mm > xdim_mm + taller = image_height_mm > zdim_mm + + if taller or wider: + pixel_scaled_image = self.change_resolution(heterogeneity_image, spacing_mm=spacing_mm, + image_pixel_spacing_mm=image_pixel_spacing_mm) + cropped_image = self.crop_image(xdim, zdim, pixel_scaled_image, crop_placement) + + if taller and wider: + area_fill_image = cropped_image + else: + area_fill_image = self.upsize_to_fill_area(xdim, zdim, cropped_image, scaling_type, constant) + + else: + pixel_scaled_image = self.change_resolution(heterogeneity_image, spacing_mm=spacing_mm, + image_pixel_spacing_mm=image_pixel_spacing_mm) + area_fill_image = self.upsize_to_fill_area(xdim, zdim, pixel_scaled_image, scaling_type, constant) + + if scaling_type == Tags.IMAGE_SCALING_STRETCH: + area_fill_image = transform.resize(heterogeneity_image, output_shape=(xdim, zdim), mode='symmetric') + + self.map = np.repeat(area_fill_image[:, np.newaxis, :], ydim, axis=1) + + def upsize_to_fill_area(self, xdim: int, zdim: int, image: np.ndarray, scaling_type: Optional[str] = None, + constant: Union[int, float] = 0) -> np.ndarray: + """ + Fills an area with an image through various methods of expansion + :param xdim: the x dimension of the area to be filled in voxels + :param zdim: the z dimension of the area to be filled in voxels + :param image: the input image + :param scaling_type: the scaling type of the heterogeneity map, with default being that no scaling occurs + OPTIONS: + TAGS.IMAGE_SCALING_SYMMETRIC: symmetric reflections of the image to span the area + TAGS.IMAGE_SCALING_STRETCH: stretch the image to span the area + TAGS.IMAGE_SCALING_WRAP: multiply the image to span the area + TAGS.IMAGE_SCALING_EDGE: continue the values at the edge of the area to fill the shape + TAGS.IMAGE_SCALING_CONSTANT: span the left-over area with a constant + :param constant: the scaling constant of the heterogeneity map, used only for scaling type 'constant' + :return:A numpy array of size (xdim, zdim) containing the filled image expanded to fill the shape + """ + if scaling_type is None or scaling_type == Tags.IMAGE_SCALING_STRETCH: + scaled_image = image + elif scaling_type == Tags.IMAGE_SCALING_CONSTANT: + pad_left = int((xdim - len(image)) / 2) + pad_height = int(zdim - len(image[0])) + pad_right = xdim - pad_left - len(image) + scaled_image = np.pad(image, ((pad_left, pad_right), (0, pad_height)), + mode=scaling_type, constant_values=constant) + else: + pad_left = int((xdim - len(image)) / 2) + pad_height = int(zdim - len(image[0])) + pad_right = xdim - pad_left - len(image) + scaled_image = np.pad(image, ((pad_left, pad_right), (0, pad_height)), + mode=scaling_type) + + self.logger.warning("The input image has filled the area by using {} scaling type".format(scaling_type)) + return scaled_image + + def crop_image(self, xdim: int, zdim: int, image: np.ndarray, + crop_placement: Union[str, tuple] = Tags.CROP_POSITION_CENTRE) -> np.ndarray: + """ + Crop the image to fit specified dimensions xdim and zdim + :param xdim: the x dimension of the area to be filled in voxels + :param zdim: the z dimension of the area to be filled in voxels + :param image: the input image + :param crop_placement: the placement of where the heterogeneity map is cropped + OPTIONS: TAGS.CROP_PLACEMENT_[TOP,BOTTOM,LEFT,RIGHT,CENTRE,RANDOM] or position of left hand corner on image + :return: cropped image + """ + (image_width_pixels, image_height_pixels) = image.shape + crop_width = min(xdim, image_width_pixels) + crop_height = min(zdim, image_height_pixels) + + if isinstance(crop_placement, tuple): + if crop_placement[0] == Tags.CROP_POSITION_LEFT: + crop_horizontal = 0 + elif crop_placement[0] == Tags.CROP_POSITION_RIGHT: + crop_horizontal = image_width_pixels-crop_width-1 + elif crop_placement[0] == Tags.CROP_POSITION_CENTRE: + crop_horizontal = round((image_width_pixels - crop_width) / 2) + elif isinstance(crop_placement[0], int): + crop_horizontal = crop_placement[0] + + if crop_placement[1] == Tags.CROP_POSITION_TOP: + crop_vertical = 0 + elif crop_placement[1] == Tags.CROP_POSITION_BOTTOM: + crop_vertical = image_height_pixels-crop_height-1 + elif crop_placement[1] == Tags.CROP_POSITION_CENTRE: + crop_vertical = round((image_height_pixels - crop_height) / 2) + elif isinstance(crop_placement[1], int): + crop_vertical = crop_placement[1] + + elif isinstance(crop_placement, str): + if crop_placement == Tags.CROP_POSITION_CENTRE: + crop_horizontal = round((image_width_pixels - crop_width) / 2) + crop_vertical = round((image_height_pixels - crop_height) / 2) + elif crop_placement == Tags.CROP_POSITION_RANDOM: + crop_horizontal = image_width_pixels - crop_width + if crop_horizontal != 0: + crop_horizontal = np.random.randint(0, crop_horizontal) + crop_vertical = image_height_pixels - crop_height + if crop_vertical != 0: + crop_vertical = np.random.randint(0, crop_vertical) + + cropped_image = image[crop_horizontal: crop_horizontal + crop_width, crop_vertical: crop_vertical + crop_height] + + self.logger.warning( + "The input image has been cropped to the dimensions of the simulation volume ({} {})".format(xdim, zdim)) + return cropped_image + + def change_resolution(self, image: np.ndarray, spacing_mm: Union[int, float], + image_pixel_spacing_mm: Union[int, float]) -> np.ndarray: + """ + Method to change the resolution of an image + :param image: input image + :param image_pixel_spacing_mm: original image pixel spacing mm + :param spacing_mm: target pixel spacing mm + :return: image with new pixel spacing + """ + (image_width_pixels, image_height_pixels) = image.shape + [image_width_mm, image_height_mm] = np.array([image_width_pixels, image_height_pixels]) * image_pixel_spacing_mm + new_image_pixel_width = round(image_width_mm / spacing_mm) + new_image_pixel_height = round(image_height_mm / spacing_mm) + + 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)) diff --git a/simpa/utils/libraries/molecule_library.py b/simpa/utils/libraries/molecule_library.py index 2da45d9e..b73a1668 100644 --- a/simpa/utils/libraries/molecule_library.py +++ b/simpa/utils/libraries/molecule_library.py @@ -3,6 +3,7 @@ # SPDX-License-Identifier: MIT import numpy as np +import torch from simpa.utils import Tags from simpa.utils.tissue_properties import TissueProperties @@ -12,7 +13,7 @@ from simpa.utils.calculate import calculate_oxygenation, calculate_gruneisen_parameter_from_temperature, calculate_bvf from simpa.utils.serializer import SerializableSIMPAClass from simpa.utils.libraries.spectrum_library import AbsorptionSpectrumLibrary - +from simpa.log import Logger from typing import Optional, Union @@ -29,14 +30,16 @@ def __init__(self, segmentation_type: Optional[str] = None, molecular_compositio """ Initialize the MolecularComposition object. - :param segmentation_type: The type of segmentation. + :param segmentation_type: The segmentation class associated with this molecular composition. :type segmentation_type: str, optional - :param molecular_composition_settings: Settings for the molecular composition. + :param molecular_composition_settings: A settings dictionary or dict containing the molecules that constitute + this composition :type molecular_composition_settings: dict, optional """ super().__init__() self.segmentation_type = segmentation_type - self.internal_properties = TissueProperties() + self.internal_properties: TissueProperties = None + self.logger = Logger() if molecular_composition_settings is None: return @@ -45,19 +48,22 @@ def __init__(self, segmentation_type: Optional[str] = None, molecular_compositio for molecule_name in _keys: self.append(molecular_composition_settings[molecule_name]) - def update_internal_properties(self): + def update_internal_properties(self, settings): """ - Update the internal tissue properties based on the molecular composition. + Re-defines the internal properties of the molecular composition. + For each data field and molecule, a linear mixing model is used to arrive at the final parameters. Raises: AssertionError: If the total volume fraction of all molecules is not exactly 100%. """ - self.internal_properties = TissueProperties() + self.internal_properties = TissueProperties(settings) self.internal_properties[Tags.DATA_FIELD_SEGMENTATION] = self.segmentation_type self.internal_properties[Tags.DATA_FIELD_OXYGENATION] = calculate_oxygenation(self) self.internal_properties[Tags.DATA_FIELD_BLOOD_VOLUME_FRACTION] = calculate_bvf(self) - for molecule in self: - self.internal_properties.volume_fraction += molecule.volume_fraction + search_list = self.copy() + + for molecule in search_list: + self.internal_properties.volume_fraction += molecule.get_volume_fraction() self.internal_properties[Tags.DATA_FIELD_GRUNEISEN_PARAMETER] += \ molecule.volume_fraction * molecule.gruneisen_parameter self.internal_properties[Tags.DATA_FIELD_DENSITY] += molecule.volume_fraction * molecule.density @@ -66,23 +72,26 @@ def update_internal_properties(self): 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-3: - raise AssertionError("Invalid Molecular composition! The volume fractions of all molecules must be" - "exactly 100%!") + if (torch.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% somewhere!") + self.logger.warning("Some of the volume has not been filled by this molecular composition. Please check" + "that this is correct") - def get_properties_for_wavelength(self, wavelength: Union[int, float]) -> TissueProperties: + def get_properties_for_wavelength(self, settings, wavelength) -> TissueProperties: """ Get the tissue properties for a specific wavelength. :param wavelength: The wavelength to get properties for. :return: The updated tissue properties. """ - self.update_internal_properties() + self.update_internal_properties(settings) self.internal_properties[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 0 self.internal_properties[Tags.DATA_FIELD_SCATTERING_PER_CM] = 0 self.internal_properties[Tags.DATA_FIELD_ANISOTROPY] = 0 - - for molecule in self: + search_list = self.copy() + for molecule in search_list: self.internal_properties[Tags.DATA_FIELD_ABSORPTION_PER_CM] += \ (molecule.volume_fraction * molecule.spectrum.get_value_for_wavelength(wavelength)) self.internal_properties[Tags.DATA_FIELD_SCATTERING_PER_CM] += \ @@ -99,6 +108,7 @@ def serialize(self) -> dict: :return: The serialized molecular composition. """ dict_items = self.__dict__ + dict_items["internal_properties"] = None list_items = [molecule for molecule in self] return {"MolecularComposition": {"dict_items": dict_items, "list_items": list_items}} @@ -175,8 +185,11 @@ def __init__(self, name: str = None, if volume_fraction is None: volume_fraction = 0.0 - if not isinstance(volume_fraction, (int, float, np.int64)): - raise TypeError(f"The given volume_fraction was not of type float instead of {type(volume_fraction)}!") + if not isinstance(volume_fraction, (int, float, np.int64, np.ndarray)): + raise TypeError(f"The given volume_fraction was not of type float or array instead of " + f"{type(volume_fraction)}!") + if isinstance(volume_fraction, np.ndarray): + volume_fraction = torch.tensor(volume_fraction, dtype=torch.float32) self.volume_fraction = volume_fraction if scattering_spectrum is None: @@ -195,29 +208,37 @@ def __init__(self, name: str = None, if gruneisen_parameter is None: gruneisen_parameter = calculate_gruneisen_parameter_from_temperature( StandardProperties.BODY_TEMPERATURE_CELCIUS) - if not isinstance(gruneisen_parameter, (int, float)): + if not isinstance(gruneisen_parameter, (np.int32, np.int64, int, float, np.ndarray)): raise TypeError(f"The given gruneisen_parameter was not of type int or float instead " f"of {type(gruneisen_parameter)}!") + if isinstance(gruneisen_parameter, np.ndarray): + gruneisen_parameter = torch.tensor(gruneisen_parameter, dtype=torch.float32) self.gruneisen_parameter = gruneisen_parameter if density is None: density = StandardProperties.DENSITY_GENERIC - if not isinstance(density, (np.int32, np.int64, int, float)): + if not isinstance(density, (np.int32, np.int64, int, float, np.ndarray)): raise TypeError(f"The given density was not of type int or float instead of {type(density)}!") + if isinstance(density, np.ndarray): + density = torch.tensor(density, dtype=torch.float32) self.density = density if speed_of_sound is None: speed_of_sound = StandardProperties.SPEED_OF_SOUND_GENERIC - if not isinstance(speed_of_sound, (np.int32, np.int64, int, float)): + if not isinstance(speed_of_sound, (np.int32, np.int64, int, float, np.ndarray)): raise TypeError("The given speed_of_sound was not of type int or float instead of {}!" .format(type(speed_of_sound))) + if isinstance(speed_of_sound, np.ndarray): + speed_of_sound = torch.tensor(speed_of_sound, dtype=torch.float32) self.speed_of_sound = speed_of_sound if alpha_coefficient is None: alpha_coefficient = StandardProperties.ALPHA_COEFF_GENERIC - if not isinstance(alpha_coefficient, (int, float)): + if not isinstance(alpha_coefficient, (np.int32, np.int64, int, float, np.ndarray)): raise TypeError("The given alpha_coefficient was not of type int or float instead of {}!" .format(type(alpha_coefficient))) + if isinstance(alpha_coefficient, np.ndarray): + alpha_coefficient = torch.tensor(alpha_coefficient, dtype=torch.float32) self.alpha_coefficient = alpha_coefficient def __eq__(self, other) -> bool: @@ -241,6 +262,9 @@ def __eq__(self, other) -> bool: else: return super().__eq__(other) + def get_volume_fraction(self): + return self.volume_fraction + def serialize(self): """ Serialize the molecule to a dictionary. @@ -280,7 +304,7 @@ class MoleculeLibrary(object): """ # Main absorbers @staticmethod - def water(volume_fraction: float = 1.0) -> Molecule: + def water(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: """ Create a water molecule with predefined properties. @@ -300,7 +324,7 @@ def water(volume_fraction: float = 1.0) -> Molecule: ) @staticmethod - def oxyhemoglobin(volume_fraction: float = 1.0) -> Molecule: + def oxyhemoglobin(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: """ Create an oxyhemoglobin molecule with predefined properties. @@ -319,7 +343,7 @@ def oxyhemoglobin(volume_fraction: float = 1.0) -> Molecule: ) @staticmethod - def deoxyhemoglobin(volume_fraction: float = 1.0) -> Molecule: + def deoxyhemoglobin(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: """ Create a deoxyhemoglobin molecule with predefined properties. @@ -338,7 +362,7 @@ def deoxyhemoglobin(volume_fraction: float = 1.0) -> Molecule: ) @staticmethod - def melanin(volume_fraction: float = 1.0) -> Molecule: + def melanin(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: """ Create a melanin molecule with predefined properties. @@ -358,7 +382,7 @@ def melanin(volume_fraction: float = 1.0) -> Molecule: ) @staticmethod - def fat(volume_fraction: float = 1.0) -> Molecule: + def fat(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: """ Create a fat molecule with predefined properties. @@ -379,7 +403,7 @@ def fat(volume_fraction: float = 1.0) -> Molecule: # Scatterers @staticmethod def constant_scatterer(scattering_coefficient: float = 100.0, anisotropy: float = 0.9, - volume_fraction: float = 1.0) -> Molecule: + volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: """ Create a constant scatterer molecule with predefined properties. @@ -400,7 +424,7 @@ def constant_scatterer(scattering_coefficient: float = 100.0, anisotropy: float ) @staticmethod - def soft_tissue_scatterer(volume_fraction: float = 1.0) -> Molecule: + def soft_tissue_scatterer(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: """ Create a soft tissue scatterer molecule with predefined properties. @@ -419,7 +443,7 @@ def soft_tissue_scatterer(volume_fraction: float = 1.0) -> Molecule: ) @staticmethod - def muscle_scatterer(volume_fraction: float = 1.0) -> Molecule: + def muscle_scatterer(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: """ Create a muscle scatterer molecule with predefined properties. @@ -438,7 +462,7 @@ def muscle_scatterer(volume_fraction: float = 1.0) -> Molecule: ) @staticmethod - def epidermal_scatterer(volume_fraction: float = 1.0) -> Molecule: + def epidermal_scatterer(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: """ Create an epidermal scatterer molecule with predefined properties. @@ -458,7 +482,7 @@ def epidermal_scatterer(volume_fraction: float = 1.0) -> Molecule: ) @staticmethod - def dermal_scatterer(volume_fraction: float = 1.0) -> Molecule: + def dermal_scatterer(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: """ Create a dermal scatterer molecule with predefined properties. @@ -479,7 +503,7 @@ def dermal_scatterer(volume_fraction: float = 1.0) -> Molecule: ) @staticmethod - def bone(volume_fraction: float = 1.0) -> Molecule: + def bone(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: """ Create a bone molecule with predefined properties. @@ -499,7 +523,7 @@ def bone(volume_fraction: float = 1.0) -> Molecule: ) @staticmethod - def mediprene(volume_fraction: float = 1.0) -> Molecule: + def mediprene(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: """ Create a mediprene molecule with predefined properties. @@ -518,7 +542,7 @@ def mediprene(volume_fraction: float = 1.0) -> Molecule: ) @staticmethod - def heavy_water(volume_fraction: float = 1.0) -> Molecule: + def heavy_water(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: """ Create a heavy water molecule with predefined properties. @@ -539,7 +563,7 @@ def heavy_water(volume_fraction: float = 1.0) -> Molecule: ) @staticmethod - def air(volume_fraction: float = 1.0) -> Molecule: + def air(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: """ Create an air molecule with predefined properties. diff --git a/simpa/utils/libraries/spectrum_library.py b/simpa/utils/libraries/spectrum_library.py index 6e5ef4c1..8ed9d196 100644 --- a/simpa/utils/libraries/spectrum_library.py +++ b/simpa/utils/libraries/spectrum_library.py @@ -7,6 +7,8 @@ import glob import numpy as np import matplotlib.pylab as plt +import torch +from scipy import interpolate from simpa.utils.libraries.literature_values import OpticalTissueProperties from simpa.utils.serializer import SerializableSIMPAClass @@ -34,18 +36,22 @@ def __init__(self, spectrum_name: str, wavelengths: np.ndarray, values: np.ndarr :raises ValueError: If the shape of wavelengths does not match the shape of values. """ + if isinstance(values, np.ndarray): + values = torch.from_numpy(values) + wavelengths = torch.from_numpy(wavelengths) self.spectrum_name = spectrum_name self.wavelengths = wavelengths - self.max_wavelength = int(np.max(wavelengths)) - self.min_wavelength = int(np.min(wavelengths)) + self.max_wavelength = int(torch.max(wavelengths)) + self.min_wavelength = int(torch.min(wavelengths)) self.values = values - if np.shape(wavelengths) != np.shape(values): + if torch.Tensor.size(wavelengths) != torch.Tensor.size(values): raise ValueError("The shape of the wavelengths and the values did not match: " + - str(np.shape(wavelengths)) + " vs " + str(np.shape(values))) + str(torch.Tensor.size(wavelengths)) + " vs " + str(torch.Tensor.size(values))) - new_wavelengths = np.arange(self.min_wavelength, self.max_wavelength+1, 1) - self.values_interp = np.interp(new_wavelengths, self.wavelengths, self.values) + new_wavelengths = torch.arange(self.min_wavelength, self.max_wavelength+1, 1) + new_absorptions_function = interpolate.interp1d(self.wavelengths, self.values) + self.values_interp = new_absorptions_function(new_wavelengths) def get_value_over_wavelength(self) -> np.ndarray: """ diff --git a/simpa/utils/libraries/structure_library/HorizontalLayerStructure.py b/simpa/utils/libraries/structure_library/HorizontalLayerStructure.py index 9a439e5d..57097731 100644 --- a/simpa/utils/libraries/structure_library/HorizontalLayerStructure.py +++ b/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, diff --git a/simpa/utils/libraries/structure_library/StructureBase.py b/simpa/utils/libraries/structure_library/StructureBase.py index 399bb988..22c39790 100644 --- a/simpa/utils/libraries/structure_library/StructureBase.py +++ b/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.update_internal_properties() + 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. @@ -80,6 +85,12 @@ def fill_internal_volume(self): indices, values = self.get_enclosed_indices() self.geometrical_volume[indices] = values + def get_volume_fractions(self): + """ + Get the volume fraction this structure takes per voxel. + """ + return self.geometrical_volume + @abstractmethod def get_enclosed_indices(self): """ @@ -89,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. @@ -97,13 +108,25 @@ def get_params_from_settings(self, single_structure_settings): """ pass - def properties_for_wavelength(self, wavelength) -> TissueProperties: + def properties_for_wavelength(self, settings, wavelength) -> TissueProperties: """ Returns the values corresponding to each optical/acoustic property used in SIMPA. + :param settings: The global settings that contains the info on the volume dimensions. :param wavelength: Wavelength of the queried properties :return: optical/acoustic properties """ - return self.molecule_composition.get_properties_for_wavelength(wavelength) + 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: diff --git a/simpa/utils/libraries/structure_library/__init__.py b/simpa/utils/libraries/structure_library/__init__.py index 7d821601..db527a4a 100644 --- a/simpa/utils/libraries/structure_library/__init__.py +++ b/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]] diff --git a/simpa/utils/libraries/tissue_library.py b/simpa/utils/libraries/tissue_library.py index 10965de1..8c4393e1 100644 --- a/simpa/utils/libraries/tissue_library.py +++ b/simpa/utils/libraries/tissue_library.py @@ -1,7 +1,8 @@ # SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT -import typing +import numpy as np +from typing import Union, List, Optional from simpa.utils import OpticalTissueProperties, SegmentationClasses, StandardProperties, MolecularCompositionGenerator from simpa.utils import Molecule @@ -12,18 +13,15 @@ from simpa.utils.calculate import randomize_uniform from simpa.utils.libraries.spectrum_library import AbsorptionSpectrumLibrary -from typing import Union, List -import torch - class TissueLibrary(object): """ A library, returning molecular compositions for various typical tissue segmentations. """ - def get_blood_volume_fractions(self, oxygenation: Union[float, int, torch.Tensor] = 1e-10, - blood_volume_fraction: Union[float, int, torch.Tensor] = 1e-10)\ - -> List[Union[int, float, torch.Tensor]]: + def get_blood_volume_fractions(self, oxygenation: Union[float, int, np.ndarray] = 1e-10, + blood_volume_fraction: Union[float, int, np.ndarray] = 1e-10)\ + -> List[Union[int, float, np.ndarray]]: """ A function that returns the volume fraction of the oxygenated and deoxygenated blood. :param oxygenation: The oxygenation level of the blood volume fraction (as a decimal). @@ -34,8 +32,8 @@ def get_blood_volume_fractions(self, oxygenation: Union[float, int, torch.Tensor """ return [blood_volume_fraction*oxygenation, blood_volume_fraction*(1-oxygenation)] - def constant(self, mua: Union[float, int, torch.Tensor] = 1e-10, mus: Union[float, int, torch.Tensor] = 1e-10, - g: Union[float, int, torch.Tensor] = 0) -> MolecularComposition: + def constant(self, mua: Union[float, int, np.ndarray] = 1e-10, mus: Union[float, int, np.ndarray] = 1e-10, + g: Union[float, int, np.ndarray] = 0) -> MolecularComposition: """ A function returning a molecular composition as specified by the user. Typically intended for the use of wanting specific mua, mus and g values. @@ -56,7 +54,7 @@ def generic_tissue(self, mua: Spectrum = AbsorptionSpectrumLibrary().CONSTANT_ABSORBER_ARBITRARY(1e-10), mus: Spectrum = AbsorptionSpectrumLibrary().CONSTANT_ABSORBER_ARBITRARY(1e-10), g: Spectrum = AbsorptionSpectrumLibrary().CONSTANT_ABSORBER_ARBITRARY(1e-10), - molecule_name: typing.Optional[str] = "generic_tissue") -> MolecularComposition: + molecule_name: Optional[str] = "generic_tissue") -> MolecularComposition: """ Returns a generic issue defined by the provided optical parameters. @@ -79,10 +77,14 @@ def generic_tissue(self, anisotropy_spectrum=g)) .get_molecular_composition(SegmentationClasses.GENERIC)) - def muscle(self, oxygenation: Union[float, int, torch.Tensor] = 0.175, - blood_volume_fraction: Union[float, int, torch.Tensor] = 0.06) -> MolecularComposition: + def muscle(self, oxygenation: Union[float, int, np.ndarray] = 0.175, + blood_volume_fraction: Union[float, int, np.ndarray] = 0.06) -> MolecularComposition: """ - + Create a molecular composition mimicking that of muscle + :param oxygenation: The oxygenation level of the blood volume fraction (as a decimal). + Default: 0.175 + :param blood_volume_fraction: The total blood volume fraction (including oxygenated and deoxygenated blood). + Default: 0.06 :return: a settings dictionary containing all min and max parameters fitting for muscle tissue. """ @@ -91,6 +93,16 @@ def muscle(self, oxygenation: Union[float, int, torch.Tensor] = 0.175, # Get the water volume fraction water_volume_fraction = OpticalTissueProperties.WATER_VOLUME_FRACTION_HUMAN_BODY + if isinstance(blood_volume_fraction, np.ndarray): + if (blood_volume_fraction + water_volume_fraction - 1 > 1e-5).any(): + raise AssertionError(f"Blood volume fraction too large, must be less than {1 - water_volume_fraction}" + f" everywhere to leave space for water") + + else: + if blood_volume_fraction + water_volume_fraction - 1 > 1e-5: + raise AssertionError(f"Blood volume fraction too large, must be less than {1 - water_volume_fraction}" + f"everywhere to leave space for water") + custom_water = MOLECULE_LIBRARY.water(water_volume_fraction) custom_water.anisotropy_spectrum = AnisotropySpectrumLibrary.CONSTANT_ANISOTROPY_ARBITRARY( OpticalTissueProperties.STANDARD_ANISOTROPY - 0.005) @@ -111,8 +123,8 @@ def muscle(self, oxygenation: Union[float, int, torch.Tensor] = 0.175, .append(custom_water) .get_molecular_composition(SegmentationClasses.MUSCLE)) - def soft_tissue(self, oxygenation: Union[float, int, torch.Tensor] = OpticalTissueProperties.BACKGROUND_OXYGENATION, - blood_volume_fraction: Union[float, int, torch.Tensor] = OpticalTissueProperties.BLOOD_VOLUME_FRACTION_MUSCLE_TISSUE) -> MolecularComposition: + def soft_tissue(self, oxygenation: Union[float, int, np.ndarray] = OpticalTissueProperties.BACKGROUND_OXYGENATION, + blood_volume_fraction: Union[float, int, np.ndarray] = OpticalTissueProperties.BLOOD_VOLUME_FRACTION_MUSCLE_TISSUE) -> MolecularComposition: """ IMPORTANT! This tissue is not tested and it is not based on a specific real tissue type. It is a mixture of muscle (mostly optical properties) and water (mostly acoustic properties). @@ -129,6 +141,16 @@ def soft_tissue(self, oxygenation: Union[float, int, torch.Tensor] = OpticalTiss # Get the water volume fraction water_volume_fraction = OpticalTissueProperties.WATER_VOLUME_FRACTION_HUMAN_BODY + if isinstance(blood_volume_fraction, np.ndarray): + if (blood_volume_fraction + water_volume_fraction - 1 > 1e-5).any(): + raise AssertionError(f"Blood volume fraction too large, must be less than {1 - water_volume_fraction}" + f"everywhere to leave space for water") + + else: + if blood_volume_fraction + water_volume_fraction - 1 > 1e-5: + raise AssertionError(f"Blood volume fraction too large, must be less than {1 - water_volume_fraction}" + f"everywhere to leave space for water") + custom_water = MOLECULE_LIBRARY.water(water_volume_fraction) custom_water.anisotropy_spectrum = AnisotropySpectrumLibrary.CONSTANT_ANISOTROPY_ARBITRARY( OpticalTissueProperties.STANDARD_ANISOTROPY - 0.005) @@ -149,20 +171,21 @@ def soft_tissue(self, oxygenation: Union[float, int, torch.Tensor] = OpticalTiss .append(custom_water) .get_molecular_composition(SegmentationClasses.SOFT_TISSUE)) - def epidermis(self, melanosom_volume_fraction: Union[float, int, torch.Tensor] = 0.014) -> MolecularComposition: + def epidermis(self, melanin_volume_fraction: Union[float, int, np.ndarray] = 0.014) -> MolecularComposition: """ - + Create a molecular composition mimicking that of dermis + :param melanin_volume_fraction: the total volume fraction of melanin :return: a settings dictionary containing all min and max parameters fitting for epidermis tissue. """ # generate the tissue dictionary return (MolecularCompositionGenerator() - .append(MOLECULE_LIBRARY.melanin(melanosom_volume_fraction)) - .append(MOLECULE_LIBRARY.epidermal_scatterer(1 - melanosom_volume_fraction)) + .append(MOLECULE_LIBRARY.melanin(melanin_volume_fraction)) + .append(MOLECULE_LIBRARY.epidermal_scatterer(1 - melanin_volume_fraction)) .get_molecular_composition(SegmentationClasses.EPIDERMIS)) - def dermis(self, oxygenation: Union[float, int, torch.Tensor] = 0.5, - blood_volume_fraction: Union[float, int, torch.Tensor] = 0.002) -> MolecularComposition: + def dermis(self, oxygenation: Union[float, int, np.ndarray] = 0.5, + blood_volume_fraction: Union[float, int, np.ndarray] = 0.002) -> MolecularComposition: """ Create a molecular composition mimicking that of dermis :param oxygenation: The oxygenation level of the blood volume fraction (as a decimal). @@ -183,8 +206,8 @@ def dermis(self, oxygenation: Union[float, int, torch.Tensor] = 0.5, .get_molecular_composition(SegmentationClasses.DERMIS)) def subcutaneous_fat(self, - oxygenation: Union[float, int, torch.Tensor] = OpticalTissueProperties.BACKGROUND_OXYGENATION, - blood_volume_fraction: Union[float, int, torch.Tensor] + oxygenation: Union[float, int, np.ndarray] = OpticalTissueProperties.BACKGROUND_OXYGENATION, + blood_volume_fraction: Union[float, int, np.ndarray] = OpticalTissueProperties.BLOOD_VOLUME_FRACTION_MUSCLE_TISSUE) -> MolecularComposition: """ Create a molecular composition mimicking that of subcutaneous fat @@ -215,7 +238,7 @@ def subcutaneous_fat(self, .append(MOLECULE_LIBRARY.water(water_volume_fraction)) .get_molecular_composition(SegmentationClasses.FAT)) - def blood(self, oxygenation: Union[float, int, torch.Tensor, None] = None) -> MolecularComposition: + def blood(self, oxygenation: Union[float, int, np.ndarray, None] = None) -> MolecularComposition: """ Create a molecular composition mimicking that of blood :param oxygenation: The oxygenation level of the blood(as a decimal). @@ -281,8 +304,9 @@ def ultrasound_gel(self) -> MolecularComposition: .append(MOLECULE_LIBRARY.water()) .get_molecular_composition(SegmentationClasses.ULTRASOUND_GEL)) - def lymph_node(self, oxygenation: Union[float, int, torch.Tensor, None] = None, - blood_volume_fraction: Union[float, int, torch.Tensor, None] = None) -> MolecularComposition: + def lymph_node(self, oxygenation: Union[float, int, np.ndarray] = OpticalTissueProperties.LYMPH_NODE_OXYGENATION, + blood_volume_fraction: Union[float, int, np.ndarray] = + OpticalTissueProperties.BLOOD_VOLUME_FRACTION_LYMPH_NODE) -> MolecularComposition: """ IMPORTANT! This tissue is not tested and it is not based on a specific real tissue type. It is a mixture of oxyhemoglobin, deoxyhemoglobin, and lymph node customized water. @@ -293,14 +317,6 @@ def lymph_node(self, oxygenation: Union[float, int, torch.Tensor, None] = None, :return: a settings dictionary fitting for generic lymph node tissue. """ - # Determine muscle oxygenation - if oxygenation is None: - oxygenation = OpticalTissueProperties.LYMPH_NODE_OXYGENATION - - # Get the blood volume fractions for oxyhemoglobin and deoxyhemoglobin - if blood_volume_fraction is None: - blood_volume_fraction = OpticalTissueProperties.BLOOD_VOLUME_FRACTION_LYMPH_NODE - [fraction_oxy, fraction_deoxy] = self.get_blood_volume_fractions(oxygenation, blood_volume_fraction) # Get the water volume fraction diff --git a/simpa/utils/matlab.py b/simpa/utils/matlab.py index 62dc1f22..1459ba0b 100644 --- a/simpa/utils/matlab.py +++ b/simpa/utils/matlab.py @@ -21,7 +21,7 @@ def generate_matlab_cmd(matlab_binary_path: str, simulation_script_path: str, da :return: list of command parts :rtype: List[str] """ - + # get path of calling script to add to matlab path base_script_path = os.path.dirname(os.path.abspath(inspect.stack()[1].filename)) # ensure data path is an absolute path diff --git a/simpa/utils/settings.py b/simpa/utils/settings.py index a6a80a30..447ff48d 100644 --- a/simpa/utils/settings.py +++ b/simpa/utils/settings.py @@ -69,6 +69,32 @@ def __delitem__(self, key): except KeyError: raise KeyError("The key '{}' is not in the Settings dictionary".format(key)) from None + def get_volume_dimensions_voxels(self): + """ + returns: tuple + the x, y, and z dimension of the volumes as a tuple + the volume dimension gets rounded after converting from a mm grid to a voxel grid of unit Tags.SPACING_MM. + """ + if Tags.SPACING_MM not in self: + raise AssertionError("The simpa.Tags.SPACING_MM tag must be set " + "to compute the volume dimensions in voxels") + if Tags.DIM_VOLUME_X_MM not in self: + raise AssertionError("The simpa.Tags.DIM_VOLUME_X_MM tag must be " + "set to compute the volume dimensions in voxels") + if Tags.DIM_VOLUME_Y_MM not in self: + raise AssertionError("The simpa.Tags.DIM_VOLUME_Y_MM tag must be " + "set to compute the volume dimensions in voxels") + if Tags.DIM_VOLUME_Z_MM not in self: + raise AssertionError("The simpa.Tags.DIM_VOLUME_Z_MM tag must be " + "set to compute the volume dimensions in voxels") + + voxel_spacing = self[Tags.SPACING_MM] + volume_x_dim = int(round(self[Tags.DIM_VOLUME_X_MM] / voxel_spacing)) + volume_y_dim = int(round(self[Tags.DIM_VOLUME_Y_MM] / voxel_spacing)) + volume_z_dim = int(round(self[Tags.DIM_VOLUME_Z_MM] / voxel_spacing)) + + return volume_x_dim, volume_y_dim, volume_z_dim + def get_optical_settings(self): """" Returns the settings for the optical forward model that are saved in this settings dictionary diff --git a/simpa/utils/tags.py b/simpa/utils/tags.py index 2e1be847..18e28778 100644 --- a/simpa/utils/tags.py +++ b/simpa/utils/tags.py @@ -1466,7 +1466,7 @@ class Tags: """ Identifier for the diffuse reflectance values at the surface of the volume (interface to 0-values voxels) Usage: simpa.core.simulation_modules.optical_simulation_module.optical_forward_model_mcx_reflectance_adapter - """ + """ DATA_FIELD_DIFFUSE_REFLECTANCE_POS = "diffuse_reflectance_pos" """ @@ -1489,6 +1489,75 @@ class Tags: Usage: simpa.core.simulation_modules.optical_simulation_module.optical_forward_model_mcx_reflectance_adapter """ + IMAGE_SCALING_SYMMETRIC = "symmetric" + """ + Flag indicating the use of reflection on edges during interpolation when rescaling an image + Usage: simpa.utils.libraries.structure_library.heterogeneity_generator + """ + + IMAGE_SCALING_STRETCH = "stretch" + """ + Flag indicating the use of reflection on edges during interpolation when rescaling an image + At the boundary, the image will reflect to fill the area + Usage: simpa.utils.libraries.structure_library.heterogeneity_generator + """ + + IMAGE_SCALING_WRAP = "wrap" + """ + Flag indicating tessellating during interpolation when rescaling an image + Usage: simpa.utils.libraries.structure_library.heterogeneity_generator + """ + + IMAGE_SCALING_CONSTANT = "constant" + """ + Flag indicating the use of a constant on edges during interpolation when rescaling an image + The rest of the area will be filled by a constant value + Usage: simpa.utils.libraries.structure_library.heterogeneity_generator + """ + + IMAGE_SCALING_EDGE = "edge" + """ + Flag indicating the expansion of the edges during interpolation when rescaling an image + The edge value will continue across the area + Usage: simpa.utils.libraries.structure_library.heterogeneity_generator + """ + + CROP_POSITION_TOP = "top" + """ + Flag indicating the crop position: along top edge of image + Usage: simpa.utils.libraries.structure_library.heterogeneity_generator + """ + + CROP_POSITION_BOTTOM = "bottom" + """ + Flag indicating the crop position: along bottom edge of image + Usage: simpa.utils.libraries.structure_library.heterogeneity_generator + """ + + CROP_POSITION_CENTRE = "centre" + """ + Flag indicating the crop position: along centre edge of image + Usage: simpa.utils.libraries.structure_library.heterogeneity_generator + """ + + CROP_POSITION_LEFT = "left" + """ + Flag indicating the crop position: along left-hand edge of image + Usage: simpa.utils.libraries.structure_library.heterogeneity_generator + """ + + CROP_POSITION_RIGHT = "right" + """ + Flag indicating the crop position: along right-hand edge of image + Usage: simpa.utils.libraries.structure_library.heterogeneity_generator + """ + + CROP_POSITION_RANDOM = "random" + """ + Flag indicating the crop position: random placement + Usage: simpa.utils.libraries.structure_library.heterogeneity_generator + """ + SIMPA_SAVE_PATH_VARNAME = "SIMPA_SAVE_PATH" """ Identifier for the environment variable that defines where the results generated with SIMPA will be sotred @@ -1511,3 +1580,8 @@ class Tags: It is assumed that if flags are specified multiple times the flag provided last is considered. This can for example be used to override predefined flags. """ + + VOLUME_FRACTION = "volume_fraction" + """ + Identifier for the volume fraction for the simulation + """ diff --git a/simpa/utils/tissue_properties.py b/simpa/utils/tissue_properties.py index 02f35013..7e84fdf1 100644 --- a/simpa/utils/tissue_properties.py +++ b/simpa/utils/tissue_properties.py @@ -3,13 +3,25 @@ # SPDX-License-Identifier: MIT from simpa.utils.constants import property_tags +from simpa.utils import Settings +from simpa.utils.processing_device import get_processing_device +import torch class TissueProperties(dict): + """ + The tissue properties contain a volumetric representation of each tissue parameter currently + modelled in the SIMPA framework. - def __init__(self): + It is a dictionary that is populated with each of the parameters. + The values of the parameters can be either numbers or numpy arrays. + It also contains a volume fraction field. + """ + + def __init__(self, settings: Settings): super().__init__() - self.volume_fraction = 0 + volume_x_dim, volume_y_dim, volume_z_dim = settings.get_volume_dimensions_voxels() + self.volume_fraction = torch.zeros((volume_x_dim, volume_y_dim, volume_z_dim), dtype=torch.float32) for key in property_tags: self[key] = 0 diff --git a/simpa_examples/__init__.py b/simpa_examples/__init__.py index 92a0960b..49f13e4d 100644 --- a/simpa_examples/__init__.py +++ b/simpa_examples/__init__.py @@ -9,3 +9,4 @@ from simpa_examples.optical_and_acoustic_simulation import run_optical_and_acoustic_simulation from simpa_examples.perform_image_reconstruction import run_perform_image_reconstruction from simpa_examples.perform_iterative_qPAI_reconstruction import run_perform_iterative_qPAI_reconstruction +from simpa_examples.segmentation_loader import run_segmentation_loader diff --git a/simpa_examples/benchmarking/performance_check.py b/simpa_examples/benchmarking/performance_check.py index 880b0d36..2f969837 100644 --- a/simpa_examples/benchmarking/performance_check.py +++ b/simpa_examples/benchmarking/performance_check.py @@ -31,13 +31,10 @@ def run_benchmarking_tests(spacing=0.4, profile: str = "TIME", savefolder: str = examples = [simpa_examples.run_linear_unmixing, simpa_examples.run_minimal_optical_simulation, simpa_examples.run_minimal_optical_simulation_uniform_cube, simpa_examples.run_msot_invision_simulation, simpa_examples.run_optical_and_acoustic_simulation, - simpa_examples.run_perform_iterative_qPAI_reconstruction] + simpa_examples.run_perform_iterative_qPAI_reconstruction, simpa_examples.run_segmentation_loader] for example in examples: - try: - example(spacing=spacing, path_manager=None, visualise=False) - except AttributeError: - print("simulation cannot be run on {} with spacing {}".format(example, spacing)) + example(spacing=spacing, path_manager=None, visualise=False) if __name__ == "__main__": diff --git a/simpa_examples/benchmarking/run_benchmarking.sh b/simpa_examples/benchmarking/run_benchmarking.sh index 96a1e904..b9757194 100644 --- a/simpa_examples/benchmarking/run_benchmarking.sh +++ b/simpa_examples/benchmarking/run_benchmarking.sh @@ -1,14 +1,16 @@ #!/bin/bash +set -e + help() { echo "Usage: calculate benchmarking for [options]" echo "For further details see readme" echo "Number of examples can be selected in performance_check.py" echo "For comparable benchmarks, please use default" echo "Options:" -echo " -i, --init First spacing to benchmark: default = 0.2mm" -echo " -c, --cease Final spacing to benchmark: default = 0.4mm" -echo " -s, --step Step between spacings: default = 0.1mm" +echo " -i, --init First spacing to benchmark: default = 0.15mm" +echo " -c, --cease Final spacing to benchmark: default = 0.25mm" +echo " -s, --step Step between spacings: default = 0.05mm" echo " -f, --file Where to store the output files: default save in current directory; 'print' prints it in console" echo " -t, --time Profile times taken: if no profile, all are set" echo " -g, --gpu Profile GPU usage: if no profile, all are set" @@ -57,15 +59,15 @@ shift 1 done if [ "$start" == 0 ]; then - start=0.2 + start=0.15 fi if [ "$stop" == 0 ]; then - stop=0.4 + stop=0.25 fi if [ "$step" == 0 ]; then - step=0.1 + step=0.05 fi if [ ${#profiles[@]} -eq 0 ]; then diff --git a/simpa_examples/minimal_optical_simulation_heterogeneous_tissue.py b/simpa_examples/minimal_optical_simulation_heterogeneous_tissue.py new file mode 100644 index 00000000..874ee093 --- /dev/null +++ b/simpa_examples/minimal_optical_simulation_heterogeneous_tissue.py @@ -0,0 +1,177 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT +# +# This file should give an example how heterogeneity can be used in the SIMPA framework to model non-homogeneous +# structures. +# The result of this simulation should be a optical simulation results with spatially varying absorption coefficients +# that is induced by a heterogeneous map of blood volume fraction and oxygenation both in the simulated vessel and +# the background. + +from simpa import Tags +import simpa as sp +import numpy as np + +# FIXME temporary workaround for newest Intel architectures +import os +os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE" + +# TODO: Please make sure that a valid path_config.env file is located in your home directory, or that you +# point to the correct file in the PathManager(). +path_manager = sp.PathManager() + +VOLUME_TRANSDUCER_DIM_IN_MM = 60 +VOLUME_PLANAR_DIM_IN_MM = 30 +VOLUME_HEIGHT_IN_MM = 60 +SPACING = 0.5 +RANDOM_SEED = 471 +VOLUME_NAME = "MyVolumeName_"+str(RANDOM_SEED) +SAVE_REFLECTANCE = False +SAVE_PHOTON_DIRECTION = False + +# If VISUALIZE is set to True, the simulation result will be plotted +VISUALIZE = True + + +def create_example_tissue(settings): + """ + This is a very simple example script of how to create a tissue definition. + It contains a muscular background, an epidermis layer on top of the muscles + and a blood vessel. + """ + dim_x, dim_y, dim_z = settings.get_volume_dimensions_voxels() + tissue_library = sp.TISSUE_LIBRARY + background_dictionary = sp.Settings() + background_dictionary[Tags.MOLECULE_COMPOSITION] = tissue_library.constant(1e-4, 1e-4, 0.9) + background_dictionary[Tags.STRUCTURE_TYPE] = Tags.BACKGROUND + + muscle_dictionary = sp.Settings() + muscle_dictionary[Tags.PRIORITY] = 1 + 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( + 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 + muscle_dictionary[Tags.ADHERE_TO_DEFORMATION] = True + muscle_dictionary[Tags.STRUCTURE_TYPE] = Tags.HORIZONTAL_LAYER_STRUCTURE + + vessel_1_dictionary = sp.Settings() + vessel_1_dictionary[Tags.PRIORITY] = 3 + vessel_1_dictionary[Tags.STRUCTURE_START_MM] = [VOLUME_TRANSDUCER_DIM_IN_MM/2, + 10, + VOLUME_HEIGHT_IN_MM/2] + vessel_1_dictionary[Tags.STRUCTURE_END_MM] = [VOLUME_TRANSDUCER_DIM_IN_MM/2, + 12, + VOLUME_HEIGHT_IN_MM/2] + vessel_1_dictionary[Tags.STRUCTURE_RADIUS_MM] = 3 + vessel_1_dictionary[Tags.MOLECULE_COMPOSITION] = tissue_library.blood( + oxygenation=sp.RandomHeterogeneity(dim_x, dim_y, dim_z, SPACING, + target_min=0.9, target_max=1.0).get_map()) + vessel_1_dictionary[Tags.CONSIDER_PARTIAL_VOLUME] = True + vessel_1_dictionary[Tags.STRUCTURE_TYPE] = Tags.CIRCULAR_TUBULAR_STRUCTURE + + epidermis_dictionary = sp.Settings() + epidermis_dictionary[Tags.PRIORITY] = 8 + 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( + 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 + + tissue_dict = sp.Settings() + tissue_dict[Tags.BACKGROUND] = background_dictionary + tissue_dict["muscle"] = muscle_dictionary + tissue_dict["epidermis"] = epidermis_dictionary + tissue_dict["vessel_1"] = vessel_1_dictionary + return tissue_dict + + +# Seed the numpy random configuration prior to creating the global_settings file in +# order to ensure that the same volume +# is generated with the same random seed every time. + +np.random.seed(RANDOM_SEED) + +general_settings = { + # These parameters set the general properties of the simulated volume + Tags.RANDOM_SEED: RANDOM_SEED, + Tags.VOLUME_NAME: VOLUME_NAME, + Tags.SIMULATION_PATH: path_manager.get_hdf5_file_save_path(), + Tags.SPACING_MM: SPACING, + Tags.DIM_VOLUME_Z_MM: VOLUME_HEIGHT_IN_MM, + Tags.DIM_VOLUME_X_MM: VOLUME_TRANSDUCER_DIM_IN_MM, + Tags.DIM_VOLUME_Y_MM: VOLUME_PLANAR_DIM_IN_MM, + Tags.WAVELENGTHS: [798], + Tags.DO_FILE_COMPRESSION: True +} + +settings = sp.Settings(general_settings) + +settings.set_volume_creation_settings({ + Tags.SIMULATE_DEFORMED_LAYERS: True, + Tags.STRUCTURES: create_example_tissue(settings) +}) +settings.set_optical_settings({ + Tags.OPTICAL_MODEL_NUMBER_PHOTONS: 5e7, + Tags.OPTICAL_MODEL_BINARY_PATH: path_manager.get_mcx_binary_path(), + Tags.COMPUTE_DIFFUSE_REFLECTANCE: SAVE_REFLECTANCE, + Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT: SAVE_PHOTON_DIRECTION +}) +settings["noise_model_1"] = { + Tags.NOISE_MEAN: 1.0, + Tags.NOISE_STD: 0.1, + Tags.NOISE_MODE: Tags.NOISE_MODE_MULTIPLICATIVE, + Tags.DATA_FIELD: Tags.DATA_FIELD_INITIAL_PRESSURE, + Tags.NOISE_NON_NEGATIVITY_CONSTRAINT: True +} + +if not SAVE_REFLECTANCE and not SAVE_PHOTON_DIRECTION: + pipeline = [ + sp.ModelBasedVolumeCreationAdapter(settings), + sp.MCXAdapter(settings), + sp.GaussianNoise(settings, "noise_model_1") + ] +else: + pipeline = [ + sp.ModelBasedVolumeCreationAdapter(settings), + sp.MCXAdapterReflectance(settings), + ] + + +class ExampleDeviceSlitIlluminationLinearDetector(sp.PhotoacousticDevice): + """ + This class represents a digital twin of a PA device with a slit as illumination next to a linear detection geometry. + + """ + + def __init__(self): + super().__init__(device_position_mm=np.asarray([VOLUME_TRANSDUCER_DIM_IN_MM/2, + VOLUME_PLANAR_DIM_IN_MM/2, 0])) + self.set_detection_geometry(sp.LinearArrayDetectionGeometry()) + self.add_illumination_geometry(sp.SlitIlluminationGeometry(slit_vector_mm=[20, 0, 0], + direction_vector_mm=[0, 0, 1])) + + +device = ExampleDeviceSlitIlluminationLinearDetector() + +sp.simulate(pipeline, settings, device) + +if Tags.WAVELENGTH in settings: + WAVELENGTH = settings[Tags.WAVELENGTH] +else: + WAVELENGTH = 700 + +if VISUALIZE: + sp.visualise_data(path_to_hdf5_file=path_manager.get_hdf5_file_save_path() + "/" + VOLUME_NAME + ".hdf5", + wavelength=WAVELENGTH, + show_initial_pressure=True, + show_oxygenation=True, + show_absorption=True, + show_diffuse_reflectance=SAVE_REFLECTANCE, + log_scale=True) diff --git a/simpa_tests/automatic_tests/test_additional_flags.py b/simpa_tests/automatic_tests/test_additional_flags.py index ca1c4731..ce0e9b33 100644 --- a/simpa_tests/automatic_tests/test_additional_flags.py +++ b/simpa_tests/automatic_tests/test_additional_flags.py @@ -1,24 +1,30 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT + import unittest import numpy as np from simpa import MCXAdapterReflectance, MCXAdapter, KWaveAdapter, TimeReversalAdapter, Tags, Settings from simpa.utils.matlab import generate_matlab_cmd + class TestAdditionalFlags(unittest.TestCase): def setUp(self) -> None: self.additional_flags = ('-l', '-a') self.settings = Settings() - + def test_get_cmd_mcx_reflectance_adapter(self): self.settings.set_optical_settings({ Tags.OPTICAL_MODEL_BINARY_PATH: '.', Tags.ADDITIONAL_FLAGS: self.additional_flags }) - mcx_reflectance_adapter = MCXAdapterReflectance(global_settings=self.settings) + mcx_reflectance_adapter = MCXAdapterReflectance(global_settings=self.settings) cmd = mcx_reflectance_adapter.get_command() for flag in self.additional_flags: - self.assertIn(flag, cmd, f"{flag} was not in command returned by mcx reflectance adapter but was defined as additional flag") - + self.assertIn( + flag, cmd, f"{flag} was not in command returned by mcx reflectance adapter but was defined as additional flag") + def test_get_cmd_mcx_adapter(self): self.settings.set_optical_settings({ Tags.OPTICAL_MODEL_BINARY_PATH: '.', @@ -27,8 +33,9 @@ def test_get_cmd_mcx_adapter(self): mcx_adapter = MCXAdapter(global_settings=self.settings) cmd = mcx_adapter.get_command() for flag in self.additional_flags: - self.assertIn(flag, cmd, f"{flag} was not in command returned by mcx adapter but was defined as additional flag") - + self.assertIn( + flag, cmd, f"{flag} was not in command returned by mcx adapter but was defined as additional flag") + def test_get_cmd_kwave_adapter(self): self.settings.set_acoustic_settings({ Tags.ADDITIONAL_FLAGS: self.additional_flags @@ -36,18 +43,20 @@ def test_get_cmd_kwave_adapter(self): kwave_adapter = KWaveAdapter(global_settings=self.settings) cmd = generate_matlab_cmd("./matlab.exe", "simulate_2D.m", "my_hdf5.mat", kwave_adapter.get_additional_flags()) for flag in self.additional_flags: - self.assertIn(flag, cmd, f"{flag} was not in command returned by kwave adapter but was defined as additional flag") + self.assertIn( + flag, cmd, f"{flag} was not in command returned by kwave adapter but was defined as additional flag") def test_get_cmd_time_reversal_adapter(self): self.settings.set_reconstruction_settings({ Tags.ADDITIONAL_FLAGS: self.additional_flags }) time_reversal_adapter = TimeReversalAdapter(global_settings=self.settings) - cmd = generate_matlab_cmd("./matlab.exe", "time_reversal_2D.m", "my_hdf5.mat", time_reversal_adapter.get_additional_flags()) + cmd = generate_matlab_cmd("./matlab.exe", "time_reversal_2D.m", "my_hdf5.mat", + time_reversal_adapter.get_additional_flags()) for flag in self.additional_flags: - self.assertIn(flag, cmd, f"{flag} was not in command returned by time reversal adapter but was defined as additional flag") - - + self.assertIn( + flag, cmd, f"{flag} was not in command returned by time reversal adapter but was defined as additional flag") + if __name__ == '__main__': - unittest.main() \ No newline at end of file + unittest.main() diff --git a/simpa_tests/automatic_tests/test_heterogeneity_generator.py b/simpa_tests/automatic_tests/test_heterogeneity_generator.py new file mode 100644 index 00000000..e186b24a --- /dev/null +++ b/simpa_tests/automatic_tests/test_heterogeneity_generator.py @@ -0,0 +1,224 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT + +import unittest +import numpy as np +import simpa as sp +from simpa.utils import Tags +import torch + + +class TestHeterogeneityGenerator(unittest.TestCase): + + def setUp(self) -> None: + self.spacing = 1.0 + self.MIN = -4.0 + self.MAX = 8.0 + self.MEAN = -332.0 + self.STD = 78.0 + self.FULL_IMAGE = np.zeros((4, 8)) + self.FULL_IMAGE[:, 1:2] = 1 + self.PARTIAL_IMAGE = np.zeros((2, 2)) + self.PARTIAL_IMAGE[:, 1] = 1 + self.TEST_SETTINGS = sp.Settings({ + # These parameters set the general properties of the simulated volume + sp.Tags.SPACING_MM: self.spacing, + sp.Tags.DIM_VOLUME_Z_MM: 8, + sp.Tags.DIM_VOLUME_X_MM: 4, + sp.Tags.DIM_VOLUME_Y_MM: 7 + }) + dimx, dimy, dimz = self.TEST_SETTINGS.get_volume_dimensions_voxels() + self.HETEROGENEITY_GENERATORS = [ + sp.RandomHeterogeneity(dimx, dimy, dimz, spacing_mm=self.spacing), + sp.RandomHeterogeneity(dimx, dimy, dimz, spacing_mm=self.spacing, gaussian_blur_size_mm=3), + sp.BlobHeterogeneity(dimx, dimy, dimz, spacing_mm=self.spacing), + sp.ImageHeterogeneity(dimx, dimy, dimz, heterogeneity_image=self.FULL_IMAGE, spacing_mm=self.spacing, + image_pixel_spacing_mm=self.spacing), + sp.ImageHeterogeneity(dimx, dimy, dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_CONSTANT, spacing_mm=self.spacing, constant=0.5), + sp.ImageHeterogeneity(dimx, dimy, dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_STRETCH, spacing_mm=self.spacing), + sp.ImageHeterogeneity(dimx, dimy, dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_SYMMETRIC, spacing_mm=self.spacing), + sp.ImageHeterogeneity(dimx, dimy, dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_WRAP, spacing_mm=self.spacing), + sp.ImageHeterogeneity(dimx, dimy, dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_EDGE, spacing_mm=self.spacing), + ] + self.HETEROGENEITY_GENERATORS_MIN_MAX = [ + sp.RandomHeterogeneity(dimx, dimy, dimz, spacing_mm=self.spacing, target_min=self.MIN, target_max=self.MAX), + sp.RandomHeterogeneity(dimx, dimy, dimz, spacing_mm=self.spacing, target_min=self.MIN, target_max=self.MAX, + gaussian_blur_size_mm=3), + sp.BlobHeterogeneity(dimx, dimy, dimz, spacing_mm=self.spacing, target_min=self.MIN, target_max=self.MAX), + sp.ImageHeterogeneity(dimx, dimy, dimz, heterogeneity_image=self.FULL_IMAGE, spacing_mm=self.spacing, + target_min=self.MIN, target_max=self.MAX), + sp.ImageHeterogeneity(dimx, dimy, dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_CONSTANT, spacing_mm=self.spacing, constant=0.5, + target_min=self.MIN, target_max=self.MAX), + sp.ImageHeterogeneity(dimx, dimy, dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_STRETCH, spacing_mm=self.spacing, + target_min=self.MIN, target_max=self.MAX), + sp.ImageHeterogeneity(dimx, dimy, dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_SYMMETRIC, spacing_mm=self.spacing, + target_min=self.MIN, target_max=self.MAX), + sp.ImageHeterogeneity(dimx, dimy, dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_WRAP, spacing_mm=self.spacing, + target_min=self.MIN, target_max=self.MAX), + sp.ImageHeterogeneity(dimx, dimy, dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_EDGE, spacing_mm=self.spacing, + target_min=self.MIN, target_max=self.MAX), + ] + self.HETEROGENEITY_GENERATORS_MEAN_STD = [ + sp.RandomHeterogeneity(dimx, dimy, dimz, spacing_mm=self.spacing, + target_mean=self.MEAN, target_std=self.STD), + sp.RandomHeterogeneity(dimx, dimy, dimz, spacing_mm=self.spacing, target_mean=self.MEAN, target_std=self.STD, + gaussian_blur_size_mm=3), + sp.BlobHeterogeneity(dimx, dimy, dimz, spacing_mm=self.spacing, target_mean=self.MEAN, target_std=self.STD), + sp.ImageHeterogeneity(dimx, dimy, dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_CONSTANT, spacing_mm=self.spacing, constant=0.5, + target_mean=self.MEAN, target_std=self.STD), + sp.ImageHeterogeneity(dimx, dimy, dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_STRETCH, spacing_mm=self.spacing, + target_mean=self.MEAN, target_std=self.STD), + sp.ImageHeterogeneity(dimx, dimy, dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_SYMMETRIC, spacing_mm=self.spacing, + target_mean=self.MEAN, target_std=self.STD), + sp.ImageHeterogeneity(dimx, dimy, dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_WRAP, spacing_mm=self.spacing, + target_mean=self.MEAN, target_std=self.STD), + sp.ImageHeterogeneity(dimx, dimy, dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_EDGE, spacing_mm=self.spacing, + target_mean=self.MEAN, target_std=self.STD), + ] + + def tearDown(self) -> None: + pass + + def assert_dimension_size(self, heterogeneity_generator, dimx, dimy, dimz): + random_map = heterogeneity_generator.get_map() + map_shape = np.shape(random_map) + self.assertAlmostEqual(dimx, map_shape[0]) + self.assertAlmostEqual(dimy, map_shape[1]) + self.assertAlmostEqual(dimz, map_shape[2]) + + def test_dimension_sizes(self): + dimx, dimy, dimz = self.TEST_SETTINGS.get_volume_dimensions_voxels() + for generator in self.HETEROGENEITY_GENERATORS: + self.assert_dimension_size(generator, dimx, dimy, dimz) + + def assert_min_max(self, heterogeneity_generator): + random_map = heterogeneity_generator.get_map() + self.assertAlmostEqual(np.min(random_map), self.MIN, 5) + self.assertAlmostEqual(np.max(random_map), self.MAX, 5) + + def test_min_max_bounds(self): + for generator in self.HETEROGENEITY_GENERATORS_MIN_MAX: + self.assert_min_max(generator) + + def assert_mean_std(self, heterogeneity_generator): + random_map = heterogeneity_generator.get_map() + self.assertAlmostEqual(np.mean(random_map), self.MEAN, 5) + self.assertAlmostEqual(np.std(random_map), self.STD, 5) + + def test_mean_std_bounds(self): + for generator in self.HETEROGENEITY_GENERATORS_MEAN_STD: + self.assert_mean_std(generator) + + +class TestImageScaling(unittest.TestCase): + """ + A set of tests for the ImageHeterogeneity class, designed to see if the scaling works. + """ + + def setUp(self): + self.spacing = 1.0 + self.MIN = -4.0 + self.MAX = 8.0 + self.PARTIAL_IMAGE = np.zeros((2, 2)) + self.PARTIAL_IMAGE[:, 1] = 1 + self.TOO_BIG_IMAGE = np.zeros((8, 8)) + self.TOO_BIG_IMAGE[:, :: 2] = 1 + self.TEST_SETTINGS = sp.Settings({ + # These parameters set the general properties of the simulated volume + sp.Tags.SPACING_MM: self.spacing, + sp.Tags.DIM_VOLUME_Z_MM: 8, + sp.Tags.DIM_VOLUME_X_MM: 4, + sp.Tags.DIM_VOLUME_Y_MM: 7 + }) + self.dimx, self.dimy, self.dimz = self.TEST_SETTINGS.get_volume_dimensions_voxels() + + def test_stretch(self): + """ + Test to see if the image can be stretched to fill th area, and then the volume. After stretched to fill the + area we should see the furthest two columns keep their values + :return: Assertion for if the image has been stretched + """ + stretched_image = sp.ImageHeterogeneity(self.dimx, self.dimy, self.dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_STRETCH, spacing_mm=self.spacing).get_map() + end_equals_1 = np.all(stretched_image[:, :, 6:] == 1) + beginning_equals_0 = np.all(stretched_image[:, :, :1] == 0) + assert end_equals_1 and beginning_equals_0 + + def test_wrap(self): + """ + Test to see if the image can be replicated to fill th area, and then the volume. Even and odd columns will keep + their values + :return: Assertion for if the image has been wrapped to fill the volume + """ + wrapped_image = sp.ImageHeterogeneity(self.dimx, self.dimy, self.dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_WRAP, spacing_mm=self.spacing).get_map() + even_equal_1 = np.all(wrapped_image[:, :, 1::2] == 1) + odd_equal_zero = np.all(wrapped_image[:, :, ::2] == 0) + assert even_equal_1 and odd_equal_zero + + def test_edge(self): + """ + Test to see if the image can fill the area by extending the edges, and then the volume. Should observe a line + of zeros and the rest ones. + :return: Assertion for if the image edges have filled the volume + """ + edge_image = sp.ImageHeterogeneity(self.dimx, self.dimy, self.dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_EDGE, spacing_mm=self.spacing).get_map() + initially_zero = np.all(edge_image[:, :, 0] == 0) + rest_ones = np.all(edge_image[:, :, 1:] == 1) + assert initially_zero and rest_ones + + def test_constant(self): + """ + Test to see if the image can fill the area with a constant, and then the volume + :return: Assertion for if the image has been filled by a constant + """ + constant_image = sp.ImageHeterogeneity(self.dimx, self.dimy, self.dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_CONSTANT, spacing_mm=self.spacing, + constant=0.5).get_map() + initial_area_same = np.all(constant_image[1:3, :, 0] == 0) and np.all(constant_image[1:3, :, 1] == 1) + rest_constant = np.all(constant_image[:, :, 2:] == 0.5) and np.all(constant_image[0, :, :] == 0.5) and \ + np.all(constant_image[3, :, :] == 0.5) + assert initial_area_same and rest_constant + + def test_symmetric(self): + """ + Test to see if the image can fill the area by symmetric reflections, and then the volume. See stripes following + two pixel 1 to 0 patterns + :return: Assertion for if the image has been reflected to fill the volume + """ + symmetric_image = sp.ImageHeterogeneity(self.dimx, self.dimy, self.dimz, heterogeneity_image=self.PARTIAL_IMAGE, + scaling_type=Tags.IMAGE_SCALING_SYMMETRIC, + spacing_mm=self.spacing).get_map() + ones_stripes_working = np.all(symmetric_image[:, :, 1:3] == 1) and np.all(symmetric_image[:, :, 5:7] == 1) + zeros_stripes_working = np.all(symmetric_image[:, :, 0] == 0) and np.all(symmetric_image[:, :, 3:5] == 0) and \ + np.all(symmetric_image[:, :, 7:] == 0) + assert ones_stripes_working and zeros_stripes_working + + def test_crop(self): + """ + Test to see if the image will crop to the desired area, thus leaving the same pattern but in a smaller shape + :return: Assertion for if the image has been cropped to the desired area + """ + crop_image = sp.ImageHeterogeneity(self.dimx, self.dimy, self.dimz, heterogeneity_image=self.TOO_BIG_IMAGE, + spacing_mm=self.spacing, image_pixel_spacing_mm=self.spacing).get_map() + odd_columns_equal_1 = np.all(crop_image[:, :, ::2] == 1) + even_columns_equal_0 = np.all(crop_image[:, :, 1::2] == 0) + size_is_right = np.all(crop_image[:, 1, :].shape == (self.dimx, self.dimz)) + assert odd_columns_equal_1 and even_columns_equal_0 and size_is_right diff --git a/simpa_tests/automatic_tests/tissue_library/test_core_assumptions.py b/simpa_tests/automatic_tests/tissue_library/test_core_assumptions.py index 1fb413ce..049d6540 100644 --- a/simpa_tests/automatic_tests/tissue_library/test_core_assumptions.py +++ b/simpa_tests/automatic_tests/tissue_library/test_core_assumptions.py @@ -7,21 +7,29 @@ import numpy as np -from simpa.utils import TISSUE_LIBRARY +from simpa.utils import Tags, Settings, TISSUE_LIBRARY from simpa.utils.calculate import calculate_bvf, calculate_oxygenation from simpa.utils.libraries.molecule_library import MolecularComposition from simpa.utils.libraries.tissue_library import TissueLibrary +TEST_SETTINGS = Settings({ + # These parameters set the general properties of the simulated volume + Tags.SPACING_MM: 1, + Tags.DIM_VOLUME_Z_MM: 4, + Tags.DIM_VOLUME_X_MM: 2, + Tags.DIM_VOLUME_Y_MM: 7 +}) + class TestCoreAssumptions(unittest.TestCase): def test_volume_fractions_sum_to_less_or_equal_one(self): for (method_name, method) in self.get_all_tissue_library_methods(): - total_volume_fraction = 0 - for molecule in method(TISSUE_LIBRARY): - total_volume_fraction += molecule.volume_fraction - self.assertAlmostEqual(total_volume_fraction, 1.0, 3, - f"Volume fraction not 1.0 +/- 0.001 for {method_name}") + molecular_composition = method(TISSUE_LIBRARY) + tissue_composition = molecular_composition.get_properties_for_wavelength(TEST_SETTINGS, 800) + total_volume_fraction = tissue_composition.volume_fraction + self.assertTrue((np.abs(total_volume_fraction-1.0) < 1e-3).all(), + f"Volume fraction not 1.0 +/- 0.001 for {method_name}") def test_bvf_and_oxygenation_consistency(self): # blood_volume_fraction (bvf) and oxygenation of tissue classes defined diff --git a/simpa_tests/manual_tests/digital_device_twins/VisualiseDevices.py b/simpa_tests/manual_tests/digital_device_twins/VisualiseDevices.py index a88dcebb..02b4ebe9 100644 --- a/simpa_tests/manual_tests/digital_device_twins/VisualiseDevices.py +++ b/simpa_tests/manual_tests/digital_device_twins/VisualiseDevices.py @@ -35,6 +35,7 @@ def visualise_result(self, show_figure_on_screen=True, save_path=None): sp.visualise_device(sp.RSOMExplorerP50(device_position_mm=np.asarray([50, 10, 0])), figure_save_path[2]) + if __name__ == "__main__": test = DeviceVisualisationTest() test.run_test(show_figure_on_screen=False) diff --git a/simpa_tests/manual_tests/executables/MATLABAdditionalFlags.py b/simpa_tests/manual_tests/executables/MATLABAdditionalFlags.py index 4e63a889..1f114949 100644 --- a/simpa_tests/manual_tests/executables/MATLABAdditionalFlags.py +++ b/simpa_tests/manual_tests/executables/MATLABAdditionalFlags.py @@ -34,7 +34,7 @@ def setup(self): """ path_manager = PathManager() - + self.settings = Settings({ Tags.WAVELENGTHS: [800], Tags.WAVELENGTH: 800, @@ -79,12 +79,12 @@ def setup(self): 0])) self.device.add_illumination_geometry(PencilBeamIlluminationGeometry()) self.device.set_detection_geometry(LinearArrayDetectionGeometry(device_position_mm=self.device.device_position_mm, - pitch_mm=0.25, - number_detector_elements=100, - field_of_view_extent_mm=np.asarray([-15, 15, 0, 0, 0, 20]))) + pitch_mm=0.25, + number_detector_elements=100, + field_of_view_extent_mm=np.asarray([-15, 15, 0, 0, 0, 20]))) output_name = f'{os.path.join(self.settings[Tags.SIMULATION_PATH], self.settings[Tags.VOLUME_NAME])}' - self.output_file_name = f'{output_name}.log' + self.output_file_name = f'{output_name}.log' def run_simulation(self): # run pipeline including volume creation and optical mcx simulation and acoustic matlab kwave simulation @@ -100,11 +100,11 @@ def test_execution_of_additional_flag(self): :raises FileNotFoundError: if log file does not exist at expected location """ - - # perform cleaning before test + + # perform cleaning before test if os.path.exists(self.output_file_name): os.remove(self.output_file_name) - + # run simulation self.settings.get_acoustic_settings()[Tags.ADDITIONAL_FLAGS] = ['-logfile', self.output_file_name] self.run_simulation() @@ -112,24 +112,26 @@ def test_execution_of_additional_flag(self): # checking if file exists afterwards if not os.path.exists(self.output_file_name): raise FileNotFoundError(f"Log file wasn't created at expected path {self.output_file_name}") - + def test_if_last_flag_is_used(self): """Tests if log file is created with correct last given name by setting multiple additional parameters :raises FileNotFoundError: if correct log file does not exist at expected location """ - + # perform cleaning before test if os.path.exists(self.output_file_name): os.remove(self.output_file_name) - + # run simulation - self.settings.get_acoustic_settings()[Tags.ADDITIONAL_FLAGS] = ['-logfile', 'temp_name', '-logfile', self.output_file_name] + self.settings.get_acoustic_settings()[Tags.ADDITIONAL_FLAGS] = [ + '-logfile', 'temp_name', '-logfile', self.output_file_name] self.run_simulation() # checking if file exists afterwards if not os.path.exists(self.output_file_name): - raise FileNotFoundError(f"Log file wasn't created with correct last given name at expected path {self.output_file_name}") + raise FileNotFoundError( + f"Log file wasn't created with correct last given name at expected path {self.output_file_name}") def perform_test(self): """ @@ -139,8 +141,9 @@ def perform_test(self): self.test_if_last_flag_is_used() def visualise_result(self, show_figure_on_screen=True, save_path=None): - pass # no figures are created that could be visualized - + pass # no figures are created that could be visualized + + if __name__ == '__main__': test = MATLABAdditionalFlags() test.run_test(show_figure_on_screen=False) diff --git a/simpa_tests/manual_tests/executables/MCXAdditionalFlags.py b/simpa_tests/manual_tests/executables/MCXAdditionalFlags.py index 4bf7da1d..b381045e 100644 --- a/simpa_tests/manual_tests/executables/MCXAdditionalFlags.py +++ b/simpa_tests/manual_tests/executables/MCXAdditionalFlags.py @@ -31,7 +31,7 @@ def setup(self): """ path_manager = PathManager() - + self.settings = Settings({ Tags.WAVELENGTHS: [800], Tags.WAVELENGTH: 800, @@ -63,7 +63,7 @@ def setup(self): self.device.add_illumination_geometry(PencilBeamIlluminationGeometry()) self.output_name = f'{os.path.join(self.settings[Tags.SIMULATION_PATH], self.settings[Tags.VOLUME_NAME])}_output' - self.output_file_name = f'{self.output_name}.log' + self.output_file_name = f'{self.output_name}.log' def run_simulation(self): # run pipeline including volume creation and optical mcx simulation @@ -78,11 +78,11 @@ def test_execution_of_additional_flag(self): :raises FileNotFoundError: if log file does not exist at expected location """ - - # perform cleaning before test + + # perform cleaning before test if os.path.exists(self. output_file_name): os.remove(self.output_file_name) - + # run simulation self.settings.get_optical_settings()[Tags.ADDITIONAL_FLAGS] = ['-l', 1, '-s', self.output_name] self.run_simulation() @@ -90,26 +90,27 @@ def test_execution_of_additional_flag(self): # checking if file exists afterwards if not os.path.exists(self.output_file_name): raise FileNotFoundError(f"Log file wasn't created at expected path {self.output_file_name}") - + def test_if_last_flag_is_used(self): """Tests if log file is created with correct last given name by setting multiple additional parameters :raises FileNotFoundError: if correct log file does not exist at expected location """ output_name = f'{os.path.join(self.settings[Tags.SIMULATION_PATH], self.settings[Tags.VOLUME_NAME])}_output' - output_file_name = f'{output_name}.log' + output_file_name = f'{output_name}.log' - # perform cleaning before test + # perform cleaning before test if os.path.exists(output_file_name): os.remove(output_file_name) - + # run simulation self.settings.get_optical_settings()[Tags.ADDITIONAL_FLAGS] = ['-l', 1, '-s', 'temp_name', '-s', output_name] self.run_simulation() # checking if file exists afterwards if not os.path.exists(output_file_name): - raise FileNotFoundError(f"Log file wasn't created with correct last given name at expected path {output_file_name}") + raise FileNotFoundError( + f"Log file wasn't created with correct last given name at expected path {output_file_name}") def perform_test(self): """ @@ -119,8 +120,9 @@ def perform_test(self): self.test_if_last_flag_is_used() def visualise_result(self, show_figure_on_screen=True, save_path=None): - pass # no figures are created that could be visualized - + pass # no figures are created that could be visualized + + if __name__ == '__main__': test = MCXAdditionalFlags() test.run_test(show_figure_on_screen=False) diff --git a/simpa_tests/manual_tests/optical_forward_models/AbsorptionAndScatteringWithinHomogenousMedium.py b/simpa_tests/manual_tests/optical_forward_models/AbsorptionAndScatteringWithinHomogenousMedium.py index 806863bd..b34937bb 100644 --- a/simpa_tests/manual_tests/optical_forward_models/AbsorptionAndScatteringWithinHomogenousMedium.py +++ b/simpa_tests/manual_tests/optical_forward_models/AbsorptionAndScatteringWithinHomogenousMedium.py @@ -38,7 +38,7 @@ os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE" -class TestAbsorptionAndScatteringWithInifinitesimalSlabExperiment(ManualIntegrationTestClass): +class TestAbsorptionAndScatteringWithinHomogeneousMedium(ManualIntegrationTestClass): def create_example_tissue(self, scattering_value=1e-30, absorption_value=1e-30, anisotropy_value=0.0): """ @@ -101,115 +101,115 @@ def test_low_scattering(self): Here, the slab is 10 mm long, mua and mus are both used with values of 0.05 mm^-1, so that mua+mus=0.1 mm^-1. We expect a decay ratio of e^1. """ - return self.test_simultion(absorption_value_1=0.01, - absorption_value_2=0.01, - scattering_value_1=1.0, - scattering_value_2=10.0, - anisotropy_value_1=0.0, - anisotropy_value_2=0.9, - title="Low Abs. Low Scat.") + return self.test_simulation(absorption_value_1=0.01, + absorption_value_2=0.01, + scattering_value_1=1.0, + scattering_value_2=10.0, + anisotropy_value_1=0.0, + anisotropy_value_2=0.9, + title="Low Abs. Low Scat.") def test_medium_scattering(self): """ Here, the slab is 10 mm long, mua and mus are both used with values of 0.05 mm^-1, so that mua+mus=0.1 mm^-1. We expect a decay ratio of e^1. """ - return self.test_simultion(absorption_value_1=0.01, - absorption_value_2=0.01, - scattering_value_1=10.0, - scattering_value_2=100.0, - anisotropy_value_1=0.0, - anisotropy_value_2=0.9, - title="Low Abs. Medium Scat.") + return self.test_simulation(absorption_value_1=0.01, + absorption_value_2=0.01, + scattering_value_1=10.0, + scattering_value_2=100.0, + anisotropy_value_1=0.0, + anisotropy_value_2=0.9, + title="Low Abs. Medium Scat.") def test_high_scattering_090(self): """ Here, the slab is 10 mm long, mua and mus are both used with values of 0.05 mm^-1, so that mua+mus=0.1 mm^-1. We expect a decay ratio of e^1. """ - return self.test_simultion(absorption_value_1=0.01, - absorption_value_2=0.01, - scattering_value_1=50.0, - scattering_value_2=500.0, - anisotropy_value_1=0.0, - anisotropy_value_2=0.9, - title="Anisotropy 0.9") + return self.test_simulation(absorption_value_1=0.01, + absorption_value_2=0.01, + scattering_value_1=50.0, + scattering_value_2=500.0, + anisotropy_value_1=0.0, + anisotropy_value_2=0.9, + title="Anisotropy 0.9") def simulate_perfect_result(self): """ Here, the slab is 10 mm long, mua and mus are both used with values of 0.05 mm^-1, so that mua+mus=0.1 mm^-1. We expect a decay ratio of e^1. """ - return self.test_simultion(absorption_value_1=0.01, - absorption_value_2=0.01, - scattering_value_1=50.0, - scattering_value_2=50.0, - anisotropy_value_1=0.0, - anisotropy_value_2=0.0, - title="Ideal Result") + return self.test_simulation(absorption_value_1=0.01, + absorption_value_2=0.01, + scattering_value_1=50.0, + scattering_value_2=50.0, + anisotropy_value_1=0.0, + anisotropy_value_2=0.0, + title="Ideal Result") def test_high_scattering_075(self): """ Here, the slab is 10 mm long, mua and mus are both used with values of 0.05 mm^-1, so that mua+mus=0.1 mm^-1. We expect a decay ratio of e^1. """ - return self.test_simultion(absorption_value_1=0.01, - absorption_value_2=0.01, - scattering_value_1=50.0, - scattering_value_2=200.0, - anisotropy_value_1=0.0, - anisotropy_value_2=0.75, - title="Anisotropy 0.75") + return self.test_simulation(absorption_value_1=0.01, + absorption_value_2=0.01, + scattering_value_1=50.0, + scattering_value_2=200.0, + anisotropy_value_1=0.0, + anisotropy_value_2=0.75, + title="Anisotropy 0.75") def test_high_scattering_025(self): """ Here, the slab is 10 mm long, mua and mus are both used with values of 0.05 mm^-1, so that mua+mus=0.1 mm^-1. We expect a decay ratio of e^1. """ - return self.test_simultion(absorption_value_1=0.01, - absorption_value_2=0.01, - scattering_value_1=50.0, - scattering_value_2=66.666666666666667, - anisotropy_value_1=0.0, - anisotropy_value_2=0.25, - title="Anisotropy 0.25") + return self.test_simulation(absorption_value_1=0.01, + absorption_value_2=0.01, + scattering_value_1=50.0, + scattering_value_2=66.666666666666667, + anisotropy_value_1=0.0, + anisotropy_value_2=0.25, + title="Anisotropy 0.25") def test_ignore_mcx_anisotropy_025(self): """ Here, the slab is 10 mm long, mua and mus are both used with values of 0.05 mm^-1, so that mua+mus=0.1 mm^-1. We expect a decay ratio of e^1. """ - return self.test_simultion(absorption_value_1=0.01, - absorption_value_2=0.01, - scattering_value_1=50.0, - scattering_value_2=66.666666666666667, - anisotropy_value_1=0.0, - anisotropy_value_2=0.25, - title="Ignore MCX Anisotropy 0.25", - use_mcx_anisotropy=False) + return self.test_simulation(absorption_value_1=0.01, + absorption_value_2=0.01, + scattering_value_1=50.0, + scattering_value_2=66.666666666666667, + anisotropy_value_1=0.0, + anisotropy_value_2=0.25, + title="Ignore MCX Anisotropy 0.25", + use_mcx_anisotropy=False) def test_ignore_mcx_anisotropy_075(self): """ Here, the slab is 10 mm long, mua and mus are both used with values of 0.05 mm^-1, so that mua+mus=0.1 mm^-1. We expect a decay ratio of e^1. """ - return self.test_simultion(absorption_value_1=0.01, - absorption_value_2=0.01, - scattering_value_1=50.0, - scattering_value_2=200.0, - anisotropy_value_1=0.0, - anisotropy_value_2=0.75, - title="Ignore MCX Anisotropy 0.75", - use_mcx_anisotropy=False) - - def test_simultion(self, scattering_value_1=1e-30, - absorption_value_1=1e-30, - anisotropy_value_1=1.0, - scattering_value_2=1e-30, - absorption_value_2=1e-30, - anisotropy_value_2=1.0, - title="Medium Abs. High Scat.", - use_mcx_anisotropy=True): + return self.test_simulation(absorption_value_1=0.01, + absorption_value_2=0.01, + scattering_value_1=50.0, + scattering_value_2=200.0, + anisotropy_value_1=0.0, + anisotropy_value_2=0.75, + title="Ignore MCX Anisotropy 0.75", + use_mcx_anisotropy=False) + + def test_simulation(self, scattering_value_1=1e-30, + absorption_value_1=1e-30, + anisotropy_value_1=1.0, + scattering_value_2=1e-30, + absorption_value_2=1e-30, + anisotropy_value_2=1.0, + title="Medium Abs. High Scat.", + use_mcx_anisotropy=True): # RUN SIMULATION 1 @@ -220,7 +220,10 @@ def test_simultion(self, scattering_value_1=1e-30, anisotropy_value=anisotropy_value_1) }) - self.settings.get_optical_settings()[Tags.MCX_ASSUMED_ANISOTROPY] = anisotropy_value_1 + if use_mcx_anisotropy: + self.settings.get_optical_settings()[Tags.MCX_ASSUMED_ANISOTROPY] = anisotropy_value_1 + else: + self.settings.get_optical_settings()[Tags.MCX_ASSUMED_ANISOTROPY] = anisotropy_value_2 pipeline = [ ModelBasedVolumeCreationAdapter(self.settings), @@ -241,10 +244,7 @@ def test_simultion(self, scattering_value_1=1e-30, anisotropy_value=anisotropy_value_2) }) - if use_mcx_anisotropy: - self.settings.get_optical_settings()[Tags.MCX_ASSUMED_ANISOTROPY] = anisotropy_value_2 - else: - self.settings.get_optical_settings()[Tags.MCX_ASSUMED_ANISOTROPY] = 0.9 + self.settings.get_optical_settings()[Tags.MCX_ASSUMED_ANISOTROPY] = anisotropy_value_2 pipeline = [ ModelBasedVolumeCreationAdapter(self.settings), @@ -309,5 +309,5 @@ def perform_test(self): if __name__ == '__main__': - test = TestAbsorptionAndScatteringWithInifinitesimalSlabExperiment() + test = TestAbsorptionAndScatteringWithinHomogeneousMedium() test.run_test(show_figure_on_screen=False) diff --git a/simpa_tests/test_utils/tissue_composition_tests.py b/simpa_tests/test_utils/tissue_composition_tests.py index bbf558c2..7adc8ca0 100644 --- a/simpa_tests/test_utils/tissue_composition_tests.py +++ b/simpa_tests/test_utils/tissue_composition_tests.py @@ -2,7 +2,7 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT -from simpa.utils import Tags, SegmentationClasses, calculate_gruneisen_parameter_from_temperature +from simpa.utils import Tags, Settings, SegmentationClasses, calculate_gruneisen_parameter_from_temperature from simpa.utils.libraries.molecule_library import MolecularComposition from simpa.utils.tissue_properties import TissueProperties from simpa.utils.constants import property_tags @@ -11,6 +11,14 @@ import matplotlib.patches as patches import matplotlib.pyplot as plt +TEST_SETTINGS = Settings({ + # These parameters set the general properties of the simulated volume + Tags.SPACING_MM: 1, + Tags.DIM_VOLUME_Z_MM: 1, + Tags.DIM_VOLUME_X_MM: 1, + Tags.DIM_VOLUME_Y_MM: 1 +}) + def validate_expected_values_dictionary(expected_values: dict): @@ -41,8 +49,9 @@ def compare_molecular_composition_against_expected_values(molecular_composition: num_subplots = len(property_tags) for wavelength in expected_values.keys(): - molecular_composition.update_internal_properties() - composition_properties = molecular_composition.get_properties_for_wavelength(wavelength=wavelength) + molecular_composition.update_internal_properties(TEST_SETTINGS) + composition_properties = molecular_composition.get_properties_for_wavelength( + TEST_SETTINGS, wavelength=wavelength) expected_properties = expected_values[wavelength] if visualise_values: @@ -92,7 +101,7 @@ def get_epidermis_reference_dictionary(): """ reference_dict = dict() - values450nm = TissueProperties() + values450nm = TissueProperties(TEST_SETTINGS) values450nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 13.5 values450nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 121.6 values450nm[Tags.DATA_FIELD_ANISOTROPY] = 0.728 @@ -104,7 +113,7 @@ def get_epidermis_reference_dictionary(): values450nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1624.0 values450nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.35 - values500nm = TissueProperties() + values500nm = TissueProperties(TEST_SETTINGS) values500nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 9.77 values500nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 93.01 values500nm[Tags.DATA_FIELD_ANISOTROPY] = 0.745 @@ -116,7 +125,7 @@ def get_epidermis_reference_dictionary(): values500nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1624.0 values500nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.35 - values550nm = TissueProperties() + values550nm = TissueProperties(TEST_SETTINGS) values550nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 6.85 values550nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 74.7 values550nm[Tags.DATA_FIELD_ANISOTROPY] = 0.759 @@ -128,7 +137,7 @@ def get_epidermis_reference_dictionary(): values550nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1624.0 values550nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.35 - values600nm = TissueProperties() + values600nm = TissueProperties(TEST_SETTINGS) values600nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 5.22 values600nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 63.76 values600nm[Tags.DATA_FIELD_ANISOTROPY] = 0.774 @@ -140,7 +149,7 @@ def get_epidermis_reference_dictionary(): values600nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1624.0 values600nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.35 - values650nm = TissueProperties() + values650nm = TissueProperties(TEST_SETTINGS) values650nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 3.68 values650nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 55.48 values650nm[Tags.DATA_FIELD_ANISOTROPY] = 0.7887 @@ -152,7 +161,7 @@ def get_epidermis_reference_dictionary(): values650nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1624.0 values650nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.35 - values700nm = TissueProperties() + values700nm = TissueProperties(TEST_SETTINGS) values700nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 3.07 values700nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 54.66 values700nm[Tags.DATA_FIELD_ANISOTROPY] = 0.804 @@ -195,7 +204,7 @@ def get_dermis_reference_dictionary(): """ reference_dict = dict() - values450nm = TissueProperties() + values450nm = TissueProperties(TEST_SETTINGS) values450nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 2.105749981 values450nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 244.6 values450nm[Tags.DATA_FIELD_ANISOTROPY] = 0.715 @@ -207,7 +216,7 @@ def get_dermis_reference_dictionary(): values450nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1624 values450nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.35 - values500nm = TissueProperties() + values500nm = TissueProperties(TEST_SETTINGS) values500nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 0.924812913 values500nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 175.0 values500nm[Tags.DATA_FIELD_ANISOTROPY] = 0.715 @@ -219,7 +228,7 @@ def get_dermis_reference_dictionary(): values500nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1624 values500nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.35 - values550nm = TissueProperties() + values550nm = TissueProperties(TEST_SETTINGS) values550nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 0.974386604 values550nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 131.1 values550nm[Tags.DATA_FIELD_ANISOTROPY] = 0.715 @@ -231,7 +240,7 @@ def get_dermis_reference_dictionary(): values550nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1624 values550nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.35 - values600nm = TissueProperties() + values600nm = TissueProperties(TEST_SETTINGS) values600nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 0.440476363 values600nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 101.9 values600nm[Tags.DATA_FIELD_ANISOTROPY] = 0.715 @@ -243,7 +252,7 @@ def get_dermis_reference_dictionary(): values600nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1624 values600nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.35 - values650nm = TissueProperties() + values650nm = TissueProperties(TEST_SETTINGS) values650nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 0.313052704 values650nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 81.7 values650nm[Tags.DATA_FIELD_ANISOTROPY] = 0.715 @@ -255,7 +264,7 @@ def get_dermis_reference_dictionary(): values650nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1624 values650nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.35 - values700nm = TissueProperties() + values700nm = TissueProperties(TEST_SETTINGS) values700nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 0.277003236 values700nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 67.1 values700nm[Tags.DATA_FIELD_ANISOTROPY] = 0.715 @@ -267,7 +276,7 @@ def get_dermis_reference_dictionary(): values700nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1624 values700nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.35 - values750nm = TissueProperties() + values750nm = TissueProperties(TEST_SETTINGS) values750nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 0.264286111 values750nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 56.3 values750nm[Tags.DATA_FIELD_ANISOTROPY] = 0.715 @@ -279,7 +288,7 @@ def get_dermis_reference_dictionary(): values750nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1624 values750nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.35 - values800nm = TissueProperties() + values800nm = TissueProperties(TEST_SETTINGS) values800nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 0.256933531 values800nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 48.1 values800nm[Tags.DATA_FIELD_ANISOTROPY] = 0.715 @@ -291,7 +300,7 @@ def get_dermis_reference_dictionary(): values800nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1624 values800nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.35 - values850nm = TissueProperties() + values850nm = TissueProperties(TEST_SETTINGS) values850nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 0.255224508 values850nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 41.8 values850nm[Tags.DATA_FIELD_ANISOTROPY] = 0.715 @@ -303,7 +312,7 @@ def get_dermis_reference_dictionary(): values850nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1624 values850nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.35 - values900nm = TissueProperties() + values900nm = TissueProperties(TEST_SETTINGS) values900nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 0.254198591 values900nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 36.7 values900nm[Tags.DATA_FIELD_ANISOTROPY] = 0.715 @@ -315,7 +324,7 @@ def get_dermis_reference_dictionary(): values900nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1624 values900nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.35 - values950nm = TissueProperties() + values950nm = TissueProperties(TEST_SETTINGS) values950nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 0.254522563 values950nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 32.6 values950nm[Tags.DATA_FIELD_ANISOTROPY] = 0.715 @@ -359,7 +368,7 @@ def get_muscle_reference_dictionary(): """ reference_dict = dict() - values650nm = TissueProperties() + values650nm = TissueProperties(TEST_SETTINGS) values650nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 1.04 values650nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 87.5 values650nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9 @@ -371,7 +380,7 @@ def get_muscle_reference_dictionary(): values650nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1588.4 values650nm[Tags.DATA_FIELD_ALPHA_COEFF] = 1.09 - values700nm = TissueProperties() + values700nm = TissueProperties(TEST_SETTINGS) values700nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 0.48 values700nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 81.8 values700nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9 @@ -383,7 +392,7 @@ def get_muscle_reference_dictionary(): values700nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1588.4 values700nm[Tags.DATA_FIELD_ALPHA_COEFF] = 1.09 - values750nm = TissueProperties() + values750nm = TissueProperties(TEST_SETTINGS) values750nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 0.41 values750nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 77.1 values750nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9 @@ -395,7 +404,7 @@ def get_muscle_reference_dictionary(): values750nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1588.4 values750nm[Tags.DATA_FIELD_ALPHA_COEFF] = 1.09 - values800nm = TissueProperties() + values800nm = TissueProperties(TEST_SETTINGS) values800nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 0.28 values800nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 70.4 values800nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9 @@ -407,7 +416,7 @@ def get_muscle_reference_dictionary(): values800nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1588.4 values800nm[Tags.DATA_FIELD_ALPHA_COEFF] = 1.09 - values850nm = TissueProperties() + values850nm = TissueProperties(TEST_SETTINGS) values850nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 0.3 values850nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 66.7 values850nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9 @@ -419,7 +428,7 @@ def get_muscle_reference_dictionary(): values850nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1588.4 values850nm[Tags.DATA_FIELD_ALPHA_COEFF] = 1.09 - values900nm = TissueProperties() + values900nm = TissueProperties(TEST_SETTINGS) values900nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 0.32 values900nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 62.1 values900nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9 @@ -431,7 +440,7 @@ def get_muscle_reference_dictionary(): values900nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1588.4 values900nm[Tags.DATA_FIELD_ALPHA_COEFF] = 1.09 - values950nm = TissueProperties() + values950nm = TissueProperties(TEST_SETTINGS) values950nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 0.46 values950nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 59.0 values950nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9 @@ -475,7 +484,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): """ reference_dict = dict() - values450nm = TissueProperties() + values450nm = TissueProperties(TEST_SETTINGS) values450nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 336 values450nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 772 values450nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9447 @@ -487,7 +496,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values450nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values450nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values500nm = TissueProperties() + values500nm = TissueProperties(TEST_SETTINGS) values500nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 112 values500nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 868.3 values500nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9761 @@ -499,7 +508,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values500nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values500nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values550nm = TissueProperties() + values550nm = TissueProperties(TEST_SETTINGS) values550nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 230 values550nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 714.9 values550nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9642 @@ -511,7 +520,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values550nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values550nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values600nm = TissueProperties() + values600nm = TissueProperties(TEST_SETTINGS) values600nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 17 values600nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 868.8 values600nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9794 @@ -523,7 +532,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values600nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values600nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values650nm = TissueProperties() + values650nm = TissueProperties(TEST_SETTINGS) values650nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 2 values650nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 880.1 values650nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9825 @@ -535,7 +544,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values650nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values650nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values700nm = TissueProperties() + values700nm = TissueProperties(TEST_SETTINGS) values700nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 1.6 values700nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 857.0 values700nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9836 @@ -547,7 +556,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values700nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values700nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values750nm = TissueProperties() + values750nm = TissueProperties(TEST_SETTINGS) values750nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 2.8 values750nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 802.2 values750nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9837 @@ -559,7 +568,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values750nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values750nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values800nm = TissueProperties() + values800nm = TissueProperties(TEST_SETTINGS) values800nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 4.4 values800nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 767.3 values800nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9833 @@ -571,7 +580,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values800nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values800nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values850nm = TissueProperties() + values850nm = TissueProperties(TEST_SETTINGS) values850nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 5.7 values850nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 742.0 values850nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9832 @@ -583,7 +592,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values850nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values850nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values900nm = TissueProperties() + values900nm = TissueProperties(TEST_SETTINGS) values900nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 6.4 values900nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 688.6 values900nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9824 @@ -595,7 +604,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values900nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values900nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values950nm = TissueProperties() + values950nm = TissueProperties(TEST_SETTINGS) values950nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 6.4 values950nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 652.1 values950nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9808 @@ -644,7 +653,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) """ reference_dict = dict() - values450nm = TissueProperties() + values450nm = TissueProperties(TEST_SETTINGS) values450nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 553 values450nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 772 values450nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9447 @@ -656,7 +665,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values450nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values450nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values500nm = TissueProperties() + values500nm = TissueProperties(TEST_SETTINGS) values500nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 112 values500nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 868.3 values500nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9761 @@ -668,7 +677,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values500nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values500nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values550nm = TissueProperties() + values550nm = TissueProperties(TEST_SETTINGS) values550nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 286 values550nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 714.9 values550nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9642 @@ -680,7 +689,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values550nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values550nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values600nm = TissueProperties() + values600nm = TissueProperties(TEST_SETTINGS) values600nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 79 values600nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 868.8 values600nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9794 @@ -692,7 +701,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values600nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values600nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values650nm = TissueProperties() + values650nm = TissueProperties(TEST_SETTINGS) values650nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 20.1 values650nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 880.1 values650nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9825 @@ -704,7 +713,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values650nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values650nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values700nm = TissueProperties() + values700nm = TissueProperties(TEST_SETTINGS) values700nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 9.6 values700nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 857.0 values700nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9836 @@ -716,7 +725,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values700nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values700nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values750nm = TissueProperties() + values750nm = TissueProperties(TEST_SETTINGS) values750nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 7.5 values750nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 802.2 values750nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9837 @@ -728,7 +737,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values750nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values750nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values800nm = TissueProperties() + values800nm = TissueProperties(TEST_SETTINGS) values800nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 4.1 values800nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 767.3 values800nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9833 @@ -740,7 +749,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values800nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values800nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values850nm = TissueProperties() + values850nm = TissueProperties(TEST_SETTINGS) values850nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 3.7 values850nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 742.0 values850nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9832 @@ -752,7 +761,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values850nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values850nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values900nm = TissueProperties() + values900nm = TissueProperties(TEST_SETTINGS) values900nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 4.1 values900nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 688.6 values900nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9824 @@ -764,7 +773,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values900nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values900nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 - values950nm = TissueProperties() + values950nm = TissueProperties(TEST_SETTINGS) values950nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 3.2 values950nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 652.1 values950nm[Tags.DATA_FIELD_ANISOTROPY] = 0.9808 @@ -799,7 +808,7 @@ def get_lymph_node_reference_dictionary(only_use_NIR_values=False): """ reference_dict = dict() - values450nm = TissueProperties() + values450nm = TissueProperties(TEST_SETTINGS) values450nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = None values450nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = None values450nm[Tags.DATA_FIELD_ANISOTROPY] = None @@ -811,7 +820,7 @@ def get_lymph_node_reference_dictionary(only_use_NIR_values=False): values450nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1586 values450nm[Tags.DATA_FIELD_ALPHA_COEFF] = 2.50 - values500nm = TissueProperties() + values500nm = TissueProperties(TEST_SETTINGS) values500nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = None values500nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = None values500nm[Tags.DATA_FIELD_ANISOTROPY] = None @@ -823,7 +832,7 @@ def get_lymph_node_reference_dictionary(only_use_NIR_values=False): values500nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1586 values500nm[Tags.DATA_FIELD_ALPHA_COEFF] = 2.50 - values550nm = TissueProperties() + values550nm = TissueProperties(TEST_SETTINGS) values550nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = None values550nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = None values550nm[Tags.DATA_FIELD_ANISOTROPY] = None @@ -835,7 +844,7 @@ def get_lymph_node_reference_dictionary(only_use_NIR_values=False): values550nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1586 values550nm[Tags.DATA_FIELD_ALPHA_COEFF] = 2.50 - values600nm = TissueProperties() + values600nm = TissueProperties(TEST_SETTINGS) values600nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = None values600nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = None values600nm[Tags.DATA_FIELD_ANISOTROPY] = None @@ -847,7 +856,7 @@ def get_lymph_node_reference_dictionary(only_use_NIR_values=False): values600nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1586 values600nm[Tags.DATA_FIELD_ALPHA_COEFF] = 2.50 - values650nm = TissueProperties() + values650nm = TissueProperties(TEST_SETTINGS) values650nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = None values650nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = None values650nm[Tags.DATA_FIELD_ANISOTROPY] = None @@ -859,7 +868,7 @@ def get_lymph_node_reference_dictionary(only_use_NIR_values=False): values650nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1586 values650nm[Tags.DATA_FIELD_ALPHA_COEFF] = 2.50 - values700nm = TissueProperties() + values700nm = TissueProperties(TEST_SETTINGS) values700nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = None values700nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = None values700nm[Tags.DATA_FIELD_ANISOTROPY] = None @@ -871,7 +880,7 @@ def get_lymph_node_reference_dictionary(only_use_NIR_values=False): values700nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1586 values700nm[Tags.DATA_FIELD_ALPHA_COEFF] = 2.50 - values750nm = TissueProperties() + values750nm = TissueProperties(TEST_SETTINGS) values750nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = None values750nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = None values750nm[Tags.DATA_FIELD_ANISOTROPY] = None @@ -883,7 +892,7 @@ def get_lymph_node_reference_dictionary(only_use_NIR_values=False): values750nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1586 values750nm[Tags.DATA_FIELD_ALPHA_COEFF] = 2.50 - values800nm = TissueProperties() + values800nm = TissueProperties(TEST_SETTINGS) values800nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = None values800nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = None values800nm[Tags.DATA_FIELD_ANISOTROPY] = None @@ -895,7 +904,7 @@ def get_lymph_node_reference_dictionary(only_use_NIR_values=False): values800nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1586 values800nm[Tags.DATA_FIELD_ALPHA_COEFF] = 2.50 - values850nm = TissueProperties() + values850nm = TissueProperties(TEST_SETTINGS) values850nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = None values850nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = None values850nm[Tags.DATA_FIELD_ANISOTROPY] = None @@ -907,7 +916,7 @@ def get_lymph_node_reference_dictionary(only_use_NIR_values=False): values850nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1586 values850nm[Tags.DATA_FIELD_ALPHA_COEFF] = 2.50 - values900nm = TissueProperties() + values900nm = TissueProperties(TEST_SETTINGS) values900nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = None values900nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = None values900nm[Tags.DATA_FIELD_ANISOTROPY] = None @@ -919,7 +928,7 @@ def get_lymph_node_reference_dictionary(only_use_NIR_values=False): values900nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1586 values900nm[Tags.DATA_FIELD_ALPHA_COEFF] = 2.50 - values950nm = TissueProperties() + values950nm = TissueProperties(TEST_SETTINGS) values950nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = None values950nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = None values950nm[Tags.DATA_FIELD_ANISOTROPY] = None diff --git a/simpa_tests/test_utils/tissue_models.py b/simpa_tests/test_utils/tissue_models.py index 77052e2e..4a234804 100644 --- a/simpa_tests/test_utils/tissue_models.py +++ b/simpa_tests/test_utils/tissue_models.py @@ -47,3 +47,46 @@ def create_simple_tissue_model(transducer_dim_in_mm: float, planar_dim_in_mm: fl adhere_to_deformation=False ) return tissue_dict + + +def create_heterogeneous_tissue_model(settings): + """ + This is a very simple example script of how to create a tissue definition. + It contains a muscular background, an epidermis layer on top of the muscles + and a blood vessel. + """ + v1_oxy = 1.0 + v2_oxy = 0.0 + bg_oxy = sp.RandomHeterogeneity().get_map() + background_dictionary = sp.Settings() + background_dictionary[Tags.MOLECULE_COMPOSITION] = (sp.MolecularCompositionGenerator() + .append(sp.MOLECULE_LIBRARY.oxyhemoglobin(bg_oxy)) + .append(sp.MOLECULE_LIBRARY.deoxyhemoglobin(1 - bg_oxy)) + .get_molecular_composition(sp.SegmentationClasses.BLOOD)) + + background_dictionary[Tags.STRUCTURE_TYPE] = Tags.BACKGROUND + + tissue_dict = sp.Settings() + tissue_dict[Tags.BACKGROUND] = background_dictionary + + tissue_dict["vessel_1"] = sp.define_circular_tubular_structure_settings( + tube_start_mm=[transducer_dim_in_mm / 2 - 10, 0, 5], + tube_end_mm=[transducer_dim_in_mm / 2 - 10, planar_dim_in_mm, 5], + molecular_composition=(sp.MolecularCompositionGenerator() + .append(sp.MOLECULE_LIBRARY.oxyhemoglobin(v1_oxy)) + .append(sp.MOLECULE_LIBRARY.deoxyhemoglobin(1 - v1_oxy)) + .get_molecular_composition(sp.SegmentationClasses.BLOOD)), + radius_mm=2, priority=3, consider_partial_volume=True, + adhere_to_deformation=False + ) + tissue_dict["vessel_2"] = sp.define_circular_tubular_structure_settings( + tube_start_mm=[transducer_dim_in_mm / 2, 0, 10], + tube_end_mm=[transducer_dim_in_mm / 2, planar_dim_in_mm, 10], + molecular_composition=(sp.MolecularCompositionGenerator() + .append(sp.MOLECULE_LIBRARY.oxyhemoglobin(v2_oxy)) + .append(sp.MOLECULE_LIBRARY.deoxyhemoglobin(1 - v2_oxy)) + .get_molecular_composition(sp.SegmentationClasses.BLOOD)), + radius_mm=3, priority=3, consider_partial_volume=True, + adhere_to_deformation=False + ) + return tissue_dict