diff --git a/README.md b/README.md index 99de8b48..f61481dc 100755 --- a/README.md +++ b/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) diff --git a/docs/source/introduction.md b/docs/source/introduction.md index f211279d..593b74d1 100644 --- a/docs/source/introduction.md +++ b/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) diff --git a/simpa/core/device_digital_twins/digital_device_twin_base.py b/simpa/core/device_digital_twins/digital_device_twin_base.py index c9280ca2..56c90128 100644 --- a/simpa/core/device_digital_twins/digital_device_twin_base.py +++ b/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 diff --git a/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_acuity_illumination.py b/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_acuity_illumination.py index 57729c9a..f1b665d1 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_acuity_illumination.py +++ b/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 diff --git a/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_invision_illumination.py b/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_invision_illumination.py index 4dbdfbd0..5713fb54 100644 --- a/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_invision_illumination.py +++ b/simpa/core/device_digital_twins/illumination_geometries/ithera_msot_invision_illumination.py @@ -13,11 +13,7 @@ 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: @@ -25,75 +21,53 @@ def __init__(self, invision_position=None, geometry_id=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: diff --git a/simpa/core/device_digital_twins/pa_devices/ithera_msot_invision.py b/simpa/core/device_digital_twins/pa_devices/ithera_msot_invision.py index 59abcdd3..f2e42945 100644 --- a/simpa/core/device_digital_twins/pa_devices/ithera_msot_invision.py +++ b/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__ 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 0e922c95..ed4617a4 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 @@ -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 diff --git a/simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_mcx_reflectance_adapter.py b/simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_mcx_reflectance_adapter.py index 4c802771..fa9a64a7 100644 --- a/simpa/core/simulation_modules/optical_simulation_module/optical_forward_model_mcx_reflectance_adapter.py +++ b/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 diff --git a/simpa/utils/quality_assurance/data_sanity_testing.py b/simpa/utils/quality_assurance/data_sanity_testing.py index 32cee0be..dd20c2f2 100644 --- a/simpa/utils/quality_assurance/data_sanity_testing.py +++ b/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):