T302 reflectance as image

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:

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
-

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]
-

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. 

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
 

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:

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

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
 

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
 

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
 
 

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)
-

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

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)
-

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())
+sys.exit(not result.wasSuccessful())

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,