diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 81292682..529ef16a 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -7,7 +7,7 @@ repos: - id: check-added-large-files # prevents adding large files, above 500kB - repo: https://github.com/pre-commit/mirrors-autopep8 # Uses MIT-License (MIT compatible) - rev: v2.0.2 # Use the sha / tag you want to point at + rev: v2.0.4 # Use the sha / tag you want to point at hooks: - id: autopep8 # formats code according to PEP8 standard. Lets commit fail if it needs to reformat code. Config for autopep8 is done in myproject.toml @@ -17,7 +17,7 @@ repos: - id: validate-cff # validates cff citation files - repo: https://github.com/Lucas-C/pre-commit-hooks # Uses MIT License (MIT compatible) - rev: v1.5.1 + rev: v1.5.4 hooks: - id: insert-license # Checks if the license header specified at license_header.txt is added in the first lines of each python file. If not, it suggests to insert them. types: [python] @@ -31,7 +31,7 @@ repos: - id: gitlint - repo: https://github.com/tcort/markdown-link-check # Uses ISC-License (MIT compatible) - rev: v3.10.3 + rev: v3.11.2 hooks: - id: markdown-link-check # checks if links in markdown files work exclude: docs/ diff --git a/pyproject.toml b/pyproject.toml index 979b6fd8..02a00faf 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -34,11 +34,12 @@ coverage = ">=6.1.2" # Uses Apache 2.0-License (MIT compatible) Deprecated = ">=1.2.13" # Uses MIT-License (MIT compatible) torch = ">=1.10.0" # Uses BSD-License (MIT compatible) python-dotenv = ">=0.19.2" # Uses BSD-License (MIT compatible) -pacfish = ">=0.4.4" # Uses BSD-License (MIT compatible) +pacfish = ">=0.4.4" # Uses BSD-License (MIT compatible) requests = ">=2.26.0" # Uses Apache 2.0-License (MIT compatible) wget = ">=3.2" # Is Public Domain (MIT compatible) jdata = ">=0.5.2" # Uses Apache 2.0-License (MIT compatible) -pre-commit = ">=3.2.2" # Uses MIT-License (MIT compatible) +pre-commit = ">=3.2.2" # Uses MIT-License (MIT compatible) +pmcx = ">=0.2.5" # Uses GNU Public License V3 or later [tool.poetry.group.docs.dependencies] sphinx-rtd-theme = "^1.0.0" diff --git a/simpa/core/simulation_modules/optical_simulation_module/__init__.py b/simpa/core/simulation_modules/optical_simulation_module/__init__.py index 2cd97ad8..0208fad0 100644 --- a/simpa/core/simulation_modules/optical_simulation_module/__init__.py +++ b/simpa/core/simulation_modules/optical_simulation_module/__init__.py @@ -3,6 +3,7 @@ # SPDX-License-Identifier: MIT from abc import abstractmethod from typing import Dict, Union +from collections.abc import Iterable import numpy as np @@ -125,27 +126,10 @@ def run_forward_model(self, :param anisotropy: Dimensionless scattering anisotropy :return: """ - if isinstance(_device, list): - # per convention this list has at least two elements - results = self.forward_model(absorption_cm=absorption, + _devices = _device if isinstance(_device, Iterable) else (_device,) + fluence = sum(self.forward_model(absorption_cm=absorption, scattering_cm=scattering, anisotropy=anisotropy, - illumination_geometry=_device[0]) - fluence = results[Tags.DATA_FIELD_FLUENCE] - for idx in range(1, len(_device)): - # we already looked at the 0th element, so go from 1 to n-1 - results = self.forward_model(absorption_cm=absorption, - scattering_cm=scattering, - anisotropy=anisotropy, - illumination_geometry=_device[idx]) - fluence += results[Tags.DATA_FIELD_FLUENCE] - - fluence = fluence / len(_device) - - else: - results = self.forward_model(absorption_cm=absorption, - scattering_cm=scattering, - anisotropy=anisotropy, - illumination_geometry=_device) - fluence = results[Tags.DATA_FIELD_FLUENCE] - return {Tags.DATA_FIELD_FLUENCE: fluence} + illumination_geometry=illumination_geometry)[Tags.DATA_FIELD_FLUENCE] + for illumination_geometry in _devices) + return {Tags.DATA_FIELD_FLUENCE: fluence / len(_devices)} 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 11aabb75..2c4c5669 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 @@ -3,20 +3,18 @@ # SPDX-License-Identifier: MIT import numpy as np -import subprocess +import pmcx from simpa.utils import Tags, Settings 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 os -import gc -from typing import List, Dict, Tuple +from typing import Dict, Tuple 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` + This class implements a bridge to the mcx framework using pmcx to integrate mcx into SIMPA. + This adapter only allows for computation of fluence, for computations of diffuse reflectance, + take a look at `simpa.ReflectanceMcxAdapter` .. note:: MCX is a GPU-enabled Monte-Carlo model simulation of photon transport in tissue: @@ -33,10 +31,7 @@ def __init__(self, global_settings: Settings): :param global_settings: global settings used during simulations """ super(MCXAdapter, self).__init__(global_settings=global_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'} def forward_model(self, absorption_cm: np.ndarray, @@ -44,8 +39,7 @@ def forward_model(self, anisotropy: np.ndarray, illumination_geometry: IlluminationGeometryBase) -> Dict: """ - runs the MCX simulations. Binary file containing scattering and absorption volumes is temporarily created as - input for MCX. A JSON serializable file containing the configuration required by MCx is also generated. + runs the MCX simulations. The set of flags parsed to MCX is built based on the Tags declared in `self.component_settings`, the results from MCX are used to populate an instance of Dict and returned. @@ -61,59 +55,54 @@ def forward_model(self, else: _assumed_anisotropy = 0.9 - self.generate_mcx_bin_input(absorption_cm=absorption_cm, - scattering_cm=scattering_cm, - anisotropy=anisotropy, - assumed_anisotropy=_assumed_anisotropy) + config = self.get_mcx_config(illumination_geometry=illumination_geometry, + absorption_cm=absorption_cm, + scattering_cm=scattering_cm, + anisotropy=anisotropy, + assumed_anisotropy=_assumed_anisotropy) - settings_dict = self.get_mcx_settings(illumination_geometry=illumination_geometry, - assumed_anisotropy=_assumed_anisotropy) + pmcx_results = pmcx.run(config) - print(settings_dict) - self.generate_mcx_json_input(settings_dict=settings_dict) - # run the simulation - cmd = self.get_command() - self.run_mcx(cmd) + return self.parse_pmcx_results(pmcx_results) - # Read output - results = self.read_mcx_output() - - # clean temporary files - self.remove_mcx_output() + def parse_pmcx_results(self, pmcx_results: Dict) -> Dict: + fluence = pmcx_results["flux"] + fluence *= 100 # Convert from J/mm^2 to J/cm^2 + if np.shape(fluence)[3] == 1: + fluence = np.squeeze(fluence, 3) + fluence = self.post_process_volumes(**{'arrays': (fluence,)})[0] + results = dict() + results[Tags.DATA_FIELD_FLUENCE] = fluence return results - def generate_mcx_json_input(self, settings_dict: Dict) -> None: - """ - generates JSON serializable file with settings needed by MCX to run simulations. - - :param settings_dict: dictionary to be saved as .json - :return: None - """ - tmp_json_filename = self.global_settings[Tags.SIMULATION_PATH] + "/" + \ - self.global_settings[Tags.VOLUME_NAME] + ".json" - self.mcx_json_config_file = tmp_json_filename - self.temporary_output_files.append(tmp_json_filename) - with open(tmp_json_filename, "w") as json_file: - json.dump(settings_dict, json_file, indent="\t") - - def get_mcx_settings(self, - illumination_geometry: IlluminationGeometryBase, - assumed_anisotropy: np.ndarray, - **kwargs) -> Dict: + def get_mcx_config(self, + illumination_geometry: IlluminationGeometryBase, + absorption_cm: np.ndarray, + scattering_cm: np.ndarray, + anisotropy: np.ndarray, + 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. + generates a pcmx config for simulations based on Tags in `self.global_settings` and + `self.component_settings`. Among others, it defines the volume array. :param illumination_geometry: and instance of `IlluminationGeometryBase` defining the illumination geometry + :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 kwargs: dummy, used for class inheritance :return: dictionary with settings to be used by MCX """ - mcx_volumetric_data_file = self.global_settings[Tags.SIMULATION_PATH] + "/" + \ - self.global_settings[Tags.VOLUME_NAME] + "_output" - for name, suffix in self.mcx_output_suffixes.items(): - self.__setattr__(name, mcx_volumetric_data_file + suffix) - self.temporary_output_files.append(mcx_volumetric_data_file + suffix) + 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 (2,nx,ny,nz) + vol = np.stack([absorption_mm, scattering_mm], dtype=np.float32) + [_, self.nx, self.ny, self.nz] = np.shape(vol) + source = illumination_geometry.get_mcx_illuminator_definition(self.global_settings) + prop = np.array([[0, 0, 1, 1], [1, 1, assumed_anisotropy, 1]], dtype=np.float32) if Tags.TIME_STEP and Tags.TOTAL_TIME in self.component_settings: dt = self.component_settings[Tags.TIME_STEP] time = self.component_settings[Tags.TOTAL_TIME] @@ -121,137 +110,30 @@ def get_mcx_settings(self, time = 5e-09 dt = 5e-09 self.frames = int(time / dt) - - source = illumination_geometry.get_mcx_illuminator_definition(self.global_settings) - settings_dict = { - "Session": { - "ID": mcx_volumetric_data_file, - "DoAutoThread": 1, - "Photons": self.component_settings[Tags.OPTICAL_MODEL_NUMBER_PHOTONS], - "DoMismatch": 0 - }, - "Forward": { - "T0": 0, - "T1": time, - "Dt": dt - }, - "Optode": { - "Source": source - }, - "Domain": { - "OriginType": 0, - "LengthUnit": self.global_settings[Tags.SPACING_MM], - "Media": [ - { - "mua": 0, - "mus": 0, - "g": 1, - "n": 1 - }, - { - "mua": 1, - "mus": 1, - "g": assumed_anisotropy, - "n": 1 - } - ], - "MediaFormat": "muamus_float", - "Dim": [self.nx, self.ny, self.nz], - "VolumeFile": self.global_settings[Tags.SIMULATION_PATH] + "/" + - self.global_settings[Tags.VOLUME_NAME] + ".bin" - }} + config = { + "nphoton": self.component_settings[Tags.OPTICAL_MODEL_NUMBER_PHOTONS], + "vol": vol, + "tstart": 0, + "tend": time, + "tstep": dt, + "prop": prop, + "unitinmm": self.global_settings[Tags.SPACING_MM], + "srctype": source["Type"], + "srcpos": source["Pos"], + "srcdir": source["Dir"], + "srcparam1": source["Param1"], + "srcparam2": source["Param2"], + "isreflect": 0, + "autopilot": 1, + "outputtype": "fluence", + } if Tags.MCX_SEED not in self.component_settings: if Tags.RANDOM_SEED in self.global_settings: - settings_dict["Session"]["RNGSeed"] = self.global_settings[Tags.RANDOM_SEED] + config["seed"] = self.global_settings[Tags.RANDOM_SEED] else: - settings_dict["Session"]["RNGSeed"] = self.component_settings[Tags.MCX_SEED] - return settings_dict - - def get_command(self) -> List: - """ - generates list of commands to be parse to MCX in a subprocess - - :return: list of MCX commands - """ - cmd = list() - cmd.append(self.component_settings[Tags.OPTICAL_MODEL_BINARY_PATH]) - cmd.append("-f") - 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") - return cmd - - @staticmethod - def run_mcx(cmd: List) -> None: - """ - runs subprocess calling MCX with the flags built with `self.get_command`. Rises a `RuntimeError` if the code - exit of the subprocess is not 0. + config["seed"] = self.component_settings[Tags.MCX_SEED] - :param cmd: list defining command to parse to `subprocess.run` - :return: None - """ - results = None - try: - results = subprocess.run(cmd) - except: - raise RuntimeError(f"MCX failed to run: {cmd}, results: {results}") - - def generate_mcx_bin_input(self, - absorption_cm: np.ndarray, - scattering_cm: np.ndarray, - anisotropy: np.ndarray, - assumed_anisotropy: 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: - :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) - [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] + "/" + \ - self.global_settings[Tags.VOLUME_NAME] + ".bin" - self.temporary_output_files.append(tmp_input_path) - # write array in 'C' order to binary file - op_array.tofile(tmp_input_path) - - def read_mcx_output(self, **kwargs) -> Dict: - """ - reads the temporary output generated with MCX - - :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) - results = dict() - results[Tags.DATA_FIELD_FLUENCE] = fluence - return results - - def remove_mcx_output(self) -> None: - """ - deletes temporary MCX output files from the file system - - :return: None - """ - for f in self.temporary_output_files: - if os.path.isfile(f): - os.remove(f) + return config def pre_process_volumes(self, **kwargs) -> Tuple: """ 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..f08d7f86 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 @@ -2,10 +2,8 @@ # SPDX-FileCopyrightText: 2021 Janek Groehl # SPDX-License-Identifier: MIT import numpy as np -import struct -import jdata -import os -from typing import List, Tuple, Dict, Union +from collections.abc import Iterable +from typing import Tuple, Dict, Union from simpa.utils import Tags, Settings from simpa.core.simulation_modules.optical_simulation_module.optical_forward_model_mcx_adapter import MCXAdapter @@ -14,8 +12,8 @@ class MCXAdapterReflectance(MCXAdapter): """ - This class implements a bridge to the mcx framework to integrate mcx into SIMPA. This class targets specifically - diffuse reflectance simulations. Specifically, it implements the capability to run diffuse reflectance simulations. + This class implements a bridge to the mcx framework using pmcx to integrate mcx into SIMPA. + This class 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 @@ -36,142 +34,49 @@ def __init__(self, global_settings: Settings): :param global_settings: global settings used during simulations """ super(MCXAdapterReflectance, self).__init__(global_settings=global_settings) - self.mcx_photon_data_file = None self.padded = None - self.mcx_output_suffixes = {'mcx_volumetric_data_file': '.jnii', - 'mcx_photon_data_file': '_detp.jdat'} - def forward_model(self, - absorption_cm: np.ndarray, - scattering_cm: np.ndarray, - anisotropy: np.ndarray, - illumination_geometry: IlluminationGeometryBase) -> Dict: - """ - runs the MCX simulations. Binary file containing scattering and absorption volumes is temporarily created as - input for MCX. A JSON serializable file containing the configuration required by MCx is also generated. - The set of flags parsed to MCX is built based on the Tags declared in `self.component_settings`, the results - from MCX are used to populate an instance of Settings and returned. - - :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 illumination_geometry: and instance of `IlluminationGeometryBase` defining the illumination geometry - :param probe_position_mm: position of a probe in `mm` units. This is parsed to - `illumination_geometry.get_mcx_illuminator_definition` - :return: `Settings` 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=_assumed_anisotropy, - assumed_anisotropy=_assumed_anisotropy) - - settings_dict = self.get_mcx_settings(illumination_geometry=illumination_geometry, - assumed_anisotropy=_assumed_anisotropy, - ) - - print(settings_dict) - self.generate_mcx_json_input(settings_dict=settings_dict) - # run the simulation - cmd = self.get_command() - self.run_mcx(cmd) - - # Read output - results = self.read_mcx_output() - struct._clearcache() - - # clean temporary files - self.remove_mcx_output() - return results - - def get_command(self) -> List: - """ - generates list of commands to be parse to MCX in a subprocess - - :return: list of MCX commands - """ - cmd = list() - cmd.append(self.component_settings[Tags.OPTICAL_MODEL_BINARY_PATH]) - cmd.append("-f") - cmd.append(self.mcx_json_config_file) - cmd.append("-O") - cmd.append("F") - cmd.append("-F") - cmd.append("jnii") - if Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT in self.component_settings and \ - self.component_settings[Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT]: - cmd.append("-H") - cmd.append(f"{int(self.component_settings[Tags.OPTICAL_MODEL_NUMBER_PHOTONS])}") - cmd.append("--bc") # save photon exit position and direction - cmd.append("______000010") - cmd.append("--savedetflag") - cmd.append("XV") - if Tags.COMPUTE_DIFFUSE_REFLECTANCE in self.component_settings and \ - self.component_settings[Tags.COMPUTE_DIFFUSE_REFLECTANCE]: - cmd.append("--saveref") # save diffuse reflectance at 0 filled voxels outside of domain - return cmd - - def read_mcx_output(self, **kwargs) -> Dict: - """ - reads the temporary output generated with MCX - - :param kwargs: dummy, used for class inheritance compatibility - :return: `Settings` instance containing the MCX output - """ - results = dict() - if os.path.isfile(self.mcx_volumetric_data_file) and self.mcx_volumetric_data_file.endswith( - self.mcx_output_suffixes['mcx_volumetric_data_file']): - content = jdata.load(self.mcx_volumetric_data_file) - fluence = content['NIFTIData'] - 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 - results[Tags.DATA_FIELD_FLUENCE] = fluence - else: - raise FileNotFoundError(f"Could not find .jnii file for {self.mcx_volumetric_data_file}") + def parse_pmcx_results(self, pmcx_results: Dict) -> Dict: + results = super(MCXAdapterReflectance, self).parse_pmcx_results(pmcx_results) if Tags.COMPUTE_DIFFUSE_REFLECTANCE in self.component_settings and \ self.component_settings[Tags.COMPUTE_DIFFUSE_REFLECTANCE]: + pos = np.where(pmcx_results["dref"] > 0) + ref = pmcx_results["dref"][pos] + pos = np.array(pos).T # reformatting to aggregate results after for multi illuminant geometries results[Tags.DATA_FIELD_DIFFUSE_REFLECTANCE] = ref - results[Tags.DATA_FIELD_DIFFUSE_REFLECTANCE_POS] = ref_pos + results[Tags.DATA_FIELD_DIFFUSE_REFLECTANCE_POS] = pos if Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT in self.component_settings and \ self.component_settings[Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT]: - content = jdata.load(self.mcx_photon_data_file) - photon_pos = content['MCXData']['PhotonData']['p'] - photon_dir = content['MCXData']['PhotonData']['v'] - results[Tags.DATA_FIELD_PHOTON_EXIT_POS] = photon_pos - results[Tags.DATA_FIELD_PHOTON_EXIT_DIR] = photon_dir + results[Tags.DATA_FIELD_PHOTON_EXIT_POS] = pmcx_results["detp"][0:3, :] # p + results[Tags.DATA_FIELD_PHOTON_EXIT_DIR] = pmcx_results["detp"][3:6, :] # v return results - @staticmethod - def extract_reflectance_from_fluence(fluence: np.ndarray) -> Tuple: - """ - extracts diffuse reflectance from volumes. MCX stores diffuse reflectance as negative values in the fluence - volume. The position where the reflectance is stored is also returned. If there are no negative values in the - fluence, `None` is returned instead of reflectance and reflectance position. Negative values in fluence are - set to `0` after extraction of the reflectance. - - :param fluence: array containing fluence as generated by MCX - :return: tuple of reflectance, reflectance position and transformed fluence - """ - if np.any(fluence < 0): - pos = np.where(fluence < 0) - ref = fluence[pos] * -1 - pos = np.array(pos).T # reformatting to aggregate results after for multi illuminant geometries - fluence[fluence < 0] = 0 - return ref, pos, fluence - else: - return None, None, fluence + def get_mcx_config(self, + illumination_geometry: IlluminationGeometryBase, + absorption_cm: np.ndarray, + scattering_cm: np.ndarray, + anisotropy: np.ndarray, + assumed_anisotropy: np.ndarray, + **kwargs) -> Dict: + config = super(MCXAdapterReflectance, self).get_mcx_config(illumination_geometry, + absorption_cm, scattering_cm, anisotropy, assumed_anisotropy, **kwargs) + if Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT in self.component_settings and \ + self.component_settings[Tags.COMPUTE_PHOTON_DIRECTION_AT_EXIT]: + config["issavedet"] = 1 + config["issrcfrom0"] = 1 + config["maxdetphoton"] = int(self.component_settings[Tags.OPTICAL_MODEL_NUMBER_PHOTONS]) + config["bc"] = "______000010" # detect photons exiting in the +y direction + config["savedetflag"] = "XV" # save photon exit position and direction + if Tags.COMPUTE_DIFFUSE_REFLECTANCE in self.component_settings and \ + self.component_settings[Tags.COMPUTE_DIFFUSE_REFLECTANCE]: + config["issaveref"] = 1 # save diffuse reflectance at 0 filled voxels outside of domain + return config def pre_process_volumes(self, **kwargs) -> Tuple: """ pre-process volumes before running simulations with MCX. The volumes are transformed to `mm` units and pads a 0-values layer along the z-axis in order to save diffuse reflectance values. All 0-valued voxels are then - transformed to `np.nan`. This last transformation `0->np.nan` is a requirement form MCX. + transformed to `np.nan`. This last transformation `0->np.nan` is a requirement from MCX. :param kwargs: dictionary containing at least the keys `scattering_cm, absorption_cm, anisotropy` and `assumed_anisotropy` @@ -197,8 +102,8 @@ def pre_process_volumes(self, **kwargs) -> Tuple: def post_process_volumes(self, **kwargs): """ - post-processes volumes after MCX simulations. Dummy function implemented for compatibility with inherited - classes. Removes layer padded by `self.pre_process_volumes` if it was added and transforms `np.nan -> 0`. + post-processes volumes after MCX simulations. + Removes layer padded by `self.pre_process_volumes` if it was added and transforms `np.nan -> 0`. :param kwargs: dictionary containing at least the key `volumes` to be transformed :return: @@ -210,7 +115,7 @@ def post_process_volumes(self, **kwargs): results = tuple(arrays) for a in results: # revert nan transformation that was done while pre-processing volumes - a[np.isnan(a)] = 0. + a[np.isnan(a)] = 0 return results def run_forward_model(self, @@ -234,46 +139,26 @@ def run_forward_model(self, reflectance_position = [] photon_position = [] photon_direction = [] - if isinstance(_device, list): - # per convention this list has at least two elements + _devices = _device if isinstance(_device, Iterable) else (_device,) + fluence = None + for illumination_geometry in _devices: results = self.forward_model(absorption_cm=absorption, scattering_cm=scattering, anisotropy=anisotropy, - illumination_geometry=_device[0]) - self._append_results(results=results, - reflectance=reflectance, - reflectance_position=reflectance_position, - photon_position=photon_position, - photon_direction=photon_direction) - fluence = results[Tags.DATA_FIELD_FLUENCE] - for idx in range(1, len(_device)): - # we already looked at the 0th element, so go from 1 to n-1 - results = self.forward_model(absorption_cm=absorption, - scattering_cm=scattering, - anisotropy=anisotropy, - illumination_geometry=_device[idx]) - self._append_results(results=results, - reflectance=reflectance, - reflectance_position=reflectance_position, - photon_position=photon_position, - photon_direction=photon_direction) + illumination_geometry=illumination_geometry) + if fluence: fluence += results[Tags.DATA_FIELD_FLUENCE] + else: + fluence = results[Tags.DATA_FIELD_FLUENCE] + if Tags.DATA_FIELD_DIFFUSE_REFLECTANCE in results: + reflectance.append(results[Tags.DATA_FIELD_DIFFUSE_REFLECTANCE]) + reflectance_position.append(results[Tags.DATA_FIELD_DIFFUSE_REFLECTANCE_POS]) + if Tags.DATA_FIELD_PHOTON_EXIT_POS in results: + photon_position.append(results[Tags.DATA_FIELD_PHOTON_EXIT_POS]) + photon_direction.append(results[Tags.DATA_FIELD_PHOTON_EXIT_DIR]) - fluence = fluence / len(_device) - - else: - results = self.forward_model(absorption_cm=absorption, - scattering_cm=scattering, - anisotropy=anisotropy, - illumination_geometry=_device) - self._append_results(results=results, - reflectance=reflectance, - reflectance_position=reflectance_position, - photon_position=photon_position, - photon_direction=photon_direction) - fluence = results[Tags.DATA_FIELD_FLUENCE] aggregated_results = dict() - aggregated_results[Tags.DATA_FIELD_FLUENCE] = fluence + aggregated_results[Tags.DATA_FIELD_FLUENCE] = fluence / len(_devices) if reflectance: aggregated_results[Tags.DATA_FIELD_DIFFUSE_REFLECTANCE] = np.concatenate(reflectance, axis=0) aggregated_results[Tags.DATA_FIELD_DIFFUSE_REFLECTANCE_POS] = np.concatenate(reflectance_position, axis=0) @@ -281,16 +166,3 @@ def run_forward_model(self, aggregated_results[Tags.DATA_FIELD_PHOTON_EXIT_POS] = np.concatenate(photon_position, axis=0) aggregated_results[Tags.DATA_FIELD_PHOTON_EXIT_DIR] = np.concatenate(photon_direction, axis=0) return aggregated_results - - @staticmethod - def _append_results(results, - reflectance, - reflectance_position, - photon_position, - photon_direction): - if Tags.DATA_FIELD_DIFFUSE_REFLECTANCE in results: - reflectance.append(results[Tags.DATA_FIELD_DIFFUSE_REFLECTANCE]) - reflectance_position.append(results[Tags.DATA_FIELD_DIFFUSE_REFLECTANCE_POS]) - if Tags.DATA_FIELD_PHOTON_EXIT_POS in results: - photon_position.append(results[Tags.DATA_FIELD_PHOTON_EXIT_POS]) - photon_direction.append(results[Tags.DATA_FIELD_PHOTON_EXIT_DIR])