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..669ccfc1 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). @@ -41,8 +46,8 @@ Please put this in the conversation of the pull request, and not add it to the f 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 +- **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..c653bfe1 100644 --- a/simpa/__init__.py +++ b/simpa/__init__.py @@ -6,37 +6,38 @@ from .log import Logger from importlib.metadata import version, PackageNotFoundError + try: __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 \ - SegmentationBasedVolumeCreationAdapter -from .core.simulation_modules.optical_simulation_module.optical_forward_model_mcx_adapter import \ + +from .core.simulation_modules.volume_creation_module.model_based_adapter import \ + ModelBasedAdapter +from .core.simulation_modules.volume_creation_module.segmentation_based_adapter import \ + SegmentationBasedAdapter +from .core.simulation_modules.optical_module.mcx_adapter import \ MCXAdapter -from .core.simulation_modules.optical_simulation_module.optical_forward_model_mcx_reflectance_adapter import \ - MCXAdapterReflectance -from .core.simulation_modules.acoustic_forward_module.acoustic_forward_module_k_wave_adapter import \ +from .core.simulation_modules.optical_module.mcx_reflectance_adapter import \ + MCXReflectanceAdapter +from .core.simulation_modules.acoustic_module.k_wave_adapter import \ KWaveAdapter -from .core.simulation_modules.reconstruction_module.reconstruction_module_delay_and_sum_adapter import \ +from .core.simulation_modules.reconstruction_module.delay_and_sum_adapter import \ DelayAndSumAdapter -from .core.simulation_modules.reconstruction_module.reconstruction_module_delay_multiply_and_sum_adapter import \ +from .core.simulation_modules.reconstruction_module.delay_multiply_and_sum_adapter import \ DelayMultiplyAndSumAdapter -from .core.simulation_modules.reconstruction_module.reconstruction_module_signed_delay_multiply_and_sum_adapter import \ +from .core.simulation_modules.reconstruction_module.signed_delay_multiply_and_sum_adapter import \ SignedDelayMultiplyAndSumAdapter -from .core.simulation_modules.reconstruction_module.reconstruction_module_time_reversal_adapter import \ +from .core.simulation_modules.reconstruction_module.time_reversal_adapter import \ TimeReversalAdapter -from .core.simulation_modules.reconstruction_module.reconstruction_module_delay_and_sum_adapter import \ +from .core.simulation_modules.reconstruction_module.delay_and_sum_adapter import \ reconstruct_delay_and_sum_pytorch -from .core.simulation_modules.reconstruction_module.reconstruction_module_delay_multiply_and_sum_adapter import \ +from .core.simulation_modules.reconstruction_module.delay_multiply_and_sum_adapter import \ reconstruct_delay_multiply_and_sum_pytorch -from .core.simulation_modules.reconstruction_module.reconstruction_module_signed_delay_multiply_and_sum_adapter import \ +from .core.simulation_modules.reconstruction_module.signed_delay_multiply_and_sum_adapter import \ reconstruct_signed_delay_multiply_and_sum_pytorch -from .core.simulation_modules.acoustic_forward_module.acoustic_forward_module_k_wave_adapter import \ +from .core.simulation_modules.acoustic_module.k_wave_adapter import \ perform_k_wave_acoustic_forward_simulation from simpa.core.processing_components.monospectral.noise import GaussianNoise diff --git a/simpa/core/__init__.py b/simpa/core/__init__.py index 9a49b941..b16f1abf 100644 --- a/simpa/core/__init__.py +++ b/simpa/core/__init__.py @@ -1,33 +1,4 @@ # SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT -from abc import abstractmethod - -from simpa.core.device_digital_twins import DigitalDeviceTwinBase -from simpa.log import Logger -from simpa.utils import Settings -from simpa.utils.processing_device import get_processing_device - - -class PipelineModule: - """ - Defines a pipeline module (either simulation or processing module) that implements a run method and can be called by running the pipeline's simulate method. - """ - - def __init__(self, global_settings: Settings): - """ - :param global_settings: The SIMPA settings dictionary - :type global_settings: Settings - """ - self.logger = Logger() - self.global_settings = global_settings - self.torch_device = get_processing_device(self.global_settings) - - @abstractmethod - def run(self, digital_device_twin: DigitalDeviceTwinBase): - """ - Executes the respective simulation module - - :param digital_device_twin: The digital twin that can be used by the digital device_twin. - """ - pass +from .pipeline_element_base import PipelineElementBase diff --git a/simpa/core/device_digital_twins/__init__.py b/simpa/core/device_digital_twins/__init__.py index 13bd5bcf..0c7de4d2 100644 --- a/simpa/core/device_digital_twins/__init__.py +++ b/simpa/core/device_digital_twins/__init__.py @@ -2,154 +2,22 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT -from abc import abstractmethod -from simpa.log import Logger -import numpy as np -import hashlib -import uuid -from simpa.utils.serializer import SerializableSIMPAClass -from simpa.utils.calculate import are_equal - - -class DigitalDeviceTwinBase(SerializableSIMPAClass): - """ - This class represents a device that can be used for illumination, detection or a combined photoacoustic device - which has representations of both. - """ - - def __init__(self, device_position_mm=None, field_of_view_extent_mm=None): - """ - :param device_position_mm: Each device has an internal position which serves as origin for internal \ - representations of e.g. detector element positions or illuminator positions. - :type device_position_mm: ndarray - :param field_of_view_extent_mm: Field of view which is defined as a numpy array of the shape \ - [xs, xe, ys, ye, zs, ze], where x, y, and z denote the coordinate axes and s and e denote the start and end \ - positions. - :type field_of_view_extent_mm: ndarray - """ - if device_position_mm is None: - self.device_position_mm = np.array([0, 0, 0]) - else: - self.device_position_mm = device_position_mm - - if field_of_view_extent_mm is None: - self.field_of_view_extent_mm = np.asarray([-10, 10, -10, 10, -10, 10]) - else: - self.field_of_view_extent_mm = field_of_view_extent_mm - - self.logger = Logger() - - def __eq__(self, other): - """ - Checks each key, value pair in the devices. - """ - if isinstance(other, DigitalDeviceTwinBase): - if self.__dict__.keys() != other.__dict__.keys(): - return False - for self_key, self_value in self.__dict__.items(): - other_value = other.__dict__[self_key] - if not are_equal(self_value, other_value): - return False - return True - return False - - @abstractmethod - def check_settings_prerequisites(self, global_settings) -> bool: - """ - It might be that certain device geometries need a certain dimensionality of the simulated PAI volume, or that - it requires the existence of certain Tags in the global global_settings. - To this end, a PAI device should use this method to inform the user about a mismatch of the desired device and - throw a ValueError if that is the case. - - :param global_settings: Settings for the entire simulation pipeline. - :type global_settings: Settings - - :raises ValueError: raises a value error if the prerequisites are not matched. - :returns: True if the prerequisites are met, False if they are not met, but no exception has been raised. - :rtype: bool - - """ - pass - - @abstractmethod - def update_settings_for_use_of_model_based_volume_creator(self, global_settings): - """ - This method can be overwritten by a PA device if the device poses special constraints to the - volume that should be considered by the model-based volume creator. - - :param global_settings: Settings for the entire simulation pipeline. - :type global_settings: Settings - """ - pass - - def get_field_of_view_mm(self) -> np.ndarray: - """ - Returns the absolute field of view in mm where the probe position is already - accounted for. - It is defined as a numpy array of the shape [xs, xe, ys, ye, zs, ze], - where x, y, and z denote the coordinate axes and s and e denote the start and end - positions. - - :return: Absolute field of view in mm where the probe position is already accounted for. - :rtype: ndarray - """ - position = self.device_position_mm - field_of_view_extent = self.field_of_view_extent_mm - - field_of_view = np.asarray([position[0] + field_of_view_extent[0], - position[0] + field_of_view_extent[1], - position[1] + field_of_view_extent[2], - position[1] + field_of_view_extent[3], - position[2] + field_of_view_extent[4], - position[2] + field_of_view_extent[5] - ]) - if min(field_of_view) < 0: - self.logger.warning(f"The field of view of the chosen device is not fully within the simulated volume, " - f"field of view is: {field_of_view}") - field_of_view[field_of_view < 0] = 0 - - return field_of_view - - def generate_uuid(self): - """ - Generates a universally unique identifier (uuid) for each device. - :return: - """ - class_dict = self.__dict__ - m = hashlib.md5() - m.update(str(class_dict).encode('utf-8')) - return str(uuid.UUID(m.hexdigest())) - - def serialize(self) -> dict: - serialized_device = self.__dict__ - return {"DigitalDeviceTwinBase": serialized_device} - - @staticmethod - def deserialize(dictionary_to_deserialize): - deserialized_device = DigitalDeviceTwinBase( - device_position_mm=dictionary_to_deserialize["device_position_mm"], - field_of_view_extent_mm=dictionary_to_deserialize["field_of_view_extent_mm"]) - return deserialized_device - - -""" -It is important to have these relative imports after the definition of the DigitalDeviceTwinBase class to avoid circular imports triggered by imported child classes -""" -from .pa_devices import PhotoacousticDevice # nopep8 -from simpa.core.device_digital_twins.detection_geometries import DetectionGeometryBase # nopep8 -from simpa.core.device_digital_twins.illumination_geometries import IlluminationGeometryBase # nopep8 -from .detection_geometries.curved_array import CurvedArrayDetectionGeometry # nopep8 -from .detection_geometries.linear_array import LinearArrayDetectionGeometry # nopep8 -from .detection_geometries.planar_array import PlanarArrayDetectionGeometry # nopep8 -from .illumination_geometries.slit_illumination import SlitIlluminationGeometry # nopep8 -from .illumination_geometries.gaussian_beam_illumination import GaussianBeamIlluminationGeometry # nopep8 -from .illumination_geometries.pencil_array_illumination import PencilArrayIlluminationGeometry # nopep8 -from .illumination_geometries.pencil_beam_illumination import PencilBeamIlluminationGeometry # nopep8 -from .illumination_geometries.disk_illumination import DiskIlluminationGeometry # nopep8 -from .illumination_geometries.rectangle_illumination import RectangleIlluminationGeometry # nopep8 -from .illumination_geometries.ring_illumination import RingIlluminationGeometry # nopep8 -from .illumination_geometries.ithera_msot_acuity_illumination import MSOTAcuityIlluminationGeometry # nopep8 -from .illumination_geometries.ithera_msot_invision_illumination import MSOTInVisionIlluminationGeometry # nopep8 -from .pa_devices.ithera_msot_invision import InVision256TF # nopep8 -from .pa_devices.ithera_msot_acuity import MSOTAcuityEcho # nopep8 -from .pa_devices.ithera_rsom import RSOMExplorerP50 # nopep8 +from .digital_device_twin_base import DigitalDeviceTwinBase +from .pa_devices import PhotoacousticDevice +from .detection_geometries import DetectionGeometryBase +from .detection_geometries.curved_array import CurvedArrayDetectionGeometry +from .detection_geometries.linear_array import LinearArrayDetectionGeometry +from .detection_geometries.planar_array import PlanarArrayDetectionGeometry +from .illumination_geometries import IlluminationGeometryBase +from .illumination_geometries.slit_illumination import SlitIlluminationGeometry +from .illumination_geometries.gaussian_beam_illumination import GaussianBeamIlluminationGeometry +from .illumination_geometries.pencil_array_illumination import PencilArrayIlluminationGeometry +from .illumination_geometries.pencil_beam_illumination import PencilBeamIlluminationGeometry +from .illumination_geometries.disk_illumination import DiskIlluminationGeometry +from .illumination_geometries.rectangle_illumination import RectangleIlluminationGeometry +from .illumination_geometries.ring_illumination import RingIlluminationGeometry +from .illumination_geometries.ithera_msot_acuity_illumination import MSOTAcuityIlluminationGeometry +from .illumination_geometries.ithera_msot_invision_illumination import MSOTInVisionIlluminationGeometry +from .pa_devices.ithera_msot_invision import InVision256TF +from .pa_devices.ithera_msot_acuity import MSOTAcuityEcho +from .pa_devices.ithera_rsom import RSOMExplorerP50 diff --git a/simpa/core/device_digital_twins/detection_geometries/__init__.py b/simpa/core/device_digital_twins/detection_geometries/__init__.py index eba6bdd5..e8d6c131 100644 --- a/simpa/core/device_digital_twins/detection_geometries/__init__.py +++ b/simpa/core/device_digital_twins/detection_geometries/__init__.py @@ -2,135 +2,4 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT -from abc import abstractmethod -from simpa.core.device_digital_twins import DigitalDeviceTwinBase -import numpy as np - - -class DetectionGeometryBase(DigitalDeviceTwinBase): - """ - This class is the base class for representing all detector geometries. - """ - - def __init__(self, number_detector_elements, detector_element_width_mm, - detector_element_length_mm, center_frequency_hz, bandwidth_percent, - sampling_frequency_mhz, device_position_mm: np.ndarray = None, - field_of_view_extent_mm: np.ndarray = None): - """ - - :param number_detector_elements: Total number of detector elements. - :type number_detector_elements: int - :param detector_element_width_mm: In-plane width of one detector element (pitch - distance between two - elements) in mm. - :type detector_element_width_mm: int, float - :param detector_element_length_mm: Out-of-plane length of one detector element in mm. - :type detector_element_length_mm: int, float - :param center_frequency_hz: Center frequency of the detector with approximately gaussian frequency response in - Hz. - :type center_frequency_hz: int, float - :param bandwidth_percent: Full width at half maximum in percent of the center frequency. - :type bandwidth_percent: int, float - :param sampling_frequency_mhz: Sampling frequency of the detector in MHz. - :type sampling_frequency_mhz: int, float - :param device_position_mm: Each device has an internal position which serves as origin for internal \ - representations of detector positions. - :type device_position_mm: ndarray - :param field_of_view_extent_mm: Field of view which is defined as a numpy array of the shape \ - [xs, xe, ys, ye, zs, ze], where x, y, and z denote the coordinate axes and s and e denote the start and end \ - positions. - :type field_of_view_extent_mm: ndarray - """ - super(DetectionGeometryBase, self).__init__(device_position_mm=device_position_mm, - field_of_view_extent_mm=field_of_view_extent_mm) - self.number_detector_elements = number_detector_elements - self.detector_element_width_mm = detector_element_width_mm - self.detector_element_length_mm = detector_element_length_mm - self.center_frequency_Hz = center_frequency_hz - self.bandwidth_percent = bandwidth_percent - self.sampling_frequency_MHz = sampling_frequency_mhz - - @abstractmethod - def get_detector_element_positions_base_mm(self) -> np.ndarray: - """ - Defines the abstract positions of the detection elements in an arbitrary coordinate system. - Typically, the center of the field of view is defined as the origin. - - To obtain the positions in an interpretable coordinate system, please use the other method:: - - get_detector_element_positions_accounting_for_device_position_mm - - :returns: A numpy array containing the position vectors of the detection elements. - - """ - pass - - def get_detector_element_positions_accounting_for_device_position_mm(self) -> np.ndarray: - """ - Similar to:: - - get_detector_element_positions_base_mm - - This method returns the absolute positions of the detection elements relative to the device - position in the imaged volume, where the device position is defined by the following tag:: - - Tags.DIGITAL_DEVICE_POSITION - - :returns: A numpy array containing the coordinates of the detection elements - - """ - abstract_element_positions = self.get_detector_element_positions_base_mm() - device_position = self.device_position_mm - return np.add(abstract_element_positions, device_position) - - def get_detector_element_positions_accounting_for_field_of_view(self) -> np.ndarray: - """ - Similar to:: - - get_detector_element_positions_base_mm - - This method returns the absolute positions of the detection elements relative to the device - position in the imaged volume, where the device position is defined by the following tag:: - - Tags.DIGITAL_DEVICE_POSITION - - :returns: A numpy array containing the coordinates of the detection elements - - """ - abstract_element_positions = np.copy(self.get_detector_element_positions_base_mm()) - field_of_view = self.field_of_view_extent_mm - x_half = (field_of_view[1] - field_of_view[0]) / 2 - y_half = (field_of_view[3] - field_of_view[2]) / 2 - if np.abs(x_half) < 1e-10: - abstract_element_positions[:, 0] = 0 - if np.abs(y_half) < 1e-10: - abstract_element_positions[:, 1] = 0 - - abstract_element_positions[:, 0] += x_half - abstract_element_positions[:, 1] += y_half - abstract_element_positions[:, 2] += field_of_view[4] - return abstract_element_positions - - @abstractmethod - def get_detector_element_orientations(self) -> np.ndarray: - """ - This method yields a normalised orientation vector for each detection element. The length of - this vector is the same as the one obtained via the position methods:: - - get_detector_element_positions_base_mm - get_detector_element_positions_accounting_for_device_position_mm - - :returns: a numpy array that contains normalised orientation vectors for each detection element - - """ - pass - - def serialize(self) -> dict: - serialized_device = self.__dict__ - return {DetectionGeometryBase: serialized_device} - - @staticmethod - def deserialize(dictionary_to_deserialize): - deserialized_device = DetectionGeometryBase() - for key, value in dictionary_to_deserialize.items(): - deserialized_device.__dict__[key] = value - return deserialized_device +from .detection_geometry_base import DetectionGeometryBase diff --git a/simpa/core/device_digital_twins/detection_geometries/detection_geometry_base.py b/simpa/core/device_digital_twins/detection_geometries/detection_geometry_base.py new file mode 100644 index 00000000..eba6bdd5 --- /dev/null +++ b/simpa/core/device_digital_twins/detection_geometries/detection_geometry_base.py @@ -0,0 +1,136 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT + +from abc import abstractmethod +from simpa.core.device_digital_twins import DigitalDeviceTwinBase +import numpy as np + + +class DetectionGeometryBase(DigitalDeviceTwinBase): + """ + This class is the base class for representing all detector geometries. + """ + + def __init__(self, number_detector_elements, detector_element_width_mm, + detector_element_length_mm, center_frequency_hz, bandwidth_percent, + sampling_frequency_mhz, device_position_mm: np.ndarray = None, + field_of_view_extent_mm: np.ndarray = None): + """ + + :param number_detector_elements: Total number of detector elements. + :type number_detector_elements: int + :param detector_element_width_mm: In-plane width of one detector element (pitch - distance between two + elements) in mm. + :type detector_element_width_mm: int, float + :param detector_element_length_mm: Out-of-plane length of one detector element in mm. + :type detector_element_length_mm: int, float + :param center_frequency_hz: Center frequency of the detector with approximately gaussian frequency response in + Hz. + :type center_frequency_hz: int, float + :param bandwidth_percent: Full width at half maximum in percent of the center frequency. + :type bandwidth_percent: int, float + :param sampling_frequency_mhz: Sampling frequency of the detector in MHz. + :type sampling_frequency_mhz: int, float + :param device_position_mm: Each device has an internal position which serves as origin for internal \ + representations of detector positions. + :type device_position_mm: ndarray + :param field_of_view_extent_mm: Field of view which is defined as a numpy array of the shape \ + [xs, xe, ys, ye, zs, ze], where x, y, and z denote the coordinate axes and s and e denote the start and end \ + positions. + :type field_of_view_extent_mm: ndarray + """ + super(DetectionGeometryBase, self).__init__(device_position_mm=device_position_mm, + field_of_view_extent_mm=field_of_view_extent_mm) + self.number_detector_elements = number_detector_elements + self.detector_element_width_mm = detector_element_width_mm + self.detector_element_length_mm = detector_element_length_mm + self.center_frequency_Hz = center_frequency_hz + self.bandwidth_percent = bandwidth_percent + self.sampling_frequency_MHz = sampling_frequency_mhz + + @abstractmethod + def get_detector_element_positions_base_mm(self) -> np.ndarray: + """ + Defines the abstract positions of the detection elements in an arbitrary coordinate system. + Typically, the center of the field of view is defined as the origin. + + To obtain the positions in an interpretable coordinate system, please use the other method:: + + get_detector_element_positions_accounting_for_device_position_mm + + :returns: A numpy array containing the position vectors of the detection elements. + + """ + pass + + def get_detector_element_positions_accounting_for_device_position_mm(self) -> np.ndarray: + """ + Similar to:: + + get_detector_element_positions_base_mm + + This method returns the absolute positions of the detection elements relative to the device + position in the imaged volume, where the device position is defined by the following tag:: + + Tags.DIGITAL_DEVICE_POSITION + + :returns: A numpy array containing the coordinates of the detection elements + + """ + abstract_element_positions = self.get_detector_element_positions_base_mm() + device_position = self.device_position_mm + return np.add(abstract_element_positions, device_position) + + def get_detector_element_positions_accounting_for_field_of_view(self) -> np.ndarray: + """ + Similar to:: + + get_detector_element_positions_base_mm + + This method returns the absolute positions of the detection elements relative to the device + position in the imaged volume, where the device position is defined by the following tag:: + + Tags.DIGITAL_DEVICE_POSITION + + :returns: A numpy array containing the coordinates of the detection elements + + """ + abstract_element_positions = np.copy(self.get_detector_element_positions_base_mm()) + field_of_view = self.field_of_view_extent_mm + x_half = (field_of_view[1] - field_of_view[0]) / 2 + y_half = (field_of_view[3] - field_of_view[2]) / 2 + if np.abs(x_half) < 1e-10: + abstract_element_positions[:, 0] = 0 + if np.abs(y_half) < 1e-10: + abstract_element_positions[:, 1] = 0 + + abstract_element_positions[:, 0] += x_half + abstract_element_positions[:, 1] += y_half + abstract_element_positions[:, 2] += field_of_view[4] + return abstract_element_positions + + @abstractmethod + def get_detector_element_orientations(self) -> np.ndarray: + """ + This method yields a normalised orientation vector for each detection element. The length of + this vector is the same as the one obtained via the position methods:: + + get_detector_element_positions_base_mm + get_detector_element_positions_accounting_for_device_position_mm + + :returns: a numpy array that contains normalised orientation vectors for each detection element + + """ + pass + + def serialize(self) -> dict: + serialized_device = self.__dict__ + return {DetectionGeometryBase: serialized_device} + + @staticmethod + def deserialize(dictionary_to_deserialize): + deserialized_device = DetectionGeometryBase() + for key, value in dictionary_to_deserialize.items(): + deserialized_device.__dict__[key] = value + return deserialized_device diff --git a/simpa/core/device_digital_twins/digital_device_twin_base.py b/simpa/core/device_digital_twins/digital_device_twin_base.py new file mode 100644 index 00000000..478361a0 --- /dev/null +++ b/simpa/core/device_digital_twins/digital_device_twin_base.py @@ -0,0 +1,132 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT + +from abc import abstractmethod +from simpa.log import Logger +import numpy as np +import hashlib +import uuid +from simpa.utils.serializer import SerializableSIMPAClass +from simpa.utils.calculate import are_equal + + +class DigitalDeviceTwinBase(SerializableSIMPAClass): + """ + This class represents a device that can be used for illumination, detection or a combined photoacoustic device + which has representations of both. + """ + + def __init__(self, device_position_mm=None, field_of_view_extent_mm=None): + """ + :param device_position_mm: Each device has an internal position which serves as origin for internal \ + representations of e.g. detector element positions or illuminator positions. + :type device_position_mm: ndarray + :param field_of_view_extent_mm: Field of view which is defined as a numpy array of the shape \ + [xs, xe, ys, ye, zs, ze], where x, y, and z denote the coordinate axes and s and e denote the start and end \ + positions. + :type field_of_view_extent_mm: ndarray + """ + if device_position_mm is None: + self.device_position_mm = np.array([0, 0, 0]) + else: + self.device_position_mm = device_position_mm + + if field_of_view_extent_mm is None: + self.field_of_view_extent_mm = np.asarray([-10, 10, -10, 10, -10, 10]) + else: + self.field_of_view_extent_mm = field_of_view_extent_mm + + self.logger = Logger() + + def __eq__(self, other): + """ + Checks each key, value pair in the devices. + """ + if isinstance(other, DigitalDeviceTwinBase): + if self.__dict__.keys() != other.__dict__.keys(): + return False + for self_key, self_value in self.__dict__.items(): + other_value = other.__dict__[self_key] + if not are_equal(self_value, other_value): + return False + return True + return False + + @abstractmethod + def check_settings_prerequisites(self, global_settings) -> bool: + """ + It might be that certain device geometries need a certain dimensionality of the simulated PAI volume, or that + it requires the existence of certain Tags in the global global_settings. + To this end, a PAI device should use this method to inform the user about a mismatch of the desired device and + throw a ValueError if that is the case. + + :param global_settings: Settings for the entire simulation pipeline. + :type global_settings: Settings + + :raises ValueError: raises a value error if the prerequisites are not matched. + :returns: True if the prerequisites are met, False if they are not met, but no exception has been raised. + :rtype: bool + + """ + pass + + @abstractmethod + def update_settings_for_use_of_model_based_volume_creator(self, global_settings): + """ + This method can be overwritten by a PA device if the device poses special constraints to the + volume that should be considered by the model-based volume creator. + + :param global_settings: Settings for the entire simulation pipeline. + :type global_settings: Settings + """ + pass + + def get_field_of_view_mm(self) -> np.ndarray: + """ + Returns the absolute field of view in mm where the probe position is already + accounted for. + It is defined as a numpy array of the shape [xs, xe, ys, ye, zs, ze], + where x, y, and z denote the coordinate axes and s and e denote the start and end + positions. + + :return: Absolute field of view in mm where the probe position is already accounted for. + :rtype: ndarray + """ + position = self.device_position_mm + field_of_view_extent = self.field_of_view_extent_mm + + field_of_view = np.asarray([position[0] + field_of_view_extent[0], + position[0] + field_of_view_extent[1], + position[1] + field_of_view_extent[2], + position[1] + field_of_view_extent[3], + position[2] + field_of_view_extent[4], + position[2] + field_of_view_extent[5] + ]) + if min(field_of_view) < 0: + self.logger.warning(f"The field of view of the chosen device is not fully within the simulated volume, " + f"field of view is: {field_of_view}") + field_of_view[field_of_view < 0] = 0 + + return field_of_view + + def generate_uuid(self): + """ + Generates a universally unique identifier (uuid) for each device. + :return: + """ + class_dict = self.__dict__ + m = hashlib.md5() + m.update(str(class_dict).encode('utf-8')) + return str(uuid.UUID(m.hexdigest())) + + def serialize(self) -> dict: + serialized_device = self.__dict__ + return {"DigitalDeviceTwinBase": serialized_device} + + @staticmethod + def deserialize(dictionary_to_deserialize): + deserialized_device = DigitalDeviceTwinBase( + device_position_mm=dictionary_to_deserialize["device_position_mm"], + field_of_view_extent_mm=dictionary_to_deserialize["field_of_view_extent_mm"]) + return deserialized_device diff --git a/simpa/core/device_digital_twins/illumination_geometries/__init__.py b/simpa/core/device_digital_twins/illumination_geometries/__init__.py index 6c5780bf..cd10d903 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/__init__.py +++ b/simpa/core/device_digital_twins/illumination_geometries/__init__.py @@ -2,69 +2,4 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT -from abc import abstractmethod -from simpa.core.device_digital_twins import DigitalDeviceTwinBase -from simpa.utils import Settings -import numpy as np - - -class IlluminationGeometryBase(DigitalDeviceTwinBase): - """ - This class is the base class for representing all illumination geometries. - """ - - def __init__(self, device_position_mm=None, source_direction_vector=None, field_of_view_extent_mm=None): - """ - :param device_position_mm: Each device has an internal position which serves as origin for internal \ - representations of illuminator positions. - :type device_position_mm: ndarray - - :param source_direction_vector: Direction of the illumination source. - :type source_direction_vector: ndarray - - :param field_of_view_extent_mm: Field of view which is defined as a numpy array of the shape \ - [xs, xe, ys, ye, zs, ze], where x, y, and z denote the coordinate axes and s and e denote the start and end \ - positions. - :type field_of_view_extent_mm: ndarray - """ - super(IlluminationGeometryBase, self).__init__(device_position_mm=device_position_mm, - field_of_view_extent_mm=field_of_view_extent_mm) - - if source_direction_vector is None: - self.source_direction_vector = [0, 0, 1] - else: - self.source_direction_vector = source_direction_vector - self.normalized_source_direction_vector = self.source_direction_vector / np.linalg.norm( - self.source_direction_vector) - - @abstractmethod - def get_mcx_illuminator_definition(self, global_settings) -> dict: - """ - IMPORTANT: This method creates a dictionary that contains tags as they are expected for the - mcx simulation tool to represent the illumination geometry of this device. - - :param global_settings: The global_settings instance containing the simulation instructions. - :type global_settings: Settings - - :return: Dictionary that includes all parameters needed for mcx. - :rtype: dict - """ - pass - - def check_settings_prerequisites(self, global_settings) -> bool: - return True - - def update_settings_for_use_of_model_based_volume_creator(self, global_settings) -> Settings: - return global_settings - - def serialize(self) -> dict: - serialized_device = self.__dict__ - device_dict = {"IlluminationGeometryBase": serialized_device} - return device_dict - - @staticmethod - def deserialize(dictionary_to_deserialize): - deserialized_device = IlluminationGeometryBase() - for key, value in dictionary_to_deserialize.items(): - deserialized_device.__dict__[key] = value - return deserialized_device +from .illumination_geometry_base import IlluminationGeometryBase diff --git a/simpa/core/device_digital_twins/illumination_geometries/illumination_geometry_base.py b/simpa/core/device_digital_twins/illumination_geometries/illumination_geometry_base.py new file mode 100644 index 00000000..6c5780bf --- /dev/null +++ b/simpa/core/device_digital_twins/illumination_geometries/illumination_geometry_base.py @@ -0,0 +1,70 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT + +from abc import abstractmethod +from simpa.core.device_digital_twins import DigitalDeviceTwinBase +from simpa.utils import Settings +import numpy as np + + +class IlluminationGeometryBase(DigitalDeviceTwinBase): + """ + This class is the base class for representing all illumination geometries. + """ + + def __init__(self, device_position_mm=None, source_direction_vector=None, field_of_view_extent_mm=None): + """ + :param device_position_mm: Each device has an internal position which serves as origin for internal \ + representations of illuminator positions. + :type device_position_mm: ndarray + + :param source_direction_vector: Direction of the illumination source. + :type source_direction_vector: ndarray + + :param field_of_view_extent_mm: Field of view which is defined as a numpy array of the shape \ + [xs, xe, ys, ye, zs, ze], where x, y, and z denote the coordinate axes and s and e denote the start and end \ + positions. + :type field_of_view_extent_mm: ndarray + """ + super(IlluminationGeometryBase, self).__init__(device_position_mm=device_position_mm, + field_of_view_extent_mm=field_of_view_extent_mm) + + if source_direction_vector is None: + self.source_direction_vector = [0, 0, 1] + else: + self.source_direction_vector = source_direction_vector + self.normalized_source_direction_vector = self.source_direction_vector / np.linalg.norm( + self.source_direction_vector) + + @abstractmethod + def get_mcx_illuminator_definition(self, global_settings) -> dict: + """ + IMPORTANT: This method creates a dictionary that contains tags as they are expected for the + mcx simulation tool to represent the illumination geometry of this device. + + :param global_settings: The global_settings instance containing the simulation instructions. + :type global_settings: Settings + + :return: Dictionary that includes all parameters needed for mcx. + :rtype: dict + """ + pass + + def check_settings_prerequisites(self, global_settings) -> bool: + return True + + def update_settings_for_use_of_model_based_volume_creator(self, global_settings) -> Settings: + return global_settings + + def serialize(self) -> dict: + serialized_device = self.__dict__ + device_dict = {"IlluminationGeometryBase": serialized_device} + return device_dict + + @staticmethod + def deserialize(dictionary_to_deserialize): + deserialized_device = IlluminationGeometryBase() + for key, value in dictionary_to_deserialize.items(): + deserialized_device.__dict__[key] = value + return deserialized_device diff --git a/simpa/core/device_digital_twins/pa_devices/__init__.py b/simpa/core/device_digital_twins/pa_devices/__init__.py index 38b0d202..7b8845f7 100644 --- a/simpa/core/device_digital_twins/pa_devices/__init__.py +++ b/simpa/core/device_digital_twins/pa_devices/__init__.py @@ -2,160 +2,4 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT -import numpy as np -from abc import ABC -from simpa.core.device_digital_twins import DigitalDeviceTwinBase - - -class PhotoacousticDevice(DigitalDeviceTwinBase, ABC): - """Base class of a photoacoustic device. It consists of one detection geometry that describes the geometry of the - single detector elements and a list of illuminators. - - A Photoacoustic Device can be initialized as follows:: - - import simpa as sp - import numpy as np - - # Initialise a PhotoacousticDevice with its position and field of view - device = sp.PhotoacousticDevice(device_position_mm=np.array([10, 10, 0]), - field_of_view_extent_mm=np.array([-20, 20, 0, 0, 0, 20])) - - # Option 1) Set the detection geometry position relative to the PhotoacousticDevice - device.set_detection_geometry(sp.DetectionGeometry(), - detector_position_relative_to_pa_device=np.array([0, 0, -10])) - - # Option 2) Set the detection geometry position absolute - device.set_detection_geometry( - sp.DetectionGeometryBase(device_position_mm=np.array([10, 10, -10]))) - - # Option 1) Add the illumination geometry position relative to the PhotoacousticDevice - device.add_illumination_geometry(sp.IlluminationGeometry(), - illuminator_position_relative_to_pa_device=np.array([0, 0, 0])) - - # Option 2) Add the illumination geometry position absolute - device.add_illumination_geometry( - sp.IlluminationGeometryBase(device_position_mm=np.array([10, 10, 0])) - - Attributes: - detection_geometry (DetectionGeometryBase): Geometry of the detector elements. - illumination_geometries (list): List of illuminations defined by :py:class:`IlluminationGeometryBase`. - """ - - def __init__(self, device_position_mm=None, field_of_view_extent_mm=None): - """ - :param device_position_mm: Each device has an internal position which serves as origin for internal \ - representations of e.g. detector element positions or illuminator positions. - :type device_position_mm: ndarray - :param field_of_view_extent_mm: Field of view which is defined as a numpy array of the shape \ - [xs, xe, ys, ye, zs, ze], where x, y, and z denote the coordinate axes and s and e denote the start and end \ - positions. - :type field_of_view_extent_mm: ndarray - """ - super(PhotoacousticDevice, self).__init__(device_position_mm=device_position_mm, - field_of_view_extent_mm=field_of_view_extent_mm) - self.detection_geometry = None - self.illumination_geometries = [] - - def set_detection_geometry(self, detection_geometry, - detector_position_relative_to_pa_device=None): - """Sets the detection geometry for the PA device. The detection geometry can be instantiated with an absolute - position or it can be instantiated without the device_position_mm argument but a position relative to the - position of the PhotoacousticDevice. If both absolute and relative positions are given, the absolute position - is chosen as position of the detection geometry. - - :param detection_geometry: Detection geometry of the PA device. - :type detection_geometry: DetectionGeometryBase - :param detector_position_relative_to_pa_device: Position of the detection geometry relative to the PA device. - :type detector_position_relative_to_pa_device: ndarray - :raises ValueError: if the detection_geometry is None - - """ - if detection_geometry is None: - msg = "The given detection_geometry must not be None!" - self.logger.critical(msg) - raise ValueError(msg) - if np.linalg.norm(detection_geometry.device_position_mm) == 0 and \ - detector_position_relative_to_pa_device is not None: - detection_geometry.device_position_mm = np.add(self.device_position_mm, - detector_position_relative_to_pa_device) - self.detection_geometry = detection_geometry - - def add_illumination_geometry(self, illumination_geometry, illuminator_position_relative_to_pa_device=None): - """Adds an illuminator to the PA device. The illumination geometry can be instantiated with an absolute - position or it can be instantiated without the device_position_mm argument but a position relative to the - position of the PhotoacousticDevice. If both absolute and relative positions are given, the absolute position - is chosen as position of the illumination geometry. - - :param illumination_geometry: Geometry of the illuminator. - :type illumination_geometry: IlluminationGeometryBase - :param illuminator_position_relative_to_pa_device: Position of the illuminator relative to the PA device. - :type illuminator_position_relative_to_pa_device: ndarray - :raises ValueError: if the illumination_geometry is None - - """ - if illumination_geometry is None: - msg = "The given illumination_geometry must not be None!" - self.logger.critical(msg) - raise ValueError(msg) - if np.linalg.norm(illumination_geometry.device_position_mm) == 0: - if illuminator_position_relative_to_pa_device is not None: - illumination_geometry.device_position_mm = np.add(self.device_position_mm, - illuminator_position_relative_to_pa_device) - else: - illumination_geometry.device_position_mm = self.device_position_mm - self.illumination_geometries.append(illumination_geometry) - - def get_detection_geometry(self): - """ - :return: None if no detection geometry was set or an instance of DetectionGeometryBase. - :rtype: None, DetectionGeometryBase - """ - return self.detection_geometry - - def get_illumination_geometry(self): - """ - :return: None, if no illumination geometry was defined, - an instance of IlluminationGeometryBase if exactly one geometry was defined, - a list of IlluminationGeometryBase instances if more than one device was defined. - :rtype: None, IlluminationGeometryBase - """ - if len(self.illumination_geometries) == 0: - return None - - if len(self.illumination_geometries) == 1: - return self.illumination_geometries[0] - - return self.illumination_geometries - - def check_settings_prerequisites(self, global_settings) -> bool: - _result = True - if self.detection_geometry is not None \ - and not self.detection_geometry.check_settings_prerequisites(global_settings): - _result = False - for illumination_geometry in self.illumination_geometries: - if illumination_geometry is not None \ - and not illumination_geometry.check_settings_prerequisites(global_settings): - _result = False - return _result - - def update_settings_for_use_of_model_based_volume_creator(self, global_settings): - pass - - def serialize(self) -> dict: - serialized_device = self.__dict__ - device_dict = {"PhotoacousticDevice": serialized_device} - return device_dict - - @staticmethod - def deserialize(dictionary_to_deserialize): - deserialized_device = PhotoacousticDevice( - device_position_mm=dictionary_to_deserialize["device_position_mm"], - field_of_view_extent_mm=dictionary_to_deserialize["field_of_view_extent_mm"]) - det_geometry = dictionary_to_deserialize["detection_geometry"] - if det_geometry != "None": - deserialized_device.set_detection_geometry(dictionary_to_deserialize["detection_geometry"]) - if "illumination_geometries" in dictionary_to_deserialize: - for illumination_geometry in dictionary_to_deserialize["illumination_geometries"]: - deserialized_device.illumination_geometries.append(illumination_geometry) - - return deserialized_device +from .photoacoustic_device import PhotoacousticDevice 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/device_digital_twins/pa_devices/photoacoustic_device.py b/simpa/core/device_digital_twins/pa_devices/photoacoustic_device.py new file mode 100644 index 00000000..38b0d202 --- /dev/null +++ b/simpa/core/device_digital_twins/pa_devices/photoacoustic_device.py @@ -0,0 +1,161 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT + +import numpy as np +from abc import ABC +from simpa.core.device_digital_twins import DigitalDeviceTwinBase + + +class PhotoacousticDevice(DigitalDeviceTwinBase, ABC): + """Base class of a photoacoustic device. It consists of one detection geometry that describes the geometry of the + single detector elements and a list of illuminators. + + A Photoacoustic Device can be initialized as follows:: + + import simpa as sp + import numpy as np + + # Initialise a PhotoacousticDevice with its position and field of view + device = sp.PhotoacousticDevice(device_position_mm=np.array([10, 10, 0]), + field_of_view_extent_mm=np.array([-20, 20, 0, 0, 0, 20])) + + # Option 1) Set the detection geometry position relative to the PhotoacousticDevice + device.set_detection_geometry(sp.DetectionGeometry(), + detector_position_relative_to_pa_device=np.array([0, 0, -10])) + + # Option 2) Set the detection geometry position absolute + device.set_detection_geometry( + sp.DetectionGeometryBase(device_position_mm=np.array([10, 10, -10]))) + + # Option 1) Add the illumination geometry position relative to the PhotoacousticDevice + device.add_illumination_geometry(sp.IlluminationGeometry(), + illuminator_position_relative_to_pa_device=np.array([0, 0, 0])) + + # Option 2) Add the illumination geometry position absolute + device.add_illumination_geometry( + sp.IlluminationGeometryBase(device_position_mm=np.array([10, 10, 0])) + + Attributes: + detection_geometry (DetectionGeometryBase): Geometry of the detector elements. + illumination_geometries (list): List of illuminations defined by :py:class:`IlluminationGeometryBase`. + """ + + def __init__(self, device_position_mm=None, field_of_view_extent_mm=None): + """ + :param device_position_mm: Each device has an internal position which serves as origin for internal \ + representations of e.g. detector element positions or illuminator positions. + :type device_position_mm: ndarray + :param field_of_view_extent_mm: Field of view which is defined as a numpy array of the shape \ + [xs, xe, ys, ye, zs, ze], where x, y, and z denote the coordinate axes and s and e denote the start and end \ + positions. + :type field_of_view_extent_mm: ndarray + """ + super(PhotoacousticDevice, self).__init__(device_position_mm=device_position_mm, + field_of_view_extent_mm=field_of_view_extent_mm) + self.detection_geometry = None + self.illumination_geometries = [] + + def set_detection_geometry(self, detection_geometry, + detector_position_relative_to_pa_device=None): + """Sets the detection geometry for the PA device. The detection geometry can be instantiated with an absolute + position or it can be instantiated without the device_position_mm argument but a position relative to the + position of the PhotoacousticDevice. If both absolute and relative positions are given, the absolute position + is chosen as position of the detection geometry. + + :param detection_geometry: Detection geometry of the PA device. + :type detection_geometry: DetectionGeometryBase + :param detector_position_relative_to_pa_device: Position of the detection geometry relative to the PA device. + :type detector_position_relative_to_pa_device: ndarray + :raises ValueError: if the detection_geometry is None + + """ + if detection_geometry is None: + msg = "The given detection_geometry must not be None!" + self.logger.critical(msg) + raise ValueError(msg) + if np.linalg.norm(detection_geometry.device_position_mm) == 0 and \ + detector_position_relative_to_pa_device is not None: + detection_geometry.device_position_mm = np.add(self.device_position_mm, + detector_position_relative_to_pa_device) + self.detection_geometry = detection_geometry + + def add_illumination_geometry(self, illumination_geometry, illuminator_position_relative_to_pa_device=None): + """Adds an illuminator to the PA device. The illumination geometry can be instantiated with an absolute + position or it can be instantiated without the device_position_mm argument but a position relative to the + position of the PhotoacousticDevice. If both absolute and relative positions are given, the absolute position + is chosen as position of the illumination geometry. + + :param illumination_geometry: Geometry of the illuminator. + :type illumination_geometry: IlluminationGeometryBase + :param illuminator_position_relative_to_pa_device: Position of the illuminator relative to the PA device. + :type illuminator_position_relative_to_pa_device: ndarray + :raises ValueError: if the illumination_geometry is None + + """ + if illumination_geometry is None: + msg = "The given illumination_geometry must not be None!" + self.logger.critical(msg) + raise ValueError(msg) + if np.linalg.norm(illumination_geometry.device_position_mm) == 0: + if illuminator_position_relative_to_pa_device is not None: + illumination_geometry.device_position_mm = np.add(self.device_position_mm, + illuminator_position_relative_to_pa_device) + else: + illumination_geometry.device_position_mm = self.device_position_mm + self.illumination_geometries.append(illumination_geometry) + + def get_detection_geometry(self): + """ + :return: None if no detection geometry was set or an instance of DetectionGeometryBase. + :rtype: None, DetectionGeometryBase + """ + return self.detection_geometry + + def get_illumination_geometry(self): + """ + :return: None, if no illumination geometry was defined, + an instance of IlluminationGeometryBase if exactly one geometry was defined, + a list of IlluminationGeometryBase instances if more than one device was defined. + :rtype: None, IlluminationGeometryBase + """ + if len(self.illumination_geometries) == 0: + return None + + if len(self.illumination_geometries) == 1: + return self.illumination_geometries[0] + + return self.illumination_geometries + + def check_settings_prerequisites(self, global_settings) -> bool: + _result = True + if self.detection_geometry is not None \ + and not self.detection_geometry.check_settings_prerequisites(global_settings): + _result = False + for illumination_geometry in self.illumination_geometries: + if illumination_geometry is not None \ + and not illumination_geometry.check_settings_prerequisites(global_settings): + _result = False + return _result + + def update_settings_for_use_of_model_based_volume_creator(self, global_settings): + pass + + def serialize(self) -> dict: + serialized_device = self.__dict__ + device_dict = {"PhotoacousticDevice": serialized_device} + return device_dict + + @staticmethod + def deserialize(dictionary_to_deserialize): + deserialized_device = PhotoacousticDevice( + device_position_mm=dictionary_to_deserialize["device_position_mm"], + field_of_view_extent_mm=dictionary_to_deserialize["field_of_view_extent_mm"]) + det_geometry = dictionary_to_deserialize["detection_geometry"] + if det_geometry != "None": + deserialized_device.set_detection_geometry(dictionary_to_deserialize["detection_geometry"]) + if "illumination_geometries" in dictionary_to_deserialize: + for illumination_geometry in dictionary_to_deserialize["illumination_geometries"]: + deserialized_device.illumination_geometries.append(illumination_geometry) + + return deserialized_device diff --git a/simpa/core/pipeline_element_base.py b/simpa/core/pipeline_element_base.py new file mode 100644 index 00000000..2fdcfb10 --- /dev/null +++ b/simpa/core/pipeline_element_base.py @@ -0,0 +1,33 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT +from abc import abstractmethod + +from simpa.core.device_digital_twins import DigitalDeviceTwinBase +from simpa.log import Logger +from simpa.utils import Settings +from simpa.utils.processing_device import get_processing_device + + +class PipelineElementBase: + """ + Defines a pipeline element (either simulation or processing module) that implements a run method and can be called by running the pipeline's simulate method. + """ + + def __init__(self, global_settings: Settings): + """ + :param global_settings: The SIMPA settings dictionary + :type global_settings: Settings + """ + self.logger = Logger() + self.global_settings = global_settings + self.torch_device = get_processing_device(self.global_settings) + + @abstractmethod + def run(self, digital_device_twin: DigitalDeviceTwinBase): + """ + Executes the respective simulation module + + :param digital_device_twin: The digital twin that can be used by the digital device_twin. + """ + pass diff --git a/simpa/core/processing_components/__init__.py b/simpa/core/processing_components/__init__.py index 881a8f2d..2c0e2695 100644 --- a/simpa/core/processing_components/__init__.py +++ b/simpa/core/processing_components/__init__.py @@ -2,20 +2,4 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT -from abc import ABC -from simpa.core import PipelineModule - - -class ProcessingComponent(PipelineModule, ABC): - """ - Defines a pipeline processing component, which can be used to pre- or post-process simulation data. - """ - - def __init__(self, global_settings, component_settings_key: str): - """ - Initialises the ProcessingComponent. - - :param component_settings_key: The key where the component settings are stored in - """ - super(ProcessingComponent, self).__init__(global_settings=global_settings) - self.component_settings = global_settings[component_settings_key] +from .processing_component_base import ProcessingComponentBase diff --git a/simpa/core/processing_components/monospectral/field_of_view_cropping.py b/simpa/core/processing_components/monospectral/field_of_view_cropping.py index b950dc8f..68ae16a2 100644 --- a/simpa/core/processing_components/monospectral/field_of_view_cropping.py +++ b/simpa/core/processing_components/monospectral/field_of_view_cropping.py @@ -5,12 +5,12 @@ from simpa.utils import Tags, Settings from simpa.utils.constants import property_tags, wavelength_independent_properties, toolkit_tags from simpa.io_handling import load_data_field, save_data_field -from simpa.core.processing_components import ProcessingComponent +from simpa.core.processing_components import ProcessingComponentBase from simpa.core.device_digital_twins import DigitalDeviceTwinBase, PhotoacousticDevice import numpy as np -class FieldOfViewCropping(ProcessingComponent): +class FieldOfViewCropping(ProcessingComponentBase): def __init__(self, global_settings, settings_key=None): if settings_key is None: diff --git a/simpa/core/processing_components/monospectral/iterative_qPAI_algorithm.py b/simpa/core/processing_components/monospectral/iterative_qPAI_algorithm.py index 3a942564..df4c61da 100644 --- a/simpa/core/processing_components/monospectral/iterative_qPAI_algorithm.py +++ b/simpa/core/processing_components/monospectral/iterative_qPAI_algorithm.py @@ -11,16 +11,16 @@ from simpa.utils.libraries.literature_values import OpticalTissueProperties, StandardProperties from simpa.utils.libraries.molecule_library import MolecularComposition from simpa.utils.calculate import calculate_gruneisen_parameter_from_temperature -from simpa.core.simulation_modules.optical_simulation_module.optical_forward_model_mcx_adapter import \ +from simpa.core.simulation_modules.optical_module.mcx_adapter import \ MCXAdapter from simpa.utils import Settings from simpa.io_handling import save_data_field, load_data_field from simpa.utils import TISSUE_LIBRARY -from simpa.core.processing_components import ProcessingComponent +from simpa.core.processing_components import ProcessingComponentBase import os -class IterativeqPAI(ProcessingComponent): +class IterativeqPAI(ProcessingComponentBase): """ Applies iterative qPAI Algorithm [1] on simulated initial pressure map and saves the reconstruction result in the hdf5 output file. If a 2-d map of initial_pressure is passed the algorithm saves @@ -45,7 +45,7 @@ class IterativeqPAI(ProcessingComponent): """ def __init__(self, global_settings, component_settings_key: str): - super(ProcessingComponent, self).__init__(global_settings=global_settings) + super(ProcessingComponentBase, self).__init__(global_settings=global_settings) self.global_settings = global_settings self.optical_settings = global_settings.get_optical_settings() @@ -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/processing_components/monospectral/noise/gamma_noise.py b/simpa/core/processing_components/monospectral/noise/gamma_noise.py index 011ddffa..5d87ca6c 100644 --- a/simpa/core/processing_components/monospectral/noise/gamma_noise.py +++ b/simpa/core/processing_components/monospectral/noise/gamma_noise.py @@ -5,13 +5,13 @@ import numpy as np import torch -from simpa.core.processing_components import ProcessingComponent +from simpa.core.processing_components import ProcessingComponentBase from simpa.io_handling import load_data_field, save_data_field from simpa.utils import Tags from simpa.utils.quality_assurance.data_sanity_testing import assert_array_well_defined -class GammaNoise(ProcessingComponent): +class GammaNoise(ProcessingComponentBase): """ Applies Gamma noise to the defined data field. The noise will be applied to all wavelengths. diff --git a/simpa/core/processing_components/monospectral/noise/gaussian_noise.py b/simpa/core/processing_components/monospectral/noise/gaussian_noise.py index 5c31a50d..98dba54f 100644 --- a/simpa/core/processing_components/monospectral/noise/gaussian_noise.py +++ b/simpa/core/processing_components/monospectral/noise/gaussian_noise.py @@ -5,13 +5,13 @@ from simpa.utils import Tags from simpa.utils import EPS from simpa.io_handling import load_data_field, save_data_field -from simpa.core.processing_components import ProcessingComponent +from simpa.core.processing_components import ProcessingComponentBase from simpa.utils.quality_assurance.data_sanity_testing import assert_array_well_defined import numpy as np import torch -class GaussianNoise(ProcessingComponent): +class GaussianNoise(ProcessingComponentBase): """ Applies Gaussian noise to the defined data field. The noise will be applied to all wavelengths. diff --git a/simpa/core/processing_components/monospectral/noise/poisson_noise.py b/simpa/core/processing_components/monospectral/noise/poisson_noise.py index 0fd544e6..0f847f17 100644 --- a/simpa/core/processing_components/monospectral/noise/poisson_noise.py +++ b/simpa/core/processing_components/monospectral/noise/poisson_noise.py @@ -4,13 +4,13 @@ from simpa.utils import Tags from simpa.io_handling import load_data_field, save_data_field -from simpa.core.processing_components import ProcessingComponent +from simpa.core.processing_components import ProcessingComponentBase from simpa.utils.quality_assurance.data_sanity_testing import assert_array_well_defined import numpy as np import torch -class PoissonNoise(ProcessingComponent): +class PoissonNoise(ProcessingComponentBase): """ Applies Poisson noise to the defined data field. The noise will be applied to all wavelengths. diff --git a/simpa/core/processing_components/monospectral/noise/salt_and_pepper_noise.py b/simpa/core/processing_components/monospectral/noise/salt_and_pepper_noise.py index 4204b32e..b6c70fba 100644 --- a/simpa/core/processing_components/monospectral/noise/salt_and_pepper_noise.py +++ b/simpa/core/processing_components/monospectral/noise/salt_and_pepper_noise.py @@ -4,13 +4,13 @@ from simpa.utils import Tags from simpa.io_handling import load_data_field, save_data_field -from simpa.core.processing_components import ProcessingComponent +from simpa.core.processing_components import ProcessingComponentBase from simpa.utils.quality_assurance.data_sanity_testing import assert_array_well_defined import numpy as np import torch -class SaltAndPepperNoise(ProcessingComponent): +class SaltAndPepperNoise(ProcessingComponentBase): """ Applies salt and pepper noise to the defined data field. The noise will be applied to all wavelengths. diff --git a/simpa/core/processing_components/monospectral/noise/uniform_noise.py b/simpa/core/processing_components/monospectral/noise/uniform_noise.py index 5610dc43..befc0de4 100644 --- a/simpa/core/processing_components/monospectral/noise/uniform_noise.py +++ b/simpa/core/processing_components/monospectral/noise/uniform_noise.py @@ -4,13 +4,13 @@ from simpa.utils import Tags from simpa.io_handling import load_data_field, save_data_field -from simpa.core.processing_components import ProcessingComponent +from simpa.core.processing_components import ProcessingComponentBase from simpa.utils.quality_assurance.data_sanity_testing import assert_array_well_defined import numpy as np import torch -class UniformNoise(ProcessingComponent): +class UniformNoise(ProcessingComponentBase): """ Applies uniform noise to the defined data field. The noise will be applied to all wavelengths. diff --git a/simpa/core/processing_components/multispectral/__init__.py b/simpa/core/processing_components/multispectral/__init__.py index d179c455..bf4ed09d 100644 --- a/simpa/core/processing_components/multispectral/__init__.py +++ b/simpa/core/processing_components/multispectral/__init__.py @@ -1,57 +1,4 @@ # SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT -from simpa.io_handling import load_data_field -from simpa.utils import Tags -from simpa.log import Logger -import numpy as np -from abc import ABC, abstractmethod - - -class MultispectralProcessingAlgorithm(ABC): - """ - A MultispectralProcessingAlgorithm class represents an algorithm that works with multispectral input data. - """ - - def __init__(self, global_settings, component_settings_key: str): - """ - Instantiates a multispectral processing algorithm. - - Per default, this methods loads all data from a certain - Tags.DATA_FIELD into a data array for all - Tags.WAVELENGTHS. - - """ - if component_settings_key is None: - raise KeyError("The component settings must be set for a multispectral" - "processing algorithm!") - self.component_settings = global_settings[component_settings_key] - - if Tags.WAVELENGTHS not in self.component_settings: - raise KeyError("Tags.WAVELENGTHS must be in the component_settings of a multispectral processing algorithm") - - if Tags.DATA_FIELD not in self.component_settings: - raise KeyError("Tags.DATA_FIELD must be in the component_settings of a multispectral processing algorithm") - - self.logger = Logger() - self.global_settings = global_settings - self.wavelengths = self.component_settings[Tags.WAVELENGTHS] - self.data_field = self.component_settings[Tags.DATA_FIELD] - - self.data = list() - for i in range(len(self.wavelengths)): - self.data.append(load_data_field(self.global_settings[Tags.SIMPA_OUTPUT_PATH], - self.data_field, - self.wavelengths[i])) - - self.data = np.asarray(self.data) - if Tags.SIGNAL_THRESHOLD in self.component_settings: - self.data[self.data < self.component_settings[Tags.SIGNAL_THRESHOLD]*np.max(self.data)] = 0 - - @abstractmethod - def run(self): - """ - This method must be implemented by the multispectral algorithm, such that - any multispectral algorithm can be executed by invoking the run method. - """ - pass +from .multispectral_processing_algorithm import MultispectralProcessingAlgorithm diff --git a/simpa/core/processing_components/multispectral/multispectral_processing_algorithm.py b/simpa/core/processing_components/multispectral/multispectral_processing_algorithm.py new file mode 100644 index 00000000..d179c455 --- /dev/null +++ b/simpa/core/processing_components/multispectral/multispectral_processing_algorithm.py @@ -0,0 +1,57 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT +from simpa.io_handling import load_data_field +from simpa.utils import Tags +from simpa.log import Logger +import numpy as np +from abc import ABC, abstractmethod + + +class MultispectralProcessingAlgorithm(ABC): + """ + A MultispectralProcessingAlgorithm class represents an algorithm that works with multispectral input data. + """ + + def __init__(self, global_settings, component_settings_key: str): + """ + Instantiates a multispectral processing algorithm. + + Per default, this methods loads all data from a certain + Tags.DATA_FIELD into a data array for all + Tags.WAVELENGTHS. + + """ + if component_settings_key is None: + raise KeyError("The component settings must be set for a multispectral" + "processing algorithm!") + self.component_settings = global_settings[component_settings_key] + + if Tags.WAVELENGTHS not in self.component_settings: + raise KeyError("Tags.WAVELENGTHS must be in the component_settings of a multispectral processing algorithm") + + if Tags.DATA_FIELD not in self.component_settings: + raise KeyError("Tags.DATA_FIELD must be in the component_settings of a multispectral processing algorithm") + + self.logger = Logger() + self.global_settings = global_settings + self.wavelengths = self.component_settings[Tags.WAVELENGTHS] + self.data_field = self.component_settings[Tags.DATA_FIELD] + + self.data = list() + for i in range(len(self.wavelengths)): + self.data.append(load_data_field(self.global_settings[Tags.SIMPA_OUTPUT_PATH], + self.data_field, + self.wavelengths[i])) + + self.data = np.asarray(self.data) + if Tags.SIGNAL_THRESHOLD in self.component_settings: + self.data[self.data < self.component_settings[Tags.SIGNAL_THRESHOLD]*np.max(self.data)] = 0 + + @abstractmethod + def run(self): + """ + This method must be implemented by the multispectral algorithm, such that + any multispectral algorithm can be executed by invoking the run method. + """ + pass diff --git a/simpa/core/processing_components/processing_component_base.py b/simpa/core/processing_components/processing_component_base.py new file mode 100644 index 00000000..e74f6686 --- /dev/null +++ b/simpa/core/processing_components/processing_component_base.py @@ -0,0 +1,21 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT + +from abc import ABC +from simpa.core import PipelineElementBase + + +class ProcessingComponentBase(PipelineElementBase, ABC): + """ + Defines a pipeline processing component, which can be used to pre- or post-process simulation data. + """ + + def __init__(self, global_settings, component_settings_key: str): + """ + Initialises the ProcessingComponent object. + + :param component_settings_key: The key where the component settings are stored in + """ + super(ProcessingComponentBase, self).__init__(global_settings=global_settings) + self.component_settings = global_settings[component_settings_key] diff --git a/simpa/core/simulation_modules/__init__.py b/simpa/core/simulation_modules/__init__.py index 9624110c..2da4c6d2 100644 --- a/simpa/core/simulation_modules/__init__.py +++ b/simpa/core/simulation_modules/__init__.py @@ -2,43 +2,4 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT -from abc import abstractmethod - -from simpa.core import PipelineModule -from simpa.utils import Settings, Tags -from typing import List - - -class SimulationModule(PipelineModule): - """ - Defines a simulation module that is a step in the simulation pipeline. - Each simulation module can only be one of Volume Creation, Light Propagation Modeling, Acoustic Wave Propagation Modeling, Image Reconstruction. - """ - - def __init__(self, global_settings: Settings): - """ - :param global_settings: The SIMPA settings dictionary - :type global_settings: Settings - """ - super(SimulationModule, self).__init__(global_settings=global_settings) - self.component_settings = self.load_component_settings() - if self.component_settings is None: - raise ValueError("The component settings should not be None at this point") - - @abstractmethod - def load_component_settings(self) -> Settings: - """ - :return: Loads component settings corresponding to this simulation component - """ - pass - - def get_additional_flags(self) -> List[str]: - """Reads the list of additional flags from the corresponding component settings Tags.ADDITIONAL_FLAGS - - :return: List[str]: list of additional flags - """ - cmd = [] - 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 +from .simulation_module_base import SimulationModuleBase diff --git a/simpa/core/simulation_modules/acoustic_module/__init__.py b/simpa/core/simulation_modules/acoustic_module/__init__.py new file mode 100644 index 00000000..9ad0c1ea --- /dev/null +++ b/simpa/core/simulation_modules/acoustic_module/__init__.py @@ -0,0 +1,4 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT +from .acoustic_adapter_base import AcousticAdapterBase diff --git a/simpa/core/simulation_modules/acoustic_forward_module/__init__.py b/simpa/core/simulation_modules/acoustic_module/acoustic_adapter_base.py similarity index 94% rename from simpa/core/simulation_modules/acoustic_forward_module/__init__.py rename to simpa/core/simulation_modules/acoustic_module/acoustic_adapter_base.py index e9f33d10..85c4dd07 100644 --- a/simpa/core/simulation_modules/acoustic_forward_module/__init__.py +++ b/simpa/core/simulation_modules/acoustic_module/acoustic_adapter_base.py @@ -4,7 +4,7 @@ from abc import abstractmethod import numpy as np -from simpa.core.simulation_modules import SimulationModule +from simpa.core.simulation_modules import SimulationModuleBase from simpa.utils import Tags, Settings from simpa.io_handling.io_hdf5 import save_hdf5 from simpa.utils.dict_path_manager import generate_dict_path @@ -12,7 +12,7 @@ from simpa.utils.quality_assurance.data_sanity_testing import assert_array_well_defined -class AcousticForwardModelBaseAdapter(SimulationModule): +class AcousticAdapterBase(SimulationModuleBase): """ This method is the entry method for running an acoustic forward model. It is invoked in the *simpa.core.simulation.simulate* method, but can also be called @@ -32,7 +32,7 @@ class AcousticForwardModelBaseAdapter(SimulationModule): """ def __init__(self, global_settings: Settings): - super(AcousticForwardModelBaseAdapter, self).__init__(global_settings=global_settings) + super(AcousticAdapterBase, self).__init__(global_settings=global_settings) def load_component_settings(self) -> Settings: """Implements abstract method to serve acoustic settings as component settings diff --git a/simpa/core/simulation_modules/acoustic_forward_module/acoustic_forward_model_test_adapter.py b/simpa/core/simulation_modules/acoustic_module/acoustic_test_adapter.py similarity index 75% rename from simpa/core/simulation_modules/acoustic_forward_module/acoustic_forward_model_test_adapter.py rename to simpa/core/simulation_modules/acoustic_module/acoustic_test_adapter.py index 9350da60..9477cb3e 100644 --- a/simpa/core/simulation_modules/acoustic_forward_module/acoustic_forward_model_test_adapter.py +++ b/simpa/core/simulation_modules/acoustic_module/acoustic_test_adapter.py @@ -4,10 +4,10 @@ import numpy as np from simpa.utils import Tags -from simpa.core.simulation_modules.acoustic_forward_module import AcousticForwardModelBaseAdapter +from simpa.core.simulation_modules.acoustic_module import AcousticAdapterBase -class AcousticForwardModelTestAdapter(AcousticForwardModelBaseAdapter): +class AcousticTestAdapter(AcousticAdapterBase): def forward_model(self, device) -> np.ndarray: diff --git a/simpa/core/simulation_modules/acoustic_forward_module/acoustic_forward_module_k_wave_adapter.py b/simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py similarity index 96% rename from simpa/core/simulation_modules/acoustic_forward_module/acoustic_forward_module_k_wave_adapter.py rename to simpa/core/simulation_modules/acoustic_module/k_wave_adapter.py index 17114651..9e4db45e 100644 --- a/simpa/core/simulation_modules/acoustic_forward_module/acoustic_forward_module_k_wave_adapter.py +++ b/simpa/core/simulation_modules/acoustic_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 @@ -13,8 +13,8 @@ from simpa.core.device_digital_twins import (CurvedArrayDetectionGeometry, DetectionGeometryBase) -from simpa.core.simulation_modules.acoustic_forward_module import \ - AcousticForwardModelBaseAdapter +from simpa.core.simulation_modules.acoustic_module import \ + AcousticAdapterBase from simpa.io_handling.io_hdf5 import load_data_field, save_hdf5 from simpa.utils import Tags from simpa.utils.matlab import generate_matlab_cmd @@ -24,9 +24,9 @@ from simpa.utils.settings import Settings -class KWaveAdapter(AcousticForwardModelBaseAdapter): +class KWaveAdapter(AcousticAdapterBase): """ - The KwaveAcousticForwardModel adapter enables acoustic simulations to be run with the + The KwaveAdapter enables acoustic simulations to be run with the k-wave MATLAB toolbox. k-Wave is a free toolbox (http://www.k-wave.org/) developed by Bradley Treeby and Ben Cox (University College London) and Jiri Jaros (Brno University of Technology). @@ -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/acoustic_forward_module/simulate_2D.m b/simpa/core/simulation_modules/acoustic_module/simulate_2D.m similarity index 100% rename from simpa/core/simulation_modules/acoustic_forward_module/simulate_2D.m rename to simpa/core/simulation_modules/acoustic_module/simulate_2D.m diff --git a/simpa/core/simulation_modules/acoustic_forward_module/simulate_3D.m b/simpa/core/simulation_modules/acoustic_module/simulate_3D.m similarity index 100% rename from simpa/core/simulation_modules/acoustic_forward_module/simulate_3D.m rename to simpa/core/simulation_modules/acoustic_module/simulate_3D.m diff --git a/simpa/core/simulation_modules/optical_module/__init__.py b/simpa/core/simulation_modules/optical_module/__init__.py new file mode 100644 index 00000000..7666c9bf --- /dev/null +++ b/simpa/core/simulation_modules/optical_module/__init__.py @@ -0,0 +1,4 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT +from .optical_adapter_base import OpticalAdapterBase diff --git a/simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_mcx_adapter.py b/simpa/core/simulation_modules/optical_module/mcx_adapter.py similarity index 98% rename from simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_mcx_adapter.py rename to simpa/core/simulation_modules/optical_module/mcx_adapter.py index b2bc652e..551787ec 100644 --- a/simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_mcx_adapter.py +++ b/simpa/core/simulation_modules/optical_module/mcx_adapter.py @@ -5,7 +5,7 @@ import numpy as np import subprocess from simpa.utils import Tags, Settings -from simpa.core.simulation_modules.optical_simulation_module import OpticalForwardModuleBase +from simpa.core.simulation_modules.optical_module import OpticalAdapterBase from simpa.core.device_digital_twins.illumination_geometries import IlluminationGeometryBase import json import jdata @@ -13,7 +13,7 @@ from typing import List, Dict, Tuple -class MCXAdapter(OpticalForwardModuleBase): +class MCXAdapter(OpticalAdapterBase): """ This class implements a bridge to the mcx framework to integrate mcx into SIMPA. This adapter only allows for computation of fluence, for computations of diffuse reflectance, take a look at `simpa.ReflectanceMcxAdapter` diff --git a/simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_mcx_reflectance_adapter.py b/simpa/core/simulation_modules/optical_module/mcx_reflectance_adapter.py similarity index 98% rename from simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_mcx_reflectance_adapter.py rename to simpa/core/simulation_modules/optical_module/mcx_reflectance_adapter.py index 73c92b7a..4473c965 100644 --- a/simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_mcx_reflectance_adapter.py +++ b/simpa/core/simulation_modules/optical_module/mcx_reflectance_adapter.py @@ -8,11 +8,11 @@ from typing import List, Tuple, Dict, Union from simpa.utils import Tags, Settings -from simpa.core.simulation_modules.optical_simulation_module.optical_forward_model_mcx_adapter import MCXAdapter +from simpa.core.simulation_modules.optical_module.mcx_adapter import MCXAdapter from simpa.core.device_digital_twins import IlluminationGeometryBase, PhotoacousticDevice -class MCXAdapterReflectance(MCXAdapter): +class MCXReflectanceAdapter(MCXAdapter): """ This class implements a bridge to the mcx framework to integrate mcx into SIMPA. This class targets specifically diffuse reflectance simulations. Specifically, it implements the capability to run diffuse reflectance simulations. @@ -35,7 +35,7 @@ def __init__(self, global_settings: Settings): :param global_settings: global settings used during simulations """ - super(MCXAdapterReflectance, self).__init__(global_settings=global_settings) + super(MCXReflectanceAdapter, self).__init__(global_settings=global_settings) self.mcx_photon_data_file = None self.padded = None self.mcx_output_suffixes = {'mcx_volumetric_data_file': '.jnii', diff --git a/simpa/core/simulation_modules/optical_simulation_module/__init__.py b/simpa/core/simulation_modules/optical_module/optical_adapter_base.py similarity index 97% rename from simpa/core/simulation_modules/optical_simulation_module/__init__.py rename to simpa/core/simulation_modules/optical_module/optical_adapter_base.py index 9e2a3b6a..fe7e3e78 100644 --- a/simpa/core/simulation_modules/optical_simulation_module/__init__.py +++ b/simpa/core/simulation_modules/optical_module/optical_adapter_base.py @@ -6,7 +6,7 @@ import numpy as np -from simpa.core.simulation_modules import SimulationModule +from simpa.core.simulation_modules import SimulationModuleBase from simpa.core.device_digital_twins import (IlluminationGeometryBase, PhotoacousticDevice) from simpa.io_handling.io_hdf5 import load_data_field, save_hdf5 @@ -16,7 +16,7 @@ assert_array_well_defined -class OpticalForwardModuleBase(SimulationModule): +class OpticalAdapterBase(SimulationModuleBase): """ Use this class as a base for implementations of optical forward models. This class has the attributes `self.temporary_output_files` which stores file paths that are temporarily created as @@ -24,7 +24,7 @@ class OpticalForwardModuleBase(SimulationModule): """ def __init__(self, global_settings: Settings): - super(OpticalForwardModuleBase, self).__init__(global_settings=global_settings) + super(OpticalAdapterBase, self).__init__(global_settings=global_settings) self.nx = None self.ny = None self.nz = None diff --git a/simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_test_adapter.py b/simpa/core/simulation_modules/optical_module/optical_test_adapter.py similarity index 74% rename from simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_test_adapter.py rename to simpa/core/simulation_modules/optical_module/optical_test_adapter.py index fb891954..bb1a7b49 100644 --- a/simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_test_adapter.py +++ b/simpa/core/simulation_modules/optical_module/optical_test_adapter.py @@ -2,11 +2,11 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT -from simpa.core.simulation_modules.optical_simulation_module import OpticalForwardModuleBase +from simpa.core.simulation_modules.optical_module import OpticalAdapterBase from simpa import Tags -class OpticalForwardModelTestAdapter(OpticalForwardModuleBase): +class OpticalTestAdapter(OpticalAdapterBase): """ This Adapter was created for testing purposes and only """ diff --git a/simpa/core/simulation_modules/reconstruction_module/__init__.py b/simpa/core/simulation_modules/reconstruction_module/__init__.py index 05f8b3ef..7c084f66 100644 --- a/simpa/core/simulation_modules/reconstruction_module/__init__.py +++ b/simpa/core/simulation_modules/reconstruction_module/__init__.py @@ -1,98 +1,9 @@ # SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT +from .reconstruction_adapter_base import ReconstructionAdapterBase -from simpa.utils import Tags -from simpa.core.device_digital_twins import DetectionGeometryBase -from simpa.core.device_digital_twins import PhotoacousticDevice -from simpa.io_handling.io_hdf5 import load_data_field -from abc import abstractmethod -from simpa.core.simulation_modules import SimulationModule -from simpa.utils.dict_path_manager import generate_dict_path -from simpa.io_handling.io_hdf5 import save_hdf5 -import numpy as np -from simpa.utils import Settings -from simpa.core.simulation_modules.reconstruction_module.reconstruction_utils import bandpass_filter_with_settings, apply_b_mode -from simpa.utils.quality_assurance.data_sanity_testing import assert_array_well_defined - - -class ReconstructionAdapterBase(SimulationModule): - """ - This class is the main entry point to perform image reconstruction using the SIMPA toolkit. - All information necessary for the respective reconstruction method must be contained in the - respective settings dictionary. - """ - - def __init__(self, global_settings: Settings): - super(ReconstructionAdapterBase, self).__init__(global_settings=global_settings) - - def load_component_settings(self) -> Settings: - """Implements abstract method to serve reconstruction settings as component settings - - :return: Settings: reconstruction component settings - """ - return self.global_settings.get_reconstruction_settings() - - @abstractmethod - def reconstruction_algorithm(self, time_series_sensor_data, - detection_geometry: DetectionGeometryBase) -> np.ndarray: - """ - A deriving class needs to implement this method according to its model. - - :param time_series_sensor_data: the time series sensor data - :param detection_geometry: - :return: a reconstructed photoacoustic image - """ - pass - - def run(self, device): - self.logger.info("Performing reconstruction...") - - time_series_sensor_data = load_data_field(self.global_settings[Tags.SIMPA_OUTPUT_PATH], - Tags.DATA_FIELD_TIME_SERIES_DATA, self.global_settings[Tags.WAVELENGTH]) - - _device = None - if isinstance(device, DetectionGeometryBase): - _device = device - elif isinstance(device, PhotoacousticDevice): - _device = device.get_detection_geometry() - else: - raise TypeError(f"Type {type(device)} is not supported for performing image reconstruction.") - - if Tags.RECONSTRUCTION_PERFORM_BANDPASS_FILTERING in self.component_settings and \ - self.component_settings[Tags.RECONSTRUCTION_PERFORM_BANDPASS_FILTERING]: - - time_series_sensor_data = bandpass_filter_with_settings(time_series_sensor_data, - self.global_settings, - self.component_settings, - _device) - - # check for B-mode methods and perform envelope detection on time series data if specified - if Tags.RECONSTRUCTION_BMODE_BEFORE_RECONSTRUCTION in self.component_settings \ - and self.component_settings[Tags.RECONSTRUCTION_BMODE_BEFORE_RECONSTRUCTION] \ - and Tags.RECONSTRUCTION_BMODE_METHOD in self.component_settings: - time_series_sensor_data = apply_b_mode( - time_series_sensor_data, method=self.component_settings[Tags.RECONSTRUCTION_BMODE_METHOD]) - - reconstruction = self.reconstruction_algorithm(time_series_sensor_data, _device) - - # check for B-mode methods and perform envelope detection on time series data if specified - if Tags.RECONSTRUCTION_BMODE_AFTER_RECONSTRUCTION in self.component_settings \ - and self.component_settings[Tags.RECONSTRUCTION_BMODE_AFTER_RECONSTRUCTION] \ - and Tags.RECONSTRUCTION_BMODE_METHOD in self.component_settings: - reconstruction = apply_b_mode( - reconstruction, method=self.component_settings[Tags.RECONSTRUCTION_BMODE_METHOD]) - - if not (Tags.IGNORE_QA_ASSERTIONS in self.global_settings and Tags.IGNORE_QA_ASSERTIONS): - assert_array_well_defined(reconstruction, array_name="reconstruction") - - reconstruction_output_path = generate_dict_path( - Tags.DATA_FIELD_RECONSTRUCTED_DATA, self.global_settings[Tags.WAVELENGTH]) - - save_hdf5(reconstruction, self.global_settings[Tags.SIMPA_OUTPUT_PATH], - reconstruction_output_path) - - self.logger.info("Performing reconstruction...[Done]") +from simpa.utils import Tags, Settings def create_reconstruction_settings(speed_of_sound_in_m_per_s: int = 1540, time_spacing_in_s: float = 2.5e-8, diff --git a/simpa/core/simulation_modules/reconstruction_module/reconstruction_module_delay_and_sum_adapter.py b/simpa/core/simulation_modules/reconstruction_module/delay_and_sum_adapter.py similarity index 100% rename from simpa/core/simulation_modules/reconstruction_module/reconstruction_module_delay_and_sum_adapter.py rename to simpa/core/simulation_modules/reconstruction_module/delay_and_sum_adapter.py diff --git a/simpa/core/simulation_modules/reconstruction_module/reconstruction_module_delay_multiply_and_sum_adapter.py b/simpa/core/simulation_modules/reconstruction_module/delay_multiply_and_sum_adapter.py similarity index 100% rename from simpa/core/simulation_modules/reconstruction_module/reconstruction_module_delay_multiply_and_sum_adapter.py rename to simpa/core/simulation_modules/reconstruction_module/delay_multiply_and_sum_adapter.py diff --git a/simpa/core/simulation_modules/reconstruction_module/reconstruction_adapter_base.py b/simpa/core/simulation_modules/reconstruction_module/reconstruction_adapter_base.py new file mode 100644 index 00000000..90d17099 --- /dev/null +++ b/simpa/core/simulation_modules/reconstruction_module/reconstruction_adapter_base.py @@ -0,0 +1,95 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT + +from simpa.utils import Tags +from simpa.core.device_digital_twins import DetectionGeometryBase +from simpa.core.device_digital_twins import PhotoacousticDevice +from simpa.io_handling.io_hdf5 import load_data_field +from abc import abstractmethod +from simpa.core.simulation_modules import SimulationModuleBase +from simpa.utils.dict_path_manager import generate_dict_path +from simpa.io_handling.io_hdf5 import save_hdf5 +import numpy as np +from simpa.utils import Settings +from simpa.core.simulation_modules.reconstruction_module.reconstruction_utils import bandpass_filter_with_settings, apply_b_mode +from simpa.utils.quality_assurance.data_sanity_testing import assert_array_well_defined + + +class ReconstructionAdapterBase(SimulationModuleBase): + """ + This class is the main entry point to perform image reconstruction using the SIMPA toolkit. + All information necessary for the respective reconstruction method must be contained in the + respective settings dictionary. + """ + + def __init__(self, global_settings: Settings): + super(ReconstructionAdapterBase, self).__init__(global_settings=global_settings) + + def load_component_settings(self) -> Settings: + """Implements abstract method to serve reconstruction settings as component settings + + :return: Settings: reconstruction component settings + """ + return self.global_settings.get_reconstruction_settings() + + @abstractmethod + def reconstruction_algorithm(self, time_series_sensor_data, + detection_geometry: DetectionGeometryBase) -> np.ndarray: + """ + A deriving class needs to implement this method according to its model. + + :param time_series_sensor_data: the time series sensor data + :param detection_geometry: + :return: a reconstructed photoacoustic image + """ + pass + + def run(self, device): + self.logger.info("Performing reconstruction...") + + time_series_sensor_data = load_data_field(self.global_settings[Tags.SIMPA_OUTPUT_PATH], + Tags.DATA_FIELD_TIME_SERIES_DATA, self.global_settings[Tags.WAVELENGTH]) + + _device = None + if isinstance(device, DetectionGeometryBase): + _device = device + elif isinstance(device, PhotoacousticDevice): + _device = device.get_detection_geometry() + else: + raise TypeError(f"Type {type(device)} is not supported for performing image reconstruction.") + + if Tags.RECONSTRUCTION_PERFORM_BANDPASS_FILTERING in self.component_settings and \ + self.component_settings[Tags.RECONSTRUCTION_PERFORM_BANDPASS_FILTERING]: + + time_series_sensor_data = bandpass_filter_with_settings(time_series_sensor_data, + self.global_settings, + self.component_settings, + _device) + + # check for B-mode methods and perform envelope detection on time series data if specified + if Tags.RECONSTRUCTION_BMODE_BEFORE_RECONSTRUCTION in self.component_settings \ + and self.component_settings[Tags.RECONSTRUCTION_BMODE_BEFORE_RECONSTRUCTION] \ + and Tags.RECONSTRUCTION_BMODE_METHOD in self.component_settings: + time_series_sensor_data = apply_b_mode( + time_series_sensor_data, method=self.component_settings[Tags.RECONSTRUCTION_BMODE_METHOD]) + + reconstruction = self.reconstruction_algorithm(time_series_sensor_data, _device) + + # check for B-mode methods and perform envelope detection on time series data if specified + if Tags.RECONSTRUCTION_BMODE_AFTER_RECONSTRUCTION in self.component_settings \ + and self.component_settings[Tags.RECONSTRUCTION_BMODE_AFTER_RECONSTRUCTION] \ + and Tags.RECONSTRUCTION_BMODE_METHOD in self.component_settings: + reconstruction = apply_b_mode( + reconstruction, method=self.component_settings[Tags.RECONSTRUCTION_BMODE_METHOD]) + + if not (Tags.IGNORE_QA_ASSERTIONS in self.global_settings and Tags.IGNORE_QA_ASSERTIONS): + assert_array_well_defined(reconstruction, array_name="reconstruction") + + reconstruction_output_path = generate_dict_path( + Tags.DATA_FIELD_RECONSTRUCTED_DATA, self.global_settings[Tags.WAVELENGTH]) + + save_hdf5(reconstruction, self.global_settings[Tags.SIMPA_OUTPUT_PATH], + reconstruction_output_path) + + self.logger.info("Performing reconstruction...[Done]") diff --git a/simpa/core/simulation_modules/reconstruction_module/reconstruction_module_test_adapter.py b/simpa/core/simulation_modules/reconstruction_module/reconstruction_test_adapter.py similarity index 85% rename from simpa/core/simulation_modules/reconstruction_module/reconstruction_module_test_adapter.py rename to simpa/core/simulation_modules/reconstruction_module/reconstruction_test_adapter.py index e3456c4e..4e2a7076 100644 --- a/simpa/core/simulation_modules/reconstruction_module/reconstruction_module_test_adapter.py +++ b/simpa/core/simulation_modules/reconstruction_module/reconstruction_test_adapter.py @@ -5,7 +5,7 @@ from simpa.core.simulation_modules.reconstruction_module import ReconstructionAdapterBase -class ReconstructionModuleTestAdapter(ReconstructionAdapterBase): +class ReconstructionTestAdapter(ReconstructionAdapterBase): def reconstruction_algorithm(self, time_series_sensor_data, detection_geometry): return time_series_sensor_data / 10 + 5 diff --git a/simpa/core/simulation_modules/reconstruction_module/reconstruction_module_signed_delay_multiply_and_sum_adapter.py b/simpa/core/simulation_modules/reconstruction_module/signed_delay_multiply_and_sum_adapter.py similarity index 100% rename from simpa/core/simulation_modules/reconstruction_module/reconstruction_module_signed_delay_multiply_and_sum_adapter.py rename to simpa/core/simulation_modules/reconstruction_module/signed_delay_multiply_and_sum_adapter.py diff --git a/simpa/core/simulation_modules/reconstruction_module/reconstruction_module_time_reversal_adapter.py b/simpa/core/simulation_modules/reconstruction_module/time_reversal_adapter.py similarity index 98% rename from simpa/core/simulation_modules/reconstruction_module/reconstruction_module_time_reversal_adapter.py rename to simpa/core/simulation_modules/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/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/simulation_module_base.py b/simpa/core/simulation_modules/simulation_module_base.py new file mode 100644 index 00000000..1ac50a33 --- /dev/null +++ b/simpa/core/simulation_modules/simulation_module_base.py @@ -0,0 +1,44 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT + +from abc import abstractmethod + +from simpa.core import PipelineElementBase +from simpa.utils import Settings, Tags +from typing import List + + +class SimulationModuleBase(PipelineElementBase): + """ + Defines a simulation module that is a step in the simulation pipeline. + Each simulation module can only be one of Volume Creation, Light Propagation Modeling, Acoustic Wave Propagation Modeling, Image Reconstruction. + """ + + def __init__(self, global_settings: Settings): + """ + :param global_settings: The SIMPA settings dictionary + :type global_settings: Settings + """ + super(SimulationModuleBase, self).__init__(global_settings=global_settings) + self.component_settings = self.load_component_settings() + if self.component_settings is None: + raise ValueError("The component settings should not be None at this point") + + @abstractmethod + def load_component_settings(self) -> Settings: + """ + :return: Loads component settings corresponding to this simulation component + """ + pass + + def get_additional_flags(self) -> List[str]: + """Reads the list of additional flags from the corresponding component settings Tags.ADDITIONAL_FLAGS + + :return: List[str]: list of additional flags + """ + cmd = [] + 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 diff --git a/simpa/core/simulation_modules/volume_creation_module/__init__.py b/simpa/core/simulation_modules/volume_creation_module/__init__.py index 6b9f0753..c39607ef 100644 --- a/simpa/core/simulation_modules/volume_creation_module/__init__.py +++ b/simpa/core/simulation_modules/volume_creation_module/__init__.py @@ -2,77 +2,4 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT -from abc import abstractmethod -from simpa.utils.settings import Settings -from simpa.utils import Tags -from simpa.utils.constants import wavelength_independent_properties, property_tags -import torch -from simpa.core.simulation_modules import SimulationModule -from simpa.io_handling import save_data_field -from simpa.utils.quality_assurance.data_sanity_testing import assert_equal_shapes, assert_array_well_defined - - -class VolumeCreatorModuleBase(SimulationModule): - """ - Use this class to define your own volume creation adapter. - - """ - - def __init__(self, global_settings: Settings): - super(VolumeCreatorModuleBase, self).__init__(global_settings=global_settings) - - def load_component_settings(self) -> Settings: - """Implements abstract method to serve volume creation settings as component settings - - :return: Settings: volume creation component settings - """ - return self.global_settings.get_volume_creation_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)) - sizes = (volume_x_dim, volume_y_dim, volume_z_dim) - - wavelength = self.global_settings[Tags.WAVELENGTH] - first_wavelength = self.global_settings[Tags.WAVELENGTHS][0] - - for key in property_tags: - # Create wavelength-independent properties only in the first wavelength run - if key in wavelength_independent_properties and wavelength != first_wavelength: - continue - volumes[key] = torch.zeros(sizes, dtype=torch.float, device=self.torch_device) - - return volumes, volume_x_dim, volume_y_dim, volume_z_dim - - @abstractmethod - def create_simulation_volume(self) -> dict: - """ - This method creates an in silico representation of a tissue as described in the settings file that is given. - - :return: A dictionary containing optical and acoustic properties as well as other characteristics of the - simulated volume such as oxygenation, and a segmentation mask. All of these are given as 3d numpy arrays. - :rtype: dict - """ - pass - - def run(self, device): - self.logger.info("VOLUME CREATION") - - volumes = self.create_simulation_volume() - # explicitly empty cache to free reserved GPU memory after volume creation - torch.cuda.empty_cache() - - if not (Tags.IGNORE_QA_ASSERTIONS in self.global_settings and Tags.IGNORE_QA_ASSERTIONS): - assert_equal_shapes(list(volumes.values())) - for _volume_name in volumes.keys(): - if _volume_name == Tags.DATA_FIELD_OXYGENATION: - # oxygenation can have NaN by definition - continue - assert_array_well_defined(volumes[_volume_name], array_name=_volume_name) - - for key, value in volumes.items(): - save_data_field(value, self.global_settings[Tags.SIMPA_OUTPUT_PATH], - data_field=key, wavelength=self.global_settings[Tags.WAVELENGTH]) +from .volume_creation_adapter_base import VolumeCreationAdapterBase 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/model_based_adapter.py similarity index 83% rename from simpa/core/simulation_modules/volume_creation_module/volume_creation_module_model_based_adapter.py rename to simpa/core/simulation_modules/volume_creation_module/model_based_adapter.py index 2379adf3..e4aad5bb 100644 --- a/simpa/core/simulation_modules/volume_creation_module/volume_creation_module_model_based_adapter.py +++ b/simpa/core/simulation_modules/volume_creation_module/model_based_adapter.py @@ -2,7 +2,7 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT -from simpa.core.simulation_modules.volume_creation_module import VolumeCreatorModuleBase +from simpa.core.simulation_modules.volume_creation_module import VolumeCreationAdapterBase from simpa.utils.libraries.structure_library import priority_sorted_structures from simpa.utils import Tags import numpy as np @@ -10,7 +10,7 @@ import torch -class ModelBasedVolumeCreationAdapter(VolumeCreatorModuleBase): +class ModelBasedAdapter(VolumeCreationAdapterBase): """ The model-based volume creator uses a set of rules how to generate structures to create a simulation volume. @@ -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/segmentation_based_adapter.py similarity index 63% rename from simpa/core/simulation_modules/volume_creation_module/volume_creation_module_segmentation_based_adapter.py rename to simpa/core/simulation_modules/volume_creation_module/segmentation_based_adapter.py index 490ce2eb..8800f484 100644 --- a/simpa/core/simulation_modules/volume_creation_module/volume_creation_module_segmentation_based_adapter.py +++ b/simpa/core/simulation_modules/volume_creation_module/segmentation_based_adapter.py @@ -2,7 +2,7 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT -from simpa.core.simulation_modules.volume_creation_module import VolumeCreatorModuleBase +from simpa.core.simulation_modules.volume_creation_module import VolumeCreationAdapterBase from simpa.utils import Tags from simpa.utils.constants import property_tags from simpa.io_handling import save_hdf5 @@ -10,9 +10,9 @@ import torch -class SegmentationBasedVolumeCreationAdapter(VolumeCreatorModuleBase): +class SegmentationBasedAdapter(VolumeCreationAdapterBase): """ - 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/core/simulation_modules/volume_creation_module/volume_creation_adapter_base.py b/simpa/core/simulation_modules/volume_creation_module/volume_creation_adapter_base.py new file mode 100644 index 00000000..af3f47d6 --- /dev/null +++ b/simpa/core/simulation_modules/volume_creation_module/volume_creation_adapter_base.py @@ -0,0 +1,78 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT + +from abc import abstractmethod +from simpa.utils.settings import Settings +from simpa.utils import Tags +from simpa.utils.constants import wavelength_independent_properties, property_tags +import torch +from simpa.core.simulation_modules import SimulationModuleBase +from simpa.io_handling import save_data_field +from simpa.utils.quality_assurance.data_sanity_testing import assert_equal_shapes, assert_array_well_defined + + +class VolumeCreationAdapterBase(SimulationModuleBase): + """ + Use this class to define your own volume creation adapter. + + """ + + def __init__(self, global_settings: Settings): + super(VolumeCreationAdapterBase, self).__init__(global_settings=global_settings) + + def load_component_settings(self) -> Settings: + """Implements abstract method to serve volume creation settings as component settings + + :return: Settings: volume creation component settings + """ + return self.global_settings.get_volume_creation_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)) + sizes = (volume_x_dim, volume_y_dim, volume_z_dim) + + wavelength = self.global_settings[Tags.WAVELENGTH] + first_wavelength = self.global_settings[Tags.WAVELENGTHS][0] + + for key in property_tags: + # Create wavelength-independent properties only in the first wavelength run + if key in wavelength_independent_properties and wavelength != first_wavelength: + continue + volumes[key] = torch.zeros(sizes, dtype=torch.float, device=self.torch_device) + + return volumes, volume_x_dim, volume_y_dim, volume_z_dim + + @abstractmethod + def create_simulation_volume(self) -> dict: + """ + This method creates an in silico representation of a tissue as described in the settings file that is given. + + :return: A dictionary containing optical and acoustic properties as well as other characteristics of the + simulated volume such as oxygenation, and a segmentation mask. All of these are given as 3d numpy arrays. + :rtype: dict + """ + pass + + def run(self, device): + self.logger.info("VOLUME CREATION") + + volumes = self.create_simulation_volume() + # explicitly empty cache to free reserved GPU memory after volume creation + torch.cuda.empty_cache() + + if not (Tags.IGNORE_QA_ASSERTIONS in self.global_settings and Tags.IGNORE_QA_ASSERTIONS): + assert_equal_shapes(list(volumes.values())) + for _volume_name in volumes.keys(): + if _volume_name == Tags.DATA_FIELD_OXYGENATION: + # oxygenation can have NaN by definition + continue + assert_array_well_defined(volumes[_volume_name], array_name=_volume_name) + + for key, value in volumes.items(): + save_data_field(value, self.global_settings[Tags.SIMPA_OUTPUT_PATH], + data_field=key, wavelength=self.global_settings[Tags.WAVELENGTH]) 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/EllipticalTubularStructure.py b/simpa/utils/libraries/structure_library/EllipticalTubularStructure.py index 6423ef06..b6b1cf4d 100644 --- a/simpa/utils/libraries/structure_library/EllipticalTubularStructure.py +++ b/simpa/utils/libraries/structure_library/EllipticalTubularStructure.py @@ -94,7 +94,7 @@ def get_enclosed_indices(self): main_axis_vector = main_axis_vector/torch.linalg.norm(main_axis_vector) * main_axis_length minor_axis_length = main_axis_length*torch.sqrt(1-eccentricity**2) - minor_axis_vector = torch.cross(cylinder_vector, main_axis_vector) + minor_axis_vector = torch.linalg.cross(cylinder_vector, main_axis_vector) minor_axis_vector = minor_axis_vector / torch.linalg.norm(minor_axis_vector) * minor_axis_length dot_product = torch.matmul(target_vector, cylinder_vector)/torch.linalg.norm(cylinder_vector) 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 31b652d3..395907ec 100644 --- a/simpa/utils/tags.py +++ b/simpa/utils/tags.py @@ -114,9 +114,9 @@ class Tags: Usage: module volume_creation_module, naming convention """ - VOLUME_CREATOR_SEGMENTATION_BASED = "volume_creator_segmentation_based" + VOLUME_CREATOR_SEGMENTATION_BASED = "segmentation_based_adapter" """ - Corresponds to the SegmentationBasedVolumeCreator.\n + Corresponds to the SegmentationBasedAdapter.\n Usage: module volume_creation_module, naming convention """ @@ -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 @@ -1541,3 +1610,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..f933f2de 100644 --- a/simpa_examples/benchmarking/performance_check.py +++ b/simpa_examples/benchmarking/performance_check.py @@ -28,16 +28,13 @@ def run_benchmarking_tests(spacing=0.4, profile: str = "TIME", savefolder: str = import simpa_examples - 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, + examples = [simpa_examples.run_minimal_optical_simulation, + simpa_examples.run_minimal_optical_simulation_uniform_cube, simpa_examples.run_optical_and_acoustic_simulation, - 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..0659d932 100644 --- a/simpa_examples/benchmarking/run_benchmarking.sh +++ b/simpa_examples/benchmarking/run_benchmarking.sh @@ -1,5 +1,7 @@ #!/bin/bash +set -e + help() { echo "Usage: calculate benchmarking for [options]" echo "For further details see readme" @@ -7,8 +9,8 @@ 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 " -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" @@ -61,11 +63,11 @@ if [ "$start" == 0 ]; then 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.py b/simpa_examples/minimal_optical_simulation.py index c3646f22..6f43eeb3 100644 --- a/simpa_examples/minimal_optical_simulation.py +++ b/simpa_examples/minimal_optical_simulation.py @@ -127,14 +127,14 @@ def create_example_tissue(): if not SAVE_REFLECTANCE and not SAVE_PHOTON_DIRECTION: pipeline = [ - sp.ModelBasedVolumeCreationAdapter(settings), + sp.ModelBasedAdapter(settings), sp.MCXAdapter(settings), sp.GaussianNoise(settings, "noise_model_1") ] else: pipeline = [ - sp.ModelBasedVolumeCreationAdapter(settings), - sp.MCXAdapterReflectance(settings), + sp.ModelBasedAdapter(settings), + sp.MCXReflectanceAdapter(settings), ] class ExampleDeviceSlitIlluminationLinearDetector(sp.PhotoacousticDevice): 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_examples/minimal_optical_simulation_uniform_cube.py b/simpa_examples/minimal_optical_simulation_uniform_cube.py index 675addd5..2e08796c 100644 --- a/simpa_examples/minimal_optical_simulation_uniform_cube.py +++ b/simpa_examples/minimal_optical_simulation_uniform_cube.py @@ -88,8 +88,8 @@ def create_example_tissue(): }) pipeline = [ - sp.ModelBasedVolumeCreationAdapter(settings), - sp.MCXAdapterReflectance(settings), + sp.ModelBasedAdapter(settings), + sp.MCXReflectanceAdapter(settings), ] device = sp.PencilBeamIlluminationGeometry(device_position_mm=np.asarray([VOLUME_TRANSDUCER_DIM_IN_MM/2, diff --git a/simpa_examples/msot_invision_simulation.py b/simpa_examples/msot_invision_simulation.py index 210daf4b..d7eb2cb6 100644 --- a/simpa_examples/msot_invision_simulation.py +++ b/simpa_examples/msot_invision_simulation.py @@ -28,7 +28,7 @@ def run_msot_invision_simulation(spacing: float | int = 0.5, path_manager=None, def create_pipeline(_settings: sp.Settings): return [ - sp.ModelBasedVolumeCreationAdapter(settings), + sp.ModelBasedAdapter(settings), sp.MCXAdapter(settings), sp.KWaveAdapter(settings), sp.FieldOfViewCropping(settings), @@ -152,7 +152,7 @@ def get_settings(): } acoustic_settings = { - Tags.ACOUSTIC_SIMULATION_3D: True, + Tags.ACOUSTIC_SIMULATION_3D: False, Tags.ACOUSTIC_MODEL_BINARY_PATH: path_manager.get_matlab_binary_path(), Tags.KWAVE_PROPERTY_ALPHA_POWER: 0.00, Tags.KWAVE_PROPERTY_SENSOR_RECORD: "p", diff --git a/simpa_examples/optical_and_acoustic_simulation.py b/simpa_examples/optical_and_acoustic_simulation.py index c92567f4..986823be 100644 --- a/simpa_examples/optical_and_acoustic_simulation.py +++ b/simpa_examples/optical_and_acoustic_simulation.py @@ -187,11 +187,10 @@ def create_example_tissue(): pitch_mm=0.25, number_detector_elements=100, field_of_view_extent_mm=np.asarray([-15, 15, 0, 0, 0, 20]))) - print(device.get_detection_geometry().get_detector_element_positions_base_mm()) device.add_illumination_geometry(sp.SlitIlluminationGeometry(slit_vector_mm=[100, 0, 0])) SIMULATION_PIPELINE = [ - sp.ModelBasedVolumeCreationAdapter(settings), + sp.ModelBasedAdapter(settings), sp.MCXAdapter(settings), sp.GaussianNoise(settings, "noise_initial_pressure"), sp.KWaveAdapter(settings), diff --git a/simpa_examples/perform_iterative_qPAI_reconstruction.py b/simpa_examples/perform_iterative_qPAI_reconstruction.py index f46ccfb7..784e741b 100644 --- a/simpa_examples/perform_iterative_qPAI_reconstruction.py +++ b/simpa_examples/perform_iterative_qPAI_reconstruction.py @@ -142,7 +142,7 @@ def create_example_tissue(): # run pipeline including iterative qPAI method pipeline = [ - sp.ModelBasedVolumeCreationAdapter(settings), + sp.ModelBasedAdapter(settings), sp.MCXAdapter(settings), sp.GaussianNoise(settings, "noise_model"), sp.IterativeqPAI(settings, "iterative_qpai_reconstruction") diff --git a/simpa_examples/segmentation_loader.py b/simpa_examples/segmentation_loader.py index 8114ce0f..685240a8 100644 --- a/simpa_examples/segmentation_loader.py +++ b/simpa_examples/segmentation_loader.py @@ -7,6 +7,7 @@ import numpy as np from skimage.data import shepp_logan_phantom from scipy.ndimage import zoom +from skimage.transform import resize # FIXME temporary workaround for newest Intel architectures import os @@ -19,11 +20,12 @@ @profile -def run_segmentation_loader(spacing: float | int = .1, path_manager=None, +def run_segmentation_loader(spacing: float | int = 1.0, input_spacing: float | int = 0.2, path_manager=None, visualise: bool = True): """ - :param spacing: The simulation spacing between voxels + :param spacing: The simulation spacing between voxels in mm + :param input_spacing: The input spacing between voxels in mm :param path_manager: the path manager to be used, typically sp.PathManager :param visualise: If VISUALIZE is set to True, the reconstruction result will be plotted :return: a run through of the example @@ -34,10 +36,9 @@ def run_segmentation_loader(spacing: float | int = .1, path_manager=None, label_mask = shepp_logan_phantom() label_mask = np.digitize(label_mask, bins=np.linspace(0.0, 1.0, 11), right=True) + label_mask = label_mask[100:300, 100:300] + label_mask = np.reshape(label_mask, (label_mask.shape[0], 1, label_mask.shape[1])) - label_mask = np.reshape(label_mask, (400, 1, 400)) - - input_spacing = 0.2 segmentation_volume_tiled = np.tile(label_mask, (1, 128, 1)) segmentation_volume_mask = np.round(zoom(segmentation_volume_tiled, input_spacing/spacing, order=0)).astype(int) @@ -67,9 +68,9 @@ def segmentation_class_mapping(): settings[Tags.RANDOM_SEED] = 1234 settings[Tags.WAVELENGTHS] = [700] settings[Tags.SPACING_MM] = spacing - settings[Tags.DIM_VOLUME_X_MM] = 400 / (spacing / input_spacing) - settings[Tags.DIM_VOLUME_Y_MM] = 128 / (spacing / input_spacing) - settings[Tags.DIM_VOLUME_Z_MM] = 400 / (spacing / input_spacing) + settings[Tags.DIM_VOLUME_X_MM] = segmentation_volume_mask.shape[0] * spacing + settings[Tags.DIM_VOLUME_Y_MM] = segmentation_volume_mask.shape[1] * spacing + settings[Tags.DIM_VOLUME_Z_MM] = segmentation_volume_mask.shape[2] * spacing settings.set_volume_creation_settings({ Tags.INPUT_SEGMENTATION_VOLUME: segmentation_volume_mask, @@ -85,7 +86,7 @@ def segmentation_class_mapping(): }) pipeline = [ - sp.SegmentationBasedVolumeCreationAdapter(settings), + sp.SegmentationBasedAdapter(settings), sp.MCXAdapter(settings) ] @@ -105,9 +106,11 @@ def segmentation_class_mapping(): if __name__ == "__main__": parser = ArgumentParser(description='Run the segmentation loader example') - parser.add_argument("--spacing", default=0.2, type=float, help='the voxel spacing in mm') + parser.add_argument("--spacing", default=1, type=float, help='the voxel spacing in mm') + parser.add_argument("--input_spacing", default=0.2, type=float, help='the input spacing in mm') parser.add_argument("--path_manager", default=None, help='the path manager, None uses sp.PathManager') parser.add_argument("--visualise", default=True, type=bool, help='whether to visualise the result') config = parser.parse_args() - run_segmentation_loader(spacing=config.spacing, path_manager=config.path_manager, visualise=config.visualise) + run_segmentation_loader(spacing=config.spacing, input_spacing=config.input_spacing, + path_manager=config.path_manager, visualise=config.visualise) diff --git a/simpa_tests/automatic_tests/test_IPASC_export.py b/simpa_tests/automatic_tests/test_IPASC_export.py index cc638ec7..29cbc91f 100644 --- a/simpa_tests/automatic_tests/test_IPASC_export.py +++ b/simpa_tests/automatic_tests/test_IPASC_export.py @@ -10,13 +10,13 @@ from simpa.core.device_digital_twins import RSOMExplorerP50 from simpa.utils import Tags, Settings from simpa_tests.test_utils import create_test_structure_parameters -from simpa import ModelBasedVolumeCreationAdapter -from simpa.core.simulation_modules.optical_simulation_module.optical_forward_model_test_adapter import \ - OpticalForwardModelTestAdapter -from simpa.core.simulation_modules.acoustic_forward_module.acoustic_forward_model_test_adapter import \ - AcousticForwardModelTestAdapter -from simpa.core.simulation_modules.reconstruction_module.reconstruction_module_test_adapter import \ - ReconstructionModuleTestAdapter +from simpa import ModelBasedAdapter +from simpa.core.simulation_modules.optical_module.optical_test_adapter import \ + OpticalTestAdapter +from simpa.core.simulation_modules.acoustic_module.acoustic_test_adapter import \ + AcousticTestAdapter +from simpa.core.simulation_modules.reconstruction_module.reconstruction_test_adapter import \ + ReconstructionTestAdapter from pacfish import load_data as load_ipasc from simpa.io_handling import load_hdf5 as load_simpa @@ -63,21 +63,21 @@ def setUp(self) -> None: }) self.acoustic_simulation_pipeline = [ - ModelBasedVolumeCreationAdapter(self.settings), - OpticalForwardModelTestAdapter(self.settings), - AcousticForwardModelTestAdapter(self.settings), + ModelBasedAdapter(self.settings), + OpticalTestAdapter(self.settings), + AcousticTestAdapter(self.settings), ] self.optical_simulation_pipeline = [ - ModelBasedVolumeCreationAdapter(self.settings), - OpticalForwardModelTestAdapter(self.settings) + ModelBasedAdapter(self.settings), + OpticalTestAdapter(self.settings) ] self.full_simulation_pipeline = [ - ModelBasedVolumeCreationAdapter(self.settings), - OpticalForwardModelTestAdapter(self.settings), - AcousticForwardModelTestAdapter(self.settings), - ReconstructionModuleTestAdapter(self.settings) + ModelBasedAdapter(self.settings), + OpticalTestAdapter(self.settings), + AcousticTestAdapter(self.settings), + ReconstructionTestAdapter(self.settings) ] self.device = RSOMExplorerP50(0.1, 12, 12) diff --git a/simpa_tests/automatic_tests/test_additional_flags.py b/simpa_tests/automatic_tests/test_additional_flags.py index ca1c4731..c8ad88b1 100644 --- a/simpa_tests/automatic_tests/test_additional_flags.py +++ b/simpa_tests/automatic_tests/test_additional_flags.py @@ -1,9 +1,10 @@ import unittest import numpy as np -from simpa import MCXAdapterReflectance, MCXAdapter, KWaveAdapter, TimeReversalAdapter, Tags, Settings +from simpa import MCXReflectanceAdapter, 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') @@ -14,7 +15,7 @@ def test_get_cmd_mcx_reflectance_adapter(self): Tags.OPTICAL_MODEL_BINARY_PATH: '.', Tags.ADDITIONAL_FLAGS: self.additional_flags }) - mcx_reflectance_adapter = MCXAdapterReflectance(global_settings=self.settings) + mcx_reflectance_adapter = MCXReflectanceAdapter(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") @@ -46,8 +47,7 @@ def test_get_cmd_time_reversal_adapter(self): 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") - - + if __name__ == '__main__': - unittest.main() \ No newline at end of file + unittest.main() diff --git a/simpa_tests/automatic_tests/test_create_a_volume.py b/simpa_tests/automatic_tests/test_create_a_volume.py index c67a8fc8..05652702 100644 --- a/simpa_tests/automatic_tests/test_create_a_volume.py +++ b/simpa_tests/automatic_tests/test_create_a_volume.py @@ -8,7 +8,7 @@ from simpa.core.simulation import simulate import os from simpa_tests.test_utils import create_test_structure_parameters -from simpa import ModelBasedVolumeCreationAdapter +from simpa import ModelBasedAdapter from simpa.core.device_digital_twins import RSOMExplorerP50 @@ -32,7 +32,7 @@ def test_create_volume(self): settings.set_volume_creation_settings({Tags.STRUCTURES: create_test_structure_parameters()}) simulation_pipeline = [ - ModelBasedVolumeCreationAdapter(settings) + ModelBasedAdapter(settings) ] simulate(simulation_pipeline, settings, RSOMExplorerP50(0.1, 1, 1)) 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/test_linear_unmixing.py b/simpa_tests/automatic_tests/test_linear_unmixing.py index fe61e885..dc6e955e 100644 --- a/simpa_tests/automatic_tests/test_linear_unmixing.py +++ b/simpa_tests/automatic_tests/test_linear_unmixing.py @@ -50,7 +50,7 @@ def setUp(self): self.device = sp.PencilBeamIlluminationGeometry() # Run empty pipeline simulation to "fill" hdf5 file following usual procedure - pipeline = [sp.ModelBasedVolumeCreationAdapter(self.settings)] + pipeline = [sp.ModelBasedAdapter(self.settings)] sp.simulate(pipeline, self.settings, self.device) def test(self): diff --git a/simpa_tests/automatic_tests/test_noise_models.py b/simpa_tests/automatic_tests/test_noise_models.py index b23ac50d..d4dd256f 100644 --- a/simpa_tests/automatic_tests/test_noise_models.py +++ b/simpa_tests/automatic_tests/test_noise_models.py @@ -13,7 +13,7 @@ from simpa.core.simulation import simulate from simpa.utils import TISSUE_LIBRARY from simpa.io_handling import load_data_field -from simpa import ModelBasedVolumeCreationAdapter +from simpa import ModelBasedAdapter class TestNoiseModels(unittest.TestCase): @@ -63,7 +63,7 @@ def validate_noise_model_results(self, noise_model, noise_model_settings, settings["noise_model_settings"] = noise_model_settings simulation_pipeline = [ - ModelBasedVolumeCreationAdapter(settings), + ModelBasedAdapter(settings), noise_model(settings, "noise_model_settings") ] diff --git a/simpa_tests/automatic_tests/test_path_manager.py b/simpa_tests/automatic_tests/test_path_manager.py index 11f4e565..eb462943 100644 --- a/simpa_tests/automatic_tests/test_path_manager.py +++ b/simpa_tests/automatic_tests/test_path_manager.py @@ -13,6 +13,7 @@ class TestPathManager(unittest.TestCase): + def setUp(self): self.path = '/path_config.env' self.save_path = "/workplace/data/" @@ -33,8 +34,10 @@ def setUp(self): self.simpa_home_exists = os.path.exists(self.simpa_home) @unittest.expectedFailure - def test_variables_not_set(): - path_manager = PathManager() + @patch.dict(os.environ, {Tags.SIMPA_SAVE_PATH_VARNAME: None, + Tags.MCX_BINARY_PATH_VARNAME: None}) + def test_variables_not_set(self): + path_manager = PathManager("/path/to/nowhere/") _ = path_manager.get_mcx_binary_path() _ = path_manager.get_hdf5_file_save_path() _ = path_manager.get_matlab_binary_path() @@ -42,7 +45,7 @@ def test_variables_not_set(): @patch.dict(os.environ, {Tags.SIMPA_SAVE_PATH_VARNAME: "test_simpa_save_path", Tags.MCX_BINARY_PATH_VARNAME: "test_mcx_path"}) def test_instantiate_without_file(self): - path_manager = PathManager() + path_manager = PathManager("/path/to/nowhere/") self.assertEqual(path_manager.get_mcx_binary_path(), "test_mcx_path") self.assertEqual(path_manager.get_hdf5_file_save_path(), "test_simpa_save_path") diff --git a/simpa_tests/automatic_tests/test_pipeline.py b/simpa_tests/automatic_tests/test_pipeline.py index 1b290fe9..e251a178 100644 --- a/simpa_tests/automatic_tests/test_pipeline.py +++ b/simpa_tests/automatic_tests/test_pipeline.py @@ -9,11 +9,11 @@ import numpy as np from simpa_tests.test_utils import create_test_structure_parameters import os -from simpa import ModelBasedVolumeCreationAdapter -from simpa.core.simulation_modules.optical_simulation_module.optical_forward_model_test_adapter import \ - OpticalForwardModelTestAdapter -from simpa.core.simulation_modules.acoustic_forward_module.acoustic_forward_model_test_adapter import \ - AcousticForwardModelTestAdapter +from simpa import ModelBasedAdapter +from simpa.core.simulation_modules.optical_module.optical_test_adapter import \ + OpticalTestAdapter +from simpa.core.simulation_modules.acoustic_module.acoustic_test_adapter import \ + AcousticTestAdapter from simpa.core.device_digital_twins import RSOMExplorerP50 @@ -72,9 +72,9 @@ def test_pipeline(self): }) simulation_pipeline = [ - ModelBasedVolumeCreationAdapter(settings), - OpticalForwardModelTestAdapter(settings), - AcousticForwardModelTestAdapter(settings), + ModelBasedAdapter(settings), + OpticalTestAdapter(settings), + AcousticTestAdapter(settings), ] simulate(simulation_pipeline, settings, RSOMExplorerP50(0.1, 1, 1)) 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/full_integration_test.bat b/simpa_tests/full_integration_test.bat old mode 100755 new mode 100644 diff --git a/simpa_tests/manual_tests/acoustic_forward_models/KWaveAcousticForwardConvenienceFunction.py b/simpa_tests/manual_tests/acoustic_forward_models/KWaveAcousticForwardConvenienceFunction.py index a9157c33..bf9f4e5f 100644 --- a/simpa_tests/manual_tests/acoustic_forward_models/KWaveAcousticForwardConvenienceFunction.py +++ b/simpa_tests/manual_tests/acoustic_forward_models/KWaveAcousticForwardConvenienceFunction.py @@ -5,9 +5,9 @@ from simpa.core.device_digital_twins import SlitIlluminationGeometry, LinearArrayDetectionGeometry, PhotoacousticDevice from simpa import perform_k_wave_acoustic_forward_simulation -from simpa.core.simulation_modules.reconstruction_module.reconstruction_module_delay_and_sum_adapter import \ +from simpa.core.simulation_modules.reconstruction_module.delay_and_sum_adapter import \ reconstruct_delay_and_sum_pytorch -from simpa import MCXAdapter, ModelBasedVolumeCreationAdapter, \ +from simpa import MCXAdapter, ModelBasedAdapter, \ GaussianNoise from simpa.utils import Tags, Settings, TISSUE_LIBRARY from simpa.core.simulation import simulate @@ -88,7 +88,7 @@ def setup(self): # run pipeline including volume creation and optical mcx simulation self.pipeline = [ - ModelBasedVolumeCreationAdapter(self.settings), + ModelBasedAdapter(self.settings), MCXAdapter(self.settings), GaussianNoise(self.settings, "noise_model") ] @@ -128,7 +128,7 @@ def test_convenience_function(self): self.reconstructed = reconstruct_delay_and_sum_pytorch( time_series_data.copy(), self.device.get_detection_geometry(), speed_of_sound_in_m_per_s=1540, - time_spacing_in_s=1/40_000_000_000, + time_spacing_in_s=1/40_000_000, sensor_spacing_in_mm=self.device.get_detection_geometry().pitch_mm, recon_mode=Tags.RECONSTRUCTION_MODE_PRESSURE, ) diff --git a/simpa_tests/manual_tests/acoustic_forward_models/MinimalKWaveTest.py b/simpa_tests/manual_tests/acoustic_forward_models/MinimalKWaveTest.py index 7f781d5b..35e966f8 100644 --- a/simpa_tests/manual_tests/acoustic_forward_models/MinimalKWaveTest.py +++ b/simpa_tests/manual_tests/acoustic_forward_models/MinimalKWaveTest.py @@ -27,7 +27,7 @@ def setup(self): if os.path.exists(p0_path): self.initial_pressure = np.load(p0_path)["initial_pressure"] else: - self.initial_pressure = np.zeros((100, 100, 100)) + self.initial_pressure = np.zeros((100, 30, 100)) self.initial_pressure[50, :, 50] = 1 self.speed_of_sound = np.ones((100, 30, 100)) * self.SPEED_OF_SOUND self.density = np.ones((100, 30, 100)) * 1000 diff --git a/simpa_tests/manual_tests/digital_device_twins/SimulationWithMSOTInvision.py b/simpa_tests/manual_tests/digital_device_twins/SimulationWithMSOTInvision.py index 397d2da5..81f5017b 100644 --- a/simpa_tests/manual_tests/digital_device_twins/SimulationWithMSOTInvision.py +++ b/simpa_tests/manual_tests/digital_device_twins/SimulationWithMSOTInvision.py @@ -12,7 +12,7 @@ from simpa.visualisation.matplotlib_data_visualisation import visualise_data import numpy as np from simpa.utils.path_manager import PathManager -from simpa import DelayAndSumAdapter, MCXAdapter, KWaveAdapter, ModelBasedVolumeCreationAdapter, FieldOfViewCropping +from simpa import DelayAndSumAdapter, MCXAdapter, KWaveAdapter, ModelBasedAdapter, FieldOfViewCropping from simpa.core.device_digital_twins import * from simpa_tests.manual_tests import ManualIntegrationTestClass import os @@ -26,7 +26,7 @@ def setup(self): def create_pipeline(self, _settings: Settings): return [ - ModelBasedVolumeCreationAdapter(_settings), + ModelBasedAdapter(_settings), MCXAdapter(_settings), KWaveAdapter(_settings), FieldOfViewCropping(_settings), 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/image_reconstruction/DelayAndSumReconstruction.py b/simpa_tests/manual_tests/image_reconstruction/DelayAndSumReconstruction.py index 0f29f65d..c20d2901 100644 --- a/simpa_tests/manual_tests/image_reconstruction/DelayAndSumReconstruction.py +++ b/simpa_tests/manual_tests/image_reconstruction/DelayAndSumReconstruction.py @@ -8,7 +8,7 @@ from simpa.io_handling import load_data_field from simpa.core.simulation import simulate from simpa import KWaveAdapter, MCXAdapter, \ - DelayAndSumAdapter, ModelBasedVolumeCreationAdapter, GaussianNoise + DelayAndSumAdapter, ModelBasedAdapter, GaussianNoise from simpa import reconstruct_delay_and_sum_pytorch from simpa_tests.manual_tests import ReconstructionAlgorithmTestBaseClass @@ -28,7 +28,7 @@ def test_reconstruction_of_simulation(self): self.device.update_settings_for_use_of_model_based_volume_creator(self.settings) SIMUATION_PIPELINE = [ - ModelBasedVolumeCreationAdapter(self.settings), + ModelBasedAdapter(self.settings), MCXAdapter(self.settings), GaussianNoise(self.settings, "noise_initial_pressure"), KWaveAdapter(self.settings), @@ -49,7 +49,7 @@ def test_convenience_function(self): # reconstruct via convenience function self.reconstructed_image_convenience = reconstruct_delay_and_sum_pytorch(time_series_sensor_data, self.device.get_detection_geometry(), reconstruction_settings[Tags.DATA_FIELD_SPEED_OF_SOUND], 1.0 / ( - self.device.get_detection_geometry().sampling_frequency_MHz * 1000), self.settings[Tags.SPACING_MM], reconstruction_settings[Tags.RECONSTRUCTION_MODE], reconstruction_settings[Tags.RECONSTRUCTION_APODIZATION_METHOD]) + self.device.get_detection_geometry().sampling_frequency_MHz * 1_000_000), self.settings[Tags.SPACING_MM], reconstruction_settings[Tags.RECONSTRUCTION_MODE], reconstruction_settings[Tags.RECONSTRUCTION_APODIZATION_METHOD]) # apply envelope detection method if set if reconstruction_settings[Tags.RECONSTRUCTION_BMODE_AFTER_RECONSTRUCTION]: diff --git a/simpa_tests/manual_tests/image_reconstruction/DelayMultiplyAndSumReconstruction.py b/simpa_tests/manual_tests/image_reconstruction/DelayMultiplyAndSumReconstruction.py index 7a0de0d5..e9a458a6 100644 --- a/simpa_tests/manual_tests/image_reconstruction/DelayMultiplyAndSumReconstruction.py +++ b/simpa_tests/manual_tests/image_reconstruction/DelayMultiplyAndSumReconstruction.py @@ -7,7 +7,7 @@ from simpa.io_handling import load_data_field from simpa.core.simulation import simulate from simpa import KWaveAdapter, MCXAdapter, \ - DelayMultiplyAndSumAdapter, ModelBasedVolumeCreationAdapter, GaussianNoise + DelayMultiplyAndSumAdapter, ModelBasedAdapter, GaussianNoise from simpa import reconstruct_delay_multiply_and_sum_pytorch from simpa_tests.manual_tests import ReconstructionAlgorithmTestBaseClass @@ -27,7 +27,7 @@ def test_reconstruction_of_simulation(self): self.device.update_settings_for_use_of_model_based_volume_creator(self.settings) SIMUATION_PIPELINE = [ - ModelBasedVolumeCreationAdapter(self.settings), + ModelBasedAdapter(self.settings), MCXAdapter(self.settings), GaussianNoise(self.settings, "noise_initial_pressure"), KWaveAdapter(self.settings), @@ -48,7 +48,7 @@ def test_convenience_function(self): # reconstruct via convenience function self.reconstructed_image_convenience = reconstruct_delay_multiply_and_sum_pytorch(time_series_sensor_data, self.device.get_detection_geometry(), reconstruction_settings[Tags.DATA_FIELD_SPEED_OF_SOUND], 1.0 / ( - self.device.get_detection_geometry().sampling_frequency_MHz * 1000), self.settings[Tags.SPACING_MM], reconstruction_settings[Tags.RECONSTRUCTION_MODE], reconstruction_settings[Tags.RECONSTRUCTION_APODIZATION_METHOD]) + self.device.get_detection_geometry().sampling_frequency_MHz * 1_000_000), self.settings[Tags.SPACING_MM], reconstruction_settings[Tags.RECONSTRUCTION_MODE], reconstruction_settings[Tags.RECONSTRUCTION_APODIZATION_METHOD]) # apply envelope detection method if set if reconstruction_settings[Tags.RECONSTRUCTION_BMODE_AFTER_RECONSTRUCTION]: diff --git a/simpa_tests/manual_tests/image_reconstruction/PointSourceReconstruction.py b/simpa_tests/manual_tests/image_reconstruction/PointSourceReconstruction.py index e59f35e3..60e39425 100644 --- a/simpa_tests/manual_tests/image_reconstruction/PointSourceReconstruction.py +++ b/simpa_tests/manual_tests/image_reconstruction/PointSourceReconstruction.py @@ -12,7 +12,7 @@ from simpa.visualisation.matplotlib_data_visualisation import visualise_data import numpy as np from simpa.utils.path_manager import PathManager -from simpa import DelayAndSumAdapter, MCXAdapter, KWaveAdapter, ModelBasedVolumeCreationAdapter, FieldOfViewCropping +from simpa import DelayAndSumAdapter, MCXAdapter, KWaveAdapter, ModelBasedAdapter, FieldOfViewCropping from simpa.core.device_digital_twins import * from simpa.io_handling import load_data_field @@ -155,7 +155,7 @@ def create_point_source(): }) SIMUATION_PIPELINE = [ - ModelBasedVolumeCreationAdapter(settings), + ModelBasedAdapter(settings), MCXAdapter(settings), KWaveAdapter(settings), FieldOfViewCropping(settings), diff --git a/simpa_tests/manual_tests/image_reconstruction/SignedDelayMultiplyAndSumReconstruction.py b/simpa_tests/manual_tests/image_reconstruction/SignedDelayMultiplyAndSumReconstruction.py index 754cdd01..ff61cc4e 100644 --- a/simpa_tests/manual_tests/image_reconstruction/SignedDelayMultiplyAndSumReconstruction.py +++ b/simpa_tests/manual_tests/image_reconstruction/SignedDelayMultiplyAndSumReconstruction.py @@ -7,7 +7,7 @@ from simpa.io_handling import load_data_field from simpa.core.simulation import simulate from simpa import KWaveAdapter, MCXAdapter, \ - SignedDelayMultiplyAndSumAdapter, ModelBasedVolumeCreationAdapter + SignedDelayMultiplyAndSumAdapter, ModelBasedAdapter from simpa.core.processing_components.monospectral.noise import GaussianNoise from simpa import reconstruct_signed_delay_multiply_and_sum_pytorch from simpa_tests.manual_tests import ReconstructionAlgorithmTestBaseClass @@ -28,7 +28,7 @@ def test_reconstruction_of_simulation(self): self.device.update_settings_for_use_of_model_based_volume_creator(self.settings) SIMUATION_PIPELINE = [ - ModelBasedVolumeCreationAdapter(self.settings), + ModelBasedAdapter(self.settings), MCXAdapter(self.settings), GaussianNoise(self.settings, "noise_initial_pressure"), KWaveAdapter(self.settings), @@ -49,7 +49,7 @@ def test_convenience_function(self): # reconstruct via convenience function self.reconstructed_image_convenience = reconstruct_signed_delay_multiply_and_sum_pytorch(time_series_sensor_data, self.device.get_detection_geometry(), reconstruction_settings[Tags.DATA_FIELD_SPEED_OF_SOUND], 1.0 / ( - self.device.get_detection_geometry().sampling_frequency_MHz * 1000), self.settings[Tags.SPACING_MM], reconstruction_settings[Tags.RECONSTRUCTION_MODE], reconstruction_settings[Tags.RECONSTRUCTION_APODIZATION_METHOD]) + self.device.get_detection_geometry().sampling_frequency_MHz * 1_000_000), self.settings[Tags.SPACING_MM], reconstruction_settings[Tags.RECONSTRUCTION_MODE], reconstruction_settings[Tags.RECONSTRUCTION_APODIZATION_METHOD]) # apply envelope detection method if set if reconstruction_settings[Tags.RECONSTRUCTION_BMODE_AFTER_RECONSTRUCTION]: diff --git a/simpa_tests/manual_tests/image_reconstruction/TimeReversalReconstruction.py b/simpa_tests/manual_tests/image_reconstruction/TimeReversalReconstruction.py index 72a8c49f..cad47201 100644 --- a/simpa_tests/manual_tests/image_reconstruction/TimeReversalReconstruction.py +++ b/simpa_tests/manual_tests/image_reconstruction/TimeReversalReconstruction.py @@ -7,7 +7,7 @@ from simpa.io_handling import load_data_field from simpa.core.simulation import simulate from simpa import KWaveAdapter, MCXAdapter, \ - TimeReversalAdapter, ModelBasedVolumeCreationAdapter, GaussianNoise + TimeReversalAdapter, ModelBasedAdapter, GaussianNoise from simpa import reconstruct_delay_and_sum_pytorch from simpa_tests.manual_tests import ReconstructionAlgorithmTestBaseClass @@ -27,7 +27,7 @@ def test_reconstruction_of_simulation(self): self.device.update_settings_for_use_of_model_based_volume_creator(self.settings) SIMULATION_PIPELINE = [ - ModelBasedVolumeCreationAdapter(self.settings), + ModelBasedAdapter(self.settings), MCXAdapter(self.settings), GaussianNoise(self.settings, "noise_initial_pressure"), KWaveAdapter(self.settings), diff --git a/simpa_tests/manual_tests/optical_forward_models/AbsorptionAndScatteringWithInifinitesimalSlabExperiment.py b/simpa_tests/manual_tests/optical_forward_models/AbsorptionAndScatteringWithInifinitesimalSlabExperiment.py index 5ac62e06..04199a5e 100644 --- a/simpa_tests/manual_tests/optical_forward_models/AbsorptionAndScatteringWithInifinitesimalSlabExperiment.py +++ b/simpa_tests/manual_tests/optical_forward_models/AbsorptionAndScatteringWithInifinitesimalSlabExperiment.py @@ -28,7 +28,7 @@ import matplotlib.pyplot as plt import numpy as np -from simpa import MCXAdapter, ModelBasedVolumeCreationAdapter +from simpa import MCXAdapter, ModelBasedAdapter from simpa.core.device_digital_twins import PhotoacousticDevice, PencilBeamIlluminationGeometry from simpa.core.simulation import simulate from simpa.io_handling import load_data_field @@ -208,7 +208,7 @@ def test_simulation(self, distance=10, expected_decay_ratio=np.e, scattering_val self.settings.get_optical_settings()[Tags.MCX_ASSUMED_ANISOTROPY] = anisotropy_value pipeline = [ - ModelBasedVolumeCreationAdapter(self.settings), + ModelBasedAdapter(self.settings), MCXAdapter(self.settings) ] diff --git a/simpa_tests/manual_tests/optical_forward_models/AbsorptionAndScatteringWithinHomogenousMedium.py b/simpa_tests/manual_tests/optical_forward_models/AbsorptionAndScatteringWithinHomogenousMedium.py index 806863bd..1d62a943 100644 --- a/simpa_tests/manual_tests/optical_forward_models/AbsorptionAndScatteringWithinHomogenousMedium.py +++ b/simpa_tests/manual_tests/optical_forward_models/AbsorptionAndScatteringWithinHomogenousMedium.py @@ -28,7 +28,7 @@ import matplotlib.pyplot as plt import numpy as np -from simpa import MCXAdapter, ModelBasedVolumeCreationAdapter +from simpa import MCXAdapter, ModelBasedAdapter from simpa.core.device_digital_twins import PhotoacousticDevice, PencilBeamIlluminationGeometry from simpa.core.simulation import simulate from simpa.io_handling import load_data_field @@ -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,10 +220,13 @@ 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), + ModelBasedAdapter(self.settings), MCXAdapter(self.settings) ] @@ -241,13 +244,10 @@ 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), + ModelBasedAdapter(self.settings), MCXAdapter(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/manual_tests/optical_forward_models/CompareMCXResultsWithDiffusionTheory.py b/simpa_tests/manual_tests/optical_forward_models/CompareMCXResultsWithDiffusionTheory.py index 898fd1f6..939fc870 100644 --- a/simpa_tests/manual_tests/optical_forward_models/CompareMCXResultsWithDiffusionTheory.py +++ b/simpa_tests/manual_tests/optical_forward_models/CompareMCXResultsWithDiffusionTheory.py @@ -4,7 +4,7 @@ from simpa.utils import Tags, PathManager, Settings, TISSUE_LIBRARY from simpa.core.simulation import simulate -from simpa import ModelBasedVolumeCreationAdapter, MCXAdapter +from simpa import ModelBasedAdapter, MCXAdapter from simpa.core.device_digital_twins import PhotoacousticDevice, PencilBeamIlluminationGeometry from simpa.io_handling import load_data_field import numpy as np @@ -130,7 +130,7 @@ def run_simulation(self, distance, spacing): # run pipeline including volume creation and optical mcx simulation pipeline = [ - ModelBasedVolumeCreationAdapter(self.settings), + ModelBasedAdapter(self.settings), MCXAdapter(self.settings), ] simulate(pipeline, self.settings, self.device) diff --git a/simpa_tests/manual_tests/optical_forward_models/ComputeDiffuseReflectance.py b/simpa_tests/manual_tests/optical_forward_models/ComputeDiffuseReflectance.py index 116e55b7..45304974 100644 --- a/simpa_tests/manual_tests/optical_forward_models/ComputeDiffuseReflectance.py +++ b/simpa_tests/manual_tests/optical_forward_models/ComputeDiffuseReflectance.py @@ -4,7 +4,7 @@ from simpa.utils import Tags, PathManager, Settings, TISSUE_LIBRARY from simpa.core.simulation import simulate -from simpa import ModelBasedVolumeCreationAdapter, MCXAdapterReflectance +from simpa import ModelBasedAdapter, MCXReflectanceAdapter from simpa.core.device_digital_twins import PhotoacousticDevice, PencilBeamIlluminationGeometry from simpa.io_handling import load_data_field import numpy as np @@ -136,8 +136,8 @@ def run_simulation(self, distance, spacing): # run pipeline including volume creation and optical mcx simulation pipeline = [ - ModelBasedVolumeCreationAdapter(self.settings), - MCXAdapterReflectance(self.settings), + ModelBasedAdapter(self.settings), + MCXReflectanceAdapter(self.settings), ] simulate(pipeline, self.settings, self.device) diff --git a/simpa_tests/manual_tests/processing_components/QPAIReconstruction.py b/simpa_tests/manual_tests/processing_components/QPAIReconstruction.py index c262bcbb..21bf4e56 100644 --- a/simpa_tests/manual_tests/processing_components/QPAIReconstruction.py +++ b/simpa_tests/manual_tests/processing_components/QPAIReconstruction.py @@ -7,7 +7,7 @@ from simpa_tests.manual_tests import ManualIntegrationTestClass from simpa.core.device_digital_twins import RSOMExplorerP50 from simpa.core.processing_components.monospectral.iterative_qPAI_algorithm import IterativeqPAI -from simpa import MCXAdapter, ModelBasedVolumeCreationAdapter, \ +from simpa import MCXAdapter, ModelBasedAdapter, \ GaussianNoise from simpa.utils import Tags, Settings, TISSUE_LIBRARY from simpa.core.simulation import simulate @@ -85,7 +85,7 @@ def setup(self): # run pipeline including volume creation and optical mcx simulation pipeline = [ - ModelBasedVolumeCreationAdapter(self.settings), + ModelBasedAdapter(self.settings), MCXAdapter(self.settings), GaussianNoise(self.settings, "noise_model") ] diff --git a/simpa_tests/manual_tests/processing_components/TestLinearUnmixingVisual.py b/simpa_tests/manual_tests/processing_components/TestLinearUnmixingVisual.py index 38896d64..7f76e1ac 100644 --- a/simpa_tests/manual_tests/processing_components/TestLinearUnmixingVisual.py +++ b/simpa_tests/manual_tests/processing_components/TestLinearUnmixingVisual.py @@ -65,7 +65,7 @@ def setup(self): # Run simulation pipeline for all wavelengths in Tag.WAVELENGTHS self.pipeline = [ - sp.ModelBasedVolumeCreationAdapter(self.settings) + sp.ModelBasedAdapter(self.settings) ] def perform_test(self): diff --git a/simpa_tests/manual_tests/test_with_experimental_measurements/ReproduceDISMeasurements.py b/simpa_tests/manual_tests/test_with_experimental_measurements/ReproduceDISMeasurements.py index dc464bec..9c8ab38f 100644 --- a/simpa_tests/manual_tests/test_with_experimental_measurements/ReproduceDISMeasurements.py +++ b/simpa_tests/manual_tests/test_with_experimental_measurements/ReproduceDISMeasurements.py @@ -24,7 +24,7 @@ from simpa.core.device_digital_twins import * import numpy as np from simpa.visualisation.matplotlib_data_visualisation import visualise_data -from simpa import ModelBasedVolumeCreationAdapter, MCXAdapter +from simpa import ModelBasedAdapter, MCXAdapter from simpa_tests.manual_tests.test_with_experimental_measurements.utils import read_reference_spectra, read_rxt_file from simpa_tests.manual_tests import ManualIntegrationTestClass import inspect @@ -177,7 +177,7 @@ def create_measurement_setup(sample_tickness_mm): }) self.pipeline = [ - ModelBasedVolumeCreationAdapter(self.settings), + ModelBasedAdapter(self.settings), MCXAdapter(self.settings), ] diff --git a/simpa_tests/manual_tests/volume_creation/SegmentationLoader.py b/simpa_tests/manual_tests/volume_creation/SegmentationLoader.py index 780d92a7..230a8e3c 100644 --- a/simpa_tests/manual_tests/volume_creation/SegmentationLoader.py +++ b/simpa_tests/manual_tests/volume_creation/SegmentationLoader.py @@ -72,7 +72,7 @@ def segmentation_class_mapping(): }) self.pipeline = [ - sp.SegmentationBasedVolumeCreationAdapter(self.settings), + sp.SegmentationBasedAdapter(self.settings), sp.MCXAdapter(self.settings) ] 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