Use mcx SLIT source for MSOTAcuity and MSOTInvision sources

README.md

@@ -75,25 +75,13 @@ acoustic simulations possible.
 
 ### mcx (Optical Forward Model)
 
-Either download suitable executables or build yourself from the following sources:
+Download the latest nightly build of [mcx](http://mcx.space/) for your operating system:
 
-[http://mcx.space/](http://mcx.space/)
+- [Linux](http://mcx.space/nightly/github/mcx-linux-x64-github-latest.zip)
+- [MacOS](http://mcx.space/nightly/github/mcx-macos-x64-github-latest.zip)
+- [Windows](http://mcx.space/nightly/github/mcx-windows-x64-github-latest.zip)
 
-In order to obtain access to all custom sources that we implemented, please build mcx yourself from the
-following mcx Github fork:
-[https://github.com/IMSY-DKFZ/mcx](https://github.com/IMSY-DKFZ/mcx)
-
-For the installation, please follow steps 1-4:
-1. `git clone https://github.com/IMSY-DKFZ/mcx.git`
-2. `cd mcx/src`
-3. In `MAKEFILE` adapt line 111 the sm version [according to your GPU](https://arnon.dk/matching-sm-architectures-arch-and-gencode-for-various-nvidia-cards/).
-4. `make`
-
-The built binary can be found in `src/bin`.
-Note, in case you can’t build mcx with the GPU-specific sm version you need to install a more recent NVIDIA driver and nvcc toolkit. 
-One option would be to install cuda in a conda environment via `conda install cuda -c nvidia`.
-Please note that there might be compatibility issues using mcx-cl with the MCX Adapter as this use case is not 
-being tested and supported by the SIMPA developers.
+Then extract the files and set `MCX_BINARY_PATH=/.../mcx/bin/mcx` in your path_config.env.
 
 ### k-Wave (Acoustic Forward Model)
 

docs/source/introduction.md

@@ -39,25 +39,13 @@ acoustic simulations possible.
 
 ### mcx (Optical Forward Model)
 
-Either download suitable executables or build yourself from the following sources:
+Download the latest nightly build of [mcx](http://mcx.space/) for your operating system:
 
-[http://mcx.space/](http://mcx.space/)
+- [Linux](http://mcx.space/nightly/github/mcx-linux-x64-github-latest.zip)
+- [MacOS](http://mcx.space/nightly/github/mcx-macos-x64-github-latest.zip)
+- [Windows](http://mcx.space/nightly/github/mcx-windows-x64-github-latest.zip)
 
-In order to obtain access to all custom sources that we implemented, please build mcx yourself from the
-following mcx Github fork:
-[https://github.com/IMSY-DKFZ/mcx](https://github.com/IMSY-DKFZ/mcx)
-
-For the installation, please follow steps 1-4:
-1. `git clone https://github.com/IMSY-DKFZ/mcx.git`
-2. `cd mcx/src`
-3. In `MAKEFILE` adapt line 111 the sm version [according to your GPU](https://arnon.dk/matching-sm-architectures-arch-and-gencode-for-various-nvidia-cards/).
-4. `make`
-
-The built binary can be found in `src/bin`.
-Note, in case you can’t build mcx with the GPU-specific sm version you need to install a more recent NVIDIA driver and nvcc toolkit. 
-One option would be to install cuda in a conda environment via `conda install cuda -c nvidia`.
-Please note that there might be compatibility issues using mcx-cl with the MCX Adapter as this use case is not 
-being tested and supported by the SIMPA developers.
+Then extract the files and set `MCX_BINARY_PATH=/.../mcx/bin/mcx` in your path_config.env.
 
 ### k-Wave (Acoustic Forward Model)
 

simpa/core/device_digital_twins/digital_device_twin_base.py

@@ -12,6 +12,15 @@
 from simpa.utils.serializer import SerializableSIMPAClass
 
 
+def is_equal(a, b) -> bool:
+    if isinstance(a, list):
+        return all(is_equal(ai, bi) for ai, bi in zip(a, b))
+    elif isinstance(a, np.ndarray):
+        return (a == b).all
+    else:
+        return a == b
+
+
 class DigitalDeviceTwinBase(SerializableSIMPAClass):
     """
     This class represents a device that can be used for illumination, detection or a combined photoacoustic device
@@ -49,14 +58,8 @@ def __eq__(self, other):
                 return False
             for self_key, self_value in self.__dict__.items():
                 other_value = other.__dict__[self_key]
-                if isinstance(self_value, np.ndarray):
-                    boolean = (other_value != self_value).all()
-                else:
-                    boolean = other_value != self_value
-                if boolean:
+                if not is_equal(self_value, other_value):
                     return False
-                else:
-                    continue
             return True
         return False
 

simpa/core/device_digital_twins/illumination_geometries/ithera_msot_acuity_illumination.py

@@ -29,21 +29,25 @@ def __init__(self):
         self.normalized_source_direction_vector = self.source_direction_vector / np.linalg.norm(
             self.source_direction_vector)
 
-    def get_mcx_illuminator_definition(self, global_settings: Settings):
+        divergence_angle = 8.66  # full beam divergence angle measured at Full Width at Half Maximum (FWHM)
+        full_width_at_half_maximum = 2.0 * np.tan(0.5 * np.deg2rad(divergence_angle))  # FWHM of beam divergence
+        # standard deviation of gaussian with FWHM
+        self.sigma = full_width_at_half_maximum / (2.0 * np.sqrt(2.0 * np.log(2.0)))
 
-        source_type = Tags.ILLUMINATION_TYPE_MSOT_ACUITY_ECHO
+    def get_mcx_illuminator_definition(self, global_settings: Settings):
+        source_type = Tags.ILLUMINATION_TYPE_SLIT
         spacing = global_settings[Tags.SPACING_MM]
         device_position = list(self.device_position_mm / spacing + 0.5)
-
+        device_length = 30 / spacing
+        source_pos = device_position
+        source_pos[0] -= 0.5 * device_length
         source_direction = list(self.normalized_source_direction_vector)
-
-        source_param1 = [30 / spacing, 0, 0, 0]
-
-        source_param2 = [0, 0, 0, 0]
+        source_param1 = [device_length, 0.0, 0.0, 0.0]
+        source_param2 = [self.sigma, self.sigma, 0.0, 0.0]
 
         return {
             "Type": source_type,
-            "Pos": device_position,
+            "Pos": source_pos,
             "Dir": source_direction,
             "Param1": source_param1,
             "Param2": source_param2

simpa/core/device_digital_twins/illumination_geometries/ithera_msot_invision_illumination.py

@@ -13,87 +13,61 @@ class MSOTInVisionIlluminationGeometry(IlluminationGeometryBase):
     This class represents the illumination geometry of the MSOT InVision photoacoustic device.
     """
 
-    def __init__(self, invision_position=None, geometry_id=0):
-        """
-        :param geometry_id: ID of the specific InVision illuminator.
-        :type geometry_id: int
-        """
+    def __init__(self, invision_position=None):
         super().__init__()
 
         if invision_position is None:
             self.invision_position = [0, 0, 0]
         else:
             self.invision_position = invision_position
 
-        self.geometry_id = geometry_id
-
-        angle = 0.0
         det_sep_half = 24.74 / 2
         detector_iso_distance = 74.05 / 2
-        illumination_angle = -0.41608649
-
-        if geometry_id == 0:
-            angle = 0.0
-        elif geometry_id == 1:
-            angle = 0.0
-            det_sep_half = -det_sep_half
-            illumination_angle = -illumination_angle
-        elif geometry_id == 2:
-            angle = 1.25664
-        elif geometry_id == 3:
-            angle = 1.25664
-            det_sep_half = -det_sep_half
-            illumination_angle = -illumination_angle
-        elif geometry_id == 4:
-            angle = -1.25664
-        elif geometry_id == 5:
-            angle = -1.25664
-            det_sep_half = -det_sep_half
-            illumination_angle = -illumination_angle
-        elif geometry_id == 6:
-            angle = 2.51327
-        elif geometry_id == 7:
-            angle = 2.51327
-            det_sep_half = -det_sep_half
-            illumination_angle = -illumination_angle
-        elif geometry_id == 8:
-            angle = -2.51327
-        elif geometry_id == 9:
-            angle = -2.51327
-            det_sep_half = -det_sep_half
-            illumination_angle = -illumination_angle
-
-        self.device_position_mm = [self.invision_position[0] + np.sin(angle) * detector_iso_distance,
-                                   self.invision_position[1] + det_sep_half,
-                                   self.invision_position[2] + np.cos(angle) * detector_iso_distance]
-
-        self.source_direction_vector = np.array([-np.sin(angle),
-                                                 np.sin(illumination_angle),
-                                                 np.cos(angle)])
-
-        self.normalized_source_direction_vector = self.source_direction_vector / np.linalg.norm(
-            self.source_direction_vector)
+        detector_width = 2 * 6.12
+
+        self.device_positions_mm = list()
+        self.source_direction_vectors = list()
+        self.slit_vectors_mm = list()
+
+        for index in [0, 1, 2, 3, 4]:
+            for y_offset_factor in [+1, -1]:
+                angle = index * 2.0 * np.pi / 5.0
+                illumination_angle = -0.41608649 * y_offset_factor
+                v = np.array([-np.sin(angle), np.sin(illumination_angle), -np.cos(angle)])
+                v /= np.linalg.norm(v)
+                slit_vector = np.array([np.cos(angle), 0, -np.sin(angle)]) * detector_width
+                slit_middle_on_circle = np.array([np.sin(angle), 0.0, np.cos(angle)]) * detector_iso_distance
+                y_offset = np.array([0.0, det_sep_half * y_offset_factor, 0.0])
+                pos = np.array(self.invision_position) + slit_middle_on_circle + y_offset - 0.5*slit_vector
+                self.device_positions_mm.append(pos)
+                self.source_direction_vectors.append(v)
+                self.slit_vectors_mm.append(slit_vector)
+
+        divergence_angle = 0.165806  # full beam divergence angle measured at Full Width at Half Maximum (FWHM)
+        full_width_at_half_maximum = 2.0 * np.tan(0.5 * divergence_angle)  # FWHM of beam divergence
+        # standard deviation of gaussian with FWHM
+        self.sigma = full_width_at_half_maximum / (2.0 * np.sqrt(2.0 * np.log(2.0)))
 
     def get_mcx_illuminator_definition(self, global_settings):
         self.logger.debug(self.invision_position)
-        source_type = Tags.ILLUMINATION_TYPE_MSOT_INVISION
+        source_type = Tags.ILLUMINATION_TYPE_SLIT
 
         spacing = global_settings[Tags.SPACING_MM]
 
-        source_position = list(np.array(self.device_position_mm) / spacing + 1)
+        source_positions = list(list(pos / spacing + 1) for pos in self.device_positions_mm)
 
-        source_direction = list(self.normalized_source_direction_vector)
+        source_directions = list(list(v) for v in self.source_direction_vectors)
 
-        source_param1 = [spacing, self.geometry_id, 0, 0]
+        source_param1s = list(list(slit_vector / spacing) + [0.0] for slit_vector in self.slit_vectors_mm)
 
-        source_param2 = [0, 0, 0, 0]
+        source_param2s = [[self.sigma, 0, 0, 0]]*len(self.device_positions_mm)
 
         return {
             "Type": source_type,
-            "Pos": source_position,
-            "Dir": source_direction,
-            "Param1": source_param1,
-            "Param2": source_param2
+            "Pos": source_positions,
+            "Dir": source_directions,
+            "Param1": source_param1s,
+            "Param2": source_param2s
         }
 
     def serialize(self) -> dict:

simpa/core/device_digital_twins/pa_devices/ithera_msot_invision.py

@@ -54,9 +54,7 @@ def __init__(self, device_position_mm: np.ndarray = None,
 
         self.field_of_view_extent_mm = detection_geometry.field_of_view_extent_mm
         self.set_detection_geometry(detection_geometry)
-        for i in range(10):
-            self.add_illumination_geometry(MSOTInVisionIlluminationGeometry(invision_position=self.device_position_mm,
-                                                                            geometry_id=i))
+        self.add_illumination_geometry(MSOTInVisionIlluminationGeometry(invision_position=self.device_position_mm))
 
     def serialize(self) -> dict:
         serialized_device = self.__dict__

simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_mcx_adapter.py

@@ -8,8 +8,8 @@
 from simpa.core.simulation_modules.optical_simulation_module import OpticalForwardModuleBase
 from simpa.core.device_digital_twins.illumination_geometries.illumination_geometry_base import IlluminationGeometryBase
 import json
+import jdata
 import os
-import gc
 from typing import List, Dict, Tuple
 
 
@@ -36,7 +36,7 @@ def __init__(self, global_settings: Settings):
         self.mcx_json_config_file = None
         self.mcx_volumetric_data_file = None
         self.frames = None
-        self.mcx_output_suffixes = {'mcx_volumetric_data_file': '.mc2'}
+        self.mcx_output_suffixes = {'mcx_volumetric_data_file': '.jnii'}
 
     def forward_model(self,
                       absorption_cm: np.ndarray,
@@ -182,9 +182,8 @@ def get_command(self) -> List:
         # use 'C' order array format for binary input file
         cmd.append("-a")
         cmd.append("1")
-        # use mc2 binary output file format
         cmd.append("-F")
-        cmd.append("mc2")
+        cmd.append("jnii")
         return cmd
 
     @staticmethod
@@ -237,11 +236,13 @@ def read_mcx_output(self, **kwargs) -> Dict:
         :param kwargs: dummy, used for class inheritance compatibility
         :return: `Dict` instance containing the MCX output
         """
-        shape = [self.nx, self.ny, self.nz, self.frames]
-        fluence = np.fromfile(self.mcx_volumetric_data_file, dtype=np.float32).reshape(shape, order='F')
-        fluence *= 100  # Convert from J/mm^2 to J/cm^2
-        if np.shape(fluence)[3] == 1:
-            fluence = np.squeeze(fluence, 3)
+        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

simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_mcx_reflectance_adapter.py

@@ -18,8 +18,8 @@ class MCXAdapterReflectance(MCXAdapter):
     diffuse reflectance simulations. Specifically, it implements the capability to run diffuse reflectance simulations.
 
     .. warning::
-        This MCX adapter requires a version of MCX containing the revision: `Rev::077060`, which was published in the
-        Nightly build  on `2022-01-26`.
+        This MCX adapter requires a version of MCX containing the commit 56eca8ae7e9abde309053759d6d6273ac4795fc5,
+        which was published in the Nightly build on `2024-03-10`.
 
     .. note::
         MCX is a GPU-enabled Monte-Carlo model simulation of photon transport in tissue:
@@ -101,6 +101,9 @@ def get_command(self) -> List:
         cmd.append(self.mcx_json_config_file)
         cmd.append("-O")
         cmd.append("F")
+        # 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 \
@@ -128,6 +131,9 @@ def read_mcx_output(self, **kwargs) -> Dict:
                 self.mcx_output_suffixes['mcx_volumetric_data_file']):
             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, fluence = self.extract_reflectance_from_fluence(fluence=fluence)
             fluence = self.post_process_volumes(**{'arrays': (fluence,)})[0]
             fluence *= 100  # Convert from J/mm^2 to J/cm^2

simpa/utils/quality_assurance/data_sanity_testing.py

@@ -45,7 +45,11 @@ def assert_array_well_defined(array: Union[np.ndarray, torch.Tensor], assume_non
     """
 
     error_message = None
-    array_as_tensor = torch.as_tensor(array)
+    if isinstance(array, np.ndarray) and any(stride < 0 for stride in array.strides):
+        # torch does not support tensors with negative strides so we need to make a copy of the array
+        array_as_tensor = torch.as_tensor(array.copy())
+    else:
+        array_as_tensor = torch.as_tensor(array)
     if not torch.all(torch.isfinite(array_as_tensor)):
         error_message = "nan, inf or -inf"
     if assume_positivity and torch.any(array_as_tensor <= 0):