fix: single_tif_phasics phase data ambiguity (recorded at different t…

…ime points)

CHANGELOG

@@ -1,3 +1,9 @@
+0.3.5
+ - fix: single_tif_phaiscs (SID4Bio) contains two phase images, the second
+   of which is recorded at a different time point than the intensity
+   image. The first image is recorded in wavelengths and not in
+   nanometers and thus is converted using phasics metadata first.
+   If this metadata is not available, the second image is used.
 0.3.4
  - fix: qpimage file formats: override identifiers (clean solution)
  - fix: add check for valid meta_data keys

qpformat/file_formats/single_tif_phasics.py

@@ -37,12 +37,9 @@ class SingleTifPhasics(SingleData):
     def __init__(self, path, meta_data={}, *args, **kwargs):
         if "wavelength" not in meta_data:
             # get wavelength if not given
-            wl_str = SingleTifPhasics._get_meta_data(path=path,
-                                                     section="analyse data",
-                                                     name="lambda(nm)")
-            if wl_str:
-                wavelength = float(wl_str) * 1e-9
-                meta_data["wavelength"] = wavelength
+            wl = self._get_wavelength(path)
+            if not np.isnan(wl):
+                meta_data["wavelength"] = wl
             else:
                 # We need the wavelength to convert OPD to phase
                 msg = "'wavelength' must be specified in `meta_data`!"
@@ -85,6 +82,16 @@ def _get_tif(path):
             path = fspath(path)
         return tifffile.TiffFile(path)
 
+    def _get_wavelength(self, path):
+        wl_str = SingleTifPhasics._get_meta_data(path=path,
+                                                 section="analyse data",
+                                                 name="lambda(nm)")
+        if wl_str:
+            wavelength = float(wl_str) * 1e-9
+        else:
+            wavelength = np.nan
+        return wavelength
+
     def get_qpimage_raw(self, idx=0):
         """Return QPImage without background correction"""
         # Load experimental data
@@ -100,16 +107,42 @@ def get_qpimage_raw(self, idx=0):
             blc = INTENSITY_BASELINE_CLAMP
             inten = tf.pages[0].asarray() * (imax - imin) / isamp + imin - blc
             # Phase
-            phatags = tf.pages[2].tags
+            # The SID4Bio records two phase images, one in wavelengths and
+            # one in nanometers. Surprisingly, these two phase images are
+            # not derived from the same recording. The first image (pages[1])
+            # (in wavelengths) matches the intensity image (pages[0]). The
+            # second image (pages[2]) is recorded at a different time point.
+            # Initially, I thought it would be best to compute the phase
+            # directly from the measured value in nanometers (pages[2]) using
+            # the known wavelength given by the qpformat user. However, since
+            # phase and amplitude won't match up in that case, the wavelength
+            # phase image (pages[1]) has to be used. Since phasics uses its own
+            # wavelength (set by the user in the acquisition/extraction
+            # software) which might be wrong, I decided to first compute
+            # the phase in nanometers from tf.pages[1] using the phasics
+            # wavelength and then proceed as before, computing the phase
+            # in radians using the correct, user-given wavelength.
+            wl_phasics = self._get_wavelength(self.path)
+            if not np.isnan(wl_phasics):
+                # proceed with phase in wavelengths
+                phaid = 1
+            else:
+                # proceed with phase in nanometers
+                phaid = 2
+            phatags = tf.pages[phaid].tags
             pmin = phatags["61243"].value
             pmax = phatags["61242"].value
             psamp = phatags["max_sample_value"].value
             if psamp == 0 or pmin == pmax:
                 # no phase data
                 pha = np.zeros_like(inten)
             else:
-                # optical path difference is in nm
-                opd = tf.pages[2].asarray() * (pmax - pmin) / psamp + pmin
+                # optical path difference
+                opd = tf.pages[phaid].asarray() * (pmax - pmin) / psamp + pmin
+                if phaid == 1:  # convert [wavelengths] to [nm]
+                    assert not np.isnan(wl_phasics)
+                    opd *= wl_phasics * 1e9
+                # convert from [nm] to [rad]
                 pha = opd / (self.meta_data["wavelength"] * 1e9) * 2 * np.pi
 
         meta_data = copy.copy(self.meta_data)
@@ -154,6 +187,8 @@ def verify(path):
                 "61243" in tf.pages[0].tags and
                 "61242" in tf.pages[0].tags and
                 "61238" in tf.pages[0].tags and
+                "61243" in tf.pages[1].tags and
+                "61242" in tf.pages[1].tags and
                 "max_sample_value" in tf.pages[0].tags and
                 (tf.pages[0].tags["61242"].value !=
                  tf.pages[1].tags["61242"].value)):

tests/test_series_zip_tif_phasics.py

@@ -33,7 +33,7 @@ def test_data_content():
     qpi = ds.get_qpimage(0)
     assert qpi.meta["wavelength"] == 550e-9
     assert np.allclose(qpi.amp.max(), 183.96992660669972)
-    assert np.allclose(qpi.pha.max() - qpi.pha.min(), 0.189023369599572)
+    assert np.allclose(qpi.pha.max() - qpi.pha.min(), 0.18902349472045898)
 
 
 if __name__ == "__main__":

tests/test_single_tif_phasics.py

@@ -62,7 +62,7 @@ def test_data_content():
     qpi = ds.get_qpimage()
     assert qpi.meta["wavelength"] == 550e-9
     assert np.allclose(qpi.amp.max(), 188.57930365845519)
-    assert np.allclose(qpi.pha.max() - qpi.pha.min(), 4.1683941115690617)
+    assert np.allclose(qpi.pha.max() - qpi.pha.min(), 4.168394088745117)
 
 
 if __name__ == "__main__":