Calibrate a phasor against a reference measurement

flimkit_bridge/phasor.py

@@ -14,19 +14,21 @@
     {'key': 'filter_size', 'label': 'Median window (px)', 'type': 'int',
      'min': 3, 'max': 15, 'applies_to': ('phasor',), 'advanced': True,
      'default': 3},
-    {'key': 'irf', 'label': 'IRF calibration', 'type': 'choice',
+    {'key': 'irf', 'label': 'IRF calibration', 'type': 'path',
      'applies_to': ('phasor',), 'advanced': False, 'default': 'none'},
 )
 
 _NOTHING = ('', 'none', 'None', 'null')
 
 
 def _choices(key):
-    if key == 'phasor_filter':
-        from flimkit.phasor.filters import phasor_filter_methods
-        return ['none'] + list(phasor_filter_methods())
+    from flimkit.phasor.filters import phasor_filter_methods
+    return ['none'] + list(phasor_filter_methods())
+
+
+def installed_irfs():
     from flimkit_bridge import irf as irf_module
-    return ['none'] + [entry['id'] for entry in irf_module.available()]
+    return [entry['id'] for entry in irf_module.available()]
 
 
 def settings():
@@ -42,6 +44,8 @@ def settings():
                 described[optional] = entry[optional]
         if entry['type'] == 'choice':
             described['choices'] = _choices(entry['key'])
+        if entry['key'] == 'irf':
+            described['installed'] = installed_irfs()
         values[entry['key']] = entry['default']
         schema.append(described)
     return {'values': values, 'schema': schema}
@@ -283,22 +287,75 @@ def compute(path, channel=None, binning=4, options=None):
     }
     irf_path = resolve_irf(found_options['irf'])
     if irf_path:
-        found.update(_calibrate(found, handle, irf_path, stack))
+        found.update(_calibrate(found, handle, irf_path, stack, channel=channel))
     found['real'], found['imag'] = apply_filter(
         found['real'], found['imag'], found['mean'], found_options)
     return found
 
 
-def _calibrate(found, handle, irf_path, stack):
+_WORKBOOK = ('.xlsx', '.xlsm', '.xls')
+
+
+def _irf_from_reference(reference, channel=None):
+    counts = np.asarray(reference.summed_decay(channel=channel), dtype=float)
+    if counts.sum() <= 0:
+        raise ValueError(
+            'the reference measurement has no photons, so it cannot be used as '
+            'an instrument response')
+    time_ns = np.asarray(getattr(reference, 'time_ns', None), dtype=float)
+    if time_ns is None or time_ns.size != counts.size:
+        raise ValueError(
+            'the reference measurement carries no per-bin time axis, so its '
+            'decay cannot be placed against the data')
+    return time_ns, counts
+
+
+def _pqres_irf(irf_path):
+    import ptufile
+    with ptufile.PqFile(str(irf_path)) as handle:
+        x = handle.tags.get('VarOverallDecayX')
+        y = handle.tags.get('VarOverallDecayY')
+    if x is None or y is None:
+        raise ValueError(
+            f'{irf_path} carries no overall decay, so there is no instrument '
+            f'response in it. Export a result that includes the TCSPC curve, '
+            f'or point at the reference measurement itself.')
+    time_ns = np.asarray(x, dtype=float) * 1e9
+    counts = np.asarray(y, dtype=float)
+    if counts.sum() <= 0:
+        raise ValueError(f'the decay in {irf_path} has no photons')
+    return time_ns, counts
+
+
+def _reference_irf(irf_path, channel=None):
+    if str(irf_path).lower().endswith('.pqres'):
+        return _pqres_irf(irf_path)
+    from flimkit.formats import FLIMFile
+    try:
+        reference = FLIMFile(str(irf_path), verbose=False)
+    except Exception as exc:
+        raise ValueError(
+            f'FLIMKit cannot read {irf_path} as an instrument response. Use a '
+            f'machine IRF .npy, an IRF workbook, a PicoQuant .pqres result '
+            f'carrying its overall decay, or a reference measurement in a '
+            f'format FLIMKit reads such as .ptu or .sdt: {exc}')
+    return _irf_from_reference(reference, channel=channel)
+
+
+def _calibrate(found, handle, irf_path, stack, channel=None):
     from flimkit.phasor.signal import (calibrate_signal_with_irf,
                                        calibrate_signal_with_machine_irf)
     signal = _signal_array(stack, handle, found['frequency'])
