Gn mcx based on t265

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:

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:

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:

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:

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
 

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

pre_commit_configs/link-config.json

@@ -0,0 +1,7 @@
+{
+    "ignorePatterns": [
+        {
+            "pattern": "https://doi.org/10.1117/1.JBO.27.8.083010"
+        }
+    ]
+}

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 \

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
 

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.

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: