implement phasics SID PHA tif file format

qpformat/core.py

@@ -11,7 +11,7 @@ def guess_format(path):
 
 
 def load_data(path, fmt=None, bg_path=None, bg_fmt=None,
-              h5out=None):
+              meta_data={}, h5out=None):
     """Load experimental data
 
     Parameters
@@ -27,13 +27,13 @@ def load_data(path, fmt=None, bg_path=None, bg_fmt=None,
         The file format to use (see `file_formats.formats`)
         for the background. If set to `None`, the file format
         is be guessed.
+    meta_data: dict
+        Meta data (see `qpimage.meta.VALID_META_KEYS`)
     h5out: str
         Path to an hdf5 output file where all data
         is written in the :py:mod:`qpimage` data
         file format. If set to `None`, nothing
         is written to disk.
-    meta_data: dict
-        Meta data (see `qpimage.meta.VALID_META_KEYS`)
 
     Returns
     -------
@@ -44,12 +44,12 @@ def load_data(path, fmt=None, bg_path=None, bg_fmt=None,
     if fmt is None:
         fmt = guess_format(path)
 
-    dataobj = formats_dict[fmt](path)
+    dataobj = formats_dict[fmt](path=path, meta_data=meta_data)
 
     if bg_path is not None:
         if bg_fmt is None:
             bg_fmt = guess_format(bg_path)
-        bgobj = formats_dict[bg_fmt](bg_path)
+        bgobj = formats_dict[bg_fmt](path=bg_path, meta_data=meta_data)
         dataobj.set_bg(bgobj)
 
     if h5out is not None:

qpformat/file_formats/__init__.py

@@ -5,7 +5,7 @@
 # from .group_hdf5_qpimage import GroupHdf5Qpimage
 # from .group_zip_tif_phasics import GroupZipTifPhasics
 from .single_hdf5_qimage import SingleHdf5Qpimage
-# from .single_tif_phasics import SingleTifPhasics
+from .single_tif_phasics import SingleTifPhasics
 
 
 class GroupFolder(DataSet):
@@ -91,7 +91,7 @@ class UnknownFileFormatError(BaseException):
            # GroupZipTifPhasics,
            # GroupHdf5Qpimage,
            SingleHdf5Qpimage,
-           # SingleTifPhasics,
+           SingleTifPhasics,
            ]
 
 # convenience dictionary

qpformat/file_formats/dataset.py

@@ -39,11 +39,11 @@ def __repr__(self):
     def __len__(self):
         """Return number of samples of a data set"""
 
-    def get_name(self, idx):
+    def get_name(self, idx=0):
         """Return name of data at index `idx`"""
         return "{} [{}]".format(self.path, idx)
 
-    def get_qpimage(self, idx):
+    def get_qpimage(self, idx=0):
         """Return background-corrected QPImage of data at index `idx`"""
         # raw data
         qpi = self.get_qpimage_raw(idx)
@@ -63,10 +63,10 @@ def get_qpimage(self, idx):
         return qpi
 
     @abc.abstractmethod
-    def get_qpimage_raw(self, idx):
+    def get_qpimage_raw(self, idx=0):
         """Return QPImage without background correction"""
 
-    def get_time(self, idx):
+    def get_time(self, idx=0):
         """Return time of data at in dex `idx`
 
         By default, this returns zero and must be

qpformat/file_formats/single_hdf5_qimage.py

@@ -8,8 +8,10 @@ class SingleHdf5Qpimage(DataSet):
     def __len__(self):
         return 1
 
-    def get_qpimage(self, idx):
+    def get_qpimage(self, idx=0):
         """Return background-corrected QPImage of data at index `idx`"""
+        if idx != 0:
+            raise ValueError("Single file format, only one entry (`idx!=0`)!")
         if self._bgdata:
             # The user has explicitly chosen different background data
             # using `get_qpimage_raw`.
@@ -29,6 +31,8 @@ def get_qpimage_raw(self, idx=0):
 
     def get_time(self, idx=0):
         """Return the time of the QPImage data"""
