diff --git a/README.md b/README.md index d044700..5c7cbb8 100644 --- a/README.md +++ b/README.md @@ -25,6 +25,39 @@ The server writes its address and a generated token to `~/.flimkit/bridge.json`, It writes `~/.flimkit/qupath-bridge.json` alongside it, carrying the same address under the older protocol name. QuPath extensions built before the server moved out check for that name and stop pairing without it. The second file goes when those versions are retired. +## Z-stacks + +`POST /v1/zstack` runs FLIMKit's z-stack fit over a folder of slices, one file per slice named `region_z1.ptu`, `region_z2.ptu` and so on. A folder holding several regions is fitted as several stacks in one run. + +``` +POST /v1/zstack/scan {"ptu_dir": "/path/slices"} +POST /v1/zstack {"ptu_dir": "/path/slices", "params": {"n_exp": 2, "z_step_um": 2.0}} +GET /v1/zstack/defaults +POST /v1/zstack/export {"group_dir": "/path/out/RegionA", "format": "ome-tiff"} +``` + +The scan is there so a client can say what it found before offering any settings, and refuse an empty folder without starting a job. The fit itself runs as a job, so `GET /v1/jobs/{id}` reports progress by slice and `DELETE` cancels at the next one. + +Each stack is fitted as one FOV: the decay is pooled over every slice, the lifetimes are fitted once from the pooled decay and locked, and each slice then gets a per-pixel fit with only the amplitudes free. A single slice rarely carries the photons to identify two lifetimes; the stack usually does. + +Each stack comes back as an OME-Zarr store, a `(C, Z, Y, X)` volume with a channel per map: intensity in photons, `tau_mean_int` and `tau_mean_amp` in ns, an `alpha_N` per component, and the bound fraction and chi-squared maps when they were computed. `zarr` and `ome-zarr` are dependencies of this package, not an extra, so the store is always what a client gets unless it asks for something else. `ome-tiff` is that something else, same channels and same axes, for a viewer that will not read the store; `/v1/zstack/export` rewrites a finished run in either format without fitting it again, which is also how to correct a z step entered wrong. + +The store carries one resolution level on purpose. A pyramid is built by interpolating between neighbours, and a lifetime map is mostly NaN where nothing was fitted, so the coarse levels a viewer shows when zoomed out would smear that NaN over the pixels that did fit. Chunking, one chunk per slice per channel, is what keeps a large volume readable. + +FLIMKit's own per-slice `.npy` maps, `_zseries.csv` and `_zseries.json` are left where it writes them, and the reply names them. Its per-slice PNGs and detail plots are off by default here, since each one costs a second summed fit per slice; `save_plots` turns them back on. + +### When the client cannot see the disk + +The reply names a path, so a client on the same machine opens the store itself and nothing large goes over HTTP. Over a forwarded port that path means nothing locally, and for that there is: + +``` +GET /v1/zstack/volume.ome.tif?group_dir=/path/out/RegionA&z_step_um=2.0 +``` + +which builds the volume from the same slice maps and streams it as one OME-TIFF, in blocks out of a temporary file rather than through one bytes object, since a fitted stack is routinely larger than memory. The temporary file goes whether the transfer finished or not. Shape, channels, units and voxel size ride on `X-FLIMKit-Volume-*` headers, and the store on the server is left alone. + +The slices themselves are still a server-side path: `ptu_dir` is opened by FLIMKit, not uploaded. A forwarded client picking a local folder gets a 404 naming the folder, before any job starts. + ## Versions `GET /v1/status` reports `protocol_version`, `bridge_version` and `flimkit_version`. A client checks `protocol_version`, which is what governs whether the two can talk. `bridge_version` is for display: the server and each client version independently now, so a difference between them is ordinary rather than a problem. diff --git a/flimkit_bridge/dataset_routes.py b/flimkit_bridge/dataset_routes.py index 540f0c7..82fd774 100644 --- a/flimkit_bridge/dataset_routes.py +++ b/flimkit_bridge/dataset_routes.py @@ -510,6 +510,121 @@ def work(progress, cancel): 'basename': basename, 'n_tiles': n_tiles, 'params_used': params} +def zstack_defaults(state): + from flimkit_bridge import zstack + return zstack.defaults() + + +def scan_zstack(state, payload): + from flimkit_bridge import zstack + directory = (payload or {}).get('ptu_dir') or (payload or {}).get('path') + if not directory: + raise RouteError(400, 'a ptu_dir is required, the folder holding the slices') + try: + groups = zstack.scan(directory) + except FileNotFoundError as exc: + raise RouteError(404, str(exc)) + except Exception as exc: + raise RouteError(400, str(exc)) + from flimkit_bridge import volumes + return {'ptu_dir': str(directory), 'stacks': groups, + 'n_stacks': len(groups), + 'n_slices': sum(group['n_slices'] for group in groups), + 'volume_formats': list(volumes.FORMATS)} + + +def run_zstack(state, payload): + from flimkit_bridge import zstack + jobs = _jobs(state) + payload = payload or {} + directory = payload.get('ptu_dir') or payload.get('path') + if not directory: + raise RouteError(400, 'a ptu_dir is required, the folder holding the slices') + try: + params = zstack.merge_params(payload.get('params')) + except ValueError as exc: + raise RouteError(400, str(exc)) + try: + ptu_dir, output_dir, groups = zstack.resolve( + directory, payload.get('output_dir')) + except FileNotFoundError as exc: + raise RouteError(404, str(exc)) + except ValueError as exc: + raise RouteError(400, str(exc)) + args = zstack.build_args(ptu_dir, output_dir, params) + n_slices = sum(group['n_slices'] for group in groups) + + def work(progress, cancel): + progress(0, n_slices, + f'{len(groups)} z-stack(s) over {n_slices} slices') + result = zstack.run(args, progress, cancel) + if cancel.is_set(): + return None + return zstack.summarise(result, output_dir, groups, params) + + job_id = jobs.submit('zstack_fit', work) + return {'job': job_id, 'ptu_dir': str(ptu_dir), + 'output_dir': str(output_dir), + 'stacks': [{k: v for k, v in group.items() if k != 'files'} + for group in groups], + 'n_stacks': len(groups), 'n_slices': n_slices, + 'params_used': params} + + +ZSTACK_VOLUME_PATH = '/v1/zstack/volume.ome.tif' + + +def zstack_volume(state, query): + from flimkit_bridge import volumes, zstack + options = {key: values[-1] for key, values in parse_qs(query).items()} + group_dir = options.get('group_dir') + if not group_dir: + raise RouteError(400, 'group_dir is required, one stack of fitted slices') + try: + z_step = float(options.get('z_step_um', 1.0)) + side = float(options['pixel_size_um']) if 'pixel_size_um' in options else None + except ValueError: + raise RouteError(400, 'z_step_um and pixel_size_um must be numbers') + try: + return zstack.stream_volume(group_dir, label=options.get('label'), + z_step_um=z_step, pixel_size_um=side) + except FileNotFoundError as exc: + raise RouteError(404, str(exc)) + except volumes.NothingToStack as exc: + raise RouteError(409, str(exc)) + except Exception as exc: + raise RouteError(500, str(exc)) + + +def export_zstack(state, payload): + from flimkit_bridge import zstack + jobs = _jobs(state) + payload = payload or {} + group_dir = payload.get('group_dir') + if not group_dir: + raise RouteError(400, 'a group_dir is required, one stack of fitted slices') + from flimkit_bridge import volumes + wanted = payload.get('format', 'ome-zarr') + if wanted not in volumes.FORMATS: + raise RouteError(400, f'format must be one of {list(volumes.FORMATS)}') + from pathlib import Path + if not Path(group_dir).expanduser().is_dir(): + raise RouteError(404, f'no such folder: {group_dir}') + + def work(progress, cancel): + progress(0, 1, f'writing {wanted}') + found = zstack.export( + group_dir, label=payload.get('label'), + output_dir=payload.get('output_dir'), volume_format=wanted, + z_step_um=payload.get('z_step_um', 1.0), + pixel_size_um=payload.get('pixel_size_um')) + progress(1, 1, 'written') + return {'products': [found]} + + job_id = jobs.submit('zstack_export', work) + return {'job': job_id, 'group_dir': str(group_dir), 'format': wanted} + + def fit_pixels(state, ident, payload): from flimkit_bridge import fitting registry = _registry(state) diff --git a/flimkit_bridge/server.py b/flimkit_bridge/server.py index a47e614..a8a13ec 100644 --- a/flimkit_bridge/server.py +++ b/flimkit_bridge/server.py @@ -108,6 +108,12 @@ def do_GET(self): return self._route(lambda r: r.pipeline_defaults(state)) return + if self.path == '/v1/zstack/defaults': + if not self._authorized(): + self.send_error(401) + return + self._route(lambda r: r.zstack_defaults(state)) + return if self.path == '/v1/phasor/settings': if not self._authorized(): self.send_error(401) @@ -199,6 +205,24 @@ def do_POST(self): return self._route(lambda routes: routes.run_pipeline(state, self._read_json())) return + if self.path == '/v1/zstack/scan': + if not self._authorized(): + self.send_error(401) + return + self._route(lambda routes: routes.scan_zstack(state, self._read_json())) + return + if self.path == '/v1/zstack/export': + if not self._authorized(): + self.send_error(401) + return + self._route(lambda routes: routes.export_zstack(state, self._read_json())) + return + if self.path == '/v1/zstack': + if not self._authorized(): + self.send_error(401) + return + self._route(lambda routes: routes.run_zstack(state, self._read_json())) + return from flimkit_bridge import dataset_routes as _routes stats_match = _routes.PLANE_STATS_RE.match(self.path) if stats_match: @@ -281,6 +305,12 @@ def do_DELETE(self): def _dataset_get(self): from flimkit_bridge import dataset_routes as routes parsed = urlparse(self.path) + if parsed.path == routes.ZSTACK_VOLUME_PATH: + if not self._authorized(): + self.send_error(401) + return True + self._send_volume(parsed.query) + return True if parsed.path == '/v1/datasets': if not self._authorized(): self.send_error(401) @@ -335,6 +365,41 @@ def _dataset_get(self): return True return False + def _send_volume(self, query): + import shutil + + from flimkit_bridge import dataset_routes as routes + try: + built = routes.zstack_volume(state, query) + except routes.RouteError as problem: + self.send_error(problem.status, str(problem)) + return + except Exception as exc: + self.send_error(500, str(exc)) + return + # A fitted stack can be gigabytes, so it goes out of a file in + # blocks rather than through one bytes object in memory. + try: + size = os.path.getsize(built['file']) + self.send_response(200) + self.send_header('Content-Type', 'image/tiff') + self.send_header('Content-Length', str(size)) + self.send_header('X-FLIMKit-Volume-Label', built['label']) + self.send_header('X-FLIMKit-Volume-Axes', 'ZCYX') + self.send_header('X-FLIMKit-Volume-Channels', + ','.join(built['channels'])) + self.send_header('X-FLIMKit-Volume-Units', + ','.join(built['units'])) + self.send_header('X-FLIMKit-Volume-Shape', + ','.join(str(n) for n in built['shape'])) + self.send_header('X-FLIMKit-Voxel-Size-Um', + ','.join(str(v) for v in built['voxel_size_um'])) + self.end_headers() + with open(built['file'], 'rb') as source: + shutil.copyfileobj(source, self.wfile, 1024 * 1024) + finally: + shutil.rmtree(built['holding'], ignore_errors=True) + def _send_plane(self, ident, name, query): from flimkit_bridge import dataset_routes as routes try: diff --git a/flimkit_bridge/volumes.py b/flimkit_bridge/volumes.py new file mode 100644 index 0000000..b555c2a --- /dev/null +++ b/flimkit_bridge/volumes.py @@ -0,0 +1,308 @@ +import re +from pathlib import Path + +import numpy as np + +SLICE_RE = re.compile(r'^z(\d+)$') + +CHANNELS = ( + ('intensity', 'photons'), + ('tau_mean_int', 'ns'), + ('tau_mean_amp', 'ns'), + ('tau_mean', 'ns'), + ('alpha_1', ''), + ('alpha_2', ''), + ('alpha_3', ''), + ('bound_fraction', ''), + ('chi2_r', ''), + ('calibrated_chi2_r', ''), +) + + +class NothingToStack(Exception): + pass + + +def slice_dirs(group_dir): + found = [] + for entry in sorted(Path(group_dir).iterdir()): + if not entry.is_dir(): + continue + matched = SLICE_RE.match(entry.name) + if matched: + found.append((int(matched.group(1)), entry)) + return [entry for _z, entry in sorted(found)] + + +def z_values(group_dir): + found = [] + for entry in Path(group_dir).iterdir(): + matched = SLICE_RE.match(entry.name) if entry.is_dir() else None + if matched: + found.append(int(matched.group(1))) + return sorted(found) + + +def _stacked(group_dir, prefix, name, dirs): + whole = Path(group_dir) / f'{prefix}_{name}_stack.npy' + if whole.exists(): + return np.load(str(whole)) + frames = [] + for entry in dirs: + plane = entry / f'{name}.npy' + if not plane.exists(): + return None + frames.append(np.load(str(plane))) + if not frames: + return None + shapes = {frame.shape for frame in frames} + if len(shapes) != 1: + raise NothingToStack( + f'the {name} maps under {group_dir} are not all one shape: ' + f'{sorted(shapes)}') + return np.stack(frames) + + +def read_group(group_dir, prefix): + dirs = slice_dirs(group_dir) + if not dirs: + raise NothingToStack(f'{group_dir} holds no z0000-style slice folders') + names = [] + units = [] + planes = [] + for name, unit in CHANNELS: + found = _stacked(group_dir, prefix, name, dirs) + if found is None: + continue + found = np.asarray(found, dtype=np.float32) + if found.ndim != 3: + continue + names.append(name) + units.append(unit) + planes.append(found) + if not planes: + raise NothingToStack(f'{group_dir} holds no per-slice maps to stack') + shapes = {plane.shape for plane in planes} + if len(shapes) != 1: + raise NothingToStack( + f'the slice maps in {group_dir} are not all one shape: {sorted(shapes)}') + return np.stack(planes), names, units + + +def _omero(names, units): + channels = [] + for name, unit in zip(names, units): + channels.append({ + 'label': f'{name} ({unit})' if unit else name, + 'active': True, + 'color': 'FFFFFF', + 'window': {'start': 0.0, 'end': 1.0, 'min': 0.0, 'max': 1.0}, + }) + return {'channels': channels, 'rdefs': {'model': 'greyscale'}} + + +ZARR_FORMAT = 2 + +NGFF_VERSION = '0.4' + + +def _compressor(): + """A codec the readers can actually decompress. + + zarr-python 3 defaults to zstd, and QuPath's jzarr refuses it outright: + "Compressor id:'zstd' not supported". zlib is understood by every Zarr v2 + reader and, on lifetime maps that are largely NaN, compresses better than + blosc anyway (4.3 MB against 6.7 MB on a 16.8 MB volume). + """ + import numcodecs + + return numcodecs.Zlib(level=5) + + +def _escape(text): + from xml.sax.saxutils import escape + + return escape(str(text), {'"': '"'}) + + +def _ome_xml(name, shape, names, units, voxel_size_um): + """The OME-XML a bioformats2raw store carries beside its arrays. + + The NGFF metadata alone gets the pixels across but nothing else: QuPath + reads channel names and physical pixel sizes out of this file, and + without it every channel arrives called "Channel 1" with no calibration. + """ + t, c, z, y, x = shape + z_um, y_um, x_um = voxel_size_um + channels = ''.join( + '' + ''.format(i=i, name=_escape(label)) + for i, label in enumerate(_labels(names, units))) + return ( + '' + '' + f'' + f'' + f'{channels}' + '') + + +def _labels(names, units): + labelled = [] + for i, name in enumerate(names): + unit = units[i] if i < len(units) else '' + labelled.append(f'{name} ({unit})' if unit else name) + return labelled + + +def _multiscales(name, voxel_size_um): + z_um, y_um, x_um = voxel_size_um + return [{ + 'version': NGFF_VERSION, + 'name': name, + 'axes': [ + {'name': 't', 'type': 'time'}, + {'name': 'c', 'type': 'channel'}, + {'name': 'z', 'type': 'space', 'unit': 'micrometer'}, + {'name': 'y', 'type': 'space', 'unit': 'micrometer'}, + {'name': 'x', 'type': 'space', 'unit': 'micrometer'}, + ], + # One resolution level, deliberately. A pyramid is built by + # interpolating between neighbours, and a lifetime map is mostly NaN + # where nothing was fitted, so the coarse levels a viewer shows when + # zoomed out would smear that NaN over the pixels that did fit. + # Chunking is what keeps a large volume readable, not a pyramid. + 'datasets': [{ + 'path': '0', + 'coordinateTransformations': [ + {'type': 'scale', 'scale': [1.0, 1.0, z_um, y_um, x_um]}], + }], + }] + + +def _write_json(path, payload): + import json + + path.write_text(json.dumps(payload, indent=2)) + + +def write_ome_zarr(path, volume, names, units, voxel_size_um): + """Write one fitted stack as an OME-Zarr in the bioformats2raw layout. + + The layout is not a preference: QuPath reads OME-Zarr through + Bio-Formats, which wants a root marked `bioformats2raw.layout`, the image + under a numbered series group, and an OME/METADATA.ome.xml beside it. A + plain NGFF store at the root opens with no channel names and no pixel + calibration, and older ome-zarr defaults would not open at all. + """ + import shutil + + import zarr + + path = Path(path) + # A second run over the same output folder would otherwise read back a + # mixture of the two writes. + if path.is_dir(): + shutil.rmtree(path) + elif path.exists(): + path.unlink() + + z_um, y_um, x_um = [float(v) for v in voxel_size_um] + volume = np.ascontiguousarray(volume, dtype=np.float32) + # Bio-Formats expects the five axes it knows; a fit has no time axis, so + # it gets a single timepoint rather than a different number of axes. + five = volume[np.newaxis] + _, channels, depth, height, width = five.shape + + path.mkdir(parents=True) + _write_json(path / '.zgroup', {'zarr_format': ZARR_FORMAT}) + _write_json(path / '.zattrs', {'bioformats2raw.layout': 3}) + + series = path / '0' + series.mkdir() + _write_json(series / '.zgroup', {'zarr_format': ZARR_FORMAT}) + _write_json(series / '.zattrs', { + 'multiscales': _multiscales(path.name, (z_um, y_um, x_um)), + 'omero': _omero(names, units), + 'flimkit': {'channels': list(names), 'units': list(units), + 'voxel_size_um': [z_um, y_um, x_um]}, + }) + array = _open_array(zarr, series / '0', five.shape, + (1, 1, 1, min(256, height), min(256, width))) + array[:] = five + + ome = path / 'OME' + ome.mkdir() + _write_json(ome / '.zgroup', {'zarr_format': ZARR_FORMAT}) + _write_json(ome / '.zattrs', {'series': ['0']}) + (ome / 'METADATA.ome.xml').write_text(_ome_xml( + path.name, five.shape, names, units, (z_um, y_um, x_um)), + encoding='utf-8') + return str(path) + + +def _open_array(zarr, where, shape, chunks): + options = dict(mode='w', shape=shape, chunks=chunks, dtype='float32', + compressor=_compressor()) + try: + return zarr.open_array(str(where), zarr_format=ZARR_FORMAT, **options) + except TypeError: + # zarr 2 has no zarr_format argument and writes v2 regardless. + return zarr.open_array(str(where), **options) + + +def write_ome_tiff(path, volume, names, units, voxel_size_um): + import tifffile + + path = Path(path) + z_um, y_um, x_um = [float(v) for v in voxel_size_um] + labelled = [f'{name} ({unit})' if unit else name + for name, unit in zip(names, units)] + tifffile.imwrite( + str(path), + np.ascontiguousarray(np.moveaxis(volume, 0, 1), dtype=np.float32), + ome=True, photometric='minisblack', + resolution=(1.0 / x_um, 1.0 / y_um) if x_um and y_um else None, + metadata={'axes': 'ZCYX', 'Channel': {'Name': labelled}, + 'PhysicalSizeX': x_um, 'PhysicalSizeXUnit': 'µm', + 'PhysicalSizeY': y_um, 'PhysicalSizeYUnit': 'µm', + 'PhysicalSizeZ': z_um, 'PhysicalSizeZUnit': 'µm'}) + return str(path) + + +def read_metadata(path): + """Read back what write_ome_zarr recorded, from the series group where + the bioformats2raw layout keeps it.""" + import json + + series = Path(path) / '0' / '.zattrs' + if not series.exists(): + series = Path(path) / '.zattrs' + attrs = json.loads(series.read_text()) if series.exists() else {} + return { + 'flimkit': attrs.get('flimkit') or {}, + 'omero': attrs.get('omero') or {}, + 'multiscales': attrs.get('multiscales') or [], + } + + +FORMATS = ('ome-zarr', 'ome-tiff') + + +def save_volume(output_dir, name, volume, names, units, voxel_size_um, + prefer='ome-zarr'): + if prefer not in FORMATS: + raise ValueError(f'format must be one of {list(FORMATS)}, got {prefer!r}') + output_dir = Path(output_dir) + if prefer == 'ome-zarr': + target = output_dir / f'{name}.ome.zarr' + return write_ome_zarr(target, volume, names, units, voxel_size_um), 'ome-zarr' + target = output_dir / f'{name}.ome.tif' + return write_ome_tiff(target, volume, names, units, voxel_size_um), 'ome-tiff' diff --git a/flimkit_bridge/zstack.py b/flimkit_bridge/zstack.py new file mode 100644 index 0000000..7a22a46 --- /dev/null +++ b/flimkit_bridge/zstack.py @@ -0,0 +1,388 @@ +import argparse +from pathlib import Path + +SCHEMA = ( + {'key': 'n_exp', 'label': 'Exponentials', 'type': 'int', 'min': 1, 'max': 3, + 'default': 2}, + {'key': 'tau_min_ns', 'label': 'Minimum tau (ns)', 'type': 'float', + 'min': 0.01, 'max': 100.0, 'default': 0.2}, + {'key': 'tau_max_ns', 'label': 'Maximum tau (ns)', 'type': 'float', + 'min': 0.05, 'max': 200.0, 'default': 6.0}, + {'key': 'min_photons', 'label': 'Minimum photons per pixel', 'type': 'int', + 'min': 0, 'max': 100000, 'default': 5}, + {'key': 'irf_strategy', 'label': 'Instrument response', 'type': 'choice', + 'choices': ('machine_irf', 'machine_irf_sigma_full', 'machine_irf_sigma_half', + 'gaussian', 'parametric', 'raw'), + 'default': 'machine_irf'}, + {'key': 'irf_path', 'label': 'IRF file', 'type': 'path', 'default': ''}, + {'key': 'z_step_um', 'label': 'Z step (um)', 'type': 'float', + 'min': 0.001, 'max': 1000.0, 'default': 1.0}, + {'key': 'bound_fraction', 'label': 'Bound fraction', 'type': 'bool', + 'default': False}, + {'key': 'correct_pileup', 'label': 'Correct pile-up', 'type': 'bool', + 'default': False}, + {'key': 'volume_format', 'label': 'Volume format', 'type': 'choice', + 'choices': ('ome-zarr', 'ome-tiff'), 'default': 'ome-zarr'}, + {'key': 'channel', 'label': 'Channel', 'type': 'int', 'min': 0, 'max': 16, + 'default': None, 'advanced': True}, + {'key': 'ref_tau1_ns', 'label': 'Reference tau 1 (ns)', 'type': 'float', + 'min': 0.0, 'max': 200.0, 'default': None, 'advanced': True}, + {'key': 'ref_tau2_ns', 'label': 'Reference tau 2 (ns)', 'type': 'float', + 'min': 0.0, 'max': 200.0, 'default': None, 'advanced': True}, + {'key': 'ref_tau3_ns', 'label': 'Reference tau 3 (ns)', 'type': 'float', + 'min': 0.0, 'max': 200.0, 'default': None, 'advanced': True}, + {'key': 'save_plots', 'label': 'Save per-slice plots', 'type': 'bool', + 'default': False, 'advanced': True}, +) + +OPTIONAL = ('channel', 'ref_tau1_ns', 'ref_tau2_ns', 'ref_tau3_ns') + + +def defaults(): + values = {entry['key']: entry['default'] for entry in SCHEMA} + schema = [] + for entry in SCHEMA: + described = {'key': entry['key'], 'label': entry['label'], + 'type': entry['type'], + 'advanced': bool(entry.get('advanced', False)), + 'applies_to': ['zstack']} + for optional in ('min', 'max'): + if optional in entry: + described[optional] = entry[optional] + if entry['type'] == 'choice': + described['choices'] = list(entry['choices']) + schema.append(described) + return {'values': values, 'schema': schema} + + +def _match_choice(key, value, choices): + if value in choices: + return value + for choice in choices: + if str(choice) == str(value): + return choice + raise ValueError(f'{key} must be one of {list(choices)}, got {value!r}') + + +def merge_params(supplied): + known = {entry['key']: entry for entry in SCHEMA} + merged = defaults()['values'] + for key, value in (supplied or {}).items(): + if key not in known: + raise ValueError(f'unknown z-stack parameter: {key}') + merged[key] = value + for key, entry in known.items(): + value = merged[key] + if value is None and key in OPTIONAL: + continue + if entry['type'] == 'int': + value = int(value) + elif entry['type'] == 'float': + value = float(value) + elif entry['type'] == 'bool': + value = bool(value) + elif entry['type'] == 'path': + value = str(value or '') + elif entry['type'] == 'choice': + value = _match_choice(key, value, entry['choices']) + if 'min' in entry and value < entry['min']: + raise ValueError(f'{key} must be at least {entry["min"]}, got {value}') + if 'max' in entry and value > entry['max']: + raise ValueError(f'{key} must be at most {entry["max"]}, got {value}') + merged[key] = value + if merged['tau_min_ns'] >= merged['tau_max_ns']: + raise ValueError('tau_min_ns must be below tau_max_ns') + supplied_taus = [merged[f'ref_tau{i}_ns'] for i in range(1, merged['n_exp'] + 1)] + if any(t is not None for t in supplied_taus) and not all( + t is not None for t in supplied_taus): + raise ValueError( + f'give a reference tau for all {merged["n_exp"]} components or none') + return merged + + +def scan(ptu_dir): + from flimkit.utils.batch_fit import group_zstack_files, zstack_group_label + + directory = Path(ptu_dir).expanduser() + if not directory.is_dir(): + raise FileNotFoundError(f'no such folder: {directory}') + groups = group_zstack_files(directory) + found = [] + for (region, t, s), zslices in sorted(groups.items()): + planes = sorted(zslices) + found.append({ + 'label': zstack_group_label(region, t, s), + 'region': region, + 't': t, + 's': s, + 'n_slices': len(planes), + 'z_first': planes[0], + 'z_last': planes[-1], + 'files': [str(zslices[z]) for z in planes], + }) + return found + + +def resolve(ptu_dir, output_dir=None): + directory = Path(ptu_dir).expanduser() + groups = scan(directory) + if not groups: + raise ValueError( + f'no z-stack PTU files in {directory}; FLIMKit expects ' + 'region_zN.ptu, optionally region_tN_sN_zN.ptu') + if output_dir: + output = Path(output_dir).expanduser() + else: + output = directory / f'{directory.name.replace(" ", "_")}_flimkit_zstack' + return directory, output, groups + + +def build_args(ptu_dir, output_dir, params): + from flimkit.configs import ( + MACHINE_IRF_DEFAULT_PATH, Optimizer, lm_restarts, de_population, + de_maxiter, n_workers, IRF_BINS, IRF_FIT_WIDTH, IRF_FWHM) + + args = argparse.Namespace() + args.ptu_dir = str(ptu_dir) + args.output_dir = str(output_dir) + args.nexp = int(params['n_exp']) + args.tau_min = float(params['tau_min_ns']) + args.tau_max = float(params['tau_max_ns']) + args.min_photons = int(params['min_photons']) + args.estimate_irf = params['irf_strategy'] + args.machine_irf = params['irf_path'] or str(MACHINE_IRF_DEFAULT_PATH) + args.correct_pileup = bool(params['correct_pileup']) + args.channel = params['channel'] + args.bound_fraction = bool(params['bound_fraction']) + args.ref_tau1 = params['ref_tau1_ns'] + args.ref_tau2 = params['ref_tau2_ns'] + args.ref_tau3 = params['ref_tau3_ns'] + args.pileup_in_model = False + args.bg_in_model = False + args.fit_start_ns = None + args.fit_end_ns = None + args.exclude_ns = None + args.cost_function = 'poisson' + args.optimizer = Optimizer + args.restarts = lm_restarts + args.de_population = de_population + args.de_maxiter = de_maxiter + args.workers = n_workers + args.no_polish = False + args.irf_bins = IRF_BINS + args.irf_fit_width = IRF_FIT_WIDTH + args.irf_fwhm = IRF_FWHM + args.save_stack = True + args.save_npy = True + args.no_plots = not bool(params['save_plots']) + args.save_lifetime = bool(params['save_plots']) + args.save_rgb = bool(params['save_plots']) + args.save_intensity = bool(params['save_plots']) + args.save_ind = False + return args + + +def run(args, progress, cancel): + from flimkit.FLIM.assemble import Cancelled + from flimkit.interactive import _run_zstack_fit + try: + return _run_zstack_fit(args, progress_callback=progress, cancel_event=cancel) + except Cancelled: + cancel.set() + return None + + +def pixel_size_um(files): + from flimkit.formats import FLIMFile + + for path in files: + try: + opened = FLIMFile(str(path), verbose=False) + tags = getattr(opened, 'tags', {}) or {} + found = tags.get('ImgHdr_PixResol') or tags.get('ImgHdr_PixRes') + if found: + return float(found) + except Exception: + continue + return None + + +def stream_volume(group_dir, label=None, z_step_um=1.0, pixel_size_um=None): + """Write one stack to a temporary OME-TIFF and describe it. + + This is the route out for a client that cannot see the bridge's disk, an + SSH-forwarded port being the usual reason. The caller streams the file and + is responsible for deleting it. + """ + import tempfile + + from flimkit_bridge import volumes + + group_dir = Path(group_dir).expanduser() + if not group_dir.is_dir(): + raise FileNotFoundError(f'no such folder: {group_dir}') + label = label or group_dir.name + volume, names, units = volumes.read_group(group_dir, label) + side = float(pixel_size_um or 1.0) + holding = Path(tempfile.mkdtemp(prefix='flimkit-zstack-')) + written = volumes.write_ome_tiff( + holding / f'{label}.ome.tif', volume, names, units, + (float(z_step_um), side, side)) + return { + 'file': written, + 'holding': str(holding), + 'label': label, + 'channels': names, + 'units': units, + 'shape': [int(n) for n in volume.shape], + 'n_z': int(volume.shape[1]), + 'voxel_size_um': [float(z_step_um), side, side], + } + + +def export(group_dir, label=None, output_dir=None, volume_format='ome-zarr', + z_step_um=1.0, pixel_size_um=None): + from flimkit_bridge import volumes + + group_dir = Path(group_dir).expanduser() + if not group_dir.is_dir(): + raise FileNotFoundError(f'no such folder: {group_dir}') + label = label or group_dir.name + output = Path(output_dir).expanduser() if output_dir else group_dir.parent + output.mkdir(parents=True, exist_ok=True) + volume, names, units = volumes.read_group(group_dir, label) + side = float(pixel_size_um or 1.0) + written, kind = volumes.save_volume( + output, label, volume, names, units, + (float(z_step_um), side, side), prefer=volume_format) + return { + 'file': written, + 'group_dir': str(group_dir), + 'image_id': label, + 'unit': 'ns', + 'format': kind, + 'axes': 'CZYX', + 'channels': names, + 'units': units, + 'n_z': int(volume.shape[1]), + 'voxel_size_um': [float(z_step_um), side, side], + 'z_series_csv': _existing(group_dir / f'{label}_zseries.csv'), + 'taus_ns': [], + } + + +def summarise(result, output_dir, groups, params): + from flimkit_bridge import volumes + + output_dir = Path(output_dir) + result = result or {} + by_label = {group['label']: group for group in groups} + lateral = None + stacks = [] + products = [] + for label in sorted(result): + found = result[label] or {} + group_dir = Path(found.get('group_dir') or (output_dir / label)) + described = { + 'label': label, + 'group_dir': str(group_dir), + 'n_slices': found.get('n_slices'), + 'taus_ns': [float(t) for t in (found.get('taus_ns') or [])], + 'z_series_csv': _existing(group_dir / f'{label}_zseries.csv'), + 'z_series_json': _existing(group_dir / f'{label}_zseries.json'), + 'reference_fit': _existing(group_dir / f'{label}_reference_fit.json'), + 'pooled': _pooled(group_dir, label), + 'z_series': _rows(found.get('z_series')), + } + if lateral is None: + lateral = pixel_size_um(by_label.get(label, {}).get('files') or []) + try: + volume, names, units = volumes.read_group(group_dir, label) + except Exception as problem: + described['error'] = f'could not stack the slice maps: {problem}' + stacks.append(described) + continue + side = float(lateral or 1.0) + try: + written, kind = volumes.save_volume( + output_dir, label, volume, names, units, + (float(params['z_step_um']), side, side), + prefer=params['volume_format']) + except Exception as problem: + described['error'] = f'could not write the volume: {problem}' + stacks.append(described) + continue + described.update({ + 'volume': written, + 'format': kind, + 'channels': names, + 'units': units, + 'shape': [int(n) for n in volume.shape], + 'voxel_size_um': [float(params['z_step_um']), side, side], + }) + products.append({ + 'file': written, + 'group_dir': str(group_dir), + 'image_id': label, + 'unit': 'ns', + 'format': kind, + 'axes': 'CZYX', + 'channels': names, + 'units': units, + 'n_z': int(volume.shape[1]), + 'voxel_size_um': [float(params['z_step_um']), side, side], + 'z_series_csv': described['z_series_csv'], + 'taus_ns': described['taus_ns'], + 'pooled': described['pooled'], + }) + stacks.append(described) + return { + 'output_dir': str(output_dir), + 'pixel_size_um': lateral, + 'stacks': stacks, + 'products': products, + } + + +def _existing(path): + return str(path) if Path(path).exists() else None + + +POOLED_SKIP = ('z_slices',) + + +def _pooled(group_dir, label): + """The fit of the decay pooled over the whole stack, flattened for a + client that would rather not go and read the file.""" + import json + + found = Path(group_dir) / f'{label}_reference_fit.json' + if not found.exists(): + return None + try: + described = json.loads(found.read_text()) + except (OSError, ValueError): + return None + flat = {} + for key, value in described.items(): + if key in POOLED_SKIP: + continue + if isinstance(value, list): + flat[key] = [v for v in value + if isinstance(v, (int, float)) and not isinstance(v, bool)] + elif isinstance(value, (int, float, str, bool)) or value is None: + flat[key] = value + return flat + + +def _rows(series): + if not series: + return [] + listed = [] + for z in sorted(series): + row = {'z': z} + for key, value in (series[z] or {}).items(): + if isinstance(value, (int, float, str, bool)) or value is None: + row[key] = value + listed.append(row) + return listed diff --git a/pyproject.toml b/pyproject.toml index 5b34888..f62d144 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = 'setuptools.build_meta' [project] name = 'flimkit-bridge' -version = '0.7.0' +version = '0.8.0' description = 'The local HTTP server FLIMKit image analysis clients talk to' readme = 'README.md' requires-python = '>=3.12' @@ -14,6 +14,8 @@ dependencies = [ 'flimkit>=0.13.2', 'numpy', 'tifffile', + 'zarr>=2.16', + 'numcodecs', ] [project.optional-dependencies] diff --git a/tests/test_zstack.py b/tests/test_zstack.py new file mode 100644 index 0000000..4551e96 --- /dev/null +++ b/tests/test_zstack.py @@ -0,0 +1,351 @@ +import os +from pathlib import Path + +import numpy as np +import pytest + +from flimkit_bridge import volumes, zstack + +STACK_DIR = os.environ.get('FLIMKIT_TEST_ZSTACK', '') + +needs_zstack = pytest.mark.skipif( + not STACK_DIR or not Path(STACK_DIR).is_dir(), + reason='set FLIMKIT_TEST_ZSTACK to a folder of region_zN.ptu files') + + +def test_defaults_describe_every_parameter(): + described = zstack.defaults() + assert set(described['values']) == {entry['key'] for entry in zstack.SCHEMA} + assert all(entry['applies_to'] == ['zstack'] for entry in described['schema']) + + +def test_unknown_parameters_are_refused(): + with pytest.raises(ValueError) as raised: + zstack.merge_params({'nexp': 2}) + assert 'unknown z-stack parameter' in str(raised.value) + + +def test_tau_bounds_are_checked(): + with pytest.raises(ValueError): + zstack.merge_params({'tau_min_ns': 6.0, 'tau_max_ns': 1.0}) + + +def test_a_half_supplied_reference_is_refused(): + with pytest.raises(ValueError) as raised: + zstack.merge_params({'n_exp': 2, 'ref_tau1_ns': 0.4}) + assert 'all 2 components' in str(raised.value) + + +def test_a_whole_reference_is_kept(): + merged = zstack.merge_params({'n_exp': 2, 'ref_tau1_ns': 0.4, + 'ref_tau2_ns': 2.9}) + assert merged['ref_tau1_ns'] == 0.4 + assert merged['ref_tau3_ns'] is None + + +def test_a_missing_folder_is_refused(): + with pytest.raises(FileNotFoundError): + zstack.scan('/nowhere/slices') + + +def test_a_folder_without_slices_is_refused(tmp_path): + (tmp_path / 'notes.txt').write_text('nothing here') + with pytest.raises(ValueError) as raised: + zstack.resolve(tmp_path) + assert 'region_zN.ptu' in str(raised.value) + + +def test_slices_group_by_their_names(tmp_path): + for z in (1, 2, 3): + (tmp_path / f'RegionA_z{z}.ptu').write_bytes(b'') + (tmp_path / 'RegionB_t1_s2_z1.ptu').write_bytes(b'') + groups = {group['label']: group for group in zstack.scan(tmp_path)} + assert groups['RegionA']['n_slices'] == 3 + assert groups['RegionA']['z_first'] == 1 and groups['RegionA']['z_last'] == 3 + assert 'RegionB_t0001_s2' in groups + + +def test_the_output_folder_sits_beside_the_slices(tmp_path): + (tmp_path / 'RegionA_z1.ptu').write_bytes(b'') + _ptu_dir, output_dir, groups = zstack.resolve(tmp_path) + assert output_dir.parent == tmp_path + assert output_dir.name.endswith('_flimkit_zstack') + assert len(groups) == 1 + + +def _fake_group(root, label, n_z=3, shape=(4, 5)): + group_dir = root / label + for z in range(n_z): + slice_dir = group_dir / f'z{z:04d}' + slice_dir.mkdir(parents=True) + for name in ('intensity', 'tau_mean_int', 'tau_mean_amp', 'alpha_1'): + np.save(str(slice_dir / f'{name}.npy'), + np.full(shape, float(z + 1), dtype=np.float32)) + return group_dir + + +def test_slice_maps_stack_into_one_volume(tmp_path): + group_dir = _fake_group(tmp_path, 'RegionA') + volume, names, units = volumes.read_group(group_dir, 'RegionA') + assert volume.shape == (4, 3, 4, 5) + assert names == ['intensity', 'tau_mean_int', 'tau_mean_amp', 'alpha_1'] + assert units == ['photons', 'ns', 'ns', ''] + assert volume[0, 2].max() == 3.0 + + +def test_a_saved_stack_is_preferred_over_the_slices(tmp_path): + group_dir = _fake_group(tmp_path, 'RegionA') + np.save(str(group_dir / 'RegionA_intensity_stack.npy'), + np.zeros((3, 4, 5), dtype=np.float32)) + volume, names, _units = volumes.read_group(group_dir, 'RegionA') + assert volume[names.index('intensity')].max() == 0.0 + + +def test_a_group_without_slices_is_refused(tmp_path): + (tmp_path / 'RegionA').mkdir() + with pytest.raises(volumes.NothingToStack): + volumes.read_group(tmp_path / 'RegionA', 'RegionA') + + +def test_ragged_slices_are_refused(tmp_path): + group_dir = _fake_group(tmp_path, 'RegionA') + np.save(str(group_dir / 'z0002' / 'intensity.npy'), + np.zeros((6, 7), dtype=np.float32)) + with pytest.raises(volumes.NothingToStack): + volumes.read_group(group_dir, 'RegionA') + + +def test_a_volume_round_trips_through_ome_tiff(tmp_path): + import tifffile + + group_dir = _fake_group(tmp_path, 'RegionA') + volume, names, units = volumes.read_group(group_dir, 'RegionA') + written, kind = volumes.save_volume( + tmp_path, 'RegionA', volume, names, units, (2.0, 0.5, 0.5), + prefer='ome-tiff') + assert kind == 'ome-tiff' + read_back = tifffile.imread(written) + assert read_back.shape == (3, 4, 4, 5) + assert np.allclose(np.moveaxis(read_back, 1, 0), volume) + + +def test_a_volume_round_trips_through_ome_zarr(tmp_path): + group_dir = _fake_group(tmp_path, 'RegionA') + volume, names, units = volumes.read_group(group_dir, 'RegionA') + written, kind = volumes.save_volume( + tmp_path, 'RegionA', volume, names, units, (2.0, 0.5, 0.5)) + assert kind == 'ome-zarr' + assert written.endswith('.ome.zarr') + + import zarr + + # (T, C, Z, Y, X): Bio-Formats expects the five axes it knows, and a fit + # has no time axis, so it gets a single timepoint. + stored = np.asarray(zarr.open_array(str(Path(written) / '0' / '0'), mode='r')) + assert stored.shape == (1,) + volume.shape + assert np.allclose(stored[0], volume) + described = volumes.read_metadata(written) + assert described['flimkit']['channels'] == names + assert described['flimkit']['voxel_size_um'] == [2.0, 0.5, 0.5] + labels = [c['label'] for c in described['omero']['channels']] + assert labels[0] == 'intensity (photons)' + assert len(described['multiscales'][0]['datasets']) == 1, \ + 'the fitted values should be the only resolution level' + scale = described['multiscales'][0]['datasets'][0][ + 'coordinateTransformations'][0]['scale'] + assert scale == [1.0, 1.0, 2.0, 0.5, 0.5] + + +def test_writing_twice_replaces_the_store(tmp_path): + group_dir = _fake_group(tmp_path, 'RegionA') + volume, names, units = volumes.read_group(group_dir, 'RegionA') + volumes.save_volume(tmp_path, 'RegionA', volume, names, units, (1.0, 1.0, 1.0)) + thinner = volume[:, :2] + written, _kind = volumes.save_volume( + tmp_path, 'RegionA', thinner, names, units, (1.0, 1.0, 1.0)) + + import zarr + + stored = zarr.open_array(str(Path(written) / '0' / '0'), mode='r') + assert stored.shape == (1,) + thinner.shape + + +def test_a_run_is_summarised_into_products(tmp_path): + group_dir = _fake_group(tmp_path, 'RegionA') + (group_dir / 'RegionA_zseries.csv').write_text('z,tau_mean_mean\n0,1.9\n') + params = zstack.merge_params({'z_step_um': 2.0, 'volume_format': 'ome-tiff'}) + summary = zstack.summarise( + {'RegionA': {'group_dir': str(group_dir), 'n_slices': 3, + 'taus_ns': [0.4, 2.9], + 'z_series': {0: {'tau_mean_mean': 1.9, 'path': '/tmp/a.ptu'}}}}, + tmp_path, [{'label': 'RegionA', 'files': []}], params) + product = summary['products'][0] + assert product['image_id'] == 'RegionA' + assert product['format'] == 'ome-tiff' + assert product['n_z'] == 3 + assert product['axes'] == 'CZYX' + assert product['voxel_size_um'][0] == 2.0 + assert Path(product['file']).exists() + stack = summary['stacks'][0] + assert stack['taus_ns'] == [0.4, 2.9] + assert stack['z_series'] == [{'z': 0, 'tau_mean_mean': 1.9, 'path': '/tmp/a.ptu'}] + assert stack['z_series_csv'].endswith('RegionA_zseries.csv') + + +def test_an_unknown_volume_format_is_refused(tmp_path): + group_dir = _fake_group(tmp_path, 'RegionA') + volume, names, units = volumes.read_group(group_dir, 'RegionA') + with pytest.raises(ValueError): + volumes.save_volume(tmp_path, 'RegionA', volume, names, units, + (1.0, 1.0, 1.0), prefer='nifti') + + +def test_a_group_that_cannot_be_stacked_is_reported_not_raised(tmp_path): + (tmp_path / 'RegionA').mkdir() + params = zstack.merge_params({}) + summary = zstack.summarise( + {'RegionA': {'group_dir': str(tmp_path / 'RegionA')}}, + tmp_path, [{'label': 'RegionA', 'files': []}], params) + assert summary['products'] == [] + assert 'could not stack' in summary['stacks'][0]['error'] + + +@needs_zstack +def test_the_fit_gets_the_arguments_it_needs(): + ptu_dir, output_dir, groups = zstack.resolve(STACK_DIR) + params = zstack.merge_params({'n_exp': 1, 'z_step_um': 1.5}) + args = zstack.build_args(ptu_dir, output_dir, params) + for name in ('ptu_dir', 'output_dir', 'nexp', 'tau_min', 'tau_max', + 'estimate_irf', 'machine_irf', 'min_photons', 'save_stack', + 'bound_fraction', 'optimizer', 'workers'): + assert hasattr(args, name), f'fit_zstack reads args.{name}' + assert args.nexp == 1 + assert groups and groups[0]['n_slices'] >= 1 + + +@needs_zstack +def test_a_real_folder_reports_its_stacks(): + groups = zstack.scan(STACK_DIR) + assert groups + assert all(group['n_slices'] >= 1 for group in groups) + + +def test_the_pooled_fit_travels_with_the_result(tmp_path): + import json + + group_dir = _fake_group(tmp_path, 'RegionA') + (group_dir / 'RegionA_reference_fit.json').write_text(json.dumps({ + 'taus_ns': [3.443, 0.568], 'nexp': 2, 'total_pooled_photons': 1731281.0, + 'estimate_irf': 'machine_irf', 'user_supplied_tau': False, + 'calibrated_chi2_pearson': 151.4, 'n_slices': 8, + 'z_slices': [1, 2, 3, 4, 5, 6, 7, 8]})) + params = zstack.merge_params({'volume_format': 'ome-tiff'}) + + summary = zstack.summarise( + {'RegionA': {'group_dir': str(group_dir)}}, + tmp_path, [{'label': 'RegionA', 'files': []}], params) + + pooled = summary['stacks'][0]['pooled'] + assert pooled['taus_ns'] == [3.443, 0.568] + assert pooled['total_pooled_photons'] == 1731281.0 + assert pooled['estimate_irf'] == 'machine_irf' + assert 'z_slices' not in pooled, 'one column per slice is not a summary' + + +def test_a_stack_without_a_pooled_fit_says_none(tmp_path): + group_dir = _fake_group(tmp_path, 'RegionA') + params = zstack.merge_params({'volume_format': 'ome-tiff'}) + + summary = zstack.summarise( + {'RegionA': {'group_dir': str(group_dir)}}, + tmp_path, [{'label': 'RegionA', 'files': []}], params) + + assert summary['stacks'][0]['pooled'] is None + + +def test_the_pooled_fit_rides_on_the_product_too(tmp_path): + import json + + group_dir = _fake_group(tmp_path, 'RegionA') + (group_dir / 'RegionA_reference_fit.json').write_text( + json.dumps({'taus_ns': [3.443, 0.568], 'nexp': 2, + 'calibrated_chi2_pearson': 151.4})) + params = zstack.merge_params({'volume_format': 'ome-tiff'}) + + summary = zstack.summarise( + {'RegionA': {'group_dir': str(group_dir)}}, + tmp_path, [{'label': 'RegionA', 'files': []}], params) + + pooled = summary['products'][0]['pooled'] + assert pooled['calibrated_chi2_pearson'] == 151.4 + assert pooled['nexp'] == 2 + + +def test_the_store_is_written_in_the_layout_readers_expect(tmp_path): + import json + + group_dir = _fake_group(tmp_path, 'RegionA') + volume, names, units = volumes.read_group(group_dir, 'RegionA') + written, _kind = volumes.save_volume( + tmp_path, 'RegionA', volume, names, units, (2.0, 0.65, 0.65)) + + store = Path(written) + # bioformats2raw on Zarr v2. QuPath reads OME-Zarr through Bio-Formats, + # which needs the layout marker, the image under a numbered series, and + # the OME-XML beside it; a plain NGFF store at the root opens with no + # channel names and no calibration, and Zarr v3 does not open at all. + assert json.loads((store / '.zattrs').read_text())['bioformats2raw.layout'] == 3 + assert not (store / 'zarr.json').exists(), 'zarr.json means Zarr v3' + assert (store / 'OME' / 'METADATA.ome.xml').exists() + assert json.loads((store / 'OME' / '.zattrs').read_text())['series'] == ['0'] + series = json.loads((store / '0' / '.zattrs').read_text()) + assert series['multiscales'][0]['version'] == '0.4' + assert [a['name'] for a in series['multiscales'][0]['axes']] == \ + ['t', 'c', 'z', 'y', 'x'] + assert json.loads((store / '0' / '0' / '.zarray').read_text())['zarr_format'] == 2 + + +def test_the_store_uses_a_codec_the_readers_have(tmp_path): + import json + + group_dir = _fake_group(tmp_path, 'RegionA') + volume, names, units = volumes.read_group(group_dir, 'RegionA') + written, _kind = volumes.save_volume( + tmp_path, 'RegionA', volume, names, units, (1.0, 1.0, 1.0)) + + codec = json.loads( + (Path(written) / '0' / '0' / '.zarray').read_text())['compressor'] + # zarr-python 3 defaults to zstd and QuPath's jzarr refuses it outright: + # "Compressor id:'zstd' not supported". + assert codec['id'] == 'zlib', codec + + +def test_the_ome_xml_carries_the_names_and_the_calibration(tmp_path): + group_dir = _fake_group(tmp_path, 'RegionA') + volume, names, units = volumes.read_group(group_dir, 'RegionA') + written, _kind = volumes.save_volume( + tmp_path, 'RegionA', volume, names, units, (2.0, 0.65, 0.65)) + + xml = (Path(written) / 'OME' / 'METADATA.ome.xml').read_text(encoding='utf-8') + # Without this file QuPath names every channel "Channel 1" and reports no + # pixel size at all; the NGFF metadata alone does not reach it. + assert 'Name="intensity (photons)"' in xml + assert 'Name="tau_mean_int (ns)"' in xml + assert 'Name="alpha_1"' in xml, 'a unitless map gains no empty brackets' + assert 'PhysicalSizeZ="2.0"' in xml + assert 'PhysicalSizeX="0.65"' in xml + assert f'SizeC="{len(names)}"' in xml and 'SizeZ="3"' in xml and 'SizeT="1"' in xml + + +def test_a_name_with_xml_in_it_cannot_break_the_metadata(tmp_path): + group_dir = _fake_group(tmp_path, 'RegionA') + volume, _names, _units = volumes.read_group(group_dir, 'RegionA') + written, _kind = volumes.save_volume( + tmp_path, 'RegionA', volume, ['a= 1 + assert payload['values']['z_step_um'] > 0 + assert all(entry['applies_to'] == ['zstack'] for entry in payload['schema']) + + +def test_the_defaults_need_the_token(served): + with pytest.raises(HTTPError) as raised: + urlopen(Request(f'{served}/v1/zstack/defaults')) + assert raised.value.code == 401 + + +def test_a_folder_is_scanned_before_anything_runs(served, tmp_path): + _slices(tmp_path) + (tmp_path / 'RegionB_t1_s1_z1.ptu').write_bytes(b'') + + payload = _call(served, '/v1/zstack/scan', 'POST', {'ptu_dir': str(tmp_path)}) + + assert payload['n_stacks'] == 2 + assert payload['n_slices'] == 4 + labels = {stack['label'] for stack in payload['stacks']} + assert 'RegionA' in labels + assert payload['volume_formats'] == ['ome-zarr', 'ome-tiff'] + + +def test_scanning_an_empty_folder_finds_nothing(served, tmp_path): + payload = _call(served, '/v1/zstack/scan', 'POST', {'ptu_dir': str(tmp_path)}) + + assert payload['n_stacks'] == 0 + assert payload['stacks'] == [] + + +def test_scanning_needs_a_folder(served): + with pytest.raises(HTTPError) as raised: + _call(served, '/v1/zstack/scan', 'POST', {}) + assert raised.value.code == 400 + + +def test_scanning_a_missing_folder_is_a_404(served, tmp_path): + with pytest.raises(HTTPError) as raised: + _call(served, '/v1/zstack/scan', 'POST', + {'ptu_dir': str(tmp_path / 'nowhere')}) + assert raised.value.code == 404 + + +def test_a_folder_without_slices_will_not_start(served, tmp_path): + with pytest.raises(HTTPError) as raised: + _call(served, '/v1/zstack', 'POST', {'ptu_dir': str(tmp_path)}) + assert raised.value.code == 400 + + +def test_bad_parameters_are_refused_before_the_job(served, tmp_path): + _slices(tmp_path) + with pytest.raises(HTTPError) as raised: + _call(served, '/v1/zstack', 'POST', + {'ptu_dir': str(tmp_path), 'params': {'tau_min_ns': 9.0}}) + assert raised.value.code == 400 + + +def test_starting_a_run_describes_what_it_found(served, tmp_path, monkeypatch): + from flimkit_bridge import zstack + + _slices(tmp_path, count=4) + monkeypatch.setattr(zstack, 'build_args', lambda *a, **k: object()) + monkeypatch.setattr(zstack, 'run', lambda *a, **k: {}) + monkeypatch.setattr(zstack, 'summarise', + lambda *a, **k: {'products': [], 'stacks': []}) + + started = _call(served, '/v1/zstack', 'POST', + {'ptu_dir': str(tmp_path), 'params': {'n_exp': 1}}) + + assert started['n_stacks'] == 1 + assert started['n_slices'] == 4 + assert started['params_used']['n_exp'] == 1 + assert started['output_dir'].endswith('_flimkit_zstack') + assert 'files' not in started['stacks'][0] + assert started['job'] + + +def test_a_chosen_output_folder_is_kept(served, tmp_path, monkeypatch): + from flimkit_bridge import zstack + + _slices(tmp_path) + monkeypatch.setattr(zstack, 'build_args', lambda *a, **k: object()) + monkeypatch.setattr(zstack, 'run', lambda *a, **k: {}) + monkeypatch.setattr(zstack, 'summarise', lambda *a, **k: {}) + + started = _call(served, '/v1/zstack', 'POST', + {'ptu_dir': str(tmp_path), + 'output_dir': str(tmp_path / 'elsewhere')}) + + assert started['output_dir'] == str(tmp_path / 'elsewhere') + + +def test_a_finished_run_hands_back_its_products(served, tmp_path, monkeypatch): + import time + + from flimkit_bridge import zstack + + _slices(tmp_path) + monkeypatch.setattr(zstack, 'build_args', lambda *a, **k: object()) + monkeypatch.setattr(zstack, 'run', lambda *a, **k: {'RegionA': {}}) + monkeypatch.setattr(zstack, 'summarise', lambda *a, **k: { + 'products': [{'file': '/tmp/RegionA.ome.zarr', 'image_id': 'RegionA', + 'format': 'ome-zarr', 'n_z': 3}], + 'stacks': [{'label': 'RegionA'}]}) + + started = _call(served, '/v1/zstack', 'POST', {'ptu_dir': str(tmp_path)}) + deadline = time.time() + 10 + while time.time() < deadline: + status = _call(served, f"/v1/jobs/{started['job']}") + if status['state'] in ('done', 'error', 'cancelled'): + break + time.sleep(0.05) + + assert status['state'] == 'done', status + result = _call(served, f"/v1/jobs/{started['job']}?result")['result'] + assert result['products'][0]['format'] == 'ome-zarr' + assert result['products'][0]['n_z'] == 3 + + +def test_a_bridge_without_a_job_registry_says_so(serve_state, tmp_path): + _slices(tmp_path) + url = serve_state(BridgeState(images={})) + with pytest.raises(HTTPError) as raised: + _call(url, '/v1/zstack', 'POST', {'ptu_dir': str(tmp_path)}) + assert raised.value.code == 503 + + +def test_a_finished_stack_can_be_re_exported(served, tmp_path): + import time + + import numpy as np + + group_dir = tmp_path / 'RegionA' + for z in range(2): + slice_dir = group_dir / f'z{z:04d}' + slice_dir.mkdir(parents=True) + for name in ('intensity', 'tau_mean_int'): + np.save(str(slice_dir / f'{name}.npy'), + np.zeros((4, 4), dtype=np.float32)) + + started = _call(served, '/v1/zstack/export', 'POST', + {'group_dir': str(group_dir), 'format': 'ome-tiff', + 'z_step_um': 3.0}) + assert started['format'] == 'ome-tiff' + deadline = time.time() + 10 + while time.time() < deadline: + status = _call(served, f"/v1/jobs/{started['job']}") + if status['state'] in ('done', 'error', 'cancelled'): + break + time.sleep(0.05) + assert status['state'] == 'done', status + product = _call(served, f"/v1/jobs/{started['job']}?result")['result']['products'][0] + assert product['file'].endswith('RegionA.ome.tif') + assert product['n_z'] == 2 + assert product['voxel_size_um'][0] == 3.0 + + +def test_re_exporting_a_missing_folder_is_a_404(served, tmp_path): + with pytest.raises(HTTPError) as raised: + _call(served, '/v1/zstack/export', 'POST', + {'group_dir': str(tmp_path / 'nowhere')}) + assert raised.value.code == 404 + + +def test_an_unknown_export_format_is_refused(served, tmp_path): + (tmp_path / 'RegionA').mkdir() + with pytest.raises(HTTPError) as raised: + _call(served, '/v1/zstack/export', 'POST', + {'group_dir': str(tmp_path / 'RegionA'), 'format': 'nifti'}) + assert raised.value.code == 400 + + +def _volume_slices(root, label='RegionA', n_z=2, shape=(4, 5)): + import numpy as np + + group_dir = root / label + for z in range(n_z): + slice_dir = group_dir / f'z{z:04d}' + slice_dir.mkdir(parents=True) + for name in ('intensity', 'tau_mean_int'): + np.save(str(slice_dir / f'{name}.npy'), + np.full(shape, float(z + 1), dtype=np.float32)) + return group_dir + + +def _fetch(url, path): + request = Request(f'{url}{path}', + headers={'Authorization': 'Bearer test-token'}) + with urlopen(request) as response: + return response.read(), dict(response.headers) + + +def test_a_volume_streams_to_a_client_that_cannot_see_the_disk(served, tmp_path): + import io + + import numpy as np + import tifffile + + group_dir = _volume_slices(tmp_path) + + body, headers = _fetch(served, '/v1/zstack/volume.ome.tif?group_dir=' + + str(group_dir) + '&z_step_um=2.5') + + assert headers['Content-Type'] == 'image/tiff' + assert int(headers['Content-Length']) == len(body) + assert headers['X-FLIMKit-Volume-Label'] == 'RegionA' + assert headers['X-FLIMKit-Volume-Axes'] == 'ZCYX' + assert headers['X-FLIMKit-Volume-Channels'] == 'intensity,tau_mean_int' + assert headers['X-FLIMKit-Volume-Units'] == 'photons,ns' + assert headers['X-FLIMKit-Volume-Shape'] == '2,2,4,5' + assert headers['X-FLIMKit-Voxel-Size-Um'].startswith('2.5,') + read_back = tifffile.imread(io.BytesIO(body)) + assert read_back.shape == (2, 2, 4, 5) + assert np.allclose(read_back[1, 0], 2.0) + + +def test_streaming_leaves_no_temporary_file_behind(served, tmp_path): + import tempfile + import time + from pathlib import Path as _Path + + group_dir = _volume_slices(tmp_path) + holding = _Path(tempfile.gettempdir()) + before = set(holding.glob('flimkit-zstack-*')) + + _fetch(served, '/v1/zstack/volume.ome.tif?group_dir=' + str(group_dir)) + + # The server drops the file once the body is out, which can land just + # after the client has read the last of it, so this waits rather than + # looking exactly once. + deadline = time.time() + 5 + while time.time() < deadline: + left = set(holding.glob('flimkit-zstack-*')) - before + if not left: + break + time.sleep(0.05) + assert set(holding.glob('flimkit-zstack-*')) == before + + +def test_streaming_needs_a_folder(served): + with pytest.raises(HTTPError) as raised: + _fetch(served, '/v1/zstack/volume.ome.tif') + assert raised.value.code == 400 + + +def test_streaming_a_missing_folder_is_a_404(served, tmp_path): + with pytest.raises(HTTPError) as raised: + _fetch(served, '/v1/zstack/volume.ome.tif?group_dir=' + + str(tmp_path / 'nowhere')) + assert raised.value.code == 404 + + +def test_streaming_a_folder_with_no_maps_is_a_409(served, tmp_path): + (tmp_path / 'RegionA').mkdir() + with pytest.raises(HTTPError) as raised: + _fetch(served, '/v1/zstack/volume.ome.tif?group_dir=' + + str(tmp_path / 'RegionA')) + assert raised.value.code == 409 + + +def test_streaming_needs_the_token(served, tmp_path): + group_dir = _volume_slices(tmp_path) + with pytest.raises(HTTPError) as raised: + urlopen(Request(f'{served}/v1/zstack/volume.ome.tif?group_dir=' + + str(group_dir))) + assert raised.value.code == 401