T412 scattering equation bug fix

simpa/utils/libraries/literature_values.py

@@ -149,6 +149,16 @@ class OpticalTissueProperties:
     #   year={2013},
     #   publisher={IOP Publishing}
     # }
+    # @article{bashkatov2011optical,
+    #   title={Optical properties of skin, subcutaneous, and muscle tissues: a review},
+    #   author={Bashkatov, Alexey N and Genina, Elina A and Tuchin, Valery V},
+    #   journal={Journal of innovative optical health sciences},
+    #   volume={4},
+    #   number={01},
+    #   pages={9--38},
+    #   year={2011},
+    #   publisher={World Scientific}
+    # }
     MUS500_BACKGROUND_TISSUE = 191.0  # Table 2: Average over all other soft tissue
     FRAY_BACKGROUND_TISSUE = 0.153  # Table 2: Average over all other soft tissue
     BMIE_BACKGROUND_TISSUE = 1.091  # Table 2: Average over all other soft tissue
@@ -159,10 +169,10 @@ class OpticalTissueProperties:
 
     MUS500_EPIDERMIS = 93.01  # Bashkatov et al. 2011 but adjusted for epidermis anisotropy
     FRAY_EPIDERMIS = 0.29  # Table 1; Salomatina et al 2006; One value for epidermis
-    BMIE_EPIDERMIS = 2.8  # Table 1; Salomatina et al 2006; One value for epidermis
+    BMIE_EPIDERMIS = 1.5  # Table 1; Salomatina et al 2006; Altered slightly for the slope to fit with BASHKATOV et al., 2011
     MUS500_DERMIS = 175.0  # Bashkatov et al. 2011 but adjusted for DERMIS_ANISOTROPY
     FRAY_DERMIS = 0.1  # Table 1; Salomatina et al 2006; One value for dermis
-    BMIE_DERMIS = 3.5  # Table 1; Salomatina et al 2006; One value for dermis
+    BMIE_DERMIS = 2.7  # Table 1; Salomatina et al 2006; Adjusted slightly to fit BASHKATOV et al., 2011
     MUS500_FAT = 193.0  # Table 2 average fatty tissue
     FRAY_FAT = 0.174  # Table 2 average fatty tissue
     BMIE_FAT = 0.447  # Table 2 average fatty tissue

simpa/utils/libraries/spectrum_library.py

