diff --git a/flimkit_bridge/phasor.py b/flimkit_bridge/phasor.py index 806f78a..88b4e78 100644 --- a/flimkit_bridge/phasor.py +++ b/flimkit_bridge/phasor.py @@ -14,7 +14,7 @@ {'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'}, ) @@ -22,11 +22,13 @@ 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']) diff --git a/tests/test_phasor.py b/tests/test_phasor.py index 788b87d..9a73f3c 100644 --- a/tests/test_phasor.py +++ b/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)