diff --git a/docs/source/simpa.core.simulation_modules.acoustic_module.rst b/docs/source/simpa.core.simulation_modules.acoustic_module.rst new file mode 100644 index 00000000..43ede865 --- /dev/null +++ b/docs/source/simpa.core.simulation_modules.acoustic_module.rst @@ -0,0 +1,24 @@ +acoustic\_module +======================================================= + +.. automodule:: simpa.core.simulation_modules.acoustic_module + :members: + :undoc-members: + :show-inheritance: + +.. automodule:: simpa.core.simulation_modules.acoustic_module.acoustic_adapter_base + :members: + :undoc-members: + :show-inheritance: + + +.. automodule:: simpa.core.simulation_modules.acoustic_module.acoustic_test_adapter + :members: + :undoc-members: + :show-inheritance: + + +.. automodule:: simpa.core.simulation_modules.acoustic_module.k_wave_adapter + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/simpa.core.simulation_modules.optical_module.rst b/docs/source/simpa.core.simulation_modules.optical_module.rst new file mode 100644 index 00000000..f8809ef0 --- /dev/null +++ b/docs/source/simpa.core.simulation_modules.optical_module.rst @@ -0,0 +1,36 @@ +optical\_module +====================================================== + +.. automodule:: simpa.core.simulation_modules.optical_module + :members: + :undoc-members: + :show-inheritance: + +.. automodule:: simpa.core.simulation_modules.optical_module.mcx_adapter + :members: + :undoc-members: + :show-inheritance: + + +.. automodule:: simpa.core.simulation_modules.optical_module.mcx_reflectance_adapter + :members: + :undoc-members: + :show-inheritance: + + +.. automodule:: simpa.core.simulation_modules.optical_module.optical_adapter_base + :members: + :undoc-members: + :show-inheritance: + + +.. automodule:: simpa.core.simulation_modules.optical_module.optical_test_adapter + :members: + :undoc-members: + :show-inheritance: + + +.. automodule:: simpa.core.simulation_modules.optical_module.volume_boundary_condition + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/simpa.core.simulation_modules.optical_simulation_module.rst b/docs/source/simpa.core.simulation_modules.optical_simulation_module.rst index 8d6bc5c9..558a8de6 100644 --- a/docs/source/simpa.core.simulation_modules.optical_simulation_module.rst +++ b/docs/source/simpa.core.simulation_modules.optical_simulation_module.rst @@ -22,3 +22,9 @@ optical\_simulation\_module :members: :undoc-members: :show-inheritance: + + +.. automodule:: simpa.core.simulation_modules.optical_simulation_module.volume_boundary_condition + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/simpa.utils.libraries.refractive_index_spectra_data.rst b/docs/source/simpa.utils.libraries.refractive_index_spectra_data.rst new file mode 100644 index 00000000..730b7918 --- /dev/null +++ b/docs/source/simpa.utils.libraries.refractive_index_spectra_data.rst @@ -0,0 +1,7 @@ +refractive\_index\_spectra\_data +============================================================== + +.. automodule:: simpa.utils.libraries.refractive_index_spectra_data + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/source/simpa.utils.libraries.rst b/docs/source/simpa.utils.libraries.rst index fdff5eef..b128b0c2 100644 --- a/docs/source/simpa.utils.libraries.rst +++ b/docs/source/simpa.utils.libraries.rst @@ -11,6 +11,7 @@ libraries simpa.utils.libraries.absorption_spectra_data simpa.utils.libraries.anisotropy_spectra_data + simpa.utils.libraries.refractive_index_spectra_data simpa.utils.libraries.scattering_spectra_data simpa.utils.libraries.structure_library diff --git a/docs/source/three_vs_two_dimensional_simulation_example.rst b/docs/source/three_vs_two_dimensional_simulation_example.rst new file mode 100644 index 00000000..f68198fb --- /dev/null +++ b/docs/source/three_vs_two_dimensional_simulation_example.rst @@ -0,0 +1,7 @@ +three_vs_two_dimensional_simulation_example +========================================= + +.. literalinclude:: ../../simpa_examples/three_vs_two_dimensional_simulation_example.py + :language: python + :lines: 1- + diff --git a/poetry.lock b/poetry.lock new file mode 100644 index 00000000..e69de29b diff --git a/pre_commit_configs/link-config.json b/pre_commit_configs/link-config.json new file mode 100644 index 00000000..f1dd61d8 --- /dev/null +++ b/pre_commit_configs/link-config.json @@ -0,0 +1,7 @@ +{ + "ignorePatterns": [ + { + "pattern": "https://doi.org/10.1117/1.JBO.27.8.083010" + } + ] +} \ No newline at end of file diff --git a/simpa/__init__.py b/simpa/__init__.py index c653bfe1..ff6a18d2 100644 --- a/simpa/__init__.py +++ b/simpa/__init__.py @@ -19,7 +19,7 @@ from .core.simulation_modules.optical_module.mcx_adapter import \ MCXAdapter from .core.simulation_modules.optical_module.mcx_reflectance_adapter import \ - MCXReflectanceAdapter + MCXReflectanceAdapter, FastMCXReflectanceAdapter from .core.simulation_modules.acoustic_module.k_wave_adapter import \ KWaveAdapter from .core.simulation_modules.reconstruction_module.delay_and_sum_adapter import \ diff --git a/simpa/core/device_digital_twins/illumination_geometries/illumination_geometry_base.py b/simpa/core/device_digital_twins/illumination_geometries/illumination_geometry_base.py index 6c5780bf..23dd6197 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/illumination_geometry_base.py +++ b/simpa/core/device_digital_twins/illumination_geometries/illumination_geometry_base.py @@ -51,6 +51,20 @@ def get_mcx_illuminator_definition(self, global_settings) -> dict: """ pass + @abstractmethod + def get_mmc_illuminator_definition(self, global_settings) -> dict: + """ + IMPORTANT: This method creates a dictionary that contains tags as they are expected for the + MMC simulation tool to represent the illumination geometry of this device. + + :param global_settings: The global_settings instance containing the simulation instructions. + :type global_settings: Settings + + :return: Dictionary that includes all parameters needed for mcx. + :rtype: dict + """ + pass + def check_settings_prerequisites(self, global_settings) -> bool: return True diff --git a/simpa/core/device_digital_twins/illumination_geometries/rectangle_illumination.py b/simpa/core/device_digital_twins/illumination_geometries/rectangle_illumination.py index 69bd3588..a0c83c23 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/rectangle_illumination.py +++ b/simpa/core/device_digital_twins/illumination_geometries/rectangle_illumination.py @@ -24,6 +24,7 @@ class RectangleIlluminationGeometry(IlluminationGeometryBase): def __init__(self, length_mm: int = 10, width_mm: int = 10, + focal_length_in_mm: float | str | None = None, device_position_mm: typing.Optional[np.ndarray] = None, source_direction_vector: typing.Optional[np.ndarray] = None, field_of_view_extent_mm: typing.Optional[np.ndarray] = None): @@ -51,6 +52,12 @@ def __init__(self, assert length_mm > 0 assert width_mm > 0 + if isinstance(focal_length_in_mm, str): + assert ((focal_length_in_mm == "_NaN_" + or focal_length_in_mm == "-_Inf_") + or focal_length_in_mm == "_Inf_"), f"{focal_length_in_mm} is not supported yet" + + self.focal_length_in_mm = focal_length_in_mm self.length_mm = length_mm self.width_mm = width_mm @@ -68,14 +75,20 @@ def get_mcx_illuminator_definition(self, global_settings: Settings) -> dict: spacing = global_settings[Tags.SPACING_MM] - device_position = list(np.rint(self.device_position_mm / spacing)) + device_position = list(self.device_position_mm / spacing) self.logger.debug(device_position) source_direction = list(self.normalized_source_direction_vector) - source_param1 = [np.rint(self.width_mm / spacing) + 1, 0, 0] - source_param2 = [0, np.rint(self.length_mm / spacing) + 1, 0] + if self.focal_length_in_mm is not None: + param4 = self.focal_length_in_mm if isinstance( + self.focal_length_in_mm, str) else self.focal_length_in_mm / spacing + source_direction.append(param4) + + source_param1 = [self.width_mm / spacing + 2, 0., 0.] + # Legit ka warum +2 but only then the whole area is illuminated correctly + source_param2 = [0., self.length_mm / spacing + 2, 0.] # If Pos=[10, 12, 0], Param1=[10, 0, 0], Param2=[0, 20, 0], # then illumination covers: x in [10, 20], y in [12, 32] @@ -88,6 +101,20 @@ def get_mcx_illuminator_definition(self, global_settings: Settings) -> dict: "Param2": source_param2 } + def get_mmc_illuminator_definition(self, global_settings: Settings) -> dict: + """ + Returns the illumination parameters for MMC simulations. + + :param global_settings: The global settings. + + :return: The illumination parameters as a dictionary. + """ + + mcx_illuminator_definition = self.get_mcx_illuminator_definition(global_settings) + mcx_illuminator_definition["Param1"] += [0.0] + mcx_illuminator_definition["Param2"] += [0.0] + return mcx_illuminator_definition + def serialize(self) -> dict: """ Serializes the object into a dictionary. diff --git a/simpa/core/simulation_modules/optical_module/mcx_adapter.py b/simpa/core/simulation_modules/optical_module/mcx_adapter.py index 551787ec..a9d1200b 100644 --- a/simpa/core/simulation_modules/optical_module/mcx_adapter.py +++ b/simpa/core/simulation_modules/optical_module/mcx_adapter.py @@ -42,6 +42,7 @@ def forward_model(self, absorption_cm: np.ndarray, scattering_cm: np.ndarray, anisotropy: np.ndarray, + refractive_index: np.ndarray, illumination_geometry: IlluminationGeometryBase) -> Dict: """ runs the MCX simulations. Binary file containing scattering and absorption volumes is temporarily created as @@ -52,24 +53,18 @@ def forward_model(self, :param absorption_cm: array containing the absorption of the tissue in `cm` units :param scattering_cm: array containing the scattering of the tissue in `cm` units :param anisotropy: array containing the anisotropy of the volume defined by `absorption_cm` and `scattering_cm` + :param refractive_index: array containing the refractive index of the volume defined by `absorption_cm` and `scattering_cm` :param illumination_geometry: and instance of `IlluminationGeometryBase` defining the illumination geometry :return: `Dict` containing the results of optical simulations, the keys in this dictionary-like object depend on the Tags defined in `self.component_settings` """ - if Tags.MCX_ASSUMED_ANISOTROPY in self.component_settings: - _assumed_anisotropy = self.component_settings[Tags.MCX_ASSUMED_ANISOTROPY] - else: - _assumed_anisotropy = 0.9 - self.generate_mcx_bin_input(absorption_cm=absorption_cm, scattering_cm=scattering_cm, anisotropy=anisotropy, - assumed_anisotropy=_assumed_anisotropy) + refractive_index=refractive_index) - settings_dict = self.get_mcx_settings(illumination_geometry=illumination_geometry, - assumed_anisotropy=_assumed_anisotropy) + settings_dict = self.get_mcx_settings(illumination_geometry=illumination_geometry) - print(settings_dict) self.generate_mcx_json_input(settings_dict=settings_dict) # run the simulation cmd = self.get_command() @@ -99,14 +94,12 @@ def generate_mcx_json_input(self, settings_dict: Dict) -> None: def get_mcx_settings(self, illumination_geometry: IlluminationGeometryBase, - assumed_anisotropy: np.ndarray, **kwargs) -> Dict: """ generates MCX-specific settings for simulations based on Tags in `self.global_settings` and `self.component_settings` . Among others, it defines the volume type, dimensions and path to binary file. :param illumination_geometry: and instance of `IlluminationGeometryBase` defining the illumination geometry - :param assumed_anisotropy: :param kwargs: dummy, used for class inheritance :return: dictionary with settings to be used by MCX """ @@ -124,6 +117,7 @@ def get_mcx_settings(self, self.frames = int(time / dt) source = illumination_geometry.get_mcx_illuminator_definition(self.global_settings) + mcx_length_unit = self.global_settings[Tags.SPACING_MM] if Tags.TRUE_SPACING_MM not in self.global_settings else self.global_settings[Tags.TRUE_SPACING_MM] settings_dict = { "Session": { "ID": mcx_volumetric_data_file, @@ -141,22 +135,16 @@ def get_mcx_settings(self, }, "Domain": { "OriginType": 0, - "LengthUnit": self.global_settings[Tags.SPACING_MM], + "LengthUnit": mcx_length_unit, "Media": [ { "mua": 0, "mus": 0, "g": 1, "n": 1 - }, - { - "mua": 1, - "mus": 1, - "g": assumed_anisotropy, - "n": 1 } ], - "MediaFormat": "muamus_float", + "MediaFormat": "asgn_float", "Dim": [self.nx, self.ny, self.nz], "VolumeFile": self.global_settings[Tags.SIMULATION_PATH] + "/" + self.global_settings[Tags.VOLUME_NAME] + ".bin" @@ -207,22 +195,23 @@ def generate_mcx_bin_input(self, absorption_cm: np.ndarray, scattering_cm: np.ndarray, anisotropy: np.ndarray, - assumed_anisotropy: np.ndarray) -> None: + refractive_index: np.ndarray) -> None: """ generates binary file containing volume scattering and absorption as input for MCX :param absorption_cm: Absorption in units of per centimeter :param scattering_cm: Scattering in units of per centimeter :param anisotropy: Dimensionless scattering anisotropy - :param assumed_anisotropy: + :param refractive_index: Refractive index :return: None """ - absorption_mm, scattering_mm = self.pre_process_volumes(**{'absorption_cm': absorption_cm, - 'scattering_cm': scattering_cm, - 'anisotropy': anisotropy, - 'assumed_anisotropy': assumed_anisotropy}) - # stack arrays to give array with shape (nx,ny,nz,2) - op_array = np.stack([absorption_mm, scattering_mm], axis=-1, dtype=np.float32) + absorption_mm = self.pre_process_volume(absorption_cm, is_mua_or_mus=True) + scattering_mm = self.pre_process_volume(scattering_cm, is_mua_or_mus=True) + anisotropy = self.pre_process_volume(anisotropy, is_mua_or_mus=False) + refractive_index = self.pre_process_volume(refractive_index, is_mua_or_mus=False) + + # stack arrays to give array with shape (nx,ny,nz,4) - where the 4 floats correspond to mua/mus/g/n + op_array = np.stack([absorption_mm, scattering_mm, anisotropy, refractive_index], axis=-1, dtype=np.float32) [self.nx, self.ny, self.nz, _] = np.shape(op_array) # # create a binary of the volume tmp_input_path = self.global_settings[Tags.SIMULATION_PATH] + "/" + \ @@ -240,11 +229,9 @@ def read_mcx_output(self, **kwargs) -> Dict: """ content = jdata.load(self.mcx_volumetric_data_file) fluence = content['NIFTIData'] - print(f"fluence.shape {fluence.shape}") if fluence.ndim > 3: # remove the 1 or 2 (for mcx >= v2024.1) additional dimensions of size 1 if present to obtain a 3d array fluence = fluence.reshape(fluence.shape[0], fluence.shape[1], -1) - print(f"fluence.shape {fluence.shape}") results = dict() results[Tags.DATA_FIELD_FLUENCE] = fluence return results @@ -259,12 +246,25 @@ def remove_mcx_output(self) -> None: if os.path.isfile(f): os.remove(f) + def pre_process_volume(self, volume: np.ndarray, is_mua_or_mus: bool) -> np.ndarray: + """ + pre-process volumes before running simulations with MCX. The volumes are transformed to `mm` units + + :param kwargs: dictionary containing at least the keys `scattering_cm, absorption_cm, anisotropy, refractive_index` + :return: `Tuple` of volumes after transformation + """ + if is_mua_or_mus: + volume /= 10 # Conversion from 1/cm to 1/mm (required by MCX) + volume[volume == 0] = np.nan + + # No preprocessing is done on anisotropy and refractive index + return volume + def pre_process_volumes(self, **kwargs) -> Tuple: """ pre-process volumes before running simulations with MCX. The volumes are transformed to `mm` units - :param kwargs: dictionary containing at least the keys `scattering_cm, absorption_cm, anisotropy` and - `assumed_anisotropy` + :param kwargs: dictionary containing at least the keys `scattering_cm, absorption_cm, anisotropy, refractive_index` :return: `Tuple` of volumes after transformation """ return self.volumes_to_mm(**kwargs) @@ -274,29 +274,16 @@ def volumes_to_mm(**kwargs) -> Tuple: """ transforms volumes into `mm` units - :param kwargs: dictionary containing at least the keys `scattering_cm, absorption_cm, anisotropy` and - `assumed_anisotropy` + :param kwargs: dictionary containing at least the keys `scattering_cm, absorption_cm, anisotropy, refractive_index` :return: `Tuple` of volumes after transformation """ scattering_cm = kwargs.get('scattering_cm') absorption_cm = kwargs.get('absorption_cm') + anisotropy = kwargs.get('anisotropy') + refractive_index = kwargs.get('refractive_index') absorption_mm = absorption_cm / 10 scattering_mm = scattering_cm / 10 - - # FIXME Currently, mcx only accepts a single value for the anisotropy. - # In order to use the correct reduced scattering coefficient throughout the simulation, - # we adjust the scattering parameter to be more accurate in the diffuse regime. - # This will lead to errors, especially in the quasi-ballistic regime. - - given_reduced_scattering = (scattering_mm * (1 - kwargs.get('anisotropy'))) - - # If the anisotropy is 1, all scattering is forward scattering which is equal to no scattering at all - if kwargs.get("assumed_anisotropy") == 1: - scattering_mm = given_reduced_scattering * 0 - else: - scattering_mm = given_reduced_scattering / (1 - kwargs.get('assumed_anisotropy')) - scattering_mm[scattering_mm < 1e-10] = 1e-10 - return absorption_mm, scattering_mm + return absorption_mm, scattering_mm, anisotropy, refractive_index @staticmethod def post_process_volumes(**kwargs) -> Tuple: diff --git a/simpa/core/simulation_modules/optical_module/mcx_reflectance_adapter.py b/simpa/core/simulation_modules/optical_module/mcx_reflectance_adapter.py index 4473c965..3ae140dd 100644 --- a/simpa/core/simulation_modules/optical_module/mcx_reflectance_adapter.py +++ b/simpa/core/simulation_modules/optical_module/mcx_reflectance_adapter.py @@ -1,12 +1,18 @@ # SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT +import pathlib +import typing + import numpy as np import struct import jdata import os -from typing import List, Tuple, Dict, Union +from typing import Tuple, Dict, Union + +from simpa.io_handling.io_hdf5 import load_data_field +from simpa.core.simulation_modules.optical_module.volume_boundary_condition import MCXVolumeBoundaryCondition from simpa.utils import Tags, Settings from simpa.core.simulation_modules.optical_module.mcx_adapter import MCXAdapter from simpa.core.device_digital_twins import IlluminationGeometryBase, PhotoacousticDevice @@ -38,13 +44,20 @@ def __init__(self, global_settings: Settings): super(MCXReflectanceAdapter, self).__init__(global_settings=global_settings) self.mcx_photon_data_file = None self.padded = None + if Tags.VOLUME_BOUNDARY_CONDITION in global_settings: + self.volume_boundary_condition_str = global_settings[Tags.VOLUME_BOUNDARY_CONDITION] + else: + self.volume_boundary_condition_str = MCXVolumeBoundaryCondition.DEFAULT.value + self.mcx_output_suffixes = {'mcx_volumetric_data_file': '.jnii', + 'mcx_volumetric_data_file_camera': '.bin', 'mcx_photon_data_file': '_detp.jdat'} def forward_model(self, absorption_cm: np.ndarray, scattering_cm: np.ndarray, anisotropy: np.ndarray, + refractive_index: np.ndarray, illumination_geometry: IlluminationGeometryBase) -> Dict: """ runs the MCX simulations. Binary file containing scattering and absorption volumes is temporarily created as @@ -55,25 +68,18 @@ def forward_model(self, :param absorption_cm: array containing the absorption of the tissue in `cm` units :param scattering_cm: array containing the scattering of the tissue in `cm` units :param anisotropy: array containing the anisotropy of the volume defined by `absorption_cm` and `scattering_cm` + :param refractive_index: array containing the refractive index of the volume defined by `absorption_cm` and `scattering_cm` :param illumination_geometry: and instance of `IlluminationGeometryBase` defining the illumination geometry - :param probe_position_mm: position of a probe in `mm` units. This is parsed to - `illumination_geometry.get_mcx_illuminator_definition` :return: `Settings` containing the results of optical simulations, the keys in this dictionary-like object depend on the Tags defined in `self.component_settings` """ - if Tags.MCX_ASSUMED_ANISOTROPY in self.component_settings: - _assumed_anisotropy = self.component_settings[Tags.MCX_ASSUMED_ANISOTROPY] - else: - _assumed_anisotropy = 0.9 self.generate_mcx_bin_input(absorption_cm=absorption_cm, scattering_cm=scattering_cm, - anisotropy=_assumed_anisotropy, - assumed_anisotropy=_assumed_anisotropy) + anisotropy=anisotropy, + refractive_index=refractive_index) - settings_dict = self.get_mcx_settings(illumination_geometry=illumination_geometry, - assumed_anisotropy=_assumed_anisotropy, - ) + settings_dict = self.get_mcx_settings(illumination_geometry=illumination_geometry) print(settings_dict) self.generate_mcx_json_input(settings_dict=settings_dict) @@ -90,9 +96,52 @@ def forward_model(self, self.remove_mcx_output() return results - def get_command(self) -> List: - """ - generates list of commands to be parse to MCX in a subprocess + def get_mcx_settings(self, + illumination_geometry: IlluminationGeometryBase, + **kwargs) -> Dict: + settings_dict = super().get_mcx_settings(illumination_geometry=illumination_geometry, **kwargs) + uses_photon_exit_data = Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT in self.component_settings and self.component_settings[ + Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT] + contains_camera_settings = Tags.MCX_CAMERA_SETTINGS in self.global_settings + contains_backtrack_settings = Tags.MCX_BACKTRACK_SETTINGS in self.global_settings + assert not contains_camera_settings or not contains_backtrack_settings + + if uses_photon_exit_data: + if Tags.MCX_DETECTOR in self.global_settings: + settings_dict["Optode"]["Detector"] = self.global_settings[Tags.MCX_DETECTOR] + else: + # For some reason, the simulation gets slower the larger the detector is + width = self.global_settings[Tags.DIM_VOLUME_X_MM] / self.global_settings[Tags.SPACING_MM] + height = self.global_settings[Tags.DIM_VOLUME_Y_MM] / self.global_settings[Tags.SPACING_MM] + position = [width / 2 + 1, height / 2 + 1, 0.0] + radius = np.sqrt(width ** 2 + height ** 2) / 2 + settings_dict["Optode"]["Detector"] = [ + { + "Pos": position, + "R": radius + } + ] + + if contains_camera_settings: + camera_settings = self.global_settings[Tags.MCX_CAMERA_SETTINGS] + settings_dict["Camera"] = { + "ObjectDistance": camera_settings[Tags.MCX_OBJECT_DISTANCE], + "ProjectionDistance": camera_settings[Tags.MCX_PROJECTION_DISTANCE], + "FocalLength": camera_settings[Tags.MCX_FOCAL_LENGTH], + "ApertureRadius": camera_settings[Tags.MCX_APERTURE_RADIUS] + } + elif contains_backtrack_settings: + backtrack_settings = self.global_settings.get_mcx_backtrack_settings() + settings_dict["Backtrack"] = { + "ObjectDistance": backtrack_settings[Tags.MCX_OBJECT_DISTANCE], + "IdealDistance": backtrack_settings[Tags.MCX_IDEAL_DISTANCE], + "ApertureRadius": backtrack_settings[Tags.MCX_APERTURE_RADIUS], + "TrueApertureRadius": backtrack_settings[Tags.MCX_TRUE_APERTURE_RADIUS] + } + return settings_dict + + def get_command(self) -> typing.List: + """Generates list of commands to be parse to MCX in a subprocess. :return: list of MCX commands """ @@ -102,22 +151,34 @@ def get_command(self) -> List: cmd.append(self.mcx_json_config_file) cmd.append("-O") cmd.append("F") + cmd.append("-b") + cmd.append("1") # use 'C' order array format for binary input file cmd.append("-a") cmd.append("1") cmd.append("-F") cmd.append("jnii") - if Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT in self.component_settings and \ - self.component_settings[Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT]: - cmd.append("-H") - cmd.append(f"{int(self.component_settings[Tags.OPTICAL_MODEL_NUMBER_PHOTONS])}") - cmd.append("--bc") # save photon exit position and direction - cmd.append("______000010") + cmd.append("-e") + cmd.append(str(1e-2)) + cmd.append("--bc") + cmd.append(self.volume_boundary_condition_str) + cmd.append("-H") + cmd.append( + f"{int(self.component_settings[Tags.OPTICAL_MODEL_NUMBER_PHOTONS])}" + ) + if ( + Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT in self.component_settings + and self.component_settings[Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT] + ): cmd.append("--savedetflag") cmd.append("XV") - if Tags.COMPUTE_DIFFUSE_REFLECTANCE in self.component_settings and \ - self.component_settings[Tags.COMPUTE_DIFFUSE_REFLECTANCE]: - cmd.append("--saveref") # save diffuse reflectance at 0 filled voxels outside of domain + + if ( + Tags.COMPUTE_DIFFUSE_REFLECTANCE in self.component_settings + and self.component_settings[Tags.COMPUTE_DIFFUSE_REFLECTANCE] + ): + cmd.append("--saveref") + cmd += self.get_additional_flags() return cmd @@ -146,6 +207,12 @@ def read_mcx_output(self, **kwargs) -> Dict: self.component_settings[Tags.COMPUTE_DIFFUSE_REFLECTANCE]: results[Tags.DATA_FIELD_DIFFUSE_REFLECTANCE] = ref results[Tags.DATA_FIELD_DIFFUSE_REFLECTANCE_POS] = ref_pos + + if (Tags.MCX_CAMERA_SETTINGS in self.global_settings) or (Tags.MCX_BACKTRACK_SETTINGS in self.global_settings): + cam_intensity_file_path = pathlib.Path(self.mcx_volumetric_data_file).with_suffix(".bin") + cam_intensity = np.fromfile(cam_intensity_file_path, dtype=np.float32) + results[Tags.DATA_FIELD_CAMERA_INTENSITY] = cam_intensity + if Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT in self.component_settings and \ self.component_settings[Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT]: content = jdata.load(self.mcx_photon_data_file) @@ -190,7 +257,9 @@ def pre_process_volumes(self, **kwargs) -> Tuple: check_padding = (Tags.COMPUTE_DIFFUSE_REFLECTANCE in self.component_settings and self.component_settings[Tags.COMPUTE_DIFFUSE_REFLECTANCE]) or \ (Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT in self.component_settings and - self.component_settings[Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT]) + self.component_settings[Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT]) or \ + Tags.MCX_CAMERA_SETTINGS in self.global_settings + # check that all volumes on first layer along z have only 0 values if np.any([np.any(a[:, :, 0] != 0)] for a in arrays) and check_padding: results = tuple(np.pad(a, ((0, 0), (0, 0), (1, 0)), "constant", constant_values=0) for a in arrays) @@ -226,7 +295,8 @@ def run_forward_model(self, device: Union[IlluminationGeometryBase, PhotoacousticDevice], absorption: np.ndarray, scattering: np.ndarray, - anisotropy: np.ndarray + anisotropy: np.ndarray, + refractive_index: np.ndarray ) -> Dict: """ runs `self.forward_model` as many times as defined by `device` and aggregates the results. @@ -236,49 +306,60 @@ def run_forward_model(self, :param absorption: Absorption volume :param scattering: Scattering volume :param anisotropy: Dimensionless scattering anisotropy + :param refractive_index: Refractive index :return: """ reflectance = [] reflectance_position = [] photon_position = [] photon_direction = [] + camera_intensity = [] + if isinstance(_device, list): # per convention this list has at least two elements results = self.forward_model(absorption_cm=absorption, scattering_cm=scattering, anisotropy=anisotropy, + refractive_index=refractive_index, illumination_geometry=_device[0]) self._append_results(results=results, reflectance=reflectance, reflectance_position=reflectance_position, photon_position=photon_position, - photon_direction=photon_direction) + photon_direction=photon_direction, + camera_intensity=camera_intensity) + fluence = results[Tags.DATA_FIELD_FLUENCE] for idx in range(1, len(_device)): # we already looked at the 0th element, so go from 1 to n-1 results = self.forward_model(absorption_cm=absorption, scattering_cm=scattering, anisotropy=anisotropy, + refractive_index=refractive_index, illumination_geometry=_device[idx]) self._append_results(results=results, reflectance=reflectance, reflectance_position=reflectance_position, photon_position=photon_position, - photon_direction=photon_direction) - fluence += results[Tags.DATA_FIELD_FLUENCE] + photon_direction=photon_direction, + camera_intensity=camera_intensity) + fluence += results[Tags.DATA_FIELD_FLUENCE] fluence = fluence / len(_device) else: results = self.forward_model(absorption_cm=absorption, scattering_cm=scattering, anisotropy=anisotropy, + refractive_index=refractive_index, illumination_geometry=_device) self._append_results(results=results, reflectance=reflectance, reflectance_position=reflectance_position, photon_position=photon_position, - photon_direction=photon_direction) + photon_direction=photon_direction, + camera_intensity=camera_intensity) + fluence = results[Tags.DATA_FIELD_FLUENCE] aggregated_results = dict() aggregated_results[Tags.DATA_FIELD_FLUENCE] = fluence @@ -288,17 +369,116 @@ def run_forward_model(self, if photon_position: aggregated_results[Tags.DATA_FIELD_PHOTON_EXIT_POS] = np.concatenate(photon_position, axis=0) aggregated_results[Tags.DATA_FIELD_PHOTON_EXIT_DIR] = np.concatenate(photon_direction, axis=0) + if camera_intensity: + aggregated_results[Tags.DATA_FIELD_CAMERA_INTENSITY] = np.concatenate(camera_intensity, axis=0) return aggregated_results @staticmethod def _append_results(results, reflectance, reflectance_position, - photon_position, - photon_direction): + photon_position: list[np.ndarray], + photon_direction: list[np.ndarray], + camera_intensity: list[np.ndarray]): if Tags.DATA_FIELD_DIFFUSE_REFLECTANCE in results: reflectance.append(results[Tags.DATA_FIELD_DIFFUSE_REFLECTANCE]) reflectance_position.append(results[Tags.DATA_FIELD_DIFFUSE_REFLECTANCE_POS]) if Tags.DATA_FIELD_PHOTON_EXIT_POS in results: photon_position.append(results[Tags.DATA_FIELD_PHOTON_EXIT_POS]) photon_direction.append(results[Tags.DATA_FIELD_PHOTON_EXIT_DIR]) + if Tags.DATA_FIELD_CAMERA_INTENSITY in results: + camera_intensity.append(results[Tags.DATA_FIELD_CAMERA_INTENSITY]) + + +class FastMCXReflectanceAdapter(MCXReflectanceAdapter): + """ + Same functionality of MCXReflectanceAdapter, but saves the simulation output through an event instead + of saving it in the HDF output file. + + Note: Does not include fluence. + """ + + def __init__(self, global_settings: Settings): + """ + initializes MCX-specific configuration and clean-up instances + + :param global_settings: global settings used during simulations + """ + super().__init__(global_settings=global_settings) + + self.simulation_finished_event_listeners = [] + + def read_mcx_output(self, **kwargs) -> Dict: + """ + reads the temporary output generated with MCX + + :param kwargs: dummy, used for class inheritance compatibility + :return: `Settings` instance containing the MCX output + """ + results = dict() + if os.path.isfile(self.mcx_volumetric_data_file) and self.mcx_volumetric_data_file.endswith( + self.mcx_output_suffixes['mcx_volumetric_data_file']): + if Tags.COMPUTE_DIFFUSE_REFLECTANCE in self.component_settings and \ + self.component_settings[Tags.COMPUTE_DIFFUSE_REFLECTANCE]: + content = jdata.load(self.mcx_volumetric_data_file) + fluence = content['NIFTIData'] + + if fluence.ndim > 3: + # remove the 1 or 2 (for mcx >= v2024.1) additional dimensions of size 1 if present to obtain a 3d array + fluence = fluence.reshape(fluence.shape[0], fluence.shape[1], -1) + + ref, ref_pos, _ = self.extract_reflectance_from_fluence(fluence=fluence) + results[Tags.DATA_FIELD_DIFFUSE_REFLECTANCE] = ref + results[Tags.DATA_FIELD_DIFFUSE_REFLECTANCE_POS] = ref_pos + else: + raise FileNotFoundError( + f"Could not find .jnii file for {self.mcx_volumetric_data_file}, something went wrong!") + + if (Tags.MCX_CAMERA_SETTINGS in self.global_settings) or (Tags.MCX_BACKTRACK_SETTINGS in self.global_settings): + cam_intensity_file_path = pathlib.Path(self.mcx_volumetric_data_file).with_suffix(".bin") + cam_intensity = np.fromfile(cam_intensity_file_path, dtype=np.float32) + results[Tags.DATA_FIELD_CAMERA_INTENSITY] = cam_intensity + + if Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT in self.component_settings and \ + self.component_settings[Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT]: + content = jdata.load(self.mcx_photon_data_file) + photon_pos = content['MCXData']['PhotonData']['p'] + photon_dir = content['MCXData']['PhotonData']['v'] + results[Tags.DATA_FIELD_PHOTON_EXIT_POS] = photon_pos + results[Tags.DATA_FIELD_PHOTON_EXIT_DIR] = photon_dir + return results + + def run(self, device: Union[IlluminationGeometryBase, PhotoacousticDevice]) -> None: + """ + runs optical simulations. Volumes are first loaded from HDF5 file and parsed to `self.forward_model`, the output + is aggregated in case multiple illuminations are defined by `device` and stored in the same HDF5 file. + + :param device: Illumination or Photoacoustic device that defines the illumination geometry + :return: None + """ + assert Tags.IGNORE_QA_ASSERTIONS not in self.global_settings, "Not implemented" + assert Tags.LASER_PULSE_ENERGY_IN_MILLIJOULE not in self.component_settings, "Not implemented" + + self.logger.info("Simulating the optical forward process...") + + file_path = self.global_settings[Tags.SIMPA_OUTPUT_FILE_PATH] + wl = str(self.global_settings[Tags.WAVELENGTH]) + + absorption = load_data_field(file_path, Tags.DATA_FIELD_ABSORPTION_PER_CM, wl) + scattering = load_data_field(file_path, Tags.DATA_FIELD_SCATTERING_PER_CM, wl) + anisotropy = load_data_field(file_path, Tags.DATA_FIELD_ANISOTROPY, wl) + refractive_index = load_data_field(file_path, Tags.DATA_FIELD_REFRACTIVE_INDEX, wl) + + assert isinstance(device, IlluminationGeometryBase), "Not implemented" + + results = self.forward_model(absorption_cm=absorption, + scattering_cm=scattering, + anisotropy=anisotropy, + refractive_index=refractive_index, + illumination_geometry=device) + + assert Tags.DATA_FIELD_FLUENCE not in results + self.logger.info("Simulating the optical forward process...[Done]") + + for event_listener in self.simulation_finished_event_listeners: + event_listener(self.global_settings[Tags.WAVELENGTH], results) diff --git a/simpa/core/simulation_modules/optical_module/optical_adapter_base.py b/simpa/core/simulation_modules/optical_module/optical_adapter_base.py index bc231a98..13473001 100644 --- a/simpa/core/simulation_modules/optical_module/optical_adapter_base.py +++ b/simpa/core/simulation_modules/optical_module/optical_adapter_base.py @@ -2,7 +2,7 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT from abc import abstractmethod -from typing import Dict, Union +from typing import Dict, Union, Iterable import numpy as np @@ -42,6 +42,7 @@ def forward_model(self, absorption_cm: np.ndarray, scattering_cm: np.ndarray, anisotropy: np.ndarray, + refractive_index: np.ndarray, illumination_geometry: IlluminationGeometryBase): """ A deriving class needs to implement this method according to its model. @@ -49,6 +50,7 @@ def forward_model(self, :param absorption_cm: Absorption in units of per centimeter :param scattering_cm: Scattering in units of per centimeter :param anisotropy: Dimensionless scattering anisotropy + :param refractive_index: Refractive index :param illumination_geometry: A device that represents a detection geometry :return: Fluence in units of J/cm^2 """ @@ -71,6 +73,7 @@ def run(self, device: Union[IlluminationGeometryBase, PhotoacousticDevice]) -> N absorption = load_data_field(file_path, Tags.DATA_FIELD_ABSORPTION_PER_CM, wl) scattering = load_data_field(file_path, Tags.DATA_FIELD_SCATTERING_PER_CM, wl) anisotropy = load_data_field(file_path, Tags.DATA_FIELD_ANISOTROPY, wl) + refractive_index = load_data_field(file_path, Tags.DATA_FIELD_REFRACTIVE_INDEX, wl) gruneisen_parameter = load_data_field(file_path, Tags.DATA_FIELD_GRUNEISEN_PARAMETER) _device = None @@ -85,7 +88,8 @@ def run(self, device: Union[IlluminationGeometryBase, PhotoacousticDevice]) -> N device=device, absorption=absorption, scattering=scattering, - anisotropy=anisotropy) + anisotropy=anisotropy, + refractive_index=refractive_index) fluence = results[Tags.DATA_FIELD_FLUENCE] if not (Tags.IGNORE_QA_ASSERTIONS in self.global_settings and Tags.IGNORE_QA_ASSERTIONS): assert_array_well_defined(fluence, assume_non_negativity=True, array_name="fluence") @@ -120,7 +124,8 @@ def run_forward_model(self, device: Union[IlluminationGeometryBase, PhotoacousticDevice], absorption: np.ndarray, scattering: np.ndarray, - anisotropy: np.ndarray) -> Dict: + anisotropy: np.ndarray, + refractive_index: np.ndarray) -> Dict: """ runs `self.forward_model` as many times as defined by `device` and aggregates the results. @@ -129,29 +134,14 @@ def run_forward_model(self, :param absorption: Absorption volume :param scattering: Scattering volume :param anisotropy: Dimensionless scattering anisotropy + :param refractive_index: Refractive index :return: """ - if isinstance(_device, list): - # per convention this list has at least two elements - results = self.forward_model(absorption_cm=absorption, - scattering_cm=scattering, - anisotropy=anisotropy, - illumination_geometry=_device[0]) - fluence = results[Tags.DATA_FIELD_FLUENCE] - for idx in range(1, len(_device)): - # we already looked at the 0th element, so go from 1 to n-1 - results = self.forward_model(absorption_cm=absorption, - scattering_cm=scattering, - anisotropy=anisotropy, - illumination_geometry=_device[idx]) - fluence += results[Tags.DATA_FIELD_FLUENCE] - - fluence = fluence / len(_device) - - else: - results = self.forward_model(absorption_cm=absorption, + _devices = _device if isinstance(_device, Iterable) else (_device,) + fluence = sum(self.forward_model(absorption_cm=absorption, scattering_cm=scattering, anisotropy=anisotropy, - illumination_geometry=_device) - fluence = results[Tags.DATA_FIELD_FLUENCE] - return {Tags.DATA_FIELD_FLUENCE: fluence} + refractive_index=refractive_index, + illumination_geometry=illumination_geometry)[Tags.DATA_FIELD_FLUENCE] + for illumination_geometry in _devices) + return {Tags.DATA_FIELD_FLUENCE: fluence / len(_devices)} diff --git a/simpa/core/simulation_modules/optical_module/optical_test_adapter.py b/simpa/core/simulation_modules/optical_module/optical_test_adapter.py index bb1a7b49..10b16678 100644 --- a/simpa/core/simulation_modules/optical_module/optical_test_adapter.py +++ b/simpa/core/simulation_modules/optical_module/optical_test_adapter.py @@ -11,6 +11,6 @@ class OpticalTestAdapter(OpticalAdapterBase): This Adapter was created for testing purposes and only """ - def forward_model(self, absorption_cm, scattering_cm, anisotropy, illumination_geometry): + def forward_model(self, absorption_cm, scattering_cm, anisotropy, refractive_index, illumination_geometry): results = {Tags.DATA_FIELD_FLUENCE: absorption_cm / ((1 - anisotropy) * scattering_cm)} return results diff --git a/simpa/core/simulation_modules/optical_module/volume_boundary_condition.py b/simpa/core/simulation_modules/optical_module/volume_boundary_condition.py new file mode 100644 index 00000000..3955fd1e --- /dev/null +++ b/simpa/core/simulation_modules/optical_module/volume_boundary_condition.py @@ -0,0 +1,29 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT + +from enum import Enum + + +# region Public enums + +class MCXVolumeBoundaryCondition(Enum): + """Defines the behavior of the photons touching the volume boundaries. + + Sets the letters of the --bc command (https://mcx.space/wiki/index.cgi?Doc/mcx_help). + Note: The behavior is only defined on the volume faces in the x- and y-axis, the behavior for the faces + on the z-axis always remains the default. + """ + + DEFAULT = "______000000" + """The default behavior.""" + MIRROR_REFLECTION = "mm_mm_000000" + """The photons are totally reflected as if the volume faces are mirrors.""" + CYCLIC = "cc_cc_000000" + """The photons reenter from the opposite volume face.""" + ABSORB = "aa_aa_000000" + """The photons are fully absorbed.""" + FRESNEL_REFLECTION = "rr_rr_000000" + """The photons are reflected based on the Fresnel equations.""" + +# endregion diff --git a/simpa/utils/__init__.py b/simpa/utils/__init__.py index c8484709..7359d37b 100644 --- a/simpa/utils/__init__.py +++ b/simpa/utils/__init__.py @@ -18,6 +18,7 @@ from .libraries.spectrum_library import Spectrum from .libraries.spectrum_library import view_saved_spectra from .libraries.spectrum_library import AnisotropySpectrumLibrary +from .libraries.spectrum_library import RefractiveIndexSpectrumLibrary from .libraries.spectrum_library import ScatteringSpectrumLibrary from .libraries.spectrum_library import get_simpa_internal_absorption_spectra_by_names diff --git a/simpa/utils/calculate.py b/simpa/utils/calculate.py index b8fd6372..b0b13840 100644 --- a/simpa/utils/calculate.py +++ b/simpa/utils/calculate.py @@ -24,7 +24,7 @@ def extract_hemoglobin_fractions(molecule_list: List) -> Dict[str, float]: } for molecule in molecule_list: - spectrum_name = molecule.spectrum.spectrum_name + spectrum_name = molecule.absorption_spectrum.spectrum_name if spectrum_name in hemoglobin: hemoglobin[spectrum_name] = molecule.volume_fraction diff --git a/simpa/utils/constants.py b/simpa/utils/constants.py index b5b5b6d2..f49f4ef8 100644 --- a/simpa/utils/constants.py +++ b/simpa/utils/constants.py @@ -34,7 +34,8 @@ class SegmentationClasses: wavelength_dependent_properties = [ Tags.DATA_FIELD_ABSORPTION_PER_CM, Tags.DATA_FIELD_SCATTERING_PER_CM, - Tags.DATA_FIELD_ANISOTROPY + Tags.DATA_FIELD_ANISOTROPY, + Tags.DATA_FIELD_REFRACTIVE_INDEX ] wavelength_independent_properties = [ @@ -59,7 +60,8 @@ class SegmentationClasses: Tags.DATA_FIELD_DIFFUSE_REFLECTANCE, Tags.DATA_FIELD_DIFFUSE_REFLECTANCE_POS, Tags.DATA_FIELD_PHOTON_EXIT_POS, - Tags.DATA_FIELD_PHOTON_EXIT_DIR] + Tags.DATA_FIELD_PHOTON_EXIT_DIR, + Tags.DATA_FIELD_CAMERA_INTENSITY] simulation_output_fields = [Tags.OPTICAL_MODEL_OUTPUT_NAME, Tags.SIMULATION_PROPERTIES] diff --git a/simpa/utils/dict_path_manager.py b/simpa/utils/dict_path_manager.py index 7a47127d..a828b287 100644 --- a/simpa/utils/dict_path_manager.py +++ b/simpa/utils/dict_path_manager.py @@ -35,7 +35,7 @@ def generate_dict_path(data_field, wavelength: (int, float) = None) -> str: elif data_field in simulation_output: if data_field in [Tags.DATA_FIELD_FLUENCE, Tags.DATA_FIELD_INITIAL_PRESSURE, Tags.OPTICAL_MODEL_UNITS, Tags.DATA_FIELD_DIFFUSE_REFLECTANCE, Tags.DATA_FIELD_DIFFUSE_REFLECTANCE_POS, - Tags.DATA_FIELD_PHOTON_EXIT_POS, Tags.DATA_FIELD_PHOTON_EXIT_DIR]: + Tags.DATA_FIELD_PHOTON_EXIT_POS, Tags.DATA_FIELD_PHOTON_EXIT_DIR, Tags.DATA_FIELD_CAMERA_INTENSITY]: dict_path = "/" + Tags.SIMULATIONS + "/" + Tags.OPTICAL_MODEL_OUTPUT_NAME + "/" + data_field + wl else: dict_path = "/" + Tags.SIMULATIONS + "/" + data_field + wl @@ -50,7 +50,7 @@ def generate_dict_path(data_field, wavelength: (int, float) = None) -> str: elif data_field in wavelength_independent_image_processing_output: dict_path = "/" + Tags.IMAGE_PROCESSING + "/" + data_field + "/" else: - raise ValueError(f"The requested data_field: '{data_field}: is not a valid argument. " + raise ValueError(f"The requested data_field: '{data_field}': is not a valid argument. " f"Please specify a valid data_field using the Tags from simpa/utils/tags.py!") return dict_path diff --git a/simpa/utils/libraries/molecule_library.py b/simpa/utils/libraries/molecule_library.py index b73a1668..cd1bc0d2 100644 --- a/simpa/utils/libraries/molecule_library.py +++ b/simpa/utils/libraries/molecule_library.py @@ -1,6 +1,7 @@ # SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT +from typing import Optional, Union import numpy as np import torch @@ -8,11 +9,11 @@ from simpa.utils.tissue_properties import TissueProperties from simpa.utils.libraries.literature_values import OpticalTissueProperties, StandardProperties -from simpa.utils.libraries.spectrum_library import AnisotropySpectrumLibrary, ScatteringSpectrumLibrary +from simpa.utils.libraries.spectrum_library import (AnisotropySpectrumLibrary, ScatteringSpectrumLibrary, + RefractiveIndexSpectrumLibrary, AbsorptionSpectrumLibrary) from simpa.utils import Spectrum 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 @@ -79,7 +80,7 @@ def update_internal_properties(self, settings): 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, settings, wavelength) -> TissueProperties: + def get_properties_for_wavelength(self, settings, wavelength: Union[int, float]) -> TissueProperties: """ Get the tissue properties for a specific wavelength. @@ -90,15 +91,21 @@ def get_properties_for_wavelength(self, settings, wavelength) -> TissuePropertie 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 + self.internal_properties[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 0 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)) + (molecule.volume_fraction * molecule.absorption_spectrum.get_value_for_wavelength(wavelength)) + self.internal_properties[Tags.DATA_FIELD_SCATTERING_PER_CM] += \ (molecule.volume_fraction * (molecule.scattering_spectrum.get_value_for_wavelength(wavelength))) + self.internal_properties[Tags.DATA_FIELD_ANISOTROPY] += \ molecule.volume_fraction * molecule.anisotropy_spectrum.get_value_for_wavelength(wavelength) + self.internal_properties[Tags.DATA_FIELD_REFRACTIVE_INDEX] += \ + molecule.volume_fraction * molecule.refractive_index.get_value_for_wavelength(wavelength) + return self.internal_properties def serialize(self) -> dict: @@ -108,7 +115,7 @@ def serialize(self) -> dict: :return: The serialized molecular composition. """ dict_items = self.__dict__ - dict_items["internal_properties"] = None + dict_items["internal_properties"] = None # Todo: Explain why. list_items = [molecule for molecule in self] return {"MolecularComposition": {"dict_items": dict_items, "list_items": list_items}} @@ -149,17 +156,18 @@ def __init__(self, name: str = None, absorption_spectrum: Spectrum = None, volume_fraction: float = None, scattering_spectrum: Spectrum = None, - anisotropy_spectrum: Spectrum = None, gruneisen_parameter: float = None, + anisotropy_spectrum: Spectrum = None, + refractive_index: Spectrum = None, + gruneisen_parameter: float = None, density: float = None, speed_of_sound: float = None, alpha_coefficient: float = None): """ - Initialize the Molecule object. - :param name: The name of the molecule. :param absorption_spectrum: The absorption spectrum of the molecule. :param volume_fraction: The volume fraction of the molecule. :param scattering_spectrum: The scattering spectrum of the molecule. :param anisotropy_spectrum: The anisotropy spectrum of the molecule. + :param refractive_index: Spectrum :param gruneisen_parameter: The Grüneisen parameter of the molecule. :param density: The density of the molecule. :param speed_of_sound: The speed of sound in the molecule. @@ -181,7 +189,7 @@ def __init__(self, name: str = None, if not isinstance(absorption_spectrum, Spectrum): raise TypeError(f"The given spectrum was not of type AbsorptionSpectrum! Instead: " f"{type(absorption_spectrum)} and reads: {absorption_spectrum}") - self.spectrum = absorption_spectrum + self.absorption_spectrum = absorption_spectrum if volume_fraction is None: volume_fraction = 0.0 @@ -200,11 +208,17 @@ def __init__(self, name: str = None, self.scattering_spectrum = scattering_spectrum if anisotropy_spectrum is None: - anisotropy_spectrum = 0.0 + anisotropy_spectrum = AnisotropySpectrumLibrary.CONSTANT_ANISOTROPY_ARBITRARY(1) if not isinstance(anisotropy_spectrum, Spectrum): raise TypeError(f"The given anisotropy was not of type Spectrum instead of {type(anisotropy_spectrum)}!") self.anisotropy_spectrum = anisotropy_spectrum + if refractive_index is None: + refractive_index = RefractiveIndexSpectrumLibrary.CONSTANT_REFRACTOR_ARBITRARY(1) + if not isinstance(refractive_index, Spectrum): + raise TypeError(f"The refractive index was not of type Spectrum instead of {type(refractive_index)}!") + self.refractive_index = refractive_index + if gruneisen_parameter is None: gruneisen_parameter = calculate_gruneisen_parameter_from_temperature( StandardProperties.BODY_TEMPERATURE_CELCIUS) @@ -250,9 +264,10 @@ def __eq__(self, other) -> bool: """ if isinstance(other, Molecule): return (self.name == other.name and - self.spectrum == other.spectrum and + self.absorption_spectrum == other.absorption_spectrum and self.volume_fraction == other.volume_fraction and self.scattering_spectrum == other.scattering_spectrum and + self.refractive_index == other.refractive_index and self.alpha_coefficient == other.alpha_coefficient and self.speed_of_sound == other.speed_of_sound and self.gruneisen_parameter == other.gruneisen_parameter and @@ -284,13 +299,14 @@ def deserialize(dictionary_to_deserialize: dict): :return: The deserialized Molecule object. """ deserialized_molecule = Molecule(name=dictionary_to_deserialize["name"], - absorption_spectrum=dictionary_to_deserialize["spectrum"], + absorption_spectrum=dictionary_to_deserialize["absorption_spectrum"], volume_fraction=dictionary_to_deserialize["volume_fraction"], scattering_spectrum=dictionary_to_deserialize["scattering_spectrum"], alpha_coefficient=dictionary_to_deserialize["alpha_coefficient"], speed_of_sound=dictionary_to_deserialize["speed_of_sound"], gruneisen_parameter=dictionary_to_deserialize["gruneisen_parameter"], anisotropy_spectrum=dictionary_to_deserialize["anisotropy_spectrum"], + refractive_index=dictionary_to_deserialize["refractive_index"], density=dictionary_to_deserialize["density"]) return deserialized_molecule @@ -318,6 +334,7 @@ def water(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: StandardProperties.WATER_MUS), anisotropy_spectrum=AnisotropySpectrumLibrary.CONSTANT_ANISOTROPY_ARBITRARY( StandardProperties.WATER_G), + refractive_index=RefractiveIndexSpectrumLibrary().get_spectrum_by_name("Water"), density=StandardProperties.DENSITY_WATER, speed_of_sound=StandardProperties.SPEED_OF_SOUND_WATER, alpha_coefficient=StandardProperties.ALPHA_COEFF_WATER @@ -337,6 +354,7 @@ def oxyhemoglobin(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: scattering_spectrum=ScatteringSpectrumLibrary().get_spectrum_by_name("blood_scattering"), anisotropy_spectrum=AnisotropySpectrumLibrary.CONSTANT_ANISOTROPY_ARBITRARY( OpticalTissueProperties.BLOOD_ANISOTROPY), + refractive_index=RefractiveIndexSpectrumLibrary().get_spectrum_by_name("Oxyhemoglobin"), density=StandardProperties.DENSITY_BLOOD, speed_of_sound=StandardProperties.SPEED_OF_SOUND_BLOOD, alpha_coefficient=StandardProperties.ALPHA_COEFF_BLOOD @@ -356,6 +374,7 @@ def deoxyhemoglobin(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: scattering_spectrum=ScatteringSpectrumLibrary().get_spectrum_by_name("blood_scattering"), anisotropy_spectrum=AnisotropySpectrumLibrary.CONSTANT_ANISOTROPY_ARBITRARY( OpticalTissueProperties.BLOOD_ANISOTROPY), + refractive_index=RefractiveIndexSpectrumLibrary().get_spectrum_by_name("Deoxyhemoglobin"), density=StandardProperties.DENSITY_BLOOD, speed_of_sound=StandardProperties.SPEED_OF_SOUND_BLOOD, alpha_coefficient=StandardProperties.ALPHA_COEFF_BLOOD @@ -376,6 +395,8 @@ def melanin(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: "epidermis", OpticalTissueProperties.MUS500_EPIDERMIS, OpticalTissueProperties.FRAY_EPIDERMIS, OpticalTissueProperties.BMIE_EPIDERMIS), anisotropy_spectrum=AnisotropySpectrumLibrary().get_spectrum_by_name("Epidermis_Anisotropy"), + # for n: DOI:10.1371/journal.pone.0150268 + refractive_index=RefractiveIndexSpectrumLibrary.CONSTANT_REFRACTOR_ARBITRARY(1.36), density=StandardProperties.DENSITY_SKIN, speed_of_sound=StandardProperties.SPEED_OF_SOUND_SKIN, alpha_coefficient=StandardProperties.ALPHA_COEFF_SKIN @@ -395,6 +416,8 @@ def fat(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: scattering_spectrum=ScatteringSpectrumLibrary().get_spectrum_by_name("fat_scattering"), anisotropy_spectrum=AnisotropySpectrumLibrary.CONSTANT_ANISOTROPY_ARBITRARY( OpticalTissueProperties.STANDARD_ANISOTROPY), + # for n: DOI:10.1371/journal.pone.0150268 + refractive_index=RefractiveIndexSpectrumLibrary.CONSTANT_REFRACTOR_ARBITRARY(1.46), density=StandardProperties.DENSITY_FAT, speed_of_sound=StandardProperties.SPEED_OF_SOUND_FAT, alpha_coefficient=StandardProperties.ALPHA_COEFF_FAT @@ -418,6 +441,7 @@ def constant_scatterer(scattering_coefficient: float = 100.0, anisotropy: float scattering_spectrum=ScatteringSpectrumLibrary.CONSTANT_SCATTERING_ARBITRARY( scattering_coefficient), anisotropy_spectrum=AnisotropySpectrumLibrary.CONSTANT_ANISOTROPY_ARBITRARY(anisotropy), + refractive_index=RefractiveIndexSpectrumLibrary.CONSTANT_REFRACTOR_ARBITRARY(refractive_index), density=StandardProperties.DENSITY_GENERIC, speed_of_sound=StandardProperties.SPEED_OF_SOUND_GENERIC, alpha_coefficient=StandardProperties.ALPHA_COEFF_GENERIC @@ -437,6 +461,8 @@ def soft_tissue_scatterer(volume_fraction: (float, torch.Tensor) = 1.0) -> Molec scattering_spectrum=ScatteringSpectrumLibrary().get_spectrum_by_name("background_scattering"), anisotropy_spectrum=AnisotropySpectrumLibrary.CONSTANT_ANISOTROPY_ARBITRARY( OpticalTissueProperties.STANDARD_ANISOTROPY), + # for n: DOI:10.1371/journal.pone.0150268 + refractive_index=RefractiveIndexSpectrumLibrary.CONSTANT_REFRACTOR_ARBITRARY(1.39), density=StandardProperties.DENSITY_GENERIC, speed_of_sound=StandardProperties.SPEED_OF_SOUND_GENERIC, alpha_coefficient=StandardProperties.ALPHA_COEFF_GENERIC @@ -456,6 +482,8 @@ def muscle_scatterer(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: scattering_spectrum=ScatteringSpectrumLibrary().get_spectrum_by_name("muscle_scattering"), anisotropy_spectrum=AnisotropySpectrumLibrary.CONSTANT_ANISOTROPY_ARBITRARY( OpticalTissueProperties.STANDARD_ANISOTROPY), + # for n: DOI:10.1371/journal.pone.0150268 + refractive_index=RefractiveIndexSpectrumLibrary.CONSTANT_REFRACTOR_ARBITRARY(1.39), density=StandardProperties.DENSITY_GENERIC, speed_of_sound=StandardProperties.SPEED_OF_SOUND_GENERIC, alpha_coefficient=StandardProperties.ALPHA_COEFF_GENERIC @@ -476,6 +504,8 @@ def epidermal_scatterer(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecul "epidermis", OpticalTissueProperties.MUS500_EPIDERMIS, OpticalTissueProperties.FRAY_EPIDERMIS, OpticalTissueProperties.BMIE_EPIDERMIS), anisotropy_spectrum=AnisotropySpectrumLibrary().get_spectrum_by_name("Epidermis_Anisotropy"), + # for n: DOI:10.1371/journal.pone.0150268 + refractive_index=RefractiveIndexSpectrumLibrary.CONSTANT_REFRACTOR_ARBITRARY(1.36), density=StandardProperties.DENSITY_SKIN, speed_of_sound=StandardProperties.SPEED_OF_SOUND_SKIN, alpha_coefficient=StandardProperties.ALPHA_COEFF_SKIN @@ -489,6 +519,8 @@ def dermal_scatterer(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: :param volume_fraction: The volume fraction of the molecule, defaults to 1.0 :return: A Molecule object representing a dermal scatterer """ + + # Todo: Add refractive index. return Molecule(name="dermal_scatterer", absorption_spectrum=AbsorptionSpectrumLibrary().get_spectrum_by_name("Skin_Baseline"), volume_fraction=volume_fraction, @@ -517,6 +549,8 @@ def bone(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: scattering_spectrum=ScatteringSpectrumLibrary().get_spectrum_by_name("bone_scattering"), anisotropy_spectrum=AnisotropySpectrumLibrary.CONSTANT_ANISOTROPY_ARBITRARY( OpticalTissueProperties.STANDARD_ANISOTROPY), + # for n: DOI:10.1371/journal.pone.0150268 + refractive_index=RefractiveIndexSpectrumLibrary.CONSTANT_REFRACTOR_ARBITRARY(1.55), density=StandardProperties.DENSITY_BONE, speed_of_sound=StandardProperties.SPEED_OF_SOUND_BONE, alpha_coefficient=StandardProperties.ALPHA_COEFF_BONE @@ -536,6 +570,8 @@ def mediprene(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: scattering_spectrum=ScatteringSpectrumLibrary.CONSTANT_SCATTERING_ARBITRARY((-np.log(0.85)) - (-np.log(0.85) / 10)), anisotropy_spectrum=AnisotropySpectrumLibrary.CONSTANT_ANISOTROPY_ARBITRARY(0.9), + # for n: This is basically just a guess + refractive_index=RefractiveIndexSpectrumLibrary.CONSTANT_REFRACTOR_ARBITRARY(1.4), density=StandardProperties.DENSITY_GEL_PAD, speed_of_sound=StandardProperties.SPEED_OF_SOUND_GEL_PAD, alpha_coefficient=StandardProperties.ALPHA_COEFF_GEL_PAD @@ -557,6 +593,7 @@ def heavy_water(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: StandardProperties.WATER_MUS), anisotropy_spectrum=AnisotropySpectrumLibrary.CONSTANT_ANISOTROPY_ARBITRARY( StandardProperties.WATER_G), + refractive_index=RefractiveIndexSpectrumLibrary().get_spectrum_by_name("Heavy_Water"), density=StandardProperties.DENSITY_HEAVY_WATER, speed_of_sound=StandardProperties.SPEED_OF_SOUND_HEAVY_WATER, alpha_coefficient=StandardProperties.ALPHA_COEFF_WATER @@ -578,6 +615,7 @@ def air(volume_fraction: (float, torch.Tensor) = 1.0) -> Molecule: StandardProperties.AIR_MUS), anisotropy_spectrum=AnisotropySpectrumLibrary.CONSTANT_ANISOTROPY_ARBITRARY( StandardProperties.AIR_G), + refractive_index=RefractiveIndexSpectrumLibrary.CONSTANT_REFRACTOR_ARBITRARY(1), density=StandardProperties.DENSITY_AIR, speed_of_sound=StandardProperties.SPEED_OF_SOUND_AIR, alpha_coefficient=StandardProperties.ALPHA_COEFF_AIR diff --git a/simpa/utils/libraries/refractive_index_spectra_data/Deoxyhemoglobin.npz b/simpa/utils/libraries/refractive_index_spectra_data/Deoxyhemoglobin.npz new file mode 100644 index 00000000..f3fb492e Binary files /dev/null and b/simpa/utils/libraries/refractive_index_spectra_data/Deoxyhemoglobin.npz differ diff --git a/simpa/utils/libraries/refractive_index_spectra_data/Heavy_Water.npz b/simpa/utils/libraries/refractive_index_spectra_data/Heavy_Water.npz new file mode 100644 index 00000000..27548220 Binary files /dev/null and b/simpa/utils/libraries/refractive_index_spectra_data/Heavy_Water.npz differ diff --git a/simpa/utils/libraries/refractive_index_spectra_data/Oxyhemoglobin.npz b/simpa/utils/libraries/refractive_index_spectra_data/Oxyhemoglobin.npz new file mode 100644 index 00000000..4f14c951 Binary files /dev/null and b/simpa/utils/libraries/refractive_index_spectra_data/Oxyhemoglobin.npz differ diff --git a/simpa/utils/libraries/refractive_index_spectra_data/Water.npz b/simpa/utils/libraries/refractive_index_spectra_data/Water.npz new file mode 100644 index 00000000..70092236 Binary files /dev/null and b/simpa/utils/libraries/refractive_index_spectra_data/Water.npz differ diff --git a/simpa/utils/libraries/refractive_index_spectra_data/__init__.py b/simpa/utils/libraries/refractive_index_spectra_data/__init__.py new file mode 100644 index 00000000..db4827b0 --- /dev/null +++ b/simpa/utils/libraries/refractive_index_spectra_data/__init__.py @@ -0,0 +1,28 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT + +""" +Sources for refractive index library: + +- Water: + M. Daimon and A. Masumura. + Measurement of the refractive index of distilled water from the near-infrared region to the ultraviolet region, + Appl. Opt. 46, 3811-3820 (2007) + +- Heavy Water + S. Kedenburg, M. Vieweg, T. Gissibl, and H. Giessen. + Linear refractive index and absorption measurements of nonlinear optical liquids + in the visible and near-infrared spectral region, + Opt. Mat. Express 2, 1588-1611 (2012) + +- Blood + Wang J, Deng Z, Wang X, Ye Q, Zhou W, Mei J, Zhang C, Tian J. + Measurement of the refractive index of hemoglobin solutions for a continuous spectral region. + Biomedical Optics Express 6, 2536-41 (2015) + +- Air + A. Börzsönyi, Z. Heiner, M. P. Kalashnikov, A. P. Kovács, and K. Osvay, + Dispersion measurement of inert gases and gas mixtures at 800 nm, + Appl. Opt. 47, 4856-4863 (2008) +""" \ No newline at end of file diff --git a/simpa/utils/libraries/spectrum_library.py b/simpa/utils/libraries/spectrum_library.py index 6a1af7ee..97a300aa 100644 --- a/simpa/utils/libraries/spectrum_library.py +++ b/simpa/utils/libraries/spectrum_library.py @@ -35,22 +35,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) + assert isinstance(wavelengths, np.ndarray), type(wavelengths) + assert isinstance(values, np.ndarray), type(values) + self.spectrum_name = spectrum_name self.wavelengths = wavelengths - self.max_wavelength = int(torch.max(wavelengths)) - self.min_wavelength = int(torch.min(wavelengths)) + self.max_wavelength = int(np.floor(np.max(wavelengths))) + self.min_wavelength = int(np.ceil(np.min(wavelengths))) self.values = values - if torch.Tensor.size(wavelengths) != torch.Tensor.size(values): + if wavelengths.shape != values.shape: raise ValueError("The shape of the wavelengths and the values did not match: " + - str(torch.Tensor.size(wavelengths)) + " vs " + str(torch.Tensor.size(values))) + str(wavelengths.shape) + " vs " + str(values.shape)) - 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) + new_wavelengths = np.arange(self.min_wavelength, self.max_wavelength + 1, 1) + values_by_wavelength_function = interpolate.interp1d(self.wavelengths, self.values) + self.values_interp = values_by_wavelength_function(new_wavelengths) def get_value_over_wavelength(self) -> np.ndarray: """ @@ -60,7 +60,7 @@ def get_value_over_wavelength(self) -> np.ndarray: """ return np.asarray([self.wavelengths, self.values]) - def get_value_for_wavelength(self, wavelength: int) -> float: + def get_value_for_wavelength(self, wavelength: int | np.ndarray) -> float: """ Retrieves the interpolated value for a given wavelength within the spectrum range. @@ -70,10 +70,10 @@ def get_value_for_wavelength(self, wavelength: int) -> float: :return: the best matching linearly interpolated values for the given wavelength. :raises ValueError: if the given wavelength is not within the range of the spectrum. """ - if wavelength < self.min_wavelength or wavelength > self.max_wavelength: + if np.min(wavelength) < self.min_wavelength or np.max(wavelength) > self.max_wavelength: raise ValueError(f"The given wavelength ({wavelength}) is not within the range of the spectrum " f"({self.min_wavelength} - {self.max_wavelength})") - return self.values_interp[wavelength-self.min_wavelength] + return self.values_interp[wavelength - self.min_wavelength] def __eq__(self, other): """ @@ -238,18 +238,33 @@ def scattering_from_rayleigh_and_mie_theory(name: str, mus_at_500_nm: float = 1. fraction_rayleigh_scattering: float = 0.0, mie_power_law_coefficient: float = 0.0) -> Spectrum: """ - Creates a scattering spectrum based on Rayleigh and Mie scattering theory. + Creates a reduced scattering spectrum based on Rayleigh and Mie scattering theory. :param name: The name of the spectrum. - :param mus_at_500_nm: Scattering coefficient at 500 nm. + :param mus_at_500_nm: Reduced scattering coefficient at 500 nm. :param fraction_rayleigh_scattering: Fraction of Rayleigh scattering. :param mie_power_law_coefficient: Power law coefficient for Mie scattering. :return: A Spectrum instance based on Rayleigh and Mie scattering theory. """ - wavelengths = np.arange(450, 1001, 1) - scattering = (mus_at_500_nm * (fraction_rayleigh_scattering * (wavelengths / 500) ** 1e-4 + - (1 - fraction_rayleigh_scattering) * (wavelengths / 500) ** -mie_power_law_coefficient)) - return Spectrum(name, wavelengths, scattering) + wavelengths = np.arange(400, 1301, 1) + reduced_scattering = (mus_at_500_nm * (fraction_rayleigh_scattering * (wavelengths / 500) ** -4 + + (1 - fraction_rayleigh_scattering) * (wavelengths / 500) ** -mie_power_law_coefficient)) + return Spectrum(name, wavelengths, reduced_scattering) + + @staticmethod + def scattering_from_scattering_power(name: str, mus_at_500_nm: float, scattering_power: float) -> Spectrum: + """ + Creates a reduced scattering spectrum based on the first reduced scattering formula of the paper. + + :param name: The name of the spectrum. + :param mus_at_500_nm: Reduced scattering coefficient at 500 nm. Corresponds to a in the formula. + :param scattering_power: The scattering power. Corresponds to b in the formula. + + :return: The reduced scattering coefficients by wavelengths as a Spectrum instance. + """ + wavelengths = np.arange(400, 1301, 1) + reduced_scattering = mus_at_500_nm * np.power(wavelengths / 500, -scattering_power) + return Spectrum(name, wavelengths, reduced_scattering) class AbsorptionSpectrumLibrary(SpectraLibrary): @@ -277,7 +292,18 @@ def CONSTANT_ABSORBER_ARBITRARY(absorption_coefficient: float = 1) -> Spectrum: np.asarray([absorption_coefficient, absorption_coefficient])) -def get_simpa_internal_absorption_spectra_by_names(absorption_spectrum_names: list) -> list: +class RefractiveIndexSpectrumLibrary(SpectraLibrary): + + def __init__(self, additional_folder_path: str = None): + super(RefractiveIndexSpectrumLibrary, self).__init__("refractive_index_spectra_data", additional_folder_path) + + @staticmethod + def CONSTANT_REFRACTOR_ARBITRARY(refractive_index: float = 1): + return Spectrum("Constant Refractor (arb)", np.asarray([450, 1000]), + np.asarray([refractive_index, refractive_index])) + + +def get_simpa_internal_absorption_spectra_by_names(absorption_spectrum_names: list): """ Retrieves SIMPA internal absorption spectra by their names. @@ -296,7 +322,7 @@ def view_saved_spectra(save_path=None, mode="absorption"): Opens a matplotlib plot and visualizes the available spectra. :param save_path: If not None, then the figure will be saved as a PNG file to the destination. - :param mode: Specifies the type of spectra to visualize ("absorption", "scattering", or "anisotropy"). + :param mode: Specifies the type of spectra to visualize ("absorption", "scattering", "anisotropy" or "refractive_index). """ plt.figure(figsize=(11, 8)) if mode == "absorption": @@ -314,8 +340,14 @@ def view_saved_spectra(save_path=None, mode="absorption"): plt.semilogy(spectrum.wavelengths, spectrum.values, label=spectrum.spectrum_name) + elif mode == "refractive_index": + for spectrum in RefractiveIndexSpectrumLibrary(): + plt.semilogy(spectrum.wavelengths, + spectrum.values, + label=spectrum.spectrum_name) else: - raise ValueError(f"Invalid mode: {mode}. Choose from 'absorption', 'scattering', or 'anisotropy'.") + raise ValueError( + f"Invalid mode: {mode}. Choose from 'absorption', 'scattering', 'anisotropy' or 'refractive_index'.") ax = plt.gca() box = ax.get_position() diff --git a/simpa/utils/libraries/tissue_library.py b/simpa/utils/libraries/tissue_library.py index 8c4393e1..4c79b0df 100644 --- a/simpa/utils/libraries/tissue_library.py +++ b/simpa/utils/libraries/tissue_library.py @@ -9,9 +9,11 @@ from simpa.utils import MOLECULE_LIBRARY from simpa.utils import Spectrum from simpa.utils.libraries.molecule_library import MolecularComposition -from simpa.utils.libraries.spectrum_library import AnisotropySpectrumLibrary, ScatteringSpectrumLibrary +from simpa.utils.libraries.spectrum_library import (AnisotropySpectrumLibrary, ScatteringSpectrumLibrary, + RefractiveIndexSpectrumLibrary, AbsorptionSpectrumLibrary) from simpa.utils.calculate import randomize_uniform -from simpa.utils.libraries.spectrum_library import AbsorptionSpectrumLibrary + +import torch class TissueLibrary(object): @@ -19,8 +21,8 @@ class TissueLibrary(object): A library, returning molecular compositions for various typical tissue segmentations. """ - def get_blood_volume_fractions(self, oxygenation: Union[float, int, np.ndarray] = 1e-10, - blood_volume_fraction: Union[float, int, np.ndarray] = 1e-10)\ + def get_blood_volume_fractions(self, oxygenation: Union[float, int, np.ndarray, torch.Tensor] = 1e-10, + blood_volume_fraction: Union[float, int, np.ndarray, torch.Tensor] = 1e-10)\ -> List[Union[int, float, np.ndarray]]: """ A function that returns the volume fraction of the oxygenated and deoxygenated blood. @@ -32,8 +34,8 @@ def get_blood_volume_fractions(self, oxygenation: Union[float, int, np.ndarray] """ return [blood_volume_fraction*oxygenation, blood_volume_fraction*(1-oxygenation)] - 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: + def constant(self, mua: Union[float, int, np.ndarray, torch.Tensor] = 1e-10, mus: Union[float, int, np.ndarray, torch.Tensor] = 1e-10, + g: Union[float, int, np.ndarray, torch.Tensor] = 0, n: float = 1.3) -> 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. @@ -48,12 +50,14 @@ def constant(self, mua: Union[float, int, np.ndarray] = 1e-10, mus: Union[float, mua_as_spectrum = AbsorptionSpectrumLibrary().CONSTANT_ABSORBER_ARBITRARY(mua) mus_as_spectrum = ScatteringSpectrumLibrary.CONSTANT_SCATTERING_ARBITRARY(mus) g_as_spectrum = AnisotropySpectrumLibrary.CONSTANT_ANISOTROPY_ARBITRARY(g) - return self.generic_tissue(mua_as_spectrum, mus_as_spectrum, g_as_spectrum, "constant_mua_mus_g") + n_as_spectrum = RefractiveIndexSpectrumLibrary.CONSTANT_REFRACTOR_ARBITRARY(n) + return self.generic_tissue(mua_as_spectrum, mus_as_spectrum, g_as_spectrum, n_as_spectrum, "constant_mua_mus_g_n") 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), + n: Spectrum = RefractiveIndexSpectrumLibrary.CONSTANT_REFRACTOR_ARBITRARY(1.3), molecule_name: Optional[str] = "generic_tissue") -> MolecularComposition: """ Returns a generic issue defined by the provided optical parameters. @@ -61,6 +65,7 @@ def generic_tissue(self, :param mua: The absorption coefficient spectrum in cm^{-1}. :param mus: The scattering coefficient spectrum in cm^{-1}. :param g: The anisotropy spectrum. + :param n: The refractive index spectrum. :param molecule_name: The molecule name. :returns: The molecular composition of the tissue. @@ -68,17 +73,19 @@ def generic_tissue(self, assert isinstance(mua, Spectrum), type(mua) assert isinstance(mus, Spectrum), type(mus) assert isinstance(g, Spectrum), type(g) + assert isinstance(n, Spectrum), type(n) assert isinstance(molecule_name, str) or molecule_name is None, type(molecule_name) return (MolecularCompositionGenerator().append(Molecule(name=molecule_name, absorption_spectrum=mua, volume_fraction=1.0, scattering_spectrum=mus, - anisotropy_spectrum=g)) + anisotropy_spectrum=g, + refractive_index=n)) .get_molecular_composition(SegmentationClasses.GENERIC)) - def muscle(self, oxygenation: Union[float, int, np.ndarray] = 0.175, - blood_volume_fraction: Union[float, int, np.ndarray] = 0.06) -> MolecularComposition: + def muscle(self, oxygenation: Union[float, int, np.ndarray, torch.Tensor] = 0.175, + blood_volume_fraction: Union[float, int, np.ndarray, torch.Tensor] = 0.06) -> MolecularComposition: """ Create a molecular composition mimicking that of muscle :param oxygenation: The oxygenation level of the blood volume fraction (as a decimal). @@ -123,8 +130,8 @@ def muscle(self, oxygenation: Union[float, int, np.ndarray] = 0.175, .append(custom_water) .get_molecular_composition(SegmentationClasses.MUSCLE)) - 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: + def soft_tissue(self, oxygenation: Union[float, int, np.ndarray, torch.Tensor] = OpticalTissueProperties.BACKGROUND_OXYGENATION, + blood_volume_fraction: Union[float, int, np.ndarray, torch.Tensor] = 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). @@ -171,7 +178,7 @@ def soft_tissue(self, oxygenation: Union[float, int, np.ndarray] = OpticalTissue .append(custom_water) .get_molecular_composition(SegmentationClasses.SOFT_TISSUE)) - def epidermis(self, melanin_volume_fraction: Union[float, int, np.ndarray] = 0.014) -> MolecularComposition: + def epidermis(self, melanin_volume_fraction: Union[float, int, np.ndarray, torch.Tensor] = 0.014) -> MolecularComposition: """ Create a molecular composition mimicking that of dermis :param melanin_volume_fraction: the total volume fraction of melanin @@ -184,8 +191,8 @@ def epidermis(self, melanin_volume_fraction: Union[float, int, np.ndarray] = 0.0 .append(MOLECULE_LIBRARY.epidermal_scatterer(1 - melanin_volume_fraction)) .get_molecular_composition(SegmentationClasses.EPIDERMIS)) - def dermis(self, oxygenation: Union[float, int, np.ndarray] = 0.5, - blood_volume_fraction: Union[float, int, np.ndarray] = 0.002) -> MolecularComposition: + def dermis(self, oxygenation: Union[float, int, np.ndarray, torch.Tensor] = 0.5, + blood_volume_fraction: Union[float, int, np.ndarray, torch.Tensor] = 0.002) -> MolecularComposition: """ Create a molecular composition mimicking that of dermis :param oxygenation: The oxygenation level of the blood volume fraction (as a decimal). @@ -206,8 +213,9 @@ def dermis(self, oxygenation: Union[float, int, np.ndarray] = 0.5, .get_molecular_composition(SegmentationClasses.DERMIS)) def subcutaneous_fat(self, - oxygenation: Union[float, int, np.ndarray] = OpticalTissueProperties.BACKGROUND_OXYGENATION, - blood_volume_fraction: Union[float, int, np.ndarray] + oxygenation: Union[float, int, np.ndarray, + torch.Tensor] = OpticalTissueProperties.BACKGROUND_OXYGENATION, + blood_volume_fraction: Union[float, int, np.ndarray, torch.Tensor] = OpticalTissueProperties.BLOOD_VOLUME_FRACTION_MUSCLE_TISSUE) -> MolecularComposition: """ Create a molecular composition mimicking that of subcutaneous fat @@ -238,7 +246,7 @@ def subcutaneous_fat(self, .append(MOLECULE_LIBRARY.water(water_volume_fraction)) .get_molecular_composition(SegmentationClasses.FAT)) - def blood(self, oxygenation: Union[float, int, np.ndarray, None] = None) -> MolecularComposition: + def blood(self, oxygenation: Union[float, int, np.ndarray, torch.Tensor, None] = None) -> MolecularComposition: """ Create a molecular composition mimicking that of blood :param oxygenation: The oxygenation level of the blood(as a decimal). diff --git a/simpa/utils/settings.py b/simpa/utils/settings.py index 447ff48d..6f61e799 100644 --- a/simpa/utils/settings.py +++ b/simpa/utils/settings.py @@ -113,6 +113,60 @@ def set_optical_settings(self, optical_settings: dict): """ self[Tags.OPTICAL_MODEL_SETTINGS] = Settings(optical_settings) + def get_optical_camera_settings(self): + """" + Returns the settings for the optical forward model that are saved in this settings dictionary + """ + optical_camera_settings = self[Tags.OPTICAL_CAMERA_SETTINGS] + if isinstance(optical_camera_settings, Settings): + return optical_camera_settings + else: + return Settings(optical_camera_settings) + + def set_optical_camera_settings(self, optical_camera_settings: dict): + """ + Replaces the currently stored optical settings with the given dictionary + + :param optical_settings: a dictionary containing the optical settings + """ + self[Tags.OPTICAL_CAMERA_SETTINGS] = Settings(optical_camera_settings) + + def get_mcx_camera_settings(self): + """" + Returns the camera settings for MCX that are saved in this settings dictionary + """ + mcx_camera_settings = self[Tags.MCX_CAMERA_SETTINGS] + if isinstance(mcx_camera_settings, Settings): + return mcx_camera_settings + else: + return Settings(mcx_camera_settings) + + def set_mcx_camera_settings(self, mcx_camera_settings: dict): + """ + Replaces the currently stored optical settings with the given dictionary + + :param mcx_camera_settings: a dictionary containing the optical settings + """ + self[Tags.MCX_CAMERA_SETTINGS] = Settings(mcx_camera_settings) + + def get_mcx_backtrack_settings(self): + """" + Returns the camera settings for MCX that are saved in this settings dictionary + """ + mcx_backtrack_settings = self[Tags.MCX_BACKTRACK_SETTINGS] + if isinstance(mcx_backtrack_settings, Settings): + return mcx_backtrack_settings + else: + return Settings(mcx_backtrack_settings) + + def set_mcx_backtrack_settings(self, mcx_backtrack_settings: dict): + """ + Replaces the currently stored optical settings with the given dictionary + + :param mcx_backtrack_settings: a dictionary containing the optical settings + """ + self[Tags.MCX_BACKTRACK_SETTINGS] = Settings(mcx_backtrack_settings) + def get_volume_creation_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 c8e29025..342a2670 100644 --- a/simpa/utils/tags.py +++ b/simpa/utils/tags.py @@ -2,11 +2,10 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT +import numpy as np from numbers import Number from typing import Iterable -import numpy as np - class Tags: """ @@ -371,6 +370,92 @@ class Tags: Usage: module optical_simulation_module """ + OPTICAL_CAMERA_SETTINGS = ("optical_camera_settings", dict) + """ + Optical camera settings + """ + + OPTICAL_CAMERA_OBJECT_TO_LENS_DISTANCE = ("optical_camera_object_to_lens_distance", float | list) + """ + The distance between the object and the camera lens. + """ + + OPTICAL_CAMERA_LENS_TO_SENSOR_DISTANCE = ("optical_camera_lens_to_sensor_distance", float | list) + """ + The distance between the camera lens and the sensor. + """ + + OPTICAL_CAMERA_FOCAL_LENGTH = ("optical_camera_focal_length", float | list) + """ + The camera lens focal length. + """ + + OPTICAL_CAMERA_F_NUMBER = ("optical_camera_f_number", float | list) + """ + The camera f number. + """ + + MCX_CAMERA_SETTINGS = ("mcx_camera_settings", dict) + """ + MCX camera settings + """ + + MCX_BACKTRACK_SETTINGS = ("mcx_backtrack_settings", dict) + """ + MCX backtrack settings + """ + + MCX_OBJECT_DISTANCE = ("mcx_object_distance", float) + """ + The distance between the object and the camera lens. + """ + + MCX_PROJECTION_DISTANCE = ("mcx_projection_distance", float) + """ + The distance between the object and the camera lens. + """ + + MCX_IDEAL_DISTANCE = ("mcx_ideal_distance", float) + """ + The ideal distance between the object and the camera lens for perfect focus. Only required in backtrack mode. + """ + + MCX_FOCAL_LENGTH = ("mcx_focal_length", float) + """ + The distance between the object and the camera lens. + """ + + MCX_APERTURE_RADIUS = ("mcx_aperture_radius", float) + """ + The aperture radius. + """ + + MCX_TRUE_APERTURE_RADIUS = ("mcx_true_aperture_radius", float) + """ + The true aperture radius. Only required in backtrack mode which uses a larger aperture radius as the radius. + """ + + MMC_IMAGE_HEIGHT = ("mmc_image_height", int) + """ + The image height. Currently only used in MMC. + """ + + MMC_IMAGE_WIDTH = ("mmc_image_width", int) + """ + The image width. Currently only used in MMC. + """ + + MMC_PIXEL_PITCH = ("mmc_pixel_pitch", float) + """ + The pixel pitch. Currently only used in MMC. + """ + + MMC_RERUNS = ("mmc_reruns", int) + """ + How often to rerun MMC simulations to effectively simulate with much more photons. + Currently only used in MMCReflectanceAdapter. + """ + LASER_PULSE_ENERGY_IN_MILLIJOULE = ("laser_pulse_energy_in_millijoule", (int, np.integer, float, list, range, tuple, np.ndarray)) """ @@ -403,13 +488,6 @@ class Tags: Usage: module optical_modelling, adapter mcx_adapter """ - MCX_ASSUMED_ANISOTROPY = ("mcx_assumed_anisotropy", (int, float)) - """ - The anisotropy that should be assumed for the mcx simulations. - If not set, a default value of 0.9 will be assumed. - Usage: module optical_modelling, adapter mcx_adapter - """ - ILLUMINATION_TYPE = ("optical_model_illumination_type", str) """ Type of the illumination geometry used in mcx.\n @@ -875,6 +953,12 @@ class Tags: Usage: SIMPA package, naming convention """ + DATA_FIELD_REFRACTIVE_INDEX = "n" + """ + Refractive index of the generated volume/structure.\n + Usage: SIMPA package, naming convention + """ + DATA_FIELD_OXYGENATION = "oxy" """ Oxygenation of the generated volume/structure.\n @@ -951,6 +1035,12 @@ class Tags: Usage: SIMPA package """ + TRUE_SPACING_MM = ("true_voxel_spacing_mm", Number) + """ + Isotropic extent of one voxels in mm used in MCX. If not set, Tags.SPACING_MM will be used instead.\n + Usage: MCXReflectanceAdapter + """ + DIM_VOLUME_X_MM = ("volume_x_dim_mm", Number) """ Extent of the x-axis of the generated volume.\n @@ -1456,6 +1546,21 @@ class Tags: Usage: simpa.core.simulation_modules.optical_simulation_module.optical_forward_model_mcx_reflectance_adapter """ + VOLUME_BOUNDARY_CONDITION = "volume_boundary_condition" + """ + FIXME + """ + + MCX_DETECTOR = ("mcx_detector", dict) + """ + The detector property list used in mcx for capturing exiting photons.\n + Only use it if you want to capture photon exit properties (COMPUTE_PHOTON_DIRECTION_AT_EXIT is enabled). + + Example usage is: + + settings[Tags.MCX_DETECTOR] = [{"Pos": [30.0, 30.0, 0.0], "R": 45.0}] + """ + COMPUTE_PHOTON_DIRECTION_AT_EXIT = "save_dir_at_exit" """ Flag that indicates if the direction of photons when they exit the volume should be stored @@ -1468,6 +1573,12 @@ class Tags: Usage: simpa.core.simulation_modules.optical_simulation_module.optical_forward_model_mcx_reflectance_adapter """ + DATA_FIELD_CAMERA_INTENSITY = "camera_intensity" + """ + Identifier for the raw camera intensity returned by MCX + Usage: simpa.core.simulation_modules.optical_simulation_module.optical_forward_model_mcx_reflectance_adapter + """ + DATA_FIELD_DIFFUSE_REFLECTANCE_POS = "diffuse_reflectance_pos" """ Identified for the position within the volumes where the diffuse reflectance was originally stored, interface to diff --git a/simpa/visualisation/matplotlib_data_visualisation.py b/simpa/visualisation/matplotlib_data_visualisation.py index 29c5e299..9c132c2b 100644 --- a/simpa/visualisation/matplotlib_data_visualisation.py +++ b/simpa/visualisation/matplotlib_data_visualisation.py @@ -20,6 +20,7 @@ def visualise_data(wavelength: int = None, show_absorption=False, show_scattering=False, show_anisotropy=False, + show_refractive_index=False, show_speed_of_sound=False, show_tissue_density=False, show_fluence=False, @@ -64,6 +65,7 @@ def visualise_data(wavelength: int = None, absorption = get_data_field_from_simpa_output(file, Tags.DATA_FIELD_ABSORPTION_PER_CM, wavelength) scattering = get_data_field_from_simpa_output(file, Tags.DATA_FIELD_SCATTERING_PER_CM, wavelength) anisotropy = get_data_field_from_simpa_output(file, Tags.DATA_FIELD_ANISOTROPY, wavelength) + refractive_index = get_data_field_from_simpa_output(file, Tags.DATA_FIELD_REFRACTIVE_INDEX, wavelength) segmentation_map = get_data_field_from_simpa_output(file, Tags.DATA_FIELD_SEGMENTATION) speed_of_sound = get_data_field_from_simpa_output(file, Tags.DATA_FIELD_SPEED_OF_SOUND) density = get_data_field_from_simpa_output(file, Tags.DATA_FIELD_DENSITY) @@ -174,6 +176,11 @@ def visualise_data(wavelength: int = None, data_item_names.append("Anisotropy") cmaps.append("gray") logscales.append(True and log_scale) + if refractive_index is not None and show_refractive_index: + data_to_show.append(refractive_index) + data_item_names.append("Refractive Index") + cmaps.append("gray") + logscales.append(True and log_scale) if speed_of_sound is not None and show_speed_of_sound: data_to_show.append(speed_of_sound) data_item_names.append("Speed of Sound") diff --git a/simpa_tests/automatic_tests/tissue_library/test_tissue_library.py b/simpa_tests/automatic_tests/tissue_library/test_tissue_library.py index b36a734b..f3b698e5 100644 --- a/simpa_tests/automatic_tests/tissue_library/test_tissue_library.py +++ b/simpa_tests/automatic_tests/tissue_library/test_tissue_library.py @@ -13,20 +13,23 @@ def test_if_optical_parameter_spectra_are_provided_correct_tissue_definition_is_ mua_sample = np.linspace(1e-6, 1e-5, wavelengths_sample.shape[0]) mus_sample = np.linspace(1e-5, 5e-6, wavelengths_sample.shape[0]) g_sample = np.linspace(0.8, 0.9, wavelengths_sample.shape[0]) + n_sample = np.linspace(1.4, 0.8, wavelengths_sample.shape[0]) mua_spectrum = Spectrum("Mua", wavelengths_sample, mua_sample) mus_spectrum = Spectrum("Mus", wavelengths_sample, mus_sample) g_spectrum = Spectrum("g", wavelengths_sample, g_sample) + n_spectrum = Spectrum("n", wavelengths_sample, n_sample) - actual_tissue = TISSUE_LIBRARY.generic_tissue(mua_spectrum, mus_spectrum, g_spectrum) + actual_tissue = TISSUE_LIBRARY.generic_tissue(mua_spectrum, mus_spectrum, g_spectrum, n_spectrum) assert actual_tissue.segmentation_type == SegmentationClasses.GENERIC assert len(actual_tissue) == 1 molecule: Molecule = actual_tissue[0] assert molecule.volume_fraction == 1.0 - assert molecule.spectrum == mua_spectrum + assert molecule.absorption_spectrum == mua_spectrum assert molecule.scattering_spectrum == mus_spectrum assert molecule.anisotropy_spectrum == g_spectrum + assert molecule.refractive_index == n_spectrum def test_if_generic_tissue_is_called_with_invalid_arguments_error_is_raised(): @@ -34,47 +37,57 @@ def test_if_generic_tissue_is_called_with_invalid_arguments_error_is_raised(): mua_sample = np.linspace(1e-6, 1e-5, wavelengths_sample.shape[0]) mus_sample = np.linspace(1e-5, 5e-6, wavelengths_sample.shape[0]) g_sample = np.linspace(0.8, 0.9, wavelengths_sample.shape[0]) + n_sample = np.linspace(1.4, 0.8, wavelengths_sample.shape[0]) mua_spectrum = Spectrum("Mua", wavelengths_sample, mua_sample) mus_spectrum = Spectrum("Mus", wavelengths_sample, mus_sample) g_spectrum = Spectrum("g", wavelengths_sample, g_sample) + n_spectrum = Spectrum("n", wavelengths_sample, n_sample) + with pytest.raises(AssertionError): + TISSUE_LIBRARY.generic_tissue(None, mus_spectrum, g_spectrum, n_spectrum) with pytest.raises(AssertionError): - TISSUE_LIBRARY.generic_tissue(None, mus_spectrum, g_spectrum) + TISSUE_LIBRARY.generic_tissue(mua_spectrum, None, g_spectrum, n_spectrum) with pytest.raises(AssertionError): - TISSUE_LIBRARY.generic_tissue(mua_spectrum, None, g_spectrum) + TISSUE_LIBRARY.generic_tissue(mua_spectrum, mus_spectrum, None, n_spectrum) with pytest.raises(AssertionError): - TISSUE_LIBRARY.generic_tissue(mua_spectrum, mus_spectrum, None) + TISSUE_LIBRARY.generic_tissue(mua_spectrum, mus_spectrum, g_spectrum, None) def test_if_optical_parameter_spectra_are_provided_correct_tissue_definition_is_returned_from_constant(): mua_sample = 1e-5 mus_sample = 3e-6 g_sample = 0.85 + n_sample = 1.13 - actual_tissue = TISSUE_LIBRARY.constant(mua_sample, mus_sample, g_sample) + actual_tissue = TISSUE_LIBRARY.constant(mua_sample, mus_sample, g_sample, n_sample) assert actual_tissue.segmentation_type == SegmentationClasses.GENERIC assert len(actual_tissue) == 1 molecule: Molecule = actual_tissue[0] assert molecule.volume_fraction == 1.0 - assert (molecule.spectrum.values == mua_sample).all() + assert (molecule.absorption_spectrum.values == mua_sample).all() assert (molecule.scattering_spectrum.values == mus_sample).all() assert (molecule.anisotropy_spectrum.values == g_sample).all() + assert (molecule.refractive_index.values == n_sample).all() def test_if_constant_is_called_with_invalid_arguments_error_is_raised(): mua_sample = 1e-5 mus_sample = 3e-6 g_sample = 0.85 + n_sample = 1.13 + + with pytest.raises(TypeError): + TISSUE_LIBRARY.constant(None, mus_sample, g_sample, n_sample) with pytest.raises(TypeError): - TISSUE_LIBRARY.constant(None, mus_sample, g_sample) + TISSUE_LIBRARY.constant(mua_sample, None, g_sample, n_sample) with pytest.raises(TypeError): - TISSUE_LIBRARY.constant(mua_sample, None, g_sample) + TISSUE_LIBRARY.constant(mua_sample, mus_sample, None, n_sample) with pytest.raises(TypeError): - TISSUE_LIBRARY.constant(mua_sample, mus_sample, None) + TISSUE_LIBRARY.constant(mua_sample, mus_sample, g_sample, None) diff --git a/simpa_tests/manual_tests/image_reconstruction/PointSourceReconstruction.py b/simpa_tests/manual_tests/image_reconstruction/PointSourceReconstruction.py index dc17dea2..83cc94b2 100644 --- a/simpa_tests/manual_tests/image_reconstruction/PointSourceReconstruction.py +++ b/simpa_tests/manual_tests/image_reconstruction/PointSourceReconstruction.py @@ -31,11 +31,12 @@ class PointSourceReconstruction(ReconstructionAlgorithmTestBaseClass): TODO """ + def __init__(self, speed_of_sound: float = 1470, volume_transducer_dim_in_mm: float = 90, volume_planar_dim_in_mm: float = 20, volume_height_in_mm: float = 90, spacing: float = 0.4): - - self.reconstructed_image_pipeline = None # TODO REMOVE + + self.reconstructed_image_pipeline = None # TODO REMOVE self.SPEED_OF_SOUND = speed_of_sound self.VOLUME_TRANSDUCER_DIM_IN_MM = volume_transducer_dim_in_mm @@ -76,11 +77,12 @@ def create_point_source(self): vessel_1_dictionary = Settings() vessel_1_dictionary[Tags.PRIORITY] = 3 vessel_1_dictionary[Tags.STRUCTURE_START_MM] = [self.VOLUME_TRANSDUCER_DIM_IN_MM/2-10, 0, 35] - vessel_1_dictionary[Tags.STRUCTURE_END_MM] = [self.VOLUME_TRANSDUCER_DIM_IN_MM/2-10, self.VOLUME_PLANAR_DIM_IN_MM, 35] + vessel_1_dictionary[Tags.STRUCTURE_END_MM] = [ + self.VOLUME_TRANSDUCER_DIM_IN_MM/2-10, self.VOLUME_PLANAR_DIM_IN_MM, 35] vessel_1_dictionary[Tags.STRUCTURE_RADIUS_MM] = self.SPACING vessel_1_dictionary[Tags.MOLECULE_COMPOSITION] = (MolecularCompositionGenerator(). - append(vessel_molecule). - get_molecular_composition(-1)) + append(vessel_molecule). + get_molecular_composition(-1)) vessel_1_dictionary[Tags.CONSIDER_PARTIAL_VOLUME] = True vessel_1_dictionary[Tags.STRUCTURE_TYPE] = Tags.CIRCULAR_TUBULAR_STRUCTURE @@ -95,7 +97,6 @@ def setup(self): # point to the correct file in the PathManager(). self.path_manager = PathManager() - # 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. @@ -173,58 +174,56 @@ def setup(self): self.settings = settings def simulate_and_evaluate_with_device(self, _device): - SIMULATION_PIPELINE = [ - ModelBasedAdapter(self.settings), - MCXAdapter(self.settings), - KWaveAdapter(self.settings), - FieldOfViewCropping(self.settings), - DelayAndSumAdapter(self.settings) - ] - - print("Simulating for device:", _device) - simulate(SIMULATION_PIPELINE, self.settings, _device) - - if Tags.WAVELENGTH in self.settings: - wavelength = self.settings[Tags.WAVELENGTH] - else: - wavelength = 700 - - initial_pressure = load_data_field(self.path_manager.get_hdf5_file_save_path() + "/" + self.VOLUME_NAME + ".hdf5", - data_field=Tags.DATA_FIELD_INITIAL_PRESSURE, - wavelength=wavelength) - reconstruction = load_data_field(self.path_manager.get_hdf5_file_save_path() + "/" + self.VOLUME_NAME + ".hdf5", - data_field=Tags.DATA_FIELD_RECONSTRUCTED_DATA, - wavelength=wavelength) - - p0_idx = np.unravel_index(np.argmax(initial_pressure), np.shape(initial_pressure)) - re_idx = np.unravel_index(np.argmax(reconstruction), np.shape(reconstruction)) - - print("x/y in initial pressure map:", p0_idx) - print("x/y in reconstruction map:", re_idx) - distance = np.sqrt((re_idx[0] - p0_idx[0]) ** 2 + (re_idx[1] - p0_idx[1]) ** 2) - print("Distance:", distance) - - if self.save_path is not None: - save_path = self.save_path + f"PointSourceReconstruction_{self.figure_number}.png" - else: - save_path = self.save_path - - visualise_data(path_to_hdf5_file=self.path_manager.get_hdf5_file_save_path() + "/" + self.VOLUME_NAME + ".hdf5", - wavelength=wavelength, - show_time_series_data=True, - show_absorption=False, - show_reconstructed_data=True, - show_xz_only=True, - show_initial_pressure=True, - show_segmentation_map=False, - log_scale=False, - save_path=save_path) - self.figure_number += 1 - return distance - - - def test_reconstruction_of_simulation(self): - + SIMULATION_PIPELINE = [ + ModelBasedAdapter(self.settings), + MCXAdapter(self.settings), + KWaveAdapter(self.settings), + FieldOfViewCropping(self.settings), + DelayAndSumAdapter(self.settings) + ] + + print("Simulating for device:", _device) + simulate(SIMULATION_PIPELINE, self.settings, _device) + + if Tags.WAVELENGTH in self.settings: + wavelength = self.settings[Tags.WAVELENGTH] + else: + wavelength = 700 + + initial_pressure = load_data_field(self.path_manager.get_hdf5_file_save_path() + "/" + self.VOLUME_NAME + ".hdf5", + data_field=Tags.DATA_FIELD_INITIAL_PRESSURE, + wavelength=wavelength) + reconstruction = load_data_field(self.path_manager.get_hdf5_file_save_path() + "/" + self.VOLUME_NAME + ".hdf5", + data_field=Tags.DATA_FIELD_RECONSTRUCTED_DATA, + wavelength=wavelength) + + p0_idx = np.unravel_index(np.argmax(initial_pressure), np.shape(initial_pressure)) + re_idx = np.unravel_index(np.argmax(reconstruction), np.shape(reconstruction)) + + print("x/y in initial pressure map:", p0_idx) + print("x/y in reconstruction map:", re_idx) + distance = np.sqrt((re_idx[0] - p0_idx[0]) ** 2 + (re_idx[1] - p0_idx[1]) ** 2) + print("Distance:", distance) + + if self.save_path is not None: + save_path = self.save_path + f"PointSourceReconstruction_{self.figure_number}.png" + else: + save_path = self.save_path + + visualise_data(path_to_hdf5_file=self.path_manager.get_hdf5_file_save_path() + "/" + self.VOLUME_NAME + ".hdf5", + wavelength=wavelength, + show_time_series_data=True, + show_absorption=False, + show_reconstructed_data=True, + show_xz_only=True, + show_initial_pressure=True, + show_segmentation_map=False, + log_scale=False, + save_path=save_path) + self.figure_number += 1 + return distance + + def test_reconstruction_of_simulation(self): dist = list() @@ -248,52 +247,52 @@ def test_reconstruction_of_simulation(self): # seed=1234, field_of_view_extent_mm=fov_e)) # device.add_illumination_geometry(PencilBeamIlluminationGeometry()) # dist.append(self.simulate_and_evaluate_with_device(device)) - + dist.append(self.simulate_and_evaluate_with_device(MSOTAcuityEcho(device_position_mm=np.array([self.VOLUME_TRANSDUCER_DIM_IN_MM/2, - self.VOLUME_PLANAR_DIM_IN_MM/2, - 35]), - field_of_view_extent_mm=np.array([-(2 * np.sin(0.34 / 40 * 128) * 40) / 2, - (2 * np.sin(0.34 / - 40 * 128) * 40) / 2, - 0, 0, -25, 25])))) + self.VOLUME_PLANAR_DIM_IN_MM/2, + 35]), + field_of_view_extent_mm=np.array([-(2 * np.sin(0.34 / 40 * 128) * 40) / 2, + (2 * np.sin(0.34 / + 40 * 128) * 40) / 2, + 0, 0, -25, 25])))) dist.append(self.simulate_and_evaluate_with_device(InVision256TF(device_position_mm=np.array([self.VOLUME_TRANSDUCER_DIM_IN_MM/2, - self.VOLUME_PLANAR_DIM_IN_MM/2, - self.VOLUME_HEIGHT_IN_MM/2])))) + self.VOLUME_PLANAR_DIM_IN_MM/2, + self.VOLUME_HEIGHT_IN_MM/2])))) device = PhotoacousticDevice(device_position_mm=np.array([self.VOLUME_TRANSDUCER_DIM_IN_MM/2, - self.VOLUME_PLANAR_DIM_IN_MM/2, - 30]), - field_of_view_extent_mm=np.asarray([-self.VOLUME_TRANSDUCER_DIM_IN_MM/2, + self.VOLUME_PLANAR_DIM_IN_MM/2, + 30]), + field_of_view_extent_mm=np.asarray([-self.VOLUME_TRANSDUCER_DIM_IN_MM/2, self.VOLUME_TRANSDUCER_DIM_IN_MM/2, 0, 0, 0, self.VOLUME_HEIGHT_IN_MM])) device.set_detection_geometry(LinearArrayDetectionGeometry(device_position_mm=device.device_position_mm, - pitch_mm=0.2, - number_detector_elements=256)) + pitch_mm=0.2, + number_detector_elements=256)) device.add_illumination_geometry(PencilBeamIlluminationGeometry(device_position_mm=device.device_position_mm)) dist.append(self.simulate_and_evaluate_with_device(device)) device = PhotoacousticDevice(device_position_mm=np.array([self.VOLUME_TRANSDUCER_DIM_IN_MM/2, - self.VOLUME_PLANAR_DIM_IN_MM/2, - 5]), - field_of_view_extent_mm=np.asarray([-self.VOLUME_TRANSDUCER_DIM_IN_MM / 2, + self.VOLUME_PLANAR_DIM_IN_MM/2, + 5]), + field_of_view_extent_mm=np.asarray([-self.VOLUME_TRANSDUCER_DIM_IN_MM / 2, self.VOLUME_TRANSDUCER_DIM_IN_MM / 2, 0, 0, 0, self.VOLUME_HEIGHT_IN_MM])) device.set_detection_geometry(LinearArrayDetectionGeometry(device_position_mm=device.device_position_mm, - pitch_mm=0.2, - number_detector_elements=256)) + pitch_mm=0.2, + number_detector_elements=256)) device.add_illumination_geometry(PencilBeamIlluminationGeometry()) dist.append(self.simulate_and_evaluate_with_device(device)) device = PhotoacousticDevice(device_position_mm=np.array([self.VOLUME_TRANSDUCER_DIM_IN_MM/2, - self.VOLUME_PLANAR_DIM_IN_MM/2, - 10]), - field_of_view_extent_mm=np.asarray([-self.VOLUME_TRANSDUCER_DIM_IN_MM / 2, + self.VOLUME_PLANAR_DIM_IN_MM/2, + 10]), + field_of_view_extent_mm=np.asarray([-self.VOLUME_TRANSDUCER_DIM_IN_MM / 2, self.VOLUME_TRANSDUCER_DIM_IN_MM / 2, 0, 0, 0, self.VOLUME_HEIGHT_IN_MM])) device.set_detection_geometry(LinearArrayDetectionGeometry(device_position_mm=device.device_position_mm, - pitch_mm=0.2, - number_detector_elements=256)) + pitch_mm=0.2, + number_detector_elements=256)) device.add_illumination_geometry(PencilBeamIlluminationGeometry()) dist.append(self.simulate_and_evaluate_with_device(device)) print("") @@ -319,6 +318,7 @@ def run_test(self, show_figure_on_screen=True, save_path=None): self.perform_test() self.tear_down() + if __name__ == '__main__': test = PointSourceReconstruction() test.run_test(show_figure_on_screen=True) diff --git a/simpa_tests/test_utils/tissue_composition_tests.py b/simpa_tests/test_utils/tissue_composition_tests.py index 7adc8ca0..9e8d080f 100644 --- a/simpa_tests/test_utils/tissue_composition_tests.py +++ b/simpa_tests/test_utils/tissue_composition_tests.py @@ -495,6 +495,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values450nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values450nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values450nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values450nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 1.373 values500nm = TissueProperties(TEST_SETTINGS) values500nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 112 @@ -507,6 +508,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values500nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values500nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values500nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values500nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 1.371 values550nm = TissueProperties(TEST_SETTINGS) values550nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 230 @@ -519,6 +521,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values550nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values550nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values550nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values550nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 1.370 values600nm = TissueProperties(TEST_SETTINGS) values600nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 17 @@ -531,6 +534,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values600nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values600nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values600nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values600nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 1.369 values650nm = TissueProperties(TEST_SETTINGS) values650nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 2 @@ -543,6 +547,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values650nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values650nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values650nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values650nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 1.368 values700nm = TissueProperties(TEST_SETTINGS) values700nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 1.6 @@ -555,6 +560,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values700nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values700nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values700nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values700nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 1.367 values750nm = TissueProperties(TEST_SETTINGS) values750nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 2.8 @@ -567,6 +573,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values750nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values750nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values750nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values750nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 1.367 values800nm = TissueProperties(TEST_SETTINGS) values800nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 4.4 @@ -579,6 +586,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values800nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values800nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values800nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values800nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 1.366 values850nm = TissueProperties(TEST_SETTINGS) values850nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 5.7 @@ -591,6 +599,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values850nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values850nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values850nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values850nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 1.365 values900nm = TissueProperties(TEST_SETTINGS) values900nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 6.4 @@ -603,6 +612,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values900nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values900nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values900nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values900nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = None values950nm = TissueProperties(TEST_SETTINGS) values950nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 6.4 @@ -615,6 +625,7 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): values950nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values950nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values950nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values950nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = None if not only_use_NIR_values: reference_dict[450] = values450nm @@ -626,8 +637,8 @@ def get_fully_oxygenated_blood_reference_dictionary(only_use_NIR_values=False): reference_dict[750] = values750nm reference_dict[800] = values800nm reference_dict[850] = values850nm - reference_dict[900] = values900nm - reference_dict[950] = values950nm + # reference_dict[900] = values900nm # FIXME: Find refractive index values for 900 and 950 nm. + # reference_dict[950] = values950nm return reference_dict @@ -664,6 +675,8 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values450nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values450nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values450nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values450nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 0.2 + values450nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 1.374 values500nm = TissueProperties(TEST_SETTINGS) values500nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 112 @@ -676,6 +689,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values500nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values500nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values500nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values500nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 1.368 values550nm = TissueProperties(TEST_SETTINGS) values550nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 286 @@ -688,6 +702,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values550nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values550nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values550nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values550nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 1.366 values600nm = TissueProperties(TEST_SETTINGS) values600nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 79 @@ -700,6 +715,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values600nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values600nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values600nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values600nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 1.364 values650nm = TissueProperties(TEST_SETTINGS) values650nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 20.1 @@ -712,6 +728,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values650nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values650nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values650nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values650nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 1.362 values700nm = TissueProperties(TEST_SETTINGS) values700nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 9.6 @@ -724,6 +741,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values700nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values700nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values700nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values700nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 1.361 values750nm = TissueProperties(TEST_SETTINGS) values750nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 7.5 @@ -736,6 +754,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values750nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values750nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values750nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values750nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 1.360 values800nm = TissueProperties(TEST_SETTINGS) values800nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 4.1 @@ -748,6 +767,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values800nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values800nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values800nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values800nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 1.358 values850nm = TissueProperties(TEST_SETTINGS) values850nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 3.7 @@ -760,6 +780,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values850nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values850nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values850nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values850nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = 1.357 values900nm = TissueProperties(TEST_SETTINGS) values900nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 4.1 @@ -772,6 +793,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values900nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values900nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values900nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values900nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = None values950nm = TissueProperties(TEST_SETTINGS) values950nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 3.2 @@ -784,6 +806,7 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) values950nm[Tags.DATA_FIELD_DENSITY] = 1049.75 values950nm[Tags.DATA_FIELD_SPEED_OF_SOUND] = 1578.2 values950nm[Tags.DATA_FIELD_ALPHA_COEFF] = 0.2 + values950nm[Tags.DATA_FIELD_REFRACTIVE_INDEX] = None if not only_use_NIR_values: reference_dict[450] = values450nm @@ -795,8 +818,8 @@ def get_fully_deoxygenated_blood_reference_dictionary(only_use_NIR_values=False) reference_dict[750] = values750nm reference_dict[800] = values800nm reference_dict[850] = values850nm - reference_dict[900] = values900nm - reference_dict[950] = values950nm + # reference_dict[900] = values900nm # FIXME: Find refractive index values for 900 and 950 nm. + # reference_dict[950] = values950nm return reference_dict