+        if idx != 0:
+            raise ValueError("Single file format, only one entry (`idx!=0`)!")
         qpi = qpimage.QPImage(h5file=self.path, h5mode="r")
         if "time" in qpi.meta:
             return qpi.meta["time"]

qpformat/file_formats/single_tif_phasics.py

@@ -0,0 +1,142 @@
+import copy
+import time
+import xml.etree.ElementTree as ET
+
+import numpy as np
+import qpimage
+from skimage.external import tifffile
+
+from .dataset import DataSet
+
+
+# baseline clamp intensity normalization for phasics tif files
+INTENSITY_BASELINE_CLAMP = 150
+
+
+class LoadTifPhasicsError(BaseException):
+    pass
+
+
+class SingleTifPhasics(DataSet):
+    def __init__(self, path, meta_data={}):
+        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
+            else:
+                # We need the wavelength to convert OPD to phase
+                msg = "'wavelength' must be specified in `meta_data`!"
+                raise LoadTifPhasicsError(msg)
+
+        super(SingleTifPhasics, self).__init__(path=path,
+                                               meta_data=meta_data)
+
+    def __len__(self):
+        return 1
+
+    def get_qpimage_raw(self, idx=0):
+        """Return QPImage without background correction"""
+        if idx != 0:
+            raise ValueError("Single file format, only one entry (`idx!=0`)!")
+        # Load experimental data
+        with tifffile.TiffFile(self.path) as tf:
+            # page 0 contains intensity
+            # page 1 contains phase in nm
+            # page 2 contains phase in wavelengths
+            # Intensity:
+            inttags = tf.pages[0].tags
+            imin = inttags["61243"].value
+            imax = inttags["61242"].value
+            isamp = inttags["max_sample_value"].value
+            blc = INTENSITY_BASELINE_CLAMP
+            inten = tf.pages[0].asarray() * (imax - imin) / isamp + imin - blc
+            # Phase
+            phatags = tf.pages[2].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
+                pha = opd / (self.meta_data["wavelength"] * 1e9) * 2 * np.pi
+
+        meta_data = copy.copy(self.meta_data)
+        if "time" not in meta_data:
+            meta_data["time"] = self.get_time(idx)
+        qpi = qpimage.QPImage(data=(pha, inten),
+                              which_data="phase,intensity",
+                              meta_data=meta_data)
+        return qpi
+
+    def get_time(self, idx=0):
+        """Return the time of the tif data since the epoch
+
+        The time is stored in the "61238" tag.
+        """
+        timestr = SingleTifPhasics._get_meta_data(path=self.path,
+                                                  section="acquisition info",
+                                                  name="date & heure")
+        if timestr is not None:
+            timestr = timestr.split(".")
+            # '2016-04-29_17h31m35s.00827'
+            structtime = time.strptime(timestr[0],
+                                       "%Y-%m-%d_%Hh%Mm%Ss")
+            fracsec = float(timestr[1]) * 1e-5
+            thetime = time.mktime(structtime) + fracsec
+        else:
+            thetime = 0
+        return thetime
+
+    @staticmethod
+    def _get_meta_data(path, section, name):
+        with tifffile.TiffFile(path) as tf:
+            meta = str(tf.pages[0].tags["61238"].value)
+
+        meta = meta.strip("'b")
+        meta = meta.replace("\\n", "\n")
+        meta = meta.replace("\\r", "")
+        root = ET.fromstring("<root>\n" + meta + "</root>")
+        for phadata in root.getchildren():
+            for cluster in phadata.getchildren():
+                sec = cluster.getchildren()[0].text
+                for child in cluster.getchildren():
+                    gchild = child.getchildren()
+                    if len(gchild) == 2:
+                        nm, val = gchild
+                        # print(sec, nm.text, val.text)
+                        if (sec.lower() == section and
+                            nm.text.lower() == name
+                            ):
+                            return val.text
+        else:
+            return None
+
+    @staticmethod
+    def verify(path):
+        """Verify that `path` has the qpimage file format
+
+        Returns `True` if the file format matches.
+        """
+        valid = False
+        if path.endswith(".tif"):
+            try:
+                tf = tifffile.TiffFile(path)
+            except:
+                pass
+            else:
+                if (len(tf) == 3 and
+                    "61243" in tf.pages[0].tags and
+                            "61242" in tf.pages[0].tags and
+                        "61238" in tf.pages[0].tags and
+                        "max_sample_value" in tf.pages[0].tags and
+                    tf.pages[0].tags["61242"].value != tf.pages[1].tags["61242"].value
+                    ):
+                    valid = True
+        return valid