-    if str(irf_path).endswith('.npy'):
+    lowered = str(irf_path).lower()
+    if lowered.endswith('.npy'):
         real_cal, imag_cal = calibrate_signal_with_machine_irf(
             signal, found['real'], found['imag'], irf_path, found['frequency'])
     else:
-        from flimkit.phasor.signal import get_phasor_irf
-        irf_time_ns, irf_counts = get_phasor_irf(irf_path)
+        if lowered.endswith(_WORKBOOK):
+            from flimkit.phasor.signal import get_phasor_irf
+            irf_time_ns, irf_counts = get_phasor_irf(irf_path)
+        else:
+            irf_time_ns, irf_counts = _reference_irf(irf_path, channel=channel)
         real_cal, imag_cal = calibrate_signal_with_irf(
             signal, found['real'], found['imag'], irf_time_ns, irf_counts,
             found['frequency'])

tests/test_phasor.py

@@ -237,15 +237,14 @@ def test_settings_lists_the_registered_filters():
     assert found['values']['phasor_filter'] == 'none'
 
 
-def test_settings_offers_every_installed_machine_irf():
+def test_settings_lists_every_installed_machine_irf():
     from flimkit_bridge import irf as irf_module
 
     found = phasor.settings()
 
     entry = next(e for e in found['schema'] if e['key'] == 'irf')
-    assert entry['choices'][0] == 'none'
     for installed in irf_module.available():
-        assert installed['id'] in entry['choices']
+        assert installed['id'] in entry['installed']
 
 
 def test_normalise_fills_the_defaults():
@@ -495,3 +494,79 @@ def test_a_polygon_reports_lifetimes_like_an_ellipse(two_populations):
     assert found[0]['n_pixels'] == 128
     assert found[0]['tau_phi_ns'] > 0
     assert found[0]['mean_g'] == pytest.approx(0.30, abs=0.01)
+
+
+def test_the_irf_setting_is_a_path_field():
+    found = phasor.settings()
+
+    entry = next(e for e in found['schema'] if e['key'] == 'irf')
+    assert entry['type'] == 'path', (
+        'a dropdown cannot point at a reference measurement on disk')
+
+
+def test_resolve_irf_still_accepts_an_installed_id():
+    from flimkit_bridge import irf as irf_module
+
+    installed = irf_module.available()
+    if not installed:
+        pytest.skip('no machine IRF installed')
+    assert phasor.resolve_irf(installed[0]['id']) == installed[0]['path']
+
+
+def test_a_reference_measurement_is_read_as_a_decay(tmp_path):
+    class Reference:
+        time_ns = np.arange(64) * 0.05
+
+        def summed_decay(self, channel=None):
+            counts = np.zeros(64)
+            counts[8] = 100.0
+            counts[9] = 40.0
+            return counts
+
+    irf_time_ns, irf_counts = phasor._irf_from_reference(Reference())
+
+    assert irf_time_ns[8] == pytest.approx(0.4)
+    assert irf_counts[8] == 100.0
+    assert irf_counts.sum() == 140.0
+
+
+def test_a_reference_with_no_photons_is_refused():
+    class Empty:
+        time_ns = np.arange(64) * 0.05
+
+        def summed_decay(self, channel=None):
+            return np.zeros(64)
+
+    with pytest.raises(ValueError, match='no photons'):
+        phasor._irf_from_reference(Empty())
+
+
+def test_an_unreadable_reference_says_what_is_supported(tmp_path):
+    bogus = tmp_path / 'reference.pqres'
+    bogus.write_bytes(b'not a flim file')
+
+    with pytest.raises(ValueError, match='pqres'):
+        phasor._reference_irf(str(bogus))
+
+
+def test_a_pqres_result_gives_its_overall_decay():
+    import os
+    sample = '/Users/as-hunt/Downloads/Picoquant/ATTO488_2_OTCSPC.pqres'
+    if not os.path.exists(sample):
+        pytest.skip('no .pqres sample on this machine')
+
+    time_ns, counts = phasor._reference_irf(sample)
+
+    assert time_ns.size == counts.size
+    assert counts.sum() > 0
+    assert time_ns[-1] > time_ns[0]
+
+
+def test_a_pqres_without_a_decay_says_so():
+    import os
+    sample = '/Users/as-hunt/Downloads/Picoquant/TCSPC_Fitting_1.pqres'
+    if not os.path.exists(sample):
+        pytest.skip('no .pqres sample on this machine')
+
+    with pytest.raises(ValueError, match='no overall decay'):
+        phasor._reference_irf(sample)