Merge branch 'refs/heads/develop' into T154_heterogeneous_tissue

# Conflicts:
#	.gitignore
#	pyproject.toml
#	simpa/utils/calculate.py
#	simpa/utils/libraries/molecule_library.py
#	simpa_tests/automatic_tests/tissue_library/test_core_assumptions.py

.gitignore

@@ -17,6 +17,9 @@ simpa_tests/figures
 simpa_tests/figures/*
 simpa_tests/**/figures
 simpa_tests/**/figures/*
+simpa_tests/manual_tests/figures
+simpa_tests/manual_tests/reference_figures
+simpa_tests/manual_tests/manual_tests_*
 
 simpa_examples/benchmarking/*.csv
 simpa_examples/benchmarking/*.txt

README.md

@@ -43,6 +43,7 @@ The SIMPA path management takes care of that.
 * [SIMPA installation instructions](#simpa-installation-instructions)
 * [External tools installation instructions](#external-tools-installation-instructions)
 * [Path Management](#path-management)
+* [Testing](#run-manual-tests)
 
 ## SIMPA installation instructions
 
@@ -114,6 +115,9 @@ one we provided in the `simpa_examples/path_config.env.example`) in the followin
 
 For this option, please follow the instructions in the `simpa_examples/path_config.env.example` file. 
 
+## Run manual tests
+To check the success of your installation ot to assess how your contributions affect the Simpa simulation outcomes, you can run the manual tests automatically. Install the testing requirements by doing `pip install .[testing]` and run the `simpa_tests/manual_tests/generate_overview.py` file. This script runs all manual tests and generates both a markdown and an HTML file that compare your results with the reference results.
+
 # Simulation examples
 
 To get started with actual simulations, SIMPA provides an [example package](simpa_examples) of simple simulation 

docs/source/introduction.md

@@ -7,6 +7,7 @@ The SIMPA path management takes care of that.
 * [SIMPA installation instructions](#simpa-installation-instructions)
 * [External tools installation instructions](#external-tools-installation-instructions)
 * [Path Management](#path-management)
+* [Testing](#run-manual-tests)
 
 ## SIMPA installation instructions
 
@@ -78,6 +79,9 @@ one we provided in the `simpa_examples/path_config.env.example`) in the followin
 
 For this option, please follow the instructions in the `simpa_examples/path_config.env.example` file. 
 
+## Run manual tests
+To check the success of your installation ot to assess how your contributions affect the Simpa simulation outcomes, you can run the manual tests automatically. Install the testing requirements by doing `pip install .[testing]` and run the `simpa_tests/manual_tests/generate_overview.py` file. This script runs all manual tests and generates both a markdown and an HTML file that compare your results with the reference results.
+
 # Simulation examples
 
 To get started with actual simulations, SIMPA provides an [example package](simpa_examples) of simple simulation 

pyproject.toml

@@ -47,6 +47,11 @@ profile = [
     "line_profiler>=4.0.0,<5.0.0",
     "memory_profiler>=0.61.0,<0.62.0"
 ]
+testing = [
+    "mdutils>=1.4.0",          # Uses MIT-License (MIT compatible)
+    "pypandoc>=1.13",          # Uses MIT-License (MIT compatible)
+    "pypandoc_binary>=1.13"    # Uses MIT-License (MIT compatible)
+    ]
 
 [project.urls]
 Homepage = "https://github.com/IMSY-DKFZ/simpa"

simpa/utils/calculate.py

@@ -3,47 +3,68 @@
 # SPDX-License-Identifier: MIT
 
 
-from typing import Union
+from typing import Union, List, Dict, Optional
 import numpy as np
 import torch
 from scipy.interpolate import interp1d
 
 
-def calculate_oxygenation(molecule_list: list) -> Union[float, int, torch.Tensor]:
+def extract_hemoglobin_fractions(molecule_list: List) -> Dict[str, float]:
     """
-    Calculate the oxygenation level based on the volume fractions of deoxyhaemoglobin and oxyhaemoglobin.
-
-    This function takes a list of molecules and returns an oxygenation value between 0 and 1 if computable,
-    otherwise returns None.
+    Extract hemoglobin volume fractions from a list of molecules.
 
     :param molecule_list: List of molecules with their spectrum information and volume fractions.
-    :return: An oxygenation value between 0 and 1 if possible, or None if not computable.
+    :return: A dictionary with hemoglobin types as keys and their volume fractions as values.
     """
-    hb = None    # Volume fraction of deoxyhaemoglobin
-    hbO2 = None  # Volume fraction of oxyhaemoglobin
+
+    # Put 0.0 as default value for both hemoglobin types in case they are not present in the molecule list.
+    hemoglobin = {
+        "Deoxyhemoglobin": 0.0,
+        "Oxyhemoglobin": 0.0
+    }
 
     for molecule in molecule_list:
-        if molecule.spectrum.spectrum_name == "Deoxyhemoglobin":
-            hb = molecule.volume_fraction
-        if molecule.spectrum.spectrum_name == "Oxyhemoglobin":
-            hbO2 = molecule.volume_fraction
+        spectrum_name = molecule.spectrum.spectrum_name
+        if spectrum_name in hemoglobin:
+            hemoglobin[spectrum_name] = molecule.volume_fraction
+
+    return hemoglobin
 
-    if hb is None and hbO2 is None:
-        return None
 
-    if hb is None:
-        hb = 0
-    elif hbO2 is None:
-        hbO2 = 0
+def calculate_oxygenation(molecule_list: List) -> Optional[float]:
+    """
+    Calculate the oxygenation level based on the volume fractions of deoxyhemoglobin and oxyhemoglobin.
 
+    :param molecule_list: List of molecules with their spectrum information and volume fractions.
+    :return: An oxygenation value between 0 and 1 if possible, or None if not computable.
+    """
+    hemoglobin = extract_hemoglobin_fractions(molecule_list)
+    hb, hbO2 = hemoglobin["Deoxyhemoglobin"], hemoglobin["Oxyhemoglobin"]
+
+    total = hb + hbO2
+
+    # Avoid division by zero. If none of the hemoglobin types are present, the oxygenation level is not computable.
     if isinstance(hb, torch.Tensor) or isinstance(hbO2, torch.Tensor):
-        return torch.where(hb + hbO2 < 1e-10, 0, hbO2 / (hb + hbO2))
+        return torch.where(total < 1e-10, 0, hbO2 / total)
 
     else:
-        if hb + hbO2 < 1e-10:
+        if total < 1e-10:
             return None
+        else:
+            return hbO2 / total
 
-    return hbO2 / (hb + hbO2)
+
+def calculate_bvf(molecule_list: List) -> Union[float, int]:
+    """
+    Calculate the blood volume fraction based on the volume fractions of deoxyhemoglobin and oxyhemoglobin.
+
+    :param molecule_list: List of molecules with their spectrum information and volume fractions.
+    :return: The blood volume fraction value between 0 and 1, or 0, if oxy and deoxy not present.
+    """
+    hemoglobin = extract_hemoglobin_fractions(molecule_list)
+    hb, hbO2 = hemoglobin["Deoxyhemoglobin"], hemoglobin["Oxyhemoglobin"]
+    # We can use the sum of hb and hb02 to compute blood volume fraction as the volume fraction of all molecules is 1.
+    return hb + hbO2
 
 
 def create_spline_for_range(xmin_mm: Union[float, int] = 0, xmax_mm: Union[float, int] = 10,

simpa/utils/constants.py

@@ -41,6 +41,7 @@ class SegmentationClasses:
     Tags.DATA_FIELD_GRUNEISEN_PARAMETER,
     Tags.DATA_FIELD_SEGMENTATION,
     Tags.DATA_FIELD_OXYGENATION,
+    Tags.DATA_FIELD_BLOOD_VOLUME_FRACTION,
     Tags.DATA_FIELD_DENSITY,
     Tags.DATA_FIELD_SPEED_OF_SOUND,
     Tags.DATA_FIELD_ALPHA_COEFF

simpa/utils/libraries/molecule_library.py

@@ -10,7 +10,7 @@
 from simpa.utils.libraries.literature_values import OpticalTissueProperties, StandardProperties
 from simpa.utils.libraries.spectrum_library import AnisotropySpectrumLibrary, ScatteringSpectrumLibrary
 from simpa.utils import Spectrum
-from simpa.utils.calculate import calculate_oxygenation, calculate_gruneisen_parameter_from_temperature
+from simpa.utils.calculate import calculate_oxygenation, calculate_gruneisen_parameter_from_temperature, calculate_bvf
 from simpa.utils.serializer import SerializableSIMPAClass
 from simpa.utils.libraries.spectrum_library import AbsorptionSpectrumLibrary
 from simpa.log import Logger
@@ -59,6 +59,7 @@ def update_internal_properties(self, settings):
         self.internal_properties = TissueProperties(settings)
         self.internal_properties[Tags.DATA_FIELD_SEGMENTATION] = self.segmentation_type
         self.internal_properties[Tags.DATA_FIELD_OXYGENATION] = calculate_oxygenation(self)
+        self.internal_properties[Tags.DATA_FIELD_BLOOD_VOLUME_FRACTION] = calculate_bvf(self)
         search_list = self.copy()
 
         for molecule in search_list:
@@ -93,10 +94,8 @@ def get_properties_for_wavelength(self, settings, wavelength) -> TissuePropertie
         for molecule in search_list:
             self.internal_properties[Tags.DATA_FIELD_ABSORPTION_PER_CM] += \
                 (molecule.volume_fraction * molecule.spectrum.get_value_for_wavelength(wavelength))
-
             self.internal_properties[Tags.DATA_FIELD_SCATTERING_PER_CM] += \
                 (molecule.volume_fraction * (molecule.scattering_spectrum.get_value_for_wavelength(wavelength)))
-
             self.internal_properties[Tags.DATA_FIELD_ANISOTROPY] += \
                 molecule.volume_fraction * molecule.anisotropy_spectrum.get_value_for_wavelength(wavelength)
 

simpa/utils/tags.py

@@ -881,6 +881,12 @@ class Tags:
     Usage: SIMPA package, naming convention
     """
 
+    DATA_FIELD_BLOOD_VOLUME_FRACTION = "bvf"
+    """
+    Blood volume fraction of the generated volume/structure.\n
+    Usage: SIMPA package, naming convention
+    """
+
     DATA_FIELD_SEGMENTATION = "seg"
     """
     Segmentation of the generated volume/structure.\n

simpa_tests/automatic_tests/test_calculation_utils.py

@@ -5,7 +5,7 @@
 import unittest
 from simpa.utils import SegmentationClasses, MolecularCompositionGenerator
 from simpa.utils.libraries.molecule_library import MOLECULE_LIBRARY
-from simpa.utils.calculate import calculate_oxygenation
+from simpa.utils.calculate import calculate_oxygenation, calculate_bvf
 from simpa.utils.calculate import randomize_uniform
 from simpa.utils.calculate import calculate_gruneisen_parameter_from_temperature
 from simpa.utils.calculate import positive_gauss
@@ -53,7 +53,7 @@ def test_oxygenation_calculation(self):
         assert abs(oxy_value) < 1e-5, ("oxy value was not 0.0 but " + str(oxy_value))
 
         # RANDOM CASES
-        for i in range(100):
+        for _ in range(100):
             oxy = np.random.random()
             deoxy = np.random.random()
             mcg = MolecularCompositionGenerator()
@@ -65,6 +65,55 @@ def test_oxygenation_calculation(self):
 
             assert abs(sO2_value - (oxy / (oxy + deoxy))) < 1e-10
 
+    def test_bvf_calculation(self):
+
+        # Neither oxy nor deoxy:
+        mcg = MolecularCompositionGenerator()
+        mcg.append(MOLECULE_LIBRARY.fat(1.0))
+        bvf_value = calculate_bvf(mcg.get_molecular_composition(segmentation_type=SegmentationClasses.GENERIC))
+        assert bvf_value == 0
+
+        mcg = MolecularCompositionGenerator()
+        mcg.append(MOLECULE_LIBRARY.fat(1.0))
+        mcg.append(MOLECULE_LIBRARY.oxyhemoglobin(0.0))
+        bvf_value = calculate_bvf(mcg.get_molecular_composition(segmentation_type=SegmentationClasses.GENERIC))
+        assert bvf_value == 0
+
+        # Just oxyhemoglobin CASES:
+        mcg = MolecularCompositionGenerator()
+        mcg.append(MOLECULE_LIBRARY.oxyhemoglobin(1.0))
+        oxy_hemo = mcg.get_molecular_composition(segmentation_type=SegmentationClasses.GENERIC)
+        bvf_value = calculate_bvf(oxy_hemo)
+        assert bvf_value == 1.0
+        mcg.append(MOLECULE_LIBRARY.deoxyhemoglobin(0.0))
+        oxy_hemo = mcg.get_molecular_composition(segmentation_type=SegmentationClasses.GENERIC)
+        bvf_value = calculate_bvf(oxy_hemo)
+        assert bvf_value == 1.0
+
+        # Just deoxyhemoglobin CASES:
+        mcg = MolecularCompositionGenerator()
+        mcg.append(MOLECULE_LIBRARY.deoxyhemoglobin(1.0))
+        deoxy_hemo = mcg.get_molecular_composition(segmentation_type=SegmentationClasses.GENERIC)
+        bvf_value = calculate_bvf(deoxy_hemo)
+        assert bvf_value == 1.0
+        mcg.append(MOLECULE_LIBRARY.oxyhemoglobin(0.0))
+        deoxy_hemo = mcg.get_molecular_composition(segmentation_type=SegmentationClasses.GENERIC)
+        bvf_value = calculate_bvf(deoxy_hemo)
+        assert bvf_value == 1.0
+
+        # RANDOM CASES
+        for _ in range(100):
+            oxy = np.random.random()
+            deoxy = np.random.random()
+            fat = np.random.random()
+            sum_oxy_deoxy_fat = oxy + deoxy + fat
+            mcg = MolecularCompositionGenerator()
+            mcg.append(MOLECULE_LIBRARY.fat(fat/sum_oxy_deoxy_fat))
+            mcg.append(MOLECULE_LIBRARY.deoxyhemoglobin(deoxy/sum_oxy_deoxy_fat))
+            mcg.append(MOLECULE_LIBRARY.oxyhemoglobin(oxy/sum_oxy_deoxy_fat))
+            bvf_value = calculate_bvf(mcg.get_molecular_composition(segmentation_type=SegmentationClasses.GENERIC))
+            assert abs(bvf_value - (oxy+deoxy)/sum_oxy_deoxy_fat) < 1e-10
+
     def test_randomize(self):
         for _ in range(1000):
             lower = np.random.randint(0, 10000000)

simpa_tests/automatic_tests/tissue_library/test_core_assumptions.py

@@ -2,11 +2,15 @@
 # SPDX-FileCopyrightText: 2021 Janek Groehl
 # SPDX-License-Identifier: MIT
 
+import inspect
 import unittest
+
+import numpy as np
+
 from simpa.utils import Tags, Settings, TISSUE_LIBRARY
+from simpa.utils.calculate import calculate_bvf, calculate_oxygenation
 from simpa.utils.libraries.molecule_library import MolecularComposition
-import inspect
-import numpy as np
+from simpa.utils.libraries.tissue_library import TissueLibrary
 
 TEST_SETTINGS = Settings({
     # These parameters set the general properties of the simulated volume
@@ -27,10 +31,39 @@ def test_volume_fractions_sum_to_less_or_equal_one(self):
             self.assertTrue((np.abs(total_volume_fraction-1.0) < 1e-3).all(),
                             f"Volume fraction not 1.0 +/- 0.001 for {method_name}")
 
+    def test_bvf_and_oxygenation_consistency(self):
+        # blood_volume_fraction (bvf) and oxygenation of tissue classes defined
+        # as input have to be the same as the calculated ones
+
+        def compare_input_with_calculations(test_tissue, oxy, bvf):
+            calculated_bvf = calculate_bvf(test_tissue)
+            calculated_sO2 = calculate_oxygenation(test_tissue)
+            if bvf < 1e-10:
+                assert calculated_sO2 is None
+                assert abs(bvf - calculated_bvf) < 1e-10
+            else:
+                assert abs(oxy - calculated_sO2) < 1e-10
+                assert abs(bvf - calculated_bvf) < 1e-10
+
+        # Test edge cases and a random one
+        oxy_values = [0., 0., 1., 1., np.random.random()]
+        bvf_values = [0., 1., 0., 1., np.random.random()]
+        for oxy in oxy_values:
+            # assert blood only with varying oxygenation_values
+            test_tissue = TISSUE_LIBRARY.blood(oxygenation=oxy)
+            compare_input_with_calculations(test_tissue, oxy, 1.)
+            # assert all other tissue classes with varying oxygenation- and bvf_values
+            for bvf in bvf_values:
+                for (_, method) in self.get_all_tissue_library_methods():
+                    args = inspect.getfullargspec(method).args
+                    if "background_oxy" in args and "blood_volume_fraction" in args:
+                        test_tissue = method(TISSUE_LIBRARY, background_oxy=oxy, blood_volume_fraction=bvf)
+                        compare_input_with_calculations(test_tissue, oxy, bvf)
+
     @staticmethod
     def get_all_tissue_library_methods():
         methods = []
-        for method in inspect.getmembers(TISSUE_LIBRARY, predicate=inspect.isfunction):
+        for method in inspect.getmembers(TissueLibrary, predicate=inspect.isfunction):
             if isinstance(method[1](TISSUE_LIBRARY), MolecularComposition):
                 methods.append(method)
         return methods

simpa_tests/manual_tests/digital_device_twins/VisualiseDevices.py

@@ -19,13 +19,21 @@ def tear_down(self):
         pass
 
     def visualise_result(self, show_figure_on_screen=True, save_path=None):
+        if show_figure_on_screen:
+            figure_save_path = [None, None, None]
+        else:
+            if save_path is None:
+                save_path = ""
+            figure_save_path = [save_path + "device_visualisation_MSOT_Acuity.png",
+                                save_path + "device_visualisation_MSOT_Invision.png",
+                                save_path + "device_visualisation_RSOM_Explorer.png"
+                                ]
         sp.visualise_device(sp.MSOTAcuityEcho(device_position_mm=np.asarray([50, 10, 0])),
-                            save_path + "device_visualisation_MSOT_Acuity.png")
+                            figure_save_path[0])
         sp.visualise_device(sp.InVision256TF(device_position_mm=np.asarray([50, 10, 50])),
-                            save_path + "device_visualisation_MSOT_Invision.png")
+                            figure_save_path[1])
         sp.visualise_device(sp.RSOMExplorerP50(device_position_mm=np.asarray([50, 10, 0])),
-                            save_path + "device_visualisation_RSOM_Explorer.png")
-
+                            figure_save_path[2])
 
 if __name__ == "__main__":
     test = DeviceVisualisationTest()