setup.py

@@ -28,9 +28,9 @@
         long_description=open('README.rst').read() if exists('README.rst') else '',
         install_requires=["nrefocus >= 0.1.5",
                           "NumPy >= 1.9.0",
+                          "qpimage",
                           "scikit-image >= 0.11.0",
                           "scipy >= 0.18.0",
-                          "qpimage",
                           ],
         setup_requires=['pytest-runner'],
         tests_require=["pytest"],

tests/test_group_folder.py

@@ -9,7 +9,7 @@
 
 
 def setup_folder_single_h5(size=2):
-    path = join(dirname(abspath(__file__)), "data/bg_ramp.h5")
+    path = join(dirname(abspath(__file__)), "data/single_qpimage.h5")
     tdir = tempfile.mkdtemp(prefix="qpformat_test_")
     files = []
     for ss in range(size):

tests/test_single_hdf5_qpimage.py

@@ -7,7 +7,7 @@
 
 
 def test_load_data():
-    path = join(dirname(abspath(__file__)), "data/bg_ramp.h5")
+    path = join(dirname(abspath(__file__)), "data/single_qpimage.h5")
     ds = qpformat.load_data(path)
     assert ds.path == path
     assert ds.get_time() == 0

tests/test_single_tif_phasics.py

@@ -0,0 +1,74 @@
+"""
+This test uses a cropped version of an original phasics tif file
+(data/single_phasics.tif), which was created using this script
+(the tags are important):
+
+    import tifffile
+
+    a = tifffile.TiffFile("in.tif")
+
+    x1 = 100
+    x2 = 180
+    y1 = 140
+    y2 = 260
+
+    for ii in range(3):
+        data = a.pages[ii].asarray()[x1:x2, y1:y2]
+
+        extratags = []
+        for tag in ["61237",
+                    "61238",
+                    "61239",
+                    "61240",
+                    "61241",
+                    "61242",
+                    "61243",
+                    "max_sample_value",
+                    "min_sample_value",
+                    ]:
+            tag = a.pages[ii].tags[tag]
+            extratags.append((tag.code,
+                              tag.dtype.strip("1234567890"),
+                              1,
+                              tag.value,
+                              False))
+
+        tifffile.imsave(file="single_phasics.tif",
+                        data=data,
+                        extratags=extratags,
+                        append=True)
+"""
+
+from os.path import abspath, dirname, join
+import sys
+
+import numpy as np
+
+# Add parent directory to beginning of path variable
+sys.path.insert(0, dirname(dirname(abspath(__file__))))
+import qpformat  # noqa: E402
+
+
+def test_load_data():
+    path = join(dirname(abspath(__file__)), "data/single_phasics.tif")
+    ds = qpformat.load_data(path)
+    assert ds.path == path
+    assert "SingleTifPhasics" in ds.__repr__()
+
+
+def test_data_content():
+    path = join(dirname(abspath(__file__)), "data/single_phasics.tif")
+    ds = qpformat.load_data(path)
+    assert ds.get_time() == 1461943895.00827
+    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)
+
+
+if __name__ == "__main__":
+    # Run all tests
+    loc = locals()
+    for key in list(loc.keys()):
+        if key.startswith("test_") and hasattr(loc[key], "__call__"):
+            loc[key]()