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/simpa/core/processing_components/__init__.py b/simpa/core/processing_components/__init__.py index ff831e5e..8c27b88c 100644 --- a/simpa/core/processing_components/__init__.py +++ b/simpa/core/processing_components/__init__.py @@ -4,6 +4,7 @@ from abc import ABC from simpa.core import SimulationModule +from simpa.utils.processing_device import get_processing_device class ProcessingComponent(SimulationModule, ABC): @@ -19,3 +20,4 @@ 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] + self.torch_device = get_processing_device(global_settings) diff --git a/simpa/core/processing_components/monospectral/noise/gamma_noise.py b/simpa/core/processing_components/monospectral/noise/gamma_noise.py index f8608cd8..011ddffa 100644 --- a/simpa/core/processing_components/monospectral/noise/gamma_noise.py +++ b/simpa/core/processing_components/monospectral/noise/gamma_noise.py @@ -3,6 +3,7 @@ # SPDX-License-Identifier: MIT import numpy as np +import torch from simpa.core.processing_components import ProcessingComponent from simpa.io_handling import load_data_field, save_data_field @@ -12,7 +13,7 @@ class GammaNoise(ProcessingComponent): """ - Applies Gaussian noise to the defined data field. + Applies Gamma noise to the defined data field. The noise will be applied to all wavelengths. Component Settings:: @@ -51,15 +52,19 @@ def run(self, device): wavelength = self.global_settings[Tags.WAVELENGTH] data_array = load_data_field(self.global_settings[Tags.SIMPA_OUTPUT_PATH], data_field, wavelength) + data_tensor = torch.as_tensor(data_array, dtype=torch.float32, device=self.torch_device) + dist = torch.distributions.gamma.Gamma(torch.tensor(shape, dtype=torch.float32, device=self.torch_device), + torch.tensor(1.0/scale, dtype=torch.float32, device=self.torch_device)) if mode == Tags.NOISE_MODE_ADDITIVE: - data_array = data_array + np.random.gamma(shape, scale, size=np.shape(data_array)) + data_tensor += dist.sample(data_tensor.shape) elif mode == Tags.NOISE_MODE_MULTIPLICATIVE: - data_array = data_array * np.random.gamma(shape, scale, size=np.shape(data_array)) + data_tensor *= dist.sample(data_tensor.shape) if not (Tags.IGNORE_QA_ASSERTIONS in self.global_settings and Tags.IGNORE_QA_ASSERTIONS): - assert_array_well_defined(data_array) + assert_array_well_defined(data_tensor) - save_data_field(data_array, self.global_settings[Tags.SIMPA_OUTPUT_PATH], data_field, wavelength) + save_data_field(data_tensor.cpu().numpy().astype(np.float64, copy=False), + self.global_settings[Tags.SIMPA_OUTPUT_PATH], data_field, wavelength) self.logger.info("Applying Gamma Noise Model...[Done]") diff --git a/simpa/core/processing_components/monospectral/noise/gaussian_noise.py b/simpa/core/processing_components/monospectral/noise/gaussian_noise.py index 92d13d0e..5c31a50d 100644 --- a/simpa/core/processing_components/monospectral/noise/gaussian_noise.py +++ b/simpa/core/processing_components/monospectral/noise/gaussian_noise.py @@ -8,6 +8,7 @@ from simpa.core.processing_components import ProcessingComponent from simpa.utils.quality_assurance.data_sanity_testing import assert_array_well_defined import numpy as np +import torch class GaussianNoise(ProcessingComponent): @@ -56,17 +57,21 @@ def run(self, device): wavelength = self.global_settings[Tags.WAVELENGTH] data_array = load_data_field(self.global_settings[Tags.SIMPA_OUTPUT_PATH], data_field, wavelength) + data_tensor = torch.as_tensor(data_array, dtype=torch.float32, device=self.torch_device) + dist = torch.distributions.normal.Normal(torch.tensor(mean, dtype=torch.float32, device=self.torch_device), + torch.tensor(std, dtype=torch.float32, device=self.torch_device)) if mode == Tags.NOISE_MODE_ADDITIVE: - data_array = data_array + np.random.normal(mean, std, size=np.shape(data_array)) + data_tensor += dist.sample(data_tensor.shape) elif mode == Tags.NOISE_MODE_MULTIPLICATIVE: - data_array = data_array * np.random.normal(mean, std, size=np.shape(data_array)) + data_tensor *= dist.sample(data_tensor.shape) if not (Tags.IGNORE_QA_ASSERTIONS in self.global_settings and Tags.IGNORE_QA_ASSERTIONS): - assert_array_well_defined(data_array) + assert_array_well_defined(data_tensor) if non_negative: - data_array[data_array < EPS] = EPS - save_data_field(data_array, self.global_settings[Tags.SIMPA_OUTPUT_PATH], data_field, wavelength) + data_tensor[data_tensor < EPS] = EPS + save_data_field(data_tensor.cpu().numpy().astype(np.float64, copy=False), + self.global_settings[Tags.SIMPA_OUTPUT_PATH], data_field, wavelength) self.logger.info("Applying Gaussian Noise Model...[Done]") diff --git a/simpa/core/processing_components/monospectral/noise/poisson_noise.py b/simpa/core/processing_components/monospectral/noise/poisson_noise.py index df440c32..0fd544e6 100644 --- a/simpa/core/processing_components/monospectral/noise/poisson_noise.py +++ b/simpa/core/processing_components/monospectral/noise/poisson_noise.py @@ -7,11 +7,12 @@ from simpa.core.processing_components import ProcessingComponent from simpa.utils.quality_assurance.data_sanity_testing import assert_array_well_defined import numpy as np +import torch class PoissonNoise(ProcessingComponent): """ - Applies Gaussian noise to the defined data field. + Applies Poisson noise to the defined data field. The noise will be applied to all wavelengths. Component Settings:: @@ -44,15 +45,18 @@ def run(self, device): wavelength = self.global_settings[Tags.WAVELENGTH] data_array = load_data_field(self.global_settings[Tags.SIMPA_OUTPUT_PATH], data_field, wavelength) + data_tensor = torch.as_tensor(data_array, dtype=torch.float32, device=self.torch_device) + dist = torch.distributions.poisson.Poisson(torch.tensor(mean, dtype=torch.float32, device=self.torch_device)) if mode == Tags.NOISE_MODE_ADDITIVE: - data_array = data_array + np.random.poisson(mean, size=np.shape(data_array)) + data_tensor += dist.sample(data_tensor.shape) elif mode == Tags.NOISE_MODE_MULTIPLICATIVE: - data_array = data_array * np.random.poisson(mean, size=np.shape(data_array)) + data_tensor *= dist.sample(data_tensor.shape) if not (Tags.IGNORE_QA_ASSERTIONS in self.global_settings and Tags.IGNORE_QA_ASSERTIONS): - assert_array_well_defined(data_array) + assert_array_well_defined(data_tensor) - save_data_field(data_array, self.global_settings[Tags.SIMPA_OUTPUT_PATH], data_field, wavelength) + save_data_field(data_tensor.cpu().numpy().astype(np.float64, copy=False), + self.global_settings[Tags.SIMPA_OUTPUT_PATH], data_field, wavelength) self.logger.info("Applying Poisson Noise Model...[Done]") diff --git a/simpa/core/processing_components/monospectral/noise/salt_and_pepper_noise.py b/simpa/core/processing_components/monospectral/noise/salt_and_pepper_noise.py index 40dc5ce6..4204b32e 100644 --- a/simpa/core/processing_components/monospectral/noise/salt_and_pepper_noise.py +++ b/simpa/core/processing_components/monospectral/noise/salt_and_pepper_noise.py @@ -7,6 +7,7 @@ from simpa.core.processing_components import ProcessingComponent from simpa.utils.quality_assurance.data_sanity_testing import assert_array_well_defined import numpy as np +import torch class SaltAndPepperNoise(ProcessingComponent): @@ -38,9 +39,10 @@ def run(self, device): wavelength = self.global_settings[Tags.WAVELENGTH] data_array = load_data_field(self.global_settings[Tags.SIMPA_OUTPUT_PATH], data_field, wavelength) + data_tensor = torch.as_tensor(data_array, dtype=torch.float32, device=self.torch_device) - min_noise = np.min(data_array) - max_noise = np.max(data_array) + min_noise = torch.min(data_tensor).item() + max_noise = torch.max(data_tensor).item() noise_frequency = 0.01 if Tags.NOISE_FREQUENCY in self.component_settings.keys(): @@ -56,18 +58,17 @@ def run(self, device): self.logger.debug(f"Noise model max: {max_noise}") self.logger.debug(f"Noise model frequency: {noise_frequency}") - shape = np.shape(data_array) - num_data_points = int(np.round(np.prod(shape) * noise_frequency / 2)) - - coords_min = tuple([np.random.randint(0, i - 1, int(num_data_points)) for i in shape]) - coords_max = tuple([np.random.randint(0, i - 1, int(num_data_points)) for i in shape]) - - data_array[coords_min] = min_noise - data_array[coords_max] = max_noise + dist = torch.distributions.uniform.Uniform(torch.tensor(-1.0, dtype=torch.float32, device=self.torch_device), + torch.tensor(1.0, dtype=torch.float32, device=self.torch_device)) + sample = dist.sample(data_tensor.shape) + sample_cutoff = 1.0 - noise_frequency + data_tensor[sample > sample_cutoff] = min_noise + data_tensor[-sample > sample_cutoff] = max_noise if not (Tags.IGNORE_QA_ASSERTIONS in self.global_settings and Tags.IGNORE_QA_ASSERTIONS): - assert_array_well_defined(data_array) + assert_array_well_defined(data_tensor) - save_data_field(data_array, self.global_settings[Tags.SIMPA_OUTPUT_PATH], data_field, wavelength) + save_data_field(data_tensor.cpu().numpy().astype(np.float64, copy=False), + self.global_settings[Tags.SIMPA_OUTPUT_PATH], data_field, wavelength) self.logger.info("Applying Salt And Pepper Noise Model...[Done]") diff --git a/simpa/core/processing_components/monospectral/noise/uniform_noise.py b/simpa/core/processing_components/monospectral/noise/uniform_noise.py index 7cc5f799..5610dc43 100644 --- a/simpa/core/processing_components/monospectral/noise/uniform_noise.py +++ b/simpa/core/processing_components/monospectral/noise/uniform_noise.py @@ -7,6 +7,7 @@ from simpa.core.processing_components import ProcessingComponent from simpa.utils.quality_assurance.data_sanity_testing import assert_array_well_defined import numpy as np +import torch class UniformNoise(ProcessingComponent): @@ -52,15 +53,19 @@ def run(self, device): wavelength = self.global_settings[Tags.WAVELENGTH] data_array = load_data_field(self.global_settings[Tags.SIMPA_OUTPUT_PATH], data_field, wavelength) + data_tensor = torch.as_tensor(data_array, dtype=torch.float32, device=self.torch_device) + dist = torch.distributions.uniform.Uniform(torch.tensor(min_noise, dtype=torch.float32, device=self.torch_device), + torch.tensor(max_noise, dtype=torch.float32, device=self.torch_device)) if mode == Tags.NOISE_MODE_ADDITIVE: - data_array = data_array + (np.random.random(size=np.shape(data_array)) * (max_noise-min_noise) + min_noise) + data_tensor += dist.sample(data_tensor.shape) elif mode == Tags.NOISE_MODE_MULTIPLICATIVE: - data_array = data_array * (np.random.random(size=np.shape(data_array)) * (max_noise-min_noise) + min_noise) + data_tensor *= dist.sample(data_tensor.shape) if not (Tags.IGNORE_QA_ASSERTIONS in self.global_settings and Tags.IGNORE_QA_ASSERTIONS): - assert_array_well_defined(data_array) + assert_array_well_defined(data_tensor) - save_data_field(data_array, self.global_settings[Tags.SIMPA_OUTPUT_PATH], data_field, wavelength) + save_data_field(data_tensor.cpu().numpy().astype(np.float64, copy=False), + self.global_settings[Tags.SIMPA_OUTPUT_PATH], data_field, wavelength) self.logger.info("Applying Uniform Noise Model...[Done]") 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..0e922c95 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 @@ -182,6 +182,9 @@ 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") return cmd @staticmethod diff --git a/simpa/utils/quality_assurance/data_sanity_testing.py b/simpa/utils/quality_assurance/data_sanity_testing.py index 0433ff8a..32cee0be 100644 --- a/simpa/utils/quality_assurance/data_sanity_testing.py +++ b/simpa/utils/quality_assurance/data_sanity_testing.py @@ -3,7 +3,9 @@ # SPDX-License-Identifier: MIT import numpy as np +import torch import inspect +from typing import Union def assert_equal_shapes(numpy_arrays: list): @@ -29,13 +31,13 @@ def assert_equal_shapes(numpy_arrays: list): f" parameters. Called from {inspect.stack()[1].function}") -def assert_array_well_defined(array: np.ndarray, assume_non_negativity: bool = False, +def assert_array_well_defined(array: Union[np.ndarray, torch.Tensor], assume_non_negativity: bool = False, assume_positivity=False, array_name: str = None): """ - This method tests if all entries of the given array are well-defined (i.e. not np.inf, np.nan, or None). + This method tests if all entries of the given array or tensor are well-defined (i.e. not np.inf, np.nan, or None). The method can be parametrised to be more strict. - :param array: The input np.ndarray + :param array: The input np.ndarray or torch.Tensor :param assume_non_negativity: bool (default: False). If true, all values must be greater than or equal to 0. :param assume_positivity: bool (default: False). If true, all values must be greater than 0. :param array_name: a string that gives more information in case of an error. @@ -43,11 +45,12 @@ def assert_array_well_defined(array: np.ndarray, assume_non_negativity: bool = F """ error_message = None - if not np.isfinite(array).all(): + 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 (array <= 0).any(): + if assume_positivity and torch.any(array_as_tensor <= 0): error_message = "not positive" - if assume_non_negativity and (array < 0).any(): + if assume_non_negativity and torch.any(array_as_tensor < 0): error_message = "negative" if error_message: if array_name is None: diff --git a/simpa_tests/automatic_tests/test_noise_models.py b/simpa_tests/automatic_tests/test_noise_models.py index 61750380..b23ac50d 100644 --- a/simpa_tests/automatic_tests/test_noise_models.py +++ b/simpa_tests/automatic_tests/test_noise_models.py @@ -5,6 +5,7 @@ import unittest import os import numpy as np +import torch from simpa.core.processing_components.monospectral.noise import * from simpa.utils import Tags, Settings @@ -32,6 +33,7 @@ def validate_noise_model_results(self, noise_model, noise_model_settings, expected_mean, expected_std, error_margin=0.05): np.random.seed(self.RANDOM_SEED) + torch.manual_seed(self.RANDOM_SEED) settings = { # These parameters set the general propeties of the simulated volume @@ -78,7 +80,7 @@ def validate_noise_model_results(self, noise_model, noise_model_settings, f"The mean was not as expected. Expected {expected_mean} but was {actual_mean}") self.assertTrue(np.abs(actual_std - expected_std) < 1e-10 or np.abs(actual_std - expected_std) / expected_std < error_margin, - f"The mean was not as expected. Expected {expected_std} but was {actual_std}") + f"The std was not as expected. Expected {expected_std} but was {actual_std}") finally: if (os.path.exists(settings[Tags.SIMPA_OUTPUT_PATH]) and os.path.isfile(settings[Tags.SIMPA_OUTPUT_PATH])): @@ -87,8 +89,8 @@ def validate_noise_model_results(self, noise_model, noise_model_settings, def setUp(self): - self.VOLUME_WIDTH_IN_MM = 10 - self.VOLUME_HEIGHT_IN_MM = 10 + self.VOLUME_WIDTH_IN_MM = 40 + self.VOLUME_HEIGHT_IN_MM = 40 self.SPACING = 1 self.RANDOM_SEED = 4711