Add support for Imspector FLIM TIFF files

phasorpy · cgohlke · Dec 16, 2024 · Dec 15, 2024 · Dec 15, 2024 · Dec 15, 2024
commit 02bcaf5f1c30f4fd130b0170ea1518490302256e
diff --git a/src/phasorpy/io.py b/src/phasorpy/io.py
@@ -16,6 +16,7 @@
 - read time-resolved and hyperspectral image data and metadata (as relevant
   to phasor analysis) from many file formats used in bio-imaging:
 
+  - :py:func:`read_imspector_tiff` - Imspector FLIM TIFF
   - :py:func:`read_lsm` - Zeiss LSM
   - :py:func:`read_ifli` - ISS IFLI
   - :py:func:`read_sdt` - Becker & Hickl SDT
@@ -123,6 +124,7 @@
     'read_fbd',
     'read_flif',
     'read_ifli',
+    'read_imspector_tiff',
     # 'read_lif',
     'read_lsm',
     # 'read_nd2',
@@ -887,6 +889,148 @@ def read_lsm(
     return DataArray(data, **metadata)
 
 
+def read_imspector_tiff(
+    filename: str | PathLike[Any],
+    /,
+) -> DataArray:
+    """Return FLIM image stack and metadata from Imspector TIFF file.
+
+    Imspector FLIM TIFF files contain TCSPC image stacks and metadata.
+
+    Parameters
+    ----------
+    filename : str or Path
+        Name of OME-TIFF file to read.
+
+    Returns
+    -------
+    xarray.DataArray
+        TCSPC image stack.
+        Usually, a 3-to-5-dimensional array of type ``uint16``.
+
+        - ``coords['H']``: times of histogram bins.
+        - ``attrs['frequency']``: repetition frequency in MHz.
+
+    Raises
+    ------
+    tifffile.TiffFileError
+        File is not a TIFF file.
+    ValueError
+        File is not an Imspector FLIM TIFF file.
+
+    Examples
+    --------
+    >>> data = read_imspector_tiff(fetch('Embryo.tif'))
+    >>> data.values
+    array(...)
+    >>> data.dtype
+    dtype('uint16')
+    >>> data.shape
+    (56, 512, 512)
+    >>> data.dims
+    ('H', 'Y', 'X')
+    >>> data.coords['H'].data  # dtime bins
+    array(...)
+    >>> data.attrs['frequency']  # doctest: +NUMBER
+    80.109
+
+    """
+    from xml.etree import ElementTree
+
+    import tifffile
+
+    with tifffile.TiffFile(filename) as tif:
+        tags = tif.pages.first.tags
+        omexml = tags.valueof(270, '')
+        make = tags.valueof(271, '')
+
+        if (
+            make != 'ImSpector'
+            or not omexml.startswith('<?xml version')
+            or len(tif.series) != 1
+            or not tif.is_ome
+        ):
+            raise ValueError(f'{tif.filename} is not an Imspector TIFF file')
+
+        series = tif.series[0]
+        ndim = series.ndim
+        axes = series.axes
+        shape = series.shape
+
+        if ndim < 3 or not axes.endswith('YX'):
+            raise ValueError(
+                f'{tif.filename} is not an Imspector FLIM TIFF file'
+            )
+
+        data = series.asarray()
+
+    attrs: dict[str, Any] = {}
+    coords = {}
+    physical_size = {}
+
+    root = ElementTree.fromstring(omexml)
+    ns = {
+        '': 'http://www.openmicroscopy.org/Schemas/OME/2008-02',
+        'ca': 'http://www.openmicroscopy.org/Schemas/CA/2008-02',
+    }
+
+    description = root.find('.//Description', ns)
+    if (
+        description is not None
+        and description.text
+        and description.text != 'not_specified'
+    ):
+        attrs['description'] = description.text
+
+    pixels = root.find('.//Image/Pixels', ns)
+    assert pixels is not None
+    for ax in 'TZYX':
+        attrib = 'TimeIncrement' if ax == 'T' else f'PhysicalSize{ax}'
+        if ax not in axes or attrib not in pixels.attrib:
+            continue
+        size = float(pixels.attrib[attrib])
+        physical_size[ax] = size
+        coords[ax] = numpy.linspace(
+            0.0,
+            size,
+            shape[axes.index(ax)],
+            endpoint=False,
+            dtype=numpy.float64,
+        )
+
+    axes_labels = root.find('.//ca:CustomAttributes/AxesLabels', ns)
+    if (
+        axes_labels is None
+        or 'X' not in axes_labels.attrib
+        or 'TCSPC' not in axes_labels.attrib['X']
+        or 'FirstAxis' not in axes_labels.attrib
+        or 'SecondAxis' not in axes_labels.attrib
+    ):
+        raise ValueError(f'{tif.filename} is not an Imspector FLIM TIFF file')
+
+    if axes_labels.attrib['FirstAxis'].endswith('TCSPC T'):
+        ax = axes[-3]
+        assert axes_labels.attrib['FirstAxis-Unit'] == 'ns'
+    elif axes_labels.attrib['SecondAxis'].endswith('TCSPC T') and ndim > 3:
+        ax = axes[-4]
+        assert axes_labels.attrib['SecondAxis-Unit'] == 'ns'
+    else:
+        raise ValueError(f'{tif.filename} is not an Imspector FLIM TIFF file')
+    axes = axes.replace(ax, 'H')
+    coords['H'] = coords[ax]
+    del coords[ax]
+
+    attrs['frequency'] = float(
+        1000.0 / (shape[axes.index('H')] * physical_size[ax])
+    )
+
+    metadata = _metadata(axes, shape, filename, attrs=attrs, **coords)
+
+    from xarray import DataArray
+
+    return DataArray(data, **metadata)
+
+
 def read_ifli(
     filename: str | PathLike[Any],
     /,
@@ -942,8 +1086,8 @@ def read_ifli(
     (256, 256, 4, 3)
     >>> data.dims
     ('Y', 'X', 'F', 'S')
-    >>> data.coords['F'].data
-    array([8.033...])
+    >>> data.coords['F'].data  # doctest: +NUMBER
+    array([8.033e+07, 1.607e+08, 2.41e+08, 4.017e+08])
     >>> data.coords['S'].data
     array(['mean', 'real', 'imag'], dtype='<U4')
     >>> data.attrs
@@ -1035,8 +1179,8 @@ def read_sdt(
     ('Y', 'X', 'H')
     >>> data.coords['H'].data
     array(...)
-    >>> data.attrs['frequency']
-    79...
+    >>> data.attrs['frequency']  # doctest: +NUMBER
+    79.99
 
     """
     import sdtfile
@@ -1148,7 +1292,7 @@ def read_ptu(
     ('T', 'Y', 'X', 'C', 'H')
     >>> data.coords['H'].data
     array(...)
-    >>> data.attrs['frequency']
+    >>> data.attrs['frequency']  # doctest: +NUMBER
     78.02
 
     """
@@ -1221,8 +1365,8 @@ def read_flif(
     ('H', 'Y', 'X')
     >>> data.coords['H'].data
     array(...)
-    >>> data.attrs['frequency']
-    80.65...
+    >>> data.attrs['frequency']  # doctest: +NUMBER
+    80.65
 
     """
     import lfdfiles

diff --git a/tests/test_io.py b/tests/test_io.py
@@ -25,6 +25,7 @@
     read_fbd,
     read_flif,
     read_ifli,
+    read_imspector_tiff,
     read_lsm,
     read_ptu,
     read_sdt,
@@ -129,6 +130,34 @@ def test_read_lsm_paramecium():
     )
 
 
+@pytest.mark.skipif(SKIP_FETCH, reason='fetch is disabled')
+def test_imspector_tiff():
+    """Test read Imspector FLIM TIFF file."""
+    data = read_imspector_tiff(fetch('Embryo.tif'))
+    assert data.values.sum(dtype=numpy.uint64) == 31348436
+    assert data.dtype == numpy.uint16
+    assert data.shape == (56, 512, 512)
+    assert data.dims == ('H', 'Y', 'X')
+    assert_almost_equal(
+        data.coords['H'][[0, -1]], [0.0, 0.218928482143], decimal=12
+    )
+    assert data.attrs['frequency'] == 80.10956424883184
+
+
+@pytest.mark.skipif(SKIP_PRIVATE, reason='file is private')
+def test_imspector_tiff_t():
+    """Test read Imspector FLIM TIFF file with TCSPC in T-axis."""
+    data = read_imspector_tiff(private_file('ZF-1100_noEF.tif'))
+    assert data.values.sum(dtype=numpy.uint64) == 18636271
+    assert data.dtype == numpy.uint16
+    assert data.shape == (56, 512, 512)
+    assert data.dims == ('H', 'Y', 'X')
+    assert_almost_equal(
+        data.coords['H'][[0, -1]], [0.0, 0.218928482143], decimal=12
+    )
+    assert data.attrs['frequency'] == 80.10956424883184
+
+
 @pytest.mark.skipif(SKIP_FETCH, reason='fetch is disabled')
 def test_read_sdt():
     """Test read Becker & Hickl SDT file."""

diff --git a/tutorials/api/phasorpy_multi-harmonic.py b/tutorials/api/phasorpy_multi-harmonic.py
@@ -15,10 +15,13 @@
 # Import required modules and functions:
 
 import numpy
-import tifffile  # TODO: from phasorpy.io import read_ometiff
 
 from phasorpy.datasets import fetch
-from phasorpy.io import phasor_from_ometiff, phasor_to_ometiff
+from phasorpy.io import (
+    phasor_from_ometiff,
+    phasor_to_ometiff,
+    read_imspector_tiff,
+)
 from phasorpy.phasor import (
     phasor_calibrate,
     phasor_filter,
@@ -34,9 +37,8 @@
 # Read a time-correlated single photon counting (TCSPC) histogram,
 # acquired at 80.11 MHz, from a file:
 
-
-signal = tifffile.imread(fetch('Embryo.tif'))
-frequency = 80.11  # MHz; from the XML metadata in the file
+signal = read_imspector_tiff(fetch('Embryo.tif'))
+frequency = signal.attrs['frequency']
 
 # %%
 # Calculate phasor coordinates
@@ -67,7 +69,8 @@
 # A homogeneous solution of Fluorescein with a fluorescence lifetime of 4.2 ns
 # was imaged as a reference for calibration:
 
-reference_signal = tifffile.imread(fetch('Fluorescein_Embryo.tif'))
+reference_signal = read_imspector_tiff(fetch('Fluorescein_Embryo.tif'))
+assert reference_signal.attrs['frequency'] == frequency
 
 # %%
 # Calculate phasor coordinates from the measured reference signal at

diff --git a/tutorials/phasorpy_introduction.py b/tutorials/phasorpy_introduction.py
@@ -61,11 +61,8 @@
 # array computing and plotting throughout this tutorial:
 
 import numpy
-import tifffile  # TODO: from phasorpy.io import read_ometiff
 from matplotlib import pyplot
 
-from phasorpy.datasets import fetch
-
 # %%
 # Read signal from file
 # ---------------------
@@ -83,11 +80,13 @@
 # files. For example, an Imspector TIFF file from the
 # `FLUTE <https://zenodo.org/records/8046636>`_  project containing a
 # time-correlated single photon counting (TCSPC) histogram
-# of a zebrafish embryo at day 3, acquired at 80 MHz:
+# of a zebrafish embryo at day 3, acquired at 80.11 MHz:
 
+from phasorpy.datasets import fetch
+from phasorpy.io import read_imspector_tiff
 
-signal = tifffile.imread(fetch('Embryo.tif'))
-frequency = 80.11  # MHz; from the XML metadata in the file
+signal = read_imspector_tiff(fetch('Embryo.tif'))
+frequency = signal.attrs['frequency']
 
 print(signal.shape, signal.dtype)
 
@@ -167,7 +166,8 @@
 #
 # Read the signal of the reference measurement from a file:
 
-reference_signal = tifffile.imread(fetch('Fluorescein_Embryo.tif'))
+reference_signal = read_imspector_tiff(fetch('Fluorescein_Embryo.tif'))
+assert reference_signal.attrs['frequency'] == frequency
 
 # %%
 # Calculate phasor coordinates from the measured reference signal:
@@ -461,3 +461,4 @@
 # sphinx_gallery_thumbnail_number = -7
 # mypy: allow-untyped-defs, allow-untyped-calls
 # mypy: disable-error-code="arg-type"
+# isort: skip_file