Merge pull request #267 from IMSY-DKFZ/main

Update develop

.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/

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)

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]")

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]")

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]")

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]")

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]")

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

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,25 +31,26 @@ 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.
     :raises AssertionError: if there are any unexpected values in the given array.
     """
 
     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:

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