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..34e5f400 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.utils + :members: + :undoc-members: + :show-inheritance: diff --git a/simpa/core/__init__.py b/simpa/core/__init__.py index 8653966c..9a49b941 100644 --- a/simpa/core/__init__.py +++ b/simpa/core/__init__.py @@ -8,10 +8,12 @@ from simpa.utils import Settings from simpa.utils.processing_device import get_processing_device + class PipelineModule: """ Defines a pipeline module (either simulation or processing module) that implements a run method and can be called by running the pipeline's simulate method. """ + def __init__(self, global_settings: Settings): """ :param global_settings: The SIMPA settings dictionary @@ -29,4 +31,3 @@ def run(self, digital_device_twin: DigitalDeviceTwinBase): :param digital_device_twin: The digital twin that can be used by the digital device_twin. """ pass - \ No newline at end of file diff --git a/simpa/core/processing_components/__init__.py b/simpa/core/processing_components/__init__.py index c8306cbf..881a8f2d 100644 --- a/simpa/core/processing_components/__init__.py +++ b/simpa/core/processing_components/__init__.py @@ -19,4 +19,3 @@ def __init__(self, global_settings, component_settings_key: str): """ super(ProcessingComponent, self).__init__(global_settings=global_settings) self.component_settings = global_settings[component_settings_key] - diff --git a/simpa/core/simulation_modules/__init__.py b/simpa/core/simulation_modules/__init__.py index 872aa0c1..c87e1e26 100644 --- a/simpa/core/simulation_modules/__init__.py +++ b/simpa/core/simulation_modules/__init__.py @@ -7,6 +7,7 @@ from simpa.core import PipelineModule from simpa.utils import Settings + class SimulationModule(PipelineModule): """ Defines a simulation module that is a step in the simulation pipeline. diff --git a/simpa/core/simulation_modules/acoustic_forward_module/__init__.py b/simpa/core/simulation_modules/acoustic_forward_module/__init__.py index 2beb0eb6..e9f33d10 100644 --- a/simpa/core/simulation_modules/acoustic_forward_module/__init__.py +++ b/simpa/core/simulation_modules/acoustic_forward_module/__init__.py @@ -33,7 +33,7 @@ class AcousticForwardModelBaseAdapter(SimulationModule): def __init__(self, global_settings: Settings): super(AcousticForwardModelBaseAdapter, self).__init__(global_settings=global_settings) - + def load_component_settings(self) -> Settings: """Implements abstract method to serve acoustic settings as component settings diff --git a/simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_mcx_adapter.py b/simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_mcx_adapter.py index da9b52d6..92a64efb 100644 --- a/simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_mcx_adapter.py +++ b/simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_mcx_adapter.py @@ -16,7 +16,7 @@ class MCXAdapter(OpticalForwardModuleBase): """ This class implements a bridge to the mcx framework to integrate mcx into SIMPA. This adapter only allows for - computation of fluence, for computations of diffuse reflectance, take a look at `simpa.ReflectanceMcxAdapter` + computation of fluence, for computations of diffuse reflectance, take a look at `simpa.MCXAdapterReflectance` .. note:: MCX is a GPU-enabled Monte-Carlo model simulation of photon transport in tissue: diff --git a/simpa/core/simulation_modules/optical_simulation_module/utils.py b/simpa/core/simulation_modules/optical_simulation_module/utils.py new file mode 100644 index 00000000..d7d8f8f7 --- /dev/null +++ b/simpa/core/simulation_modules/optical_simulation_module/utils.py @@ -0,0 +1,64 @@ +# 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 + +from simpa import load_data_field, Tags + + +def get_spectral_image_from_optical_simulation(simulation_hdf_file_path: typing.Union[str, pathlib.Path], + target_height: int, + target_width: int) -> np.ndarray: + """ + Returns the spectral image from the simulated reflectance data. + :param simulation_hdf_file_path: The path to the file containing the MCX simulation results. + :param target_height: The height of the spectral image in voxels. + :param target_width: The width of the spectral image in voxels. + :return: The spectral image in H x W x C. + """ + refl_by_wavelength = load_data_field(simulation_hdf_file_path, Tags.DATA_FIELD_DIFFUSE_REFLECTANCE) + refl_pos_by_wavelength = load_data_field(simulation_hdf_file_path, Tags.DATA_FIELD_DIFFUSE_REFLECTANCE_POS) + + assert refl_by_wavelength.keys() == refl_pos_by_wavelength.keys(), \ + (f"The reflectance values contain different wavelengths\n({refl_by_wavelength.keys()})\n" + f"than the reflectance positions\n({refl_pos_by_wavelength.keys()})\n!") + wavelengths_in_nm = [int(w) for w in refl_by_wavelength.keys()] + + return get_spectral_image_from_simulated_reflectances(wavelengths_in_nm=wavelengths_in_nm, + refl_by_wavelength=refl_by_wavelength, + refl_pos_by_wavelength=refl_pos_by_wavelength, + target_height=target_height, + target_width=target_width) + + +def get_spectral_image_from_simulated_reflectances(wavelengths_in_nm: typing.Union[np.ndarray, list[int]], + refl_by_wavelength: dict[str, np.ndarray], + refl_pos_by_wavelength: dict[str, np.ndarray], + target_height: int, target_width: int) -> np.ndarray: + """ + Returns the spectral image from the simulated reflectance data. + :param wavelengths_in_nm: The wavelengths in mm. + :param refl_by_wavelength: The reflectance values by wavelength from the simulation. + :param refl_pos_by_wavelength: The reflectance positions in [x, y, z] by wavelength from the simulation. + :param target_height: The height of the spectral image in voxels. + :param target_width: The width of the spectral image in voxels. + + :return: The spectral image in H x W x C. + """ + spectral_image = np.zeros((target_height, target_width, len(wavelengths_in_nm))) + + for j, wavelength in enumerate(wavelengths_in_nm): + reflectances = refl_by_wavelength[str(wavelength)] + reflectance_positions = refl_pos_by_wavelength[str(wavelength)] + assert reflectance_positions.shape[0] == reflectances.shape[0] and reflectance_positions.shape[ + 1] == 3, reflectance_positions.shape + reflectance_positions[:, 2] = j # Width x Height x Channel + reflectance_positions = reflectance_positions[:, [1, 0, 2]] + spectral_image[reflectance_positions[:, 0], reflectance_positions[:, 1], + reflectance_positions[:, 2]] = reflectances + + return spectral_image diff --git a/simpa/core/simulation_modules/reconstruction_module/__init__.py b/simpa/core/simulation_modules/reconstruction_module/__init__.py index 118874dc..05f8b3ef 100644 --- a/simpa/core/simulation_modules/reconstruction_module/__init__.py +++ b/simpa/core/simulation_modules/reconstruction_module/__init__.py @@ -25,7 +25,7 @@ class ReconstructionAdapterBase(SimulationModule): def __init__(self, global_settings: Settings): super(ReconstructionAdapterBase, self).__init__(global_settings=global_settings) - + def load_component_settings(self) -> Settings: """Implements abstract method to serve reconstruction settings as component settings diff --git a/simpa/core/simulation_modules/volume_creation_module/__init__.py b/simpa/core/simulation_modules/volume_creation_module/__init__.py index 7f2b3ffa..6b9f0753 100644 --- a/simpa/core/simulation_modules/volume_creation_module/__init__.py +++ b/simpa/core/simulation_modules/volume_creation_module/__init__.py @@ -20,7 +20,7 @@ class VolumeCreatorModuleBase(SimulationModule): def __init__(self, global_settings: Settings): super(VolumeCreatorModuleBase, self).__init__(global_settings=global_settings) - + def load_component_settings(self) -> Settings: """Implements abstract method to serve volume creation settings as component settings diff --git a/simpa/utils/deformation_manager.py b/simpa/utils/deformation_manager.py index 78d0a1ed..0c786a0d 100644 --- a/simpa/utils/deformation_manager.py +++ b/simpa/utils/deformation_manager.py @@ -26,7 +26,8 @@ def create_deformation_settings(bounds_mm, maximum_z_elevation_mm=1, filter_sigm # Add random permutations to the y-axis of the division knots all_scaling_value = np.multiply.outer( - np.cos(x_positions_vector / (bounds_mm[0][1] * (cosine_scaling_factor / np.pi)) - np.pi / (cosine_scaling_factor * 2)) ** 2, + np.cos(x_positions_vector / (bounds_mm[0][1] * (cosine_scaling_factor / + np.pi)) - np.pi / (cosine_scaling_factor * 2)) ** 2, np.cos(y_positions_vector / (bounds_mm[1][1] * (cosine_scaling_factor / np.pi)) - np.pi / (cosine_scaling_factor * 2)) ** 2) surface_elevations *= all_scaling_value @@ -60,7 +61,8 @@ def get_functional_from_deformation_settings(deformation_settings: dict): z_elevations_mm = deformation_settings[Tags.DEFORMATION_Z_ELEVATIONS_MM] order = "cubic" - functional_mm = RegularGridInterpolator(points=[x_coordinates_mm, y_coordinates_mm], values=z_elevations_mm, method=order) + functional_mm = RegularGridInterpolator( + points=[x_coordinates_mm, y_coordinates_mm], values=z_elevations_mm, method=order) return functional_mm diff --git a/simpa_tests/automatic_tests/device_tests/test_field_of_view.py b/simpa_tests/automatic_tests/device_tests/test_field_of_view.py index c1c79d87..cc63680c 100644 --- a/simpa_tests/automatic_tests/device_tests/test_field_of_view.py +++ b/simpa_tests/automatic_tests/device_tests/test_field_of_view.py @@ -101,4 +101,3 @@ def test_symmetric_with_odd_number_of_elements(self): self.assertAlmostEqual(ydim_end, 40) self.assertAlmostEqual(zdim_start, 0) self.assertAlmostEqual(zdim_end, 0) - diff --git a/simpa_tests/automatic_tests/simulation_tests/optical_simulation_module/test_utils.py b/simpa_tests/automatic_tests/simulation_tests/optical_simulation_module/test_utils.py new file mode 100644 index 00000000..f6397cad --- /dev/null +++ b/simpa_tests/automatic_tests/simulation_tests/optical_simulation_module/test_utils.py @@ -0,0 +1,106 @@ +# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ +# SPDX-FileCopyrightText: 2021 Janek Groehl +# SPDX-License-Identifier: MIT + +import numpy as np + +from simpa import Settings, Tags, TISSUE_LIBRARY, define_horizontal_layer_structure_settings, DiskIlluminationGeometry, \ + MCXAdapterReflectance, ModelBasedVolumeCreationAdapter, simulate, PathManager +from simpa.core.simulation_modules.optical_simulation_module.utils import \ + get_spectral_image_from_simulated_reflectances, get_spectral_image_from_optical_simulation + + +def test_if_get_spectral_image_from_simulated_reflectances_is_called_correct_spectral_image_is_returned(): + # Arrange + sample_wavelengths = np.arange(400, 600 + 1, 100) + reflectance_values = {"400": np.array([0.01, 0.34]), "500": np.array([0.04, 0.06, 0.52]), "600": np.array([0.05])} + reflectance_pos = {"400": np.array([[0, 0, 0], + [2, 1, 0]]), + "500": np.array([[0, 1, 0], + [1, 0, 0], + [1, 1, 0]]), + "600": np.array([[2, 0, 0]])} + target_height = 2 + target_width = 3 + expected_spectral_image = np.array([[[0.01, 0.0, 0.0], + [0.0, 0.06, 0.0], + [0.0, 0.0, 0.05]], + [[0.0, 0.04, 0.0], + [0.0, 0.52, 0.0], + [0.34, 0.0, 0.0]]]) + # Act + actual_spectral_image = get_spectral_image_from_simulated_reflectances(sample_wavelengths, + reflectance_values, + reflectance_pos, + target_height, + target_width) + # Assert + assert actual_spectral_image.shape == expected_spectral_image.shape + np.testing.assert_allclose(actual_spectral_image, expected_spectral_image, rtol=1e-8, atol=1e-8) + + +def test_spectral_image_is_returned_after_optical_simulation(): + dim_volume_x_y_z_mm = [4, 6, 7] + spacing = 0.5 + + target_height = round(dim_volume_x_y_z_mm[1] / spacing) + target_width = round(dim_volume_x_y_z_mm[0] / spacing) + + background_dictionary = Settings() + background_dictionary[Tags.MOLECULE_COMPOSITION] = TISSUE_LIBRARY.constant(1e-4, 1e-4, 0.9) + background_dictionary[Tags.STRUCTURE_TYPE] = Tags.BACKGROUND + + tissue_settings = Settings() + tissue_settings[Tags.BACKGROUND] = background_dictionary + + tissue_settings["tissue"] = define_horizontal_layer_structure_settings( + molecular_composition=TISSUE_LIBRARY.muscle(), + z_start_mm=0, + thickness_mm=dim_volume_x_y_z_mm[2]) + + path_manager = PathManager() + volume_name = "volume_name" + + general_settings = { + Tags.RANDOM_SEED: 0, + Tags.VOLUME_NAME: volume_name, + Tags.SIMULATION_PATH: path_manager.get_hdf5_file_save_path(), + Tags.SPACING_MM: spacing, + Tags.DIM_VOLUME_Z_MM: dim_volume_x_y_z_mm[2], + Tags.DIM_VOLUME_X_MM: dim_volume_x_y_z_mm[0], + Tags.DIM_VOLUME_Y_MM: dim_volume_x_y_z_mm[1], + Tags.WAVELENGTHS: np.arange(500, 800 + 1, 50), + Tags.DO_FILE_COMPRESSION: True + } + + expected_spectral_image_shape = (target_height, target_width, 7) + + settings = Settings(general_settings) + + settings.set_volume_creation_settings({ + Tags.SIMULATE_DEFORMED_LAYERS: True, + Tags.STRUCTURES: tissue_settings + }) + + settings.set_optical_settings({ + Tags.OPTICAL_MODEL_NUMBER_PHOTONS: 1e7, + Tags.OPTICAL_MODEL_BINARY_PATH: path_manager.get_mcx_binary_path(), + Tags.COMPUTE_DIFFUSE_REFLECTANCE: True, + Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT: True + }) + + pipeline = [ + ModelBasedVolumeCreationAdapter(settings), + MCXAdapterReflectance(settings), + ] + + device = DiskIlluminationGeometry(beam_radius_mm=2, + device_position_mm=np.array( + [dim_volume_x_y_z_mm[0] / 2., dim_volume_x_y_z_mm[1] / 2., 0])) + + simulate(pipeline, settings, device) + hdf_file_path = path_manager.get_hdf5_file_save_path() + "/" + volume_name + ".hdf5" + actual_spectral_image = get_spectral_image_from_optical_simulation(hdf_file_path, target_height, target_width) + assert isinstance(actual_spectral_image, np.ndarray) + assert actual_spectral_image.shape == expected_spectral_image_shape, actual_spectral_image.shape + assert actual_spectral_image.max() > 1e-4 diff --git a/simpa_tests/automatic_tests/test_bandpass_filter.py b/simpa_tests/automatic_tests/test_bandpass_filter.py index c499e2fd..404527d1 100644 --- a/simpa_tests/automatic_tests/test_bandpass_filter.py +++ b/simpa_tests/automatic_tests/test_bandpass_filter.py @@ -260,4 +260,3 @@ def test_butter_filter_with_random_signal(self, show_figure_on_screen=False): if show_figure_on_screen: self.visualize_filtered_spectrum(FILTERED_SIGNAL) - diff --git a/simpa_tests/do_coverage.py b/simpa_tests/do_coverage.py index cf9c4f69..88490637 100644 --- a/simpa_tests/do_coverage.py +++ b/simpa_tests/do_coverage.py @@ -25,4 +25,4 @@ cov.html_report(directory="../docs/test_coverage") # Exit with an appropriate code based on the test results -sys.exit(not result.wasSuccessful()) \ No newline at end of file +sys.exit(not result.wasSuccessful()) diff --git a/simpa_tests/manual_tests/acoustic_forward_models/KWaveAcousticForwardConvenienceFunction.py b/simpa_tests/manual_tests/acoustic_forward_models/KWaveAcousticForwardConvenienceFunction.py index 93b0a5e6..a9157c33 100644 --- a/simpa_tests/manual_tests/acoustic_forward_models/KWaveAcousticForwardConvenienceFunction.py +++ b/simpa_tests/manual_tests/acoustic_forward_models/KWaveAcousticForwardConvenienceFunction.py @@ -114,7 +114,7 @@ def test_convenience_function(self): get_detection_geometry(), speed_of_sound=1540, density=1000, alpha_coeff=0.0, spacing_mm=0.25) - + # reconstruct the time series data to compare it with initial pressure self.settings.set_reconstruction_settings({ Tags.RECONSTRUCTION_MODE: Tags.RECONSTRUCTION_MODE_PRESSURE,