@@ -247,7 +247,7 @@ def scattering_from_rayleigh_and_mie_theory(name: str, mus_at_500_nm: float = 1.
         :return: A Spectrum instance based on Rayleigh and Mie scattering theory.
         """
         wavelengths = np.arange(450, 1001, 1)
-        scattering = (mus_at_500_nm * (fraction_rayleigh_scattering * (wavelengths / 500) ** 1e-4 +
+        scattering = (mus_at_500_nm * (fraction_rayleigh_scattering * (wavelengths / 500) ** -4 +
                       (1 - fraction_rayleigh_scattering) * (wavelengths / 500) ** -mie_power_law_coefficient))
         return Spectrum(name, wavelengths, scattering)
 

simpa_tests/automatic_tests/tissue_library/test_scattering_equation.py

@@ -0,0 +1,60 @@
+import unittest
+from simpa.utils import ScatteringSpectrumLibrary
+
+
+class TestScatteringEquation(unittest.TestCase):
+    """
+    This test was written to ensure #412 does not reappear again
+    """
+    def test_scattering_equation_is_correct(self):
+        lib = ScatteringSpectrumLibrary()
+
+        # Testing both mixed
+        test_spectrum = lib.scattering_from_rayleigh_and_mie_theory("test",
+                                                    mus_at_500_nm=10,
+                                                    fraction_rayleigh_scattering=0.5,
+                                                    mie_power_law_coefficient=0.1)
+
+        self.assertAlmostEqual(test_spectrum.get_value_for_wavelength(500), 10.0, 5)
+        self.assertAlmostEqual(test_spectrum.get_value_for_wavelength(750), 5.788976824948744, 5)
+        self.assertAlmostEqual(test_spectrum.get_value_for_wavelength(1000), 4.977664957684037, 5)
+
+        # Testing Rayleigh 0
+        test_spectrum = lib.scattering_from_rayleigh_and_mie_theory("test",
+                                                                      mus_at_500_nm=100,
+                                                                      fraction_rayleigh_scattering=0.0,
+                                                                      mie_power_law_coefficient=10)
+
+        self.assertAlmostEqual(test_spectrum.get_value_for_wavelength(500), 100.0, 5)
+        self.assertAlmostEqual(test_spectrum.get_value_for_wavelength(750), 1.7341529915832612, 5)
+        self.assertAlmostEqual(test_spectrum.get_value_for_wavelength(1000), 0.09765625, 5)
+
+        # testing no decay at all
+        test_spectrum = lib.scattering_from_rayleigh_and_mie_theory("test",
+                                                                    mus_at_500_nm=100,
+                                                                    fraction_rayleigh_scattering=0.0,
+                                                                    mie_power_law_coefficient=0)
+
+        self.assertAlmostEqual(test_spectrum.get_value_for_wavelength(500), 100.0, 5)
+        self.assertAlmostEqual(test_spectrum.get_value_for_wavelength(750), 100.0, 5)
+        self.assertAlmostEqual(test_spectrum.get_value_for_wavelength(1000), 100.0, 5)
+
+        # testing only Rayleigh
+        test_spectrum = lib.scattering_from_rayleigh_and_mie_theory("test",
+                                                                    mus_at_500_nm=100,
+                                                                    fraction_rayleigh_scattering=0.75,
+                                                                    mie_power_law_coefficient=0)
+
+        self.assertAlmostEqual(test_spectrum.get_value_for_wavelength(500), 100.0, 5)
+        self.assertAlmostEqual(test_spectrum.get_value_for_wavelength(750), 39.81481481481482, 5)
+        self.assertAlmostEqual(test_spectrum.get_value_for_wavelength(1000), 29.6875, 5)
+
+        # testing scaling of zero scattering
+        test_spectrum = lib.scattering_from_rayleigh_and_mie_theory("test",
+                                                                    mus_at_500_nm=0,
+                                                                    fraction_rayleigh_scattering=0.2342345,
+                                                                    mie_power_law_coefficient=0.123123)
+
+        self.assertAlmostEqual(test_spectrum.get_value_for_wavelength(500), 0.0, 5)
+        self.assertAlmostEqual(test_spectrum.get_value_for_wavelength(750), 0.0, 5)
+        self.assertAlmostEqual(test_spectrum.get_value_for_wavelength(1000), 0.0, 5)

simpa_tests/automatic_tests/tissue_library/test_tissue_library_against_literature_values.py

@@ -10,9 +10,6 @@
     get_fully_deoxygenated_blood_reference_dictionary, \
     get_lymph_node_reference_dictionary
 
-# FIXME temporary workaround for newest Intel architectures
-import os
-os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE"
 
 VISUALISE = False
 

simpa_tests/test_utils/tissue_composition_tests.py

@@ -9,6 +9,8 @@
 from simpa.utils.libraries.tissue_library import TissueLibrary
 import numpy as np
 import matplotlib.patches as patches
+import matplotlib
+matplotlib.use("Agg")
 import matplotlib.pyplot as plt
 
 TEST_SETTINGS = Settings({
@@ -44,7 +46,7 @@ def compare_molecular_composition_against_expected_values(molecular_composition:
 
     validate_expected_values_dictionary(expected_values)
     if visualise_values:
-        plt.figure(figsize=(12, 8))
+        plt.figure(figsize=(12, 8), layout="constrained")
         plt.suptitle(title + f" [green=expected, blue=actual, red={tolerated_margin_in_percent*100}% margin]")
         num_subplots = len(property_tags)
 
@@ -79,14 +81,13 @@ def compare_molecular_composition_against_expected_values(molecular_composition:
                                              f"expected to be {expected_properties[tag]})")
 
     if visualise_values:
-        plt.tight_layout()
-        plt.show()
+        plt.savefig(f"{title}.png", dpi=300)
         plt.close()
 
 
 def get_epidermis_reference_dictionary():
     """
-    The
+    These values come from table 1; Epidermis (medium pigmented)
         @article{bashkatov2011optical,
               title={Optical properties of skin, subcutaneous, and muscle tissues: a review},
               author={Bashkatov, Alexey N and Genina, Elina A and Tuchin, Valery V},
@@ -163,7 +164,7 @@ def get_epidermis_reference_dictionary():
 
     values700nm = TissueProperties(TEST_SETTINGS)
     values700nm[Tags.DATA_FIELD_ABSORPTION_PER_CM] = 3.07
-    values700nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 54.66
+    values700nm[Tags.DATA_FIELD_SCATTERING_PER_CM] = 47.4 # altered from the original table to fit exponential scattering decay.
     values700nm[Tags.DATA_FIELD_ANISOTROPY] = 0.804
     values700nm[Tags.DATA_FIELD_GRUNEISEN_PARAMETER] = calculate_gruneisen_parameter_from_temperature(37.0)
     values700nm[Tags.DATA_FIELD_SEGMENTATION] = SegmentationClasses.EPIDERMIS