T115 include bvf

simpa/utils/calculate.py

@@ -3,43 +3,64 @@
 # 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.
+    Extract hemoglobin volume fractions from a list of molecules.
+
+    :param molecule_list: List of molecules with their spectrum information and volume fractions.
+    :return: A dictionary with hemoglobin types as keys and their volume fractions as values.
+    """
+
+    # 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:
+        spectrum_name = molecule.spectrum.spectrum_name
+        if spectrum_name in hemoglobin:
+            hemoglobin[spectrum_name] = molecule.volume_fraction
+
+    return hemoglobin
 
-    This function takes a list of molecules and returns an oxygenation value between 0 and 1 if computable,
-    otherwise returns None.
 
+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.
     """
-    hb = None    # Volume fraction of deoxyhaemoglobin
-    hbO2 = None  # Volume fraction of oxyhaemoglobin
+    hemoglobin = extract_hemoglobin_fractions(molecule_list)
+    hb, hbO2 = hemoglobin["Deoxyhemoglobin"], hemoglobin["Oxyhemoglobin"]
 
-    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
+    total = hb + hbO2
 
-    if hb is None and hbO2 is None:
+    # Avoid division by zero. If none of the hemoglobin types are present, the oxygenation level is not computable.
+    if total < 1e-10:
         return None
 
-    if hb is None:
-        hb = 0
-    elif hbO2 is None:
-        hbO2 = 0
+    return hbO2 / total
 
-    if hb + hbO2 < 1e-10:  # negative values are not allowed and division by (approx) zero
-        return None        # will lead to negative side effects.
 
-    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

@@ -9,7 +9,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
 
@@ -55,7 +55,7 @@ def update_internal_properties(self):
         self.internal_properties = TissueProperties()
         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)
         for molecule in self:
             self.internal_properties.volume_fraction += molecule.volume_fraction
             self.internal_properties[Tags.DATA_FIELD_GRUNEISEN_PARAMETER] += \

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 TISSUE_LIBRARY
-from simpa.utils.libraries.tissue_library import TissueLibrary
+from simpa.utils.calculate import calculate_bvf, calculate_oxygenation
 from simpa.utils.libraries.molecule_library import MolecularComposition
-import inspect
+from simpa.utils.libraries.tissue_library import TissueLibrary
 
 
 class TestCoreAssumptions(unittest.TestCase):
@@ -19,6 +23,35 @@ def test_volume_fractions_sum_to_less_or_equal_one(self):
             self.assertAlmostEqual(total_volume_fraction, 1.0, 3,
                                    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 = []