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