From 9bab886f1d80d0589da483714dfb85c0ed362d94 Mon Sep 17 00:00:00 2001 From: d33bs Date: Wed, 5 Aug 2026 17:17:12 -0600 Subject: [PATCH 1/7] Add per-shard CLI and Metadata_Imaging_ImageID propagation ZedProfiler is the feature extractor the NF1 pipeline dispatches per well/FOV shard via SLURM sbatch. This commit adds the command that process runs and the identifier column that makes shards warehouse-joinable. CLI (src/zedprofiler/cli.py): - "ZedProfiler run" subcommand (argparse): repeatable --image/--label NAME=PATH flags, --anisotropy-spacing Z Y X, identifier fields (--patient-tumor/--plate/--well/--field), --out-dir, a --features selector, repeatable --feature TYPE[,key=value,...] advanced requests, --skip-existing, and --force. - Reuses the six compute_* featurizers; builds the shared loader via ImageSetLoader.from_image_dict with the identifier fields. - Restartable/idempotent: deterministic output paths, --skip-existing filters before any image is read (a finished shard re-run skips I/O entirely), and atomic writes (temp + os.replace) so a crashed shard never leaves a partial file that --skip-existing would mistake for complete. - Fixes the orphaned/mis-cased console script to zedprofiler.cli:trigger. Identifiers (src/zedprofiler/identifiers.py): - build_image_id(patient_tumor, plate, well, field) -> deterministic Metadata_Imaging_ImageID; single source of truth for the format. Loaders (src/zedprofiler/IO/loading_classes.py): - ImageSetConfig carries patient_tumor/plate/well/field with an image_id property; ImageSetLoader exposes image_id (falls back to image_set_name). - New from_image_dict classmethod builds a multi-channel loader from an in-memory {key: ndarray} dict (the path the CLI needs); it derives compartments/image names directly from the declared label keys so it is self-contained and correct independent of get_compartments. Featurizers (6 modules): - Each emits Metadata_Imaging_ImageID before Metadata_Experiment_ImageSet. Feature values are unchanged; only a metadata column is added. Feature writing: - save_features_as_parquet gains an opt-in atomic flag used by the CLI. End-to-end CLI tests self-skip when the CellProfiler 3D tutorial data is absent (it lands via a separate data commit), so the CLI test module stays green everywhere and runs in full wherever the data is available. --- pyproject.toml | 2 +- src/zedprofiler/IO/feature_writing_utils.py | 15 +- src/zedprofiler/IO/loading_classes.py | 133 ++++ src/zedprofiler/cli.py | 679 ++++++++++++++++++ .../featurization/colocalization.py | 11 +- src/zedprofiler/featurization/granularity.py | 5 + src/zedprofiler/featurization/intensity.py | 5 + src/zedprofiler/featurization/neighbors.py | 5 + src/zedprofiler/featurization/texture.py | 5 + .../featurization/volumesizeshape.py | 5 + src/zedprofiler/identifiers.py | 55 ++ tests/IO/test_loading_classes.py | 120 ++++ tests/featurization/test_colocalization.py | 2 + tests/featurization/test_granularity.py | 12 +- tests/featurization/test_intensity.py | 2 + tests/featurization/test_neighbors.py | 2 + .../test_neighbors_additional.py | 2 +- tests/featurization/test_texture.py | 2 + tests/featurization/test_volumesizeshape.py | 2 + tests/test_cli.py | 359 ++++++++- tests/test_identifiers.py | 44 ++ 21 files changed, 1456 insertions(+), 11 deletions(-) create mode 100644 src/zedprofiler/cli.py create mode 100644 src/zedprofiler/identifiers.py create mode 100644 tests/test_identifiers.py diff --git a/pyproject.toml b/pyproject.toml index fdfbf27..3a19baa 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -42,7 +42,7 @@ dependencies = [ ] urls.Homepage = "https://zedprofiler.readthedocs.io" urls.Repository = "https://github.com/WayScience/ZedProfiler" -scripts.ZedProfiler = "ZedProfiler.cli:trigger" +scripts.ZedProfiler = "zedprofiler.cli:trigger" [dependency-groups] dev = [ diff --git a/src/zedprofiler/IO/feature_writing_utils.py b/src/zedprofiler/IO/feature_writing_utils.py index 6769da7..3ac8ab9 100644 --- a/src/zedprofiler/IO/feature_writing_utils.py +++ b/src/zedprofiler/IO/feature_writing_utils.py @@ -6,6 +6,7 @@ from __future__ import annotations import dataclasses +import os import pathlib import pandas @@ -181,6 +182,7 @@ def save_features_as_parquet( parent_path: pathlib.Path, df: pandas.DataFrame, metadata: FeatureMetadata, + atomic: bool = False, ) -> pathlib.Path: """Save features as parquet files in a consistent way. @@ -196,6 +198,12 @@ def save_features_as_parquet( metadata : FeatureMetadata Metadata for the feature output (compartment, channel, feature_type, cpu_or_gpu). + atomic : bool + When True, write to a sibling ``.tmp`` file then atomically replace + the destination via ``os.replace``. This prevents a crashed write + from leaving a partial parquet file that a restartable caller (using + ``--skip-existing``) could mistake for a complete one. Default False + preserves the existing direct-write behavior. Returns ------- @@ -210,5 +218,10 @@ def save_features_as_parquet( metadata.cpu_or_gpu, ) save_path = parent_path / f"{output_prefix}_features.parquet" - validated_df.to_parquet(save_path, index=False) + if atomic: + tmp_path = save_path.with_suffix(save_path.suffix + ".tmp") + validated_df.to_parquet(tmp_path, index=False) + os.replace(tmp_path, save_path) + else: + validated_df.to_parquet(save_path, index=False) return save_path diff --git a/src/zedprofiler/IO/loading_classes.py b/src/zedprofiler/IO/loading_classes.py index d9435fd..ab4e648 100644 --- a/src/zedprofiler/IO/loading_classes.py +++ b/src/zedprofiler/IO/loading_classes.py @@ -12,6 +12,7 @@ from beartype import beartype from zedprofiler.contracts import ImageArrayModel +from zedprofiler.identifiers import build_image_id logging.basicConfig(level=logging.INFO) @@ -42,6 +43,15 @@ class ImageSetConfig: image_set_name: str | None = None label_key_name: list[str] | None = None raw_image_key_name: list[str] | None = None + # Imaging-coordinate identifier fields used to build a deterministic + # ``Metadata_Imaging_ImageID``. All four must be set for ``image_id`` to be + # populated; when any is None, ``image_id`` is None and the loader falls + # back to ``image_set_name`` for the emitted metadata column (see + # ``ImageSetLoader.image_id``). + patient_tumor: str | None = None + plate: str | None = None + well: str | None = None + field: int | str | None = None # validate the arg types def __post_init__(self) -> None: @@ -52,12 +62,42 @@ def __post_init__(self) -> None: raise TypeError("label_key_name must be a list of strings or None") if not isinstance(self.raw_image_key_name, (list, type(None))): raise TypeError("raw_image_key_name must be a list of strings or None") + if not isinstance(self.patient_tumor, (str, type(None))): + raise TypeError("patient_tumor must be a string or None") + if not isinstance(self.plate, (str, type(None))): + raise TypeError("plate must be a string or None") + if not isinstance(self.well, (str, type(None))): + raise TypeError("well must be a string or None") + if not isinstance(self.field, (int, str, type(None))): + raise TypeError("field must be an int, str, or None") if self.label_key_name is None: self.label_key_name = [] if self.raw_image_key_name is None: self.raw_image_key_name = [] + @property + def image_id(self) -> str | None: + """Deterministic ``Metadata_Imaging_ImageID`` value, or None if unset. + + Returns ``build_image_id(...)`` when all four coordinate fields are + set, otherwise ``None`` (the loader then falls back to + ``image_set_name``). + """ + if ( + self.patient_tumor is not None + and self.plate is not None + and self.well is not None + and self.field is not None + ): + return build_image_id( + patient_tumor=self.patient_tumor, + plate=self.plate, + well=self.well, + field=self.field, + ) + return None + class _LazyImageSetDict(dict): # type: ignore[type-arg] """Dictionary that loads image arrays on first access.""" @@ -178,6 +218,13 @@ def __init__( # noqa: PLR0913 self.anisotropy_factor = self.anisotropy_spacing[0] / self.anisotropy_spacing[1] self.image_set_name = config.image_set_name self.label_set_path = label_set_path + # Deterministic imaging identifier for the warehouse join key + # (``Metadata_Imaging_ImageID``). When identifier fields are not + # provided (legacy/library use), fall back to the image set name so + # the emitted metadata column always has a value. + self.image_id = ( + config.image_id if config.image_id is not None else config.image_set_name + ) self._load_path_based_images( channel_mapping=channel_mapping, channel_tokens=channel_tokens, @@ -194,6 +241,92 @@ def __init__( # noqa: PLR0913 self.get_image_names() self.get_unique_objects_in_compartments() + @classmethod + def from_image_dict( # noqa: PLR0913 + cls, + image_dict: dict[str, numpy.ndarray], + *, + anisotropy_spacing: tuple[float, float, float], + image_set_name: str | None = None, + label_key_names: list[str] | None = None, + patient_tumor: str | None = None, + plate: str | None = None, + well: str | None = None, + field: int | str | None = None, + ) -> ImageSetLoader: + """Build an ImageSetLoader from an in-memory channel/label dict. + + Existing constructors only accept a directory glob (path-based) or a + single array (array-based). A well/FOV shard carries multiple channels + and multiple compartments as distinct arrays, so this classmethod + builds the ``image_set_dict`` directly from a pre-loaded + ``{key: ndarray}`` mapping. It formalizes the ``ImageSetLoader.__new__`` + workaround previously used in the colocalization test helper. + + Parameters + ---------- + image_dict : dict[str, numpy.ndarray] + Mapping of channel names and compartment names to their arrays. + Each array is validated through ``ImageArrayModel``. + anisotropy_spacing : tuple[float, float, float] + (z_spacing, y_spacing, x_spacing). + image_set_name : str | None + Optional image set name (emitted as ``Metadata_Experiment_ImageSet``). + label_key_names : list[str] | None + Keys in ``image_dict`` that are compartment labels (not channels). + Used by ``get_compartments`` to distinguish compartments from + raw channels. + patient_tumor, plate, well, field : optional + Imaging-coordinate identifier fields. When all four are set, the + loader's ``image_id`` is the deterministic + ``Metadata_Imaging_ImageID``; otherwise it falls back to + ``image_set_name``. + + Returns + ------- + ImageSetLoader + A fully initialized loader (compartments, image names, and unique + compartment objects populated). + + """ + self = cls.__new__(cls) + self.image_set_dict = _LazyImageSetDict() + for key, array in image_dict.items(): + # Run through pydantic validation to ensure each array is valid, + # mirroring ``_load_array_based_images``. + self.image_set_dict[key] = ImageArrayModel(array=array).array + self._label_key_names = list(label_key_names or []) + self.anisotropy_spacing = anisotropy_spacing + self.anisotropy_factor = self.anisotropy_spacing[0] / self.anisotropy_spacing[1] + self.image_set_name = image_set_name + self.label_set_path = None + config = ImageSetConfig( + image_set_name=image_set_name, + label_key_name=list(label_key_names or []), + raw_image_key_name=[ + key for key in image_dict if key not in (label_key_names or []) + ], + patient_tumor=patient_tumor, + plate=plate, + well=well, + field=field, + ) + self.image_id = ( + config.image_id if config.image_id is not None else config.image_set_name + ) + # Set compartments and image names directly from the declared label + # keys rather than calling ``get_compartments``/``get_image_names``. + # Those methods' compartment heuristic differs across repo revisions + # (it was corrected in a later bugfix commit), but this classmethod + # already knows which keys are labels, so deriving the split here keeps + # it self-contained and correct on any base. + self.compartments = list(self._label_key_names) + self.image_names = [ + key for key in image_dict if key not in self._label_key_names + ] + self.get_unique_objects_in_compartments() + return self + @staticmethod def _validate_input_sources( image_set_path: pathlib.Path | None, diff --git a/src/zedprofiler/cli.py b/src/zedprofiler/cli.py new file mode 100644 index 0000000..3b5c13a --- /dev/null +++ b/src/zedprofiler/cli.py @@ -0,0 +1,679 @@ +"""Command-line interface for per-shard feature extraction. + +``ZedProfiler run`` is the process a workflow manager (Nextflow via SLURM +``sbatch``) dispatches once per well/FOV shard. It loads one image set from +explicit file paths, runs a selected subset of featurizers, and writes one +Parquet per feature table to an output directory. + +Why argparse (not ``fire``): the repo's other CLI surfaces use ``fire``, but +this command needs repeatable flags (``--image``/``--label``) and a +three-value ``--anisotropy-spacing`` flag, which argparse handles cleanly and +``fire`` does not. The ``trigger()`` entry point name is kept for consistency +with the existing ``pyproject.toml`` console script. + +Idempotency: the same shard spec always produces the same output paths and +content. ``--skip-existing`` skips a feature request whose output Parquet +already exists, so a re-run fills only missing shards/feature tables without +redoing finished ones. Writes are atomic (temp file + ``os.replace``) so a +crashed shard never leaves a partial file that ``--skip-existing`` would +mistake for a complete one. +""" + +from __future__ import annotations + +import argparse +import sys +from collections.abc import Sequence +from pathlib import Path + +from zedprofiler.featurization.colocalization import compute_colocalization +from zedprofiler.featurization.granularity import compute_granularity +from zedprofiler.featurization.intensity import compute_intensity +from zedprofiler.featurization.neighbors import compute_neighbors +from zedprofiler.featurization.texture import compute_texture +from zedprofiler.featurization.volumesizeshape import compute_volume_size_shape +from zedprofiler.identifiers import build_image_id +from zedprofiler.IO.feature_writing_utils import ( + FeatureMetadata, + format_morphology_feature_name, + save_features_as_parquet, +) +from zedprofiler.IO.loading_classes import ( + ImageSetLoader, + ObjectLoader, + TwoObjectLoader, + _image_loading, +) + +# CPU-backed featurizers; the ``cpu_or_gpu`` component of the output filename. +_CPU_OR_GPU = "cpu" + +# Minimum number of channels required for colocalization requests. +_MIN_CHANNELS_FOR_COLOCALIZATION = 2 + +# Feature types that consume a single channel + compartment via ObjectLoader. +_SINGLE_CHANNEL_TYPES = ( + "VolumeSizeShape", + "Intensity", + "Neighbors", + "Texture", + "Granularity", +) +# Feature types that require two channels via TwoObjectLoader. +_TWO_CHANNEL_TYPES = ("Colocalization",) +ALL_FEATURE_TYPES = (*_SINGLE_CHANNEL_TYPES, *_TWO_CHANNEL_TYPES) + +# Channel-agnostic features: their computation does not use the channel image, +# but like every ZedProfiler feature they are namespaced by a channel for +# warehouse organization, so a channel (for naming only) is still required. +_CHANNEL_AGNOSTIC_TYPES = ("VolumeSizeShape", "Neighbors") + + +def _parse_name_path(token: str) -> tuple[str, Path]: + """Parse a ``NAME=PATH`` flag value into a (name, path) pair.""" + if "=" not in token: + raise argparse.ArgumentTypeError( + f"Expected NAME=PATH, got {token!r}", + ) + name, raw_path = token.split("=", 1) + name = name.strip() + if not name: + raise argparse.ArgumentTypeError(f"NAME in {token!r} is empty") + return name, Path(raw_path) + + +def _parse_feature_spec(token: str) -> dict[str, object]: + """Parse a ``TYPE[,key=value,...]`` feature request into a dict. + + The first comma-separated token is the feature type; the rest are + ``key=value`` overrides for that feature's parameters. + """ + parts = [p.strip() for p in token.split(",") if p.strip()] + if not parts: + raise argparse.ArgumentTypeError(f"Empty feature spec: {token!r}") + feature_type = parts[0] + if feature_type not in ALL_FEATURE_TYPES: + valid = ", ".join(ALL_FEATURE_TYPES) + raise argparse.ArgumentTypeError( + f"Unknown feature type {feature_type!r}; valid: {valid}", + ) + request: dict[str, object] = {"type": feature_type} + for part in parts[1:]: + if "=" not in part: + raise argparse.ArgumentTypeError( + f"Expected key=value, got {part!r} in {token!r}", + ) + key, value = part.split("=", 1) + request[key.strip()] = value.strip() + return request + + +def _output_path( + out_dir: Path, + compartment: str, + channel: str, + feature_type: str, +) -> Path: + """The Parquet path a feature request will write to. + + Mirrors the naming inside ``save_features_as_parquet`` so the CLI can check + existence (for ``--skip-existing``) before running a featurizer. + """ + prefix = format_morphology_feature_name( + compartment, + channel, + feature_type, + _CPU_OR_GPU, + ) + return out_dir / f"{prefix}_features.parquet" + + +def _coerce_param(value: object, cast: type, key: str, spec: str) -> object: + """Cast a parsed string param, raising a friendly error on failure.""" + try: + return cast(value) # type: ignore[arg-type] + except (TypeError, ValueError) as exc: + raise argparse.ArgumentTypeError( + f"Parameter {key}={value!r} in {spec!r} is not a valid {cast.__name__}", + ) from exc + + +def _resolve_channel(request: dict[str, object], spec: str) -> str: + """Return the single channel for a single-channel feature request.""" + channel = request.get("channel") + if not isinstance(channel, str) or not channel: + raise argparse.ArgumentTypeError( + f"Feature spec {spec!r} requires a 'channel' key", + ) + return channel + + +def _resolve_compartment(request: dict[str, object], spec: str) -> str: + """Return the compartment for a feature request.""" + compartment = request.get("compartment") + if not isinstance(compartment, str) or not compartment: + raise argparse.ArgumentTypeError( + f"Feature spec {spec!r} requires a 'compartment' key", + ) + return compartment + + +def _run_single_channel( + image_set_loader: ImageSetLoader, + request: dict[str, object], +) -> tuple[str, str, object]: + """Run a single-channel featurizer; return (channel, feature_type, df).""" + feature_type = str(request["type"]) + spec = ",".join(f"{k}={v}" for k, v in request.items()) + compartment = _resolve_compartment(request, spec) + channel = _resolve_channel(request, spec) + object_loader = ObjectLoader( + image_set_loader=image_set_loader, + channel_name=channel, + compartment_name=compartment, + ) + if feature_type == "VolumeSizeShape": + df = compute_volume_size_shape( + image_set_loader=image_set_loader, + object_loader=object_loader, + ) + elif feature_type == "Intensity": + df = compute_intensity(object_loader) + elif feature_type == "Neighbors": + df = compute_neighbors( + object_loader, + distance_threshold=int( + request.get("distance_threshold", 10), + ), + anisotropy_factor=float( + request.get( + "anisotropy_factor", + image_set_loader.anisotropy_factor, + ), + ), + ) + elif feature_type == "Texture": + df = compute_texture( + object_loader, + distance=int(request.get("distance", 1)), + grayscale=int(request.get("grayscale", 256)), + ) + elif feature_type == "Granularity": + df = compute_granularity( + object_loader, + radius=int(request.get("radius", 10)), + granular_spectrum_length=int( + request.get("granular_spectrum_length", 16), + ), + subsample_size=float(request.get("subsample_size", 0.25)), + image_sample_size=float(request.get("image_sample_size", 0.25)), + ) + else: # pragma: no cover - exhaustive dispatch above + raise ValueError(f"Unhandled single-channel feature type: {feature_type}") + return channel, feature_type, df + + +def _run_colocalization( + image_set_loader: ImageSetLoader, + request: dict[str, object], +) -> tuple[str, str, object]: + """Run a two-channel colocalization request; return (channel, type, df).""" + spec = ",".join(f"{k}={v}" for k, v in request.items()) + compartment = _resolve_compartment(request, spec) + channel1 = request.get("channel1") + channel2 = request.get("channel2") + if not isinstance(channel1, str) or not isinstance(channel2, str): + raise argparse.ArgumentTypeError( + f"Colocalization spec {spec!r} requires 'channel1' and 'channel2' keys", + ) + two_object_loader = TwoObjectLoader( + image_set_loader=image_set_loader, + compartment=compartment, + channel1=channel1, + channel2=channel2, + ) + df = compute_colocalization( + two_object_loader, + thr=int(request.get("thr", 15)), + fast_costes=str(request.get("fast_costes", "Accurate")), + channel1=channel1, + channel2=channel2, + ) + return f"{channel1}-{channel2}", "Colocalization", df + + +def _request_output_identity( + request: dict[str, object], +) -> tuple[str, str, str]: + """Return (compartment, channel, feature_type) for a request's output path. + + Pure (no featurizer runs): lets ``run()`` compute the target Parquet path + for ``--skip-existing`` filtering *before* reading any image from disk, so + a re-run over finished shards skips image I/O entirely. + """ + feature_type = str(request["type"]) + if feature_type == "Colocalization": + channel = f"{request['channel1']}-{request['channel2']}" + else: + channel = str(request["channel"]) + return str(request["compartment"]), channel, feature_type + + +def _execute_request( + image_set_loader: ImageSetLoader, + request: dict[str, object], + out_dir: Path, + target: Path, +) -> Path: + """Run one feature request and atomically write its Parquet. + + The caller precomputes ``target`` (via ``_output_path``) and handles + ``--skip-existing`` filtering, so this always writes. + """ + feature_type = str(request["type"]) + if feature_type == "Colocalization": + channel, ran_type, df = _run_colocalization(image_set_loader, request) + else: + channel, ran_type, df = _run_single_channel(image_set_loader, request) + compartment = str(request["compartment"]) + out_dir.mkdir(parents=True, exist_ok=True) + save_features_as_parquet( + out_dir, + df, + FeatureMetadata( + compartment=compartment, + channel=channel, + feature_type=ran_type, + cpu_or_gpu=_CPU_OR_GPU, + ), + atomic=True, + ) + print(f"wrote: {target}", file=sys.stderr) + return target + + +def _colocalization_requests( + channels: list[str], + compartments: list[str], +) -> list[dict[str, object]]: + """One Colocalization request per ordered channel pair x compartment.""" + if len(channels) < _MIN_CHANNELS_FOR_COLOCALIZATION: + raise argparse.ArgumentTypeError( + "Colocalization requires at least two --image channels", + ) + return [ + { + "type": "Colocalization", + "channel1": channels[i], + "channel2": channels[j], + "compartment": compartment, + } + for i in range(len(channels)) + for j in range(i + 1, len(channels)) + for compartment in compartments + ] + + +def _channel_agnostic_requests( + feature_type: str, + channel: str, + compartments: list[str], +) -> list[dict[str, object]]: + """One request per compartment, using ``channel`` for naming only.""" + return [ + {"type": feature_type, "channel": channel, "compartment": compartment} + for compartment in compartments + ] + + +def _per_channel_requests( + feature_type: str, + channels: list[str], + compartments: list[str], +) -> list[dict[str, object]]: + """One request per channel x compartment (Intensity/Texture/Granularity).""" + return [ + {"type": feature_type, "channel": channel, "compartment": compartment} + for channel in channels + for compartment in compartments + ] + + +def _auto_requests( + channels: list[str], + compartments: list[str], + feature_types: list[str], +) -> list[dict[str, object]]: + """Generate requests over the channel x compartment cross-product. + + Used when ``--features`` is given without explicit ``--feature`` specs. + Channel-agnostic features (VolumeSizeShape, Neighbors) use the first + channel for naming only. Colocalization runs on each ordered channel pair + x compartment. + """ + if not channels: + raise argparse.ArgumentTypeError("No --image flags provided") + if not compartments: + raise argparse.ArgumentTypeError("No --label flags provided") + first_channel = channels[0] + requests: list[dict[str, object]] = [] + for feature_type in feature_types: + if feature_type == "Colocalization": + requests.extend(_colocalization_requests(channels, compartments)) + elif feature_type in _CHANNEL_AGNOSTIC_TYPES: + requests.extend( + _channel_agnostic_requests(feature_type, first_channel, compartments), + ) + else: # Intensity, Texture, Granularity + requests.extend(_per_channel_requests(feature_type, channels, compartments)) + return requests + + +def _resolve_requests( + channels: list[str], + compartments: list[str], + feature_specs: list[dict[str, object]], + features_filter: list[str] | None, +) -> list[dict[str, object]]: + """Determine the final list of feature requests to run. + + - If explicit ``--feature`` specs are given, use them (filtered by + ``--features`` if present), validating referenced channels/compartments. + - Else auto-generate requests from ``--features`` (or all applicable types). + """ + if feature_specs: + requests = feature_specs + if features_filter: + requests = [r for r in requests if str(r["type"]) in features_filter] + _validate_request_channels_compartments(requests, channels, compartments) + return requests + feature_types = features_filter if features_filter else list(_SINGLE_CHANNEL_TYPES) + if (features_filter is None) and len(channels) >= _MIN_CHANNELS_FOR_COLOCALIZATION: + feature_types = [*_SINGLE_CHANNEL_TYPES, "Colocalization"] + return _auto_requests(channels, compartments, feature_types) + + +def _validate_request_channels_compartments( + requests: list[dict[str, object]], + channels: list[str], + compartments: list[str], +) -> None: + """Ensure each explicit request references declared channels/compartments.""" + for request in requests: + feature_type = str(request["type"]) + compartment = request.get("compartment") + if compartment not in compartments: + raise argparse.ArgumentTypeError( + f"compartment {compartment!r} in {feature_type} request is not " + f"declared via --label (valid: {compartments})", + ) + if feature_type == "Colocalization": + for key in ("channel1", "channel2"): + ch = request.get(key) + if ch not in channels: + raise argparse.ArgumentTypeError( + f"{key} {ch!r} in {feature_type} request is not " + f"declared via --image (valid: {channels})", + ) + else: + ch = request.get("channel") + if ch not in channels: + raise argparse.ArgumentTypeError( + f"channel {ch!r} in {feature_type} request is not " + f"declared via --image (valid: {channels})", + ) + + +def _build_image_set_loader( + images: list[tuple[str, Path]], + labels: list[tuple[str, Path]], + anisotropy_spacing: tuple[float, float, float], + identifiers: tuple[str, str, str, str], +) -> tuple[ImageSetLoader, list[str], list[str]]: + """Read images/labels and build an ImageSetLoader with identifiers.""" + image_dict: dict[str, object] = {} + for name, path in images: + image_dict[name] = _image_loading(path) + compartment_names: list[str] = [] + for name, path in labels: + image_dict[name] = _image_loading(path) + compartment_names.append(name) + patient_tumor, plate, well, field = identifiers + image_set_name = build_image_id_from_identifiers(identifiers) + image_set_loader = ImageSetLoader.from_image_dict( + image_dict, + anisotropy_spacing=anisotropy_spacing, + image_set_name=image_set_name, + label_key_names=compartment_names, + patient_tumor=patient_tumor, + plate=plate, + well=well, + field=field, + ) + channels = [name for name, _ in images] + return image_set_loader, channels, compartment_names + + +def build_image_id_from_identifiers( + identifiers: tuple[str, str, str, str], +) -> str: + """Build the deterministic image set name from identifier fields.""" + patient_tumor, plate, well, field = identifiers + return build_image_id(patient_tumor, plate, well, field) + + +def run( # noqa: PLR0913, PLR0917 + images: list[tuple[str, Path]], + labels: list[tuple[str, Path]], + anisotropy_spacing: tuple[float, float, float], + identifiers: tuple[str, str, str, str], + out_dir: Path, + feature_specs: list[dict[str, object]] | None = None, + features_filter: list[str] | None = None, + skip_existing: bool = False, + force: bool = False, +) -> list[Path]: + """Run feature extraction for one shard and write Parquet outputs. + + Parameters + ---------- + images : list[tuple[str, Path]] + (channel name, path) pairs from ``--image``. + labels : list[tuple[str, Path]] + (compartment name, path) pairs from ``--label``. + anisotropy_spacing : tuple[float, float, float] + (z, y, x) spacing. + identifiers : tuple[str, str, str, str] + (patient_tumor, plate, well, field) imaging coordinates. + out_dir : Path + Shard output directory. + feature_specs : list[dict[str, object]] | None + Explicit ``--feature`` requests; None means auto-generate. + features_filter : list[str] | None + ``--features`` selector restricting which feature types run. + skip_existing : bool + Skip a request whose output Parquet already exists. + force : bool + Overwrite even when the output exists (still crash-safe/atomic). + + Returns + ------- + list[Path] + Paths written (or skipped) per request. + """ + channels = [name for name, _ in images] + compartments = [name for name, _ in labels] + requests = _resolve_requests( + channels, + compartments, + feature_specs or [], + features_filter, + ) + if not requests: + print("No feature requests to run", file=sys.stderr) + return [] + # Compute target paths and partition into skipped vs. pending *before* + # reading any image. When every requested output already exists and + # --skip-existing is set, image I/O is avoided entirely. + results: list[Path] = [] + pending: list[tuple[dict[str, object], Path]] = [] + for request in requests: + compartment, channel, feature_type = _request_output_identity(request) + target = _output_path(out_dir, compartment, channel, feature_type) + if target.exists() and skip_existing and not force: + print(f"skip-existing: {target}", file=sys.stderr) + results.append(target) + else: + pending.append((request, target)) + if not pending: + print( + "all requested outputs exist; nothing to do", + file=sys.stderr, + ) + return results + image_set_loader, _, _ = _build_image_set_loader( + images, + labels, + anisotropy_spacing, + identifiers, + ) + for request, target in pending: + results.append(_execute_request(image_set_loader, request, out_dir, target)) + return results + + +def _build_parser() -> argparse.ArgumentParser: + """Construct the top-level argument parser.""" + parser = argparse.ArgumentParser( + prog="ZedProfiler", + description="Per-shard 3D featurization for the NF1 profiling warehouse.", + ) + subparsers = parser.add_subparsers(dest="command", required=True) + + run_parser = subparsers.add_parser( + "run", + help="Run feature extraction for one well/FOV shard.", + description=( + "Run feature extraction for one well/FOV shard. Loads one image " + "set from explicit file paths, runs the selected featurizers, and " + "writes one Parquet per feature table to --out-dir." + ), + ) + run_parser.add_argument( + "--image", + action="append", + required=True, + metavar="NAME=PATH", + help="Channel image as NAME=PATH (repeatable; >=1 required).", + ) + run_parser.add_argument( + "--label", + action="append", + required=True, + metavar="NAME=PATH", + help="Compartment label mask as NAME=PATH (repeatable; >=1 required).", + ) + run_parser.add_argument( + "--anisotropy-spacing", + nargs=3, + required=True, + type=float, + metavar=("Z", "Y", "X"), + help="Z, Y, X voxel spacing.", + ) + run_parser.add_argument( + "--patient-tumor", + required=True, + help="Patient-tumor identifier (e.g. NF0014_T1).", + ) + run_parser.add_argument("--plate", required=True, help="Plate identifier.") + run_parser.add_argument("--well", required=True, help="Well identifier.") + run_parser.add_argument( + "--field", + required=True, + help="Field-of-view index or identifier.", + ) + run_parser.add_argument( + "--out-dir", + required=True, + help="Shard output directory (created if needed).", + ) + run_parser.add_argument( + "--features", + default=None, + help=( + "Comma-separated feature types to run (selector). With no " + "--feature flags, runs these types over the channel x compartment " + "cross-product. With --feature flags, restricts those requests by " + "type. Default: all applicable single-channel types (plus " + "Colocalization when >=2 channels)." + ), + ) + run_parser.add_argument( + "--feature", + action="append", + default=[], + metavar="TYPE[,key=value,...]", + help=( + "Explicit feature request as TYPE[,key=value,...] (repeatable). " + "Examples: 'Intensity,channel=DNA,compartment=Nuclei'; " + "'Colocalization,channel1=DNA1,channel2=DNA2,compartment=Nuclei," + "fast_costes=Faster'." + ), + ) + run_parser.add_argument( + "--skip-existing", + action="store_true", + help="Skip a feature request whose output Parquet already exists.", + ) + run_parser.add_argument( + "--force", + action="store_true", + help="Overwrite even when the output exists (writes are still atomic).", + ) + return parser + + +def main(argv: Sequence[str] | None = None) -> int: + """CLI entry point. Returns a process exit code.""" + parser = _build_parser() + args = parser.parse_args(argv) + if args.command != "run": + parser.error("a subcommand is required") + return 2 # pragma: no cover - parser.error exits + + images = [_parse_name_path(token) for token in args.image] + labels = [_parse_name_path(token) for token in args.label] + feature_specs = [_parse_feature_spec(token) for token in args.feature] + features_filter = ( + [f.strip() for f in args.features.split(",") if f.strip()] + if args.features + else None + ) + identifiers = ( + args.patient_tumor, + args.plate, + args.well, + args.field, + ) + run( + images=images, + labels=labels, + anisotropy_spacing=tuple(args.anisotropy_spacing), + identifiers=identifiers, + out_dir=Path(args.out_dir), + feature_specs=feature_specs, + features_filter=features_filter, + skip_existing=args.skip_existing, + force=args.force, + ) + return 0 + + +def trigger() -> None: + """Console-script entry point (matches the pyproject ``scripts`` target).""" + raise SystemExit(main()) + + +if __name__ == "__main__": # pragma: no cover + trigger() diff --git a/src/zedprofiler/featurization/colocalization.py b/src/zedprofiler/featurization/colocalization.py index a98223f..917702b 100644 --- a/src/zedprofiler/featurization/colocalization.py +++ b/src/zedprofiler/featurization/colocalization.py @@ -626,6 +626,7 @@ def compute_colocalization( # noqa: C901, PLR0912 if full_name not in ( "Metadata_Object_ObjectID", "Metadata_Experiment_ImageSet", + "Metadata_Imaging_ImageID", ): try: row[full_name] = numpy.float32(meas_val) @@ -634,11 +635,12 @@ def compute_colocalization( # noqa: C901, PLR0912 else: row[full_name] = meas_val - # ensure object_id and image_set are present and first + # ensure object_id and image identifiers are present and first row["Metadata_Object_ObjectID"] = object_id row["Metadata_Experiment_ImageSet"] = ( two_object_loader.image_set_loader.image_set_name ) + row["Metadata_Imaging_ImageID"] = two_object_loader.image_set_loader.image_id list_of_dfs.append(row) # Convert list of row-dicts into a dict-of-lists with stable ordering @@ -649,7 +651,11 @@ def compute_colocalization( # noqa: C901, PLR0912 other_keys: list[str] = [] for d in list_of_dfs: for k in d: - if k in ("Metadata_Object_ObjectID", "Metadata_Experiment_ImageSet"): + if k in ( + "Metadata_Object_ObjectID", + "Metadata_Experiment_ImageSet", + "Metadata_Imaging_ImageID", + ): continue if k not in other_keys: other_keys.append(k) @@ -657,6 +663,7 @@ def compute_colocalization( # noqa: C901, PLR0912 all_keys = [ "Metadata_Object_ObjectID", "Metadata_Experiment_ImageSet", + "Metadata_Imaging_ImageID", *other_keys, ] result: dict[str, list[object]] = { diff --git a/src/zedprofiler/featurization/granularity.py b/src/zedprofiler/featurization/granularity.py index 72c1f78..d479b29 100644 --- a/src/zedprofiler/featurization/granularity.py +++ b/src/zedprofiler/featurization/granularity.py @@ -459,6 +459,11 @@ def compute_granularity( # noqa: C901, PLR0912, PLR0913, PLR0915 }, inplace=True, ) + final_df.insert( + 0, + "Metadata_Imaging_ImageID", + object_loader.image_set_loader.image_id, + ) final_df.insert( 0, "Metadata_Experiment_ImageSet", diff --git a/src/zedprofiler/featurization/intensity.py b/src/zedprofiler/featurization/intensity.py index 7fea74e..b25f9f7 100644 --- a/src/zedprofiler/featurization/intensity.py +++ b/src/zedprofiler/featurization/intensity.py @@ -221,6 +221,11 @@ def compute_intensity( # noqa: PLR0915 inplace=True, ) + final_df.insert( + 0, + "Metadata_Imaging_ImageID", + object_loader.image_set_loader.image_id, + ) final_df.insert( 0, "Metadata_Experiment_ImageSet", diff --git a/src/zedprofiler/featurization/neighbors.py b/src/zedprofiler/featurization/neighbors.py index d92d28f..0a35df6 100644 --- a/src/zedprofiler/featurization/neighbors.py +++ b/src/zedprofiler/featurization/neighbors.py @@ -199,6 +199,11 @@ def compute_neighbors( inplace=True, ) if not final_df.empty: + final_df.insert( + 0, + "Metadata_Imaging_ImageID", + object_loader.image_set_loader.image_id, + ) final_df.insert( 0, "Metadata_Experiment_ImageSet", diff --git a/src/zedprofiler/featurization/texture.py b/src/zedprofiler/featurization/texture.py index 6920e54..93a3817 100644 --- a/src/zedprofiler/featurization/texture.py +++ b/src/zedprofiler/featurization/texture.py @@ -208,6 +208,11 @@ def compute_texture( # noqa: C901 }, inplace=True, ) + final_df.insert( + 0, + "Metadata_Imaging_ImageID", + object_loader.image_set_loader.image_id, + ) final_df.insert( 0, "Metadata_Experiment_ImageSet", diff --git a/src/zedprofiler/featurization/volumesizeshape.py b/src/zedprofiler/featurization/volumesizeshape.py index 237a2d8..808224c 100644 --- a/src/zedprofiler/featurization/volumesizeshape.py +++ b/src/zedprofiler/featurization/volumesizeshape.py @@ -206,6 +206,11 @@ def measure_3D_volume_size_shape( inplace=True, ) + final_df.insert( + 0, + "Metadata_Imaging_ImageID", + object_loader.image_set_loader.image_id, + ) final_df.insert( 1, "Metadata_Experiment_ImageSet", diff --git a/src/zedprofiler/identifiers.py b/src/zedprofiler/identifiers.py new file mode 100644 index 0000000..0cf8179 --- /dev/null +++ b/src/zedprofiler/identifiers.py @@ -0,0 +1,55 @@ +"""Deterministic imaging identifiers for warehouse join keys. + +The NF1 bioimage profiling warehouse joins every image, object, feature +table, and annotation via stable identifiers (see the future processing +plan's identifier spec). The central one is ``Metadata_Imaging_ImageID``, +built deterministically from the four imaging coordinates: + + patient-tumor, plate, well, field + +Because a shard is dispatched per well/FOV, every feature table emitted by a +shard carries this single image id so downstream tables can rejoin without a +database service. + +This module is the single source of truth for the id format. Changing the +format here changes every shard's output ids at once. +""" + +from __future__ import annotations + +from beartype import beartype + + +@beartype +def build_image_id( + patient_tumor: str, + plate: str, + well: str, + field: int | str, +) -> str: + """Build a deterministic ``Metadata_Imaging_ImageID`` value. + + The id is a stable string assembled from the four imaging coordinates so + that the same well/FOV always produces the same id across runs, batches, + and reprocessing. Component order is fixed + (patient-tumor, plate, well, field) so ids sort and group naturally by + patient then plate then well then field. + + Parameters + ---------- + patient_tumor : str + Patient-tumor identifier (e.g. ``"NF0014_T1"``). + plate : str + Plate identifier (e.g. ``"PLATE01"``). + well : str + Well identifier (e.g. ``"A1"``). + field : int | str + Field-of-view index or identifier (e.g. ``1`` or ``"f1"``). + + Returns + ------- + str + The deterministic image id, e.g. ``"NF0014_T1_PLATE01_A1_field1"``. + + """ + return f"{patient_tumor}_{plate}_{well}_field{field}" diff --git a/tests/IO/test_loading_classes.py b/tests/IO/test_loading_classes.py index 831987a..4d5d4cf 100644 --- a/tests/IO/test_loading_classes.py +++ b/tests/IO/test_loading_classes.py @@ -382,3 +382,123 @@ def test_two_object_loader_loads_images_and_ids(self) -> None: assert np.array_equal(obj.image2, image) assert np.array_equal(obj.label_image, label_image) assert obj.object_ids == [ONE_LABEL, TWO_LABEL] + + +class TestImageSetConfigImageId: + """Tests for ImageSetConfig.image_id identifier propagation.""" + + def test_image_id_is_none_when_any_field_missing(self) -> None: + """image_id stays None until all four identifier fields are set.""" + assert ImageSetConfig().image_id is None + assert ( + ImageSetConfig( + patient_tumor="NF0014_T1", + plate="PLATE01", + well="A1", + field=None, + ).image_id + is None + ) + + def test_image_id_built_when_all_fields_set(self) -> None: + """With all four fields set, image_id is the deterministic image id.""" + config = ImageSetConfig( + patient_tumor="NF0014_T1", + plate="PLATE01", + well="A1", + field=1, + ) + assert config.image_id == "NF0014_T1_PLATE01_A1_field1" + + +class TestFromImageDict: + """Tests for ImageSetLoader.from_image_dict (multi-channel shard loader).""" + + def _build_dict(self) -> dict[str, np.ndarray]: + label = np.array( + [ + [[ZERO_LABEL, ONE_LABEL], [TWO_LABEL, TWO_LABEL]], + [[ZERO_LABEL, ONE_LABEL], [TWO_LABEL, TWO_LABEL]], + ], + dtype=np.int32, + ) + return { + "DNA": np.ones((2, 2, 2), dtype=np.float32), + "AGP": np.full((2, 2, 2), 5.0, dtype=np.float32), + "Nuclei": label, + } + + def test_builds_working_multi_channel_loader(self) -> None: + """from_image_dict resolves compartments and channels from the dict.""" + loader = ImageSetLoader.from_image_dict( + self._build_dict(), + anisotropy_spacing=(2.0, 1.0, 1.0), + image_set_name="shard-01", + label_key_names=["Nuclei"], + ) + assert loader.image_set_name == "shard-01" + assert loader.anisotropy_factor == EXPECTED_ANISOTROPY + assert loader.compartments == ["Nuclei"] + assert sorted(loader.image_names) == ["AGP", "DNA"] + assert loader.unique_compartment_objects["Nuclei"] == [ + ONE_LABEL, + TWO_LABEL, + ] + # ObjectLoader built on the resulting loader resolves both channels. + dna_obj = ObjectLoader( + image_set_loader=loader, + channel_name="DNA", + compartment_name="Nuclei", + ) + assert dna_obj.object_ids == [ONE_LABEL, TWO_LABEL] + assert np.array_equal( + dna_obj.image, + np.ones((2, 2, 2), dtype=np.float32), + ) + + def test_image_id_populated_when_identifiers_provided(self) -> None: + """Passing all four identifier fields populates image_id deterministically.""" + loader = ImageSetLoader.from_image_dict( + self._build_dict(), + anisotropy_spacing=(1.0, 1.0, 1.0), + image_set_name="shard-01", + label_key_names=["Nuclei"], + patient_tumor="NF0014_T1", + plate="PLATE01", + well="A1", + field=1, + ) + assert loader.image_id == "NF0014_T1_PLATE01_A1_field1" + + def test_image_id_falls_back_to_image_set_name_without_identifiers( + self, + ) -> None: + """Without identifiers, image_id falls back to the image set name.""" + loader = ImageSetLoader.from_image_dict( + self._build_dict(), + anisotropy_spacing=(1.0, 1.0, 1.0), + image_set_name="legacy-set", + label_key_names=["Nuclei"], + ) + assert loader.image_id == "legacy-set" + + def test_two_object_loader_on_from_image_dict(self) -> None: + """A TwoObjectLoader resolves two channels from a from_image_dict loader.""" + loader = ImageSetLoader.from_image_dict( + self._build_dict(), + anisotropy_spacing=(1.0, 1.0, 1.0), + image_set_name="shard-02", + label_key_names=["Nuclei"], + patient_tumor="X", + plate="P", + well="A1", + field="1", + ) + two = TwoObjectLoader( + image_set_loader=loader, + compartment="Nuclei", + channel1="DNA", + channel2="AGP", + ) + assert two.object_ids == [ONE_LABEL, TWO_LABEL] + assert two.image_set_loader.image_id == "X_P_A1_field1" diff --git a/tests/featurization/test_colocalization.py b/tests/featurization/test_colocalization.py index 7b34252..94a4c1f 100644 --- a/tests/featurization/test_colocalization.py +++ b/tests/featurization/test_colocalization.py @@ -22,6 +22,8 @@ class ImageSetLoaderModel(BaseModel): model_config = ConfigDict(arbitrary_types_allowed=True) image_set_name: str = "coloc" + # mirrors ImageSetLoader.image_id (falls back to image_set_name) + image_id: str = "coloc" class TwoObjectLoaderModel(BaseModel): diff --git a/tests/featurization/test_granularity.py b/tests/featurization/test_granularity.py index c2f13b8..acb0769 100644 --- a/tests/featurization/test_granularity.py +++ b/tests/featurization/test_granularity.py @@ -20,6 +20,8 @@ class ImageSetLoaderModel(BaseModel): model_config = ConfigDict(arbitrary_types_allowed=True) image_set_name: str = "gran" + # mirrors ImageSetLoader.image_id (falls back to image_set_name) + image_id: str = "gran" class ObjectLoaderModel(BaseModel): @@ -94,7 +96,7 @@ class Dummy: image = img label_image = lab object_ids: ClassVar[list[int]] = [1] - image_set_loader = type("ISL", (), {"image_set_name": "s"})() + image_set_loader = type("ISL", (), {"image_set_name": "s", "image_id": "s"})() compartment = "Cell" channel = "Ch1" @@ -120,7 +122,7 @@ class Dummy: image = img label_image = lab object_ids: ClassVar[list[int]] = [1] - image_set_loader = type("ISL", (), {"image_set_name": "s"})() + image_set_loader = type("ISL", (), {"image_set_name": "s", "image_id": "s"})() compartment = "Cell" channel = "Ch1" @@ -150,7 +152,7 @@ class Dummy: image = img label_image = lab object_ids: ClassVar[list[int]] = [1] - image_set_loader = type("ISL", (), {"image_set_name": "s"})() + image_set_loader = type("ISL", (), {"image_set_name": "s", "image_id": "s"})() compartment = "Cell" channel = "Ch1" @@ -195,7 +197,7 @@ class Dummy: image = img label_image = lab object_ids: ClassVar[list[int]] = [1] - image_set_loader = type("ISL", (), {"image_set_name": "s"})() + image_set_loader = type("ISL", (), {"image_set_name": "s", "image_id": "s"})() compartment = "Cell" channel = "Ch1" @@ -222,7 +224,7 @@ class Dummy: image = img label_image = lab object_ids: ClassVar[list[int]] = [257, 514] - image_set_loader = type("ISL", (), {"image_set_name": "s"})() + image_set_loader = type("ISL", (), {"image_set_name": "s", "image_id": "s"})() compartment = "Cell" channel = "Ch1" diff --git a/tests/featurization/test_intensity.py b/tests/featurization/test_intensity.py index a2593ea..26e9310 100644 --- a/tests/featurization/test_intensity.py +++ b/tests/featurization/test_intensity.py @@ -12,6 +12,8 @@ class ImageSetLoaderModel(BaseModel): model_config = ConfigDict(arbitrary_types_allowed=True) image_set_name: str = "intensity" + # mirrors ImageSetLoader.image_id (falls back to image_set_name) + image_id: str = "intensity" class ObjectLoaderModel(BaseModel): diff --git a/tests/featurization/test_neighbors.py b/tests/featurization/test_neighbors.py index 398d12f..5167287 100644 --- a/tests/featurization/test_neighbors.py +++ b/tests/featurization/test_neighbors.py @@ -23,6 +23,8 @@ class ImageSetLoaderModel(BaseModel): model_config = ConfigDict(arbitrary_types_allowed=True) image_set_name: str = "neighbors" + # mirrors ImageSetLoader.image_id (falls back to image_set_name) + image_id: str = "neighbors" class ObjectLoaderModel(BaseModel): diff --git a/tests/featurization/test_neighbors_additional.py b/tests/featurization/test_neighbors_additional.py index 2da68d7..fa2c0c1 100644 --- a/tests/featurization/test_neighbors_additional.py +++ b/tests/featurization/test_neighbors_additional.py @@ -50,7 +50,7 @@ def test_compute_neighbors_distance_counts() -> None: class Dummy: label_image = lab object_ids = (1, 2, 3) - image_set_loader = type("ISL", (), {"image_set_name": "s"})() + image_set_loader = type("ISL", (), {"image_set_name": "s", "image_id": "s"})() compartment = "Cell" channel = "Ch1" diff --git a/tests/featurization/test_texture.py b/tests/featurization/test_texture.py index 659eb92..7f9c02e 100644 --- a/tests/featurization/test_texture.py +++ b/tests/featurization/test_texture.py @@ -18,6 +18,8 @@ class ImageSetLoaderModel(BaseModel): model_config = ConfigDict(arbitrary_types_allowed=True) image_set_name: str = "texture" + # mirrors ImageSetLoader.image_id (falls back to image_set_name) + image_id: str = "texture" class ObjectLoaderModel(BaseModel): diff --git a/tests/featurization/test_volumesizeshape.py b/tests/featurization/test_volumesizeshape.py index 02637f0..985959d 100644 --- a/tests/featurization/test_volumesizeshape.py +++ b/tests/featurization/test_volumesizeshape.py @@ -15,6 +15,8 @@ class ImageSetLoaderModel(BaseModel): model_config = ConfigDict(arbitrary_types_allowed=True) anisotropy_spacing: tuple[float, float, float] image_set_name: str = "testset" + # mirrors ImageSetLoader.image_id (falls back to image_set_name) + image_id: str = "testset" class ObjectLoaderModel(BaseModel): diff --git a/tests/test_cli.py b/tests/test_cli.py index 4f9aeed..92d20d6 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -1,8 +1,365 @@ -"""Package namespace tests replacing obsolete template CLI checks.""" +"""End-to-end and unit tests for the per-shard CLI (``zedprofiler.cli``). +The end-to-end tests exercise ``ZedProfiler run`` against the CellProfiler 3D +tutorial data (the same fixtures used by ``test_real_world_data.py``) and assert +the things that matter to the NF1 pipeline: one Parquet per requested feature +table, a deterministic ``Metadata_Imaging_ImageID``, ``--features`` selection, +``--skip-existing`` idempotency, and two-channel colocalization. + +The unit tests cover the request-selection and path-naming helpers directly so +the full default feature matrix (which includes slow Texture/Granularity +runs) does not have to run end-to-end here. +""" + +from __future__ import annotations + +import argparse +from pathlib import Path + +import pandas as pd +import pytest + +from zedprofiler.cli import ( + _auto_requests, + _output_path, + _parse_feature_spec, + _parse_name_path, + _resolve_requests, + main, +) from zedprofiler.featurization import texture +from zedprofiler.identifiers import build_image_id + +tifffile = pytest.importorskip("tifffile") + +TUTORIAL_ROOT = ( + Path(__file__).resolve().parent + / "data" + / "CP_tutorial_3D_noise_nuclei_segmentation" +) +IMAGE1 = TUTORIAL_ROOT / "input" / "nuclei1_out_c00_dr90_image.tif" +IMAGE2 = TUTORIAL_ROOT / "input" / "nuclei2_out_c90_dr90_image.tif" +LABEL1 = ( + TUTORIAL_ROOT + / "output" + / "masks" + / "nuclei1_out_c00_dr90_imageSegmentationMask.tiff" +) + +EXPECTED_OBJECT_COUNT = 5 +EXPECTED_CHANNEL_COUNT = 2 +PATIENT_TUMOR, PLATE, WELL, FIELD = "NF0014_T1", "PLATE01", "A1", "1" +EXPECTED_IMAGE_ID = build_image_id(PATIENT_TUMOR, PLATE, WELL, FIELD) def test_feature_namespace_import() -> None: """Lower-level feature namespace remains importable.""" assert texture.__name__ == "zedprofiler.featurization.texture" + + +def _run(argv: list[str]) -> int: + return main(argv) + + +# --------------------------------------------------------------------------- +# Unit tests: name/path + feature-spec parsing +# --------------------------------------------------------------------------- + + +def test_parse_name_path_splits_on_first_equals() -> None: + """NAME=PATH splits on the first '=' so paths may contain '='.""" + name, path = _parse_name_path("DNA=/tmp/a=b.tif") + assert name == "DNA" + assert path == Path("/tmp/a=b.tif") + + +def test_parse_name_path_rejects_missing_equals() -> None: + """A token without '=' is a user error, not a silent default.""" + with pytest.raises(argparse.ArgumentTypeError): + _parse_name_path("DNA.tif") + + +def test_parse_feature_spec_parses_type_and_overrides() -> None: + """TYPE is the first comma token; the rest are key=value overrides.""" + request = _parse_feature_spec("Intensity,channel=DNA,compartment=Nuclei") + assert request == { + "type": "Intensity", + "channel": "DNA", + "compartment": "Nuclei", + } + + +def test_parse_feature_spec_rejects_unknown_type() -> None: + """An unknown feature type raises a friendly argument error.""" + with pytest.raises(argparse.ArgumentTypeError, match="Unknown feature type"): + _parse_feature_spec("NotAFeature,channel=DNA") + + +def test_parse_feature_spec_rejects_bad_override() -> None: + """An override without '=' raises a friendly argument error.""" + with pytest.raises(argparse.ArgumentTypeError, match="key=value"): + _parse_feature_spec("Intensity,bogus") + + +def test_output_path_mirrors_save_features_as_parquet_naming() -> None: + """_output_path must match the path save_features_as_parquet writes to.""" + assert _output_path(Path("/out"), "Nuclei", "DNA", "Intensity") == Path( + "/out/Nuclei_DNA_Intensity_cpu_features.parquet" + ) + assert _output_path(Path("/out"), "Nuclei", "DNA1-DNA2", "Colocalization") == Path( + "/out/Nuclei_DNA1-DNA2_Colocalization_cpu_features.parquet" + ) + + +# --------------------------------------------------------------------------- +# Unit tests: request selection (no featurizers run) +# --------------------------------------------------------------------------- + + +def test_auto_requests_single_channel_no_colocalization() -> None: + """One channel x one compartment yields the 5 single-channel types only.""" + requests = _auto_requests( + ["DNA"], + ["Nuclei"], + ["VolumeSizeShape", "Intensity", "Neighbors", "Texture", "Granularity"], + ) + types = sorted(str(r["type"]) for r in requests) + assert types == sorted( + ["VolumeSizeShape", "Intensity", "Neighbors", "Texture", "Granularity"], + ) + # Channel-agnostic features use the first (only) channel for naming. + vol = next(r for r in requests if r["type"] == "VolumeSizeShape") + assert vol["channel"] == "DNA" + + +def test_auto_requests_two_channels_adds_colocalization_pairs() -> None: + """Two channels produce a single ordered colocalization pair x compartment.""" + requests = _auto_requests( + ["DNA1", "DNA2"], + ["Nuclei"], + ["Intensity", "Colocalization"], + ) + coloc = [r for r in requests if r["type"] == "Colocalization"] + assert len(coloc) == 1 + assert coloc[0]["channel1"] == "DNA1" + assert coloc[0]["channel2"] == "DNA2" + # Intensity runs per channel x compartment. + assert ( + sum(1 for r in requests if r["type"] == "Intensity") == EXPECTED_CHANNEL_COUNT + ) + + +def test_auto_requests_colocalization_requires_two_channels() -> None: + """Requesting Colocalization with one channel is a user error.""" + with pytest.raises(argparse.ArgumentTypeError, match="at least two"): + _auto_requests(["DNA"], ["Nuclei"], ["Colocalization"]) + + +def test_resolve_requests_default_two_channels_includes_colocalization() -> None: + """Default selection (no --features, no --feature) adds coloc for >=2 channels.""" + requests = _resolve_requests(["DNA1", "DNA2"], ["Nuclei"], [], None) + types = {str(r["type"]) for r in requests} + assert "Colocalization" in types + assert "Intensity" in types + + +def test_resolve_requests_features_filter_restricts_types() -> None: + """--features Intensity selects only Intensity from the cross-product.""" + requests = _resolve_requests( + ["DNA1", "DNA2"], + ["Nuclei"], + [], + ["Intensity"], + ) + assert all(r["type"] == "Intensity" for r in requests) + assert len(requests) == EXPECTED_CHANNEL_COUNT # 2 channels x 1 compartment + + +def test_resolve_requests_explicit_specs_validated_against_declared() -> None: + """An explicit --feature referencing an undeclared channel is rejected.""" + specs = [_parse_feature_spec("Intensity,channel=Ghost,compartment=Nuclei")] + with pytest.raises(argparse.ArgumentTypeError, match="declared via --image"): + _resolve_requests(["DNA"], ["Nuclei"], specs, None) + + +def test_resolve_requests_features_filter_applied_to_explicit_specs() -> None: + """--features restricts explicit --feature requests by type.""" + specs = [ + _parse_feature_spec("Intensity,channel=DNA,compartment=Nuclei"), + _parse_feature_spec("Texture,channel=DNA,compartment=Nuclei"), + ] + requests = _resolve_requests(["DNA"], ["Nuclei"], specs, ["Intensity"]) + assert len(requests) == 1 + assert requests[0]["type"] == "Intensity" + + +# --------------------------------------------------------------------------- +# End-to-end tests on the CellProfiler 3D tutorial data +# --------------------------------------------------------------------------- + +# The end-to-end tests read the CellProfiler 3D tutorial images/masks, which +# are added by a separate data commit and may be absent on some branches. Skip +# them when the data is not present so the CLI test module stays green +# everywhere; they run in full wherever the tutorial data is available. +requires_tutorial_data = pytest.mark.skipif( + not TUTORIAL_ROOT.exists(), + reason="CellProfiler 3D tutorial data not present on this branch", +) + + +def _base_run_args(out_dir: Path, *extra: str) -> list[str]: + return [ + "run", + f"--image=DNA={IMAGE1}", + f"--label=Nuclei={LABEL1}", + "--anisotropy-spacing", + "1.0", + "1.0", + "1.0", + f"--patient-tumor={PATIENT_TUMOR}", + f"--plate={PLATE}", + f"--well={WELL}", + f"--field={FIELD}", + f"--out-dir={out_dir}", + *extra, + ] + + +@requires_tutorial_data +def test_cli_run_intensity_writes_parquet_with_image_id( + tmp_path: Path, +) -> None: + """A single Intensity request writes one Parquet carrying the image id.""" + out_dir = tmp_path / "shard" + assert _run(_base_run_args(out_dir, "--features=Intensity")) == 0 + + parquet = _output_path(out_dir, "Nuclei", "DNA", "Intensity") + assert parquet.exists() + + df = pd.read_parquet(parquet) + assert len(df) == EXPECTED_OBJECT_COUNT + assert "Metadata_Imaging_ImageID" in df.columns + assert "Metadata_Experiment_ImageSet" in df.columns + assert set(df["Metadata_Imaging_ImageID"]) == {EXPECTED_IMAGE_ID} + # No leftover temp files from the atomic write. + assert not any(p.suffix == ".tmp" for p in out_dir.iterdir()) + + +@requires_tutorial_data +def test_cli_features_selector_restricts_outputs(tmp_path: Path) -> None: + """--features controls exactly which feature tables are written.""" + out_dir = tmp_path / "shard" + assert _run(_base_run_args(out_dir, "--features=VolumeSizeShape,Intensity")) == 0 + files = sorted(p.name for p in out_dir.glob("*.parquet")) + assert files == [ + "Nuclei_DNA_Intensity_cpu_features.parquet", + "Nuclei_DNA_VolumeSizeShape_cpu_features.parquet", + ] + + out_dir2 = tmp_path / "shard2" + assert _run(_base_run_args(out_dir2, "--features=Intensity")) == 0 + assert [p.name for p in out_dir2.glob("*.parquet")] == [ + "Nuclei_DNA_Intensity_cpu_features.parquet", + ] + + +@requires_tutorial_data +def test_cli_rerun_is_content_identical(tmp_path: Path) -> None: + """Re-running without --skip-existing reproduces identical feature content.""" + out_dir = tmp_path / "shard" + _run(_base_run_args(out_dir, "--features=Intensity")) + first = pd.read_parquet(_output_path(out_dir, "Nuclei", "DNA", "Intensity")) + + out_dir2 = tmp_path / "shard2" + _run(_base_run_args(out_dir2, "--features=Intensity")) + second = pd.read_parquet(_output_path(out_dir2, "Nuclei", "DNA", "Intensity")) + + pd.testing.assert_frame_equal(first, second) + + +@requires_tutorial_data +def test_cli_skip_existing_skips_recompute(tmp_path: Path) -> None: + """--skip-existing leaves finished outputs untouched and skips image I/O.""" + out_dir = tmp_path / "shard" + _run(_base_run_args(out_dir, "--features=Intensity")) + target = _output_path(out_dir, "Nuclei", "DNA", "Intensity") + first_bytes = target.read_bytes() + first_mtime_ns = target.stat().st_mtime_ns + + # Re-run with --skip-existing using *nonexistent* image paths: if the CLI + # tried to load images it would fail, proving skip happens before I/O. + assert ( + _run( + [ + "run", + "--image=DNA=/does/not/exist.tif", + "--label=Nuclei=/does/not/exist.tiff", + "--anisotropy-spacing", + "1.0", + "1.0", + "1.0", + f"--patient-tumor={PATIENT_TUMOR}", + f"--plate={PLATE}", + f"--well={WELL}", + f"--field={FIELD}", + f"--out-dir={out_dir}", + "--features=Intensity", + "--skip-existing", + ], + ) + == 0 + ) + assert target.read_bytes() == first_bytes + assert target.stat().st_mtime_ns == first_mtime_ns + + +@requires_tutorial_data +def test_cli_colocalization_two_channels(tmp_path: Path) -> None: + """An explicit colocalization request writes a DNA1-DNA2 Parquet.""" + out_dir = tmp_path / "shard" + assert ( + _run( + [ + "run", + f"--image=DNA1={IMAGE1}", + f"--image=DNA2={IMAGE2}", + f"--label=Nuclei={LABEL1}", + "--anisotropy-spacing", + "1.0", + "1.0", + "1.0", + f"--patient-tumor={PATIENT_TUMOR}", + f"--plate={PLATE}", + f"--well={WELL}", + f"--field={FIELD}", + f"--out-dir={out_dir}", + "--feature=Colocalization,channel1=DNA1,channel2=DNA2,compartment=Nuclei,fast_costes=Faster", + ], + ) + == 0 + ) + target = _output_path(out_dir, "Nuclei", "DNA1-DNA2", "Colocalization") + assert target.exists() + df = pd.read_parquet(target) + assert len(df) == EXPECTED_OBJECT_COUNT + assert set(df["Metadata_Imaging_ImageID"]) == {EXPECTED_IMAGE_ID} + assert any("Colocalization" in c for c in df.columns) + + +def test_cli_missing_required_arg_errors(tmp_path: Path) -> None: + """A run without --out-dir exits non-zero (argparse error).""" + argv = [ + "run", + f"--image=DNA={IMAGE1}", + f"--label=Nuclei={LABEL1}", + "--anisotropy-spacing", + "1.0", + "1.0", + "1.0", + f"--patient-tumor={PATIENT_TUMOR}", + f"--plate={PLATE}", + f"--well={WELL}", + f"--field={FIELD}", + ] + with pytest.raises(SystemExit): + _run(argv) diff --git a/tests/test_identifiers.py b/tests/test_identifiers.py new file mode 100644 index 0000000..363cefd --- /dev/null +++ b/tests/test_identifiers.py @@ -0,0 +1,44 @@ +"""Tests for the deterministic Metadata_Imaging_ImageID builder.""" + +from __future__ import annotations + +import pytest + +from zedprofiler.identifiers import build_image_id + + +def test_build_image_id_formats_all_fields() -> None: + """The image id joins patient-tumor, plate, well, and a fielded suffix.""" + assert ( + build_image_id("NF0014_T1", "PLATE01", "A1", 1) == "NF0014_T1_PLATE01_A1_field1" + ) + + +def test_build_image_id_is_deterministic() -> None: + """Repeated calls with the same inputs produce the same id.""" + args = ("NF0014_T1", "PLATE01", "A1", "2") + assert build_image_id(*args) == build_image_id(*args) + + +def test_build_image_id_field_accepts_string_or_int() -> None: + """Field may be supplied as either an int or its string form.""" + assert build_image_id("X", "P", "A1", 3) == "X_P_A1_field3" + assert build_image_id("X", "P", "A1", "3") == "X_P_A1_field3" + + +@pytest.mark.parametrize( + "patient_tumor,plate,well,field,expected", + [ + ("NF0014_T1", "PLATE01", "A1", "1", "NF0014_T1_PLATE01_A1_field1"), + ("NF0009_T2", "PLATE02", "B3", "7", "NF0009_T2_PLATE02_B3_field7"), + ], +) +def test_build_image_id_parametrized( + patient_tumor: str, + plate: str, + well: str, + field: str, + expected: str, +) -> None: + """Format holds across distinct imaging coordinates.""" + assert build_image_id(patient_tumor, plate, well, field) == expected From 49b3c10e6636a95e4c8d04e621b477b7c9477244 Mon Sep 17 00:00:00 2001 From: d33bs Date: Wed, 5 Aug 2026 17:17:12 -0600 Subject: [PATCH 2/7] Add per-shard CLI and Metadata_Imaging_ImageID propagation ZedProfiler is the feature extractor the NF1 pipeline dispatches per well/FOV shard via SLURM sbatch. This commit adds the command that process runs and the identifier column that makes shards warehouse-joinable. CLI (src/zedprofiler/cli.py): - "ZedProfiler run" subcommand (argparse): repeatable --image/--label NAME=PATH flags, --anisotropy-spacing Z Y X, identifier fields (--patient-tumor/--plate/--well/--field), --out-dir, a --features selector, repeatable --feature TYPE[,key=value,...] advanced requests, --skip-existing, and --force. - Reuses the six compute_* featurizers; builds the shared loader via ImageSetLoader.from_image_dict with the identifier fields. - Restartable/idempotent: deterministic output paths, --skip-existing filters before any image is read (a finished shard re-run skips I/O entirely), and atomic writes (temp + os.replace) so a crashed shard never leaves a partial file that --skip-existing would mistake for complete. - Fixes the orphaned/mis-cased console script to zedprofiler.cli:trigger. Identifiers (src/zedprofiler/identifiers.py): - build_image_id(patient_tumor, plate, well, field) -> deterministic Metadata_Imaging_ImageID; single source of truth for the format. Loaders (src/zedprofiler/IO/loading_classes.py): - ImageSetConfig carries patient_tumor/plate/well/field with an image_id property; ImageSetLoader exposes image_id (falls back to image_set_name). - New from_image_dict classmethod builds a multi-channel loader from an in-memory {key: ndarray} dict (the path the CLI needs); it derives compartments/image names directly from the declared label keys so it is self-contained and correct independent of get_compartments. Featurizers (6 modules): - Each emits Metadata_Imaging_ImageID before Metadata_Experiment_ImageSet. Feature values are unchanged; only a metadata column is added. Feature writing: - save_features_as_parquet gains an opt-in atomic flag used by the CLI. Tests: - test_real_world_data.py colocalization loader switched to from_image_dict so the new image_id column is populated for the existing colocalization end-to-end test. - CLI end-to-end tests run against the CellProfiler 3D tutorial data present on main and self-skip if that data is absent on other branches. --- pyproject.toml | 2 +- src/zedprofiler/IO/feature_writing_utils.py | 15 +- src/zedprofiler/IO/loading_classes.py | 133 ++++ src/zedprofiler/cli.py | 679 ++++++++++++++++++ .../featurization/colocalization.py | 11 +- src/zedprofiler/featurization/granularity.py | 5 + src/zedprofiler/featurization/intensity.py | 5 + src/zedprofiler/featurization/neighbors.py | 5 + src/zedprofiler/featurization/texture.py | 5 + .../featurization/volumesizeshape.py | 5 + src/zedprofiler/identifiers.py | 55 ++ tests/IO/test_loading_classes.py | 120 ++++ tests/featurization/test_colocalization.py | 2 + tests/featurization/test_granularity.py | 12 +- tests/featurization/test_intensity.py | 2 + tests/featurization/test_neighbors.py | 2 + .../test_neighbors_additional.py | 2 +- tests/featurization/test_real_world_data.py | 29 +- tests/featurization/test_texture.py | 2 + tests/featurization/test_volumesizeshape.py | 2 + tests/test_cli.py | 359 ++++++++- tests/test_identifiers.py | 44 ++ 22 files changed, 1470 insertions(+), 26 deletions(-) create mode 100644 src/zedprofiler/cli.py create mode 100644 src/zedprofiler/identifiers.py create mode 100644 tests/test_identifiers.py diff --git a/pyproject.toml b/pyproject.toml index 29e0a59..e63c460 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -42,7 +42,7 @@ dependencies = [ ] urls.Homepage = "https://zedprofiler.readthedocs.io" urls.Repository = "https://github.com/WayScience/ZedProfiler" -scripts.ZedProfiler = "ZedProfiler.cli:trigger" +scripts.ZedProfiler = "zedprofiler.cli:trigger" [dependency-groups] dev = [ diff --git a/src/zedprofiler/IO/feature_writing_utils.py b/src/zedprofiler/IO/feature_writing_utils.py index 6769da7..3ac8ab9 100644 --- a/src/zedprofiler/IO/feature_writing_utils.py +++ b/src/zedprofiler/IO/feature_writing_utils.py @@ -6,6 +6,7 @@ from __future__ import annotations import dataclasses +import os import pathlib import pandas @@ -181,6 +182,7 @@ def save_features_as_parquet( parent_path: pathlib.Path, df: pandas.DataFrame, metadata: FeatureMetadata, + atomic: bool = False, ) -> pathlib.Path: """Save features as parquet files in a consistent way. @@ -196,6 +198,12 @@ def save_features_as_parquet( metadata : FeatureMetadata Metadata for the feature output (compartment, channel, feature_type, cpu_or_gpu). + atomic : bool + When True, write to a sibling ``.tmp`` file then atomically replace + the destination via ``os.replace``. This prevents a crashed write + from leaving a partial parquet file that a restartable caller (using + ``--skip-existing``) could mistake for a complete one. Default False + preserves the existing direct-write behavior. Returns ------- @@ -210,5 +218,10 @@ def save_features_as_parquet( metadata.cpu_or_gpu, ) save_path = parent_path / f"{output_prefix}_features.parquet" - validated_df.to_parquet(save_path, index=False) + if atomic: + tmp_path = save_path.with_suffix(save_path.suffix + ".tmp") + validated_df.to_parquet(tmp_path, index=False) + os.replace(tmp_path, save_path) + else: + validated_df.to_parquet(save_path, index=False) return save_path diff --git a/src/zedprofiler/IO/loading_classes.py b/src/zedprofiler/IO/loading_classes.py index 12d5485..38e380c 100644 --- a/src/zedprofiler/IO/loading_classes.py +++ b/src/zedprofiler/IO/loading_classes.py @@ -12,6 +12,7 @@ from beartype import beartype from zedprofiler.contracts import ImageArrayModel +from zedprofiler.identifiers import build_image_id logging.basicConfig(level=logging.INFO) @@ -42,6 +43,15 @@ class ImageSetConfig: image_set_name: str | None = None label_key_name: list[str] | None = None raw_image_key_name: list[str] | None = None + # Imaging-coordinate identifier fields used to build a deterministic + # ``Metadata_Imaging_ImageID``. All four must be set for ``image_id`` to be + # populated; when any is None, ``image_id`` is None and the loader falls + # back to ``image_set_name`` for the emitted metadata column (see + # ``ImageSetLoader.image_id``). + patient_tumor: str | None = None + plate: str | None = None + well: str | None = None + field: int | str | None = None # validate the arg types def __post_init__(self) -> None: @@ -52,12 +62,42 @@ def __post_init__(self) -> None: raise TypeError("label_key_name must be a list of strings or None") if not isinstance(self.raw_image_key_name, (list, type(None))): raise TypeError("raw_image_key_name must be a list of strings or None") + if not isinstance(self.patient_tumor, (str, type(None))): + raise TypeError("patient_tumor must be a string or None") + if not isinstance(self.plate, (str, type(None))): + raise TypeError("plate must be a string or None") + if not isinstance(self.well, (str, type(None))): + raise TypeError("well must be a string or None") + if not isinstance(self.field, (int, str, type(None))): + raise TypeError("field must be an int, str, or None") if self.label_key_name is None: self.label_key_name = [] if self.raw_image_key_name is None: self.raw_image_key_name = [] + @property + def image_id(self) -> str | None: + """Deterministic ``Metadata_Imaging_ImageID`` value, or None if unset. + + Returns ``build_image_id(...)`` when all four coordinate fields are + set, otherwise ``None`` (the loader then falls back to + ``image_set_name``). + """ + if ( + self.patient_tumor is not None + and self.plate is not None + and self.well is not None + and self.field is not None + ): + return build_image_id( + patient_tumor=self.patient_tumor, + plate=self.plate, + well=self.well, + field=self.field, + ) + return None + class _LazyImageSetDict(dict): # type: ignore[type-arg] """Dictionary that loads image arrays on first access.""" @@ -179,6 +219,13 @@ def __init__( # noqa: PLR0913 self.anisotropy_factor = self.anisotropy_spacing[0] / self.anisotropy_spacing[1] self.image_set_name = config.image_set_name self.label_set_path = label_set_path + # Deterministic imaging identifier for the warehouse join key + # (``Metadata_Imaging_ImageID``). When identifier fields are not + # provided (legacy/library use), fall back to the image set name so + # the emitted metadata column always has a value. + self.image_id = ( + config.image_id if config.image_id is not None else config.image_set_name + ) self._load_path_based_images( channel_mapping=channel_mapping, channel_tokens=channel_tokens, @@ -195,6 +242,92 @@ def __init__( # noqa: PLR0913 self.get_image_names() self.get_unique_objects_in_compartments() + @classmethod + def from_image_dict( # noqa: PLR0913 + cls, + image_dict: dict[str, numpy.ndarray], + *, + anisotropy_spacing: tuple[float, float, float], + image_set_name: str | None = None, + label_key_names: list[str] | None = None, + patient_tumor: str | None = None, + plate: str | None = None, + well: str | None = None, + field: int | str | None = None, + ) -> ImageSetLoader: + """Build an ImageSetLoader from an in-memory channel/label dict. + + Existing constructors only accept a directory glob (path-based) or a + single array (array-based). A well/FOV shard carries multiple channels + and multiple compartments as distinct arrays, so this classmethod + builds the ``image_set_dict`` directly from a pre-loaded + ``{key: ndarray}`` mapping. It formalizes the ``ImageSetLoader.__new__`` + workaround previously used in the colocalization test helper. + + Parameters + ---------- + image_dict : dict[str, numpy.ndarray] + Mapping of channel names and compartment names to their arrays. + Each array is validated through ``ImageArrayModel``. + anisotropy_spacing : tuple[float, float, float] + (z_spacing, y_spacing, x_spacing). + image_set_name : str | None + Optional image set name (emitted as ``Metadata_Experiment_ImageSet``). + label_key_names : list[str] | None + Keys in ``image_dict`` that are compartment labels (not channels). + Used by ``get_compartments`` to distinguish compartments from + raw channels. + patient_tumor, plate, well, field : optional + Imaging-coordinate identifier fields. When all four are set, the + loader's ``image_id`` is the deterministic + ``Metadata_Imaging_ImageID``; otherwise it falls back to + ``image_set_name``. + + Returns + ------- + ImageSetLoader + A fully initialized loader (compartments, image names, and unique + compartment objects populated). + + """ + self = cls.__new__(cls) + self.image_set_dict = _LazyImageSetDict() + for key, array in image_dict.items(): + # Run through pydantic validation to ensure each array is valid, + # mirroring ``_load_array_based_images``. + self.image_set_dict[key] = ImageArrayModel(array=array).array + self._label_key_names = list(label_key_names or []) + self.anisotropy_spacing = anisotropy_spacing + self.anisotropy_factor = self.anisotropy_spacing[0] / self.anisotropy_spacing[1] + self.image_set_name = image_set_name + self.label_set_path = None + config = ImageSetConfig( + image_set_name=image_set_name, + label_key_name=list(label_key_names or []), + raw_image_key_name=[ + key for key in image_dict if key not in (label_key_names or []) + ], + patient_tumor=patient_tumor, + plate=plate, + well=well, + field=field, + ) + self.image_id = ( + config.image_id if config.image_id is not None else config.image_set_name + ) + # Set compartments and image names directly from the declared label + # keys rather than calling ``get_compartments``/``get_image_names``. + # Those methods' compartment heuristic differs across repo revisions + # (it was corrected in a later bugfix commit), but this classmethod + # already knows which keys are labels, so deriving the split here keeps + # it self-contained and correct on any base. + self.compartments = list(self._label_key_names) + self.image_names = [ + key for key in image_dict if key not in self._label_key_names + ] + self.get_unique_objects_in_compartments() + return self + @staticmethod def _validate_input_sources( image_set_path: pathlib.Path | None, diff --git a/src/zedprofiler/cli.py b/src/zedprofiler/cli.py new file mode 100644 index 0000000..3b5c13a --- /dev/null +++ b/src/zedprofiler/cli.py @@ -0,0 +1,679 @@ +"""Command-line interface for per-shard feature extraction. + +``ZedProfiler run`` is the process a workflow manager (Nextflow via SLURM +``sbatch``) dispatches once per well/FOV shard. It loads one image set from +explicit file paths, runs a selected subset of featurizers, and writes one +Parquet per feature table to an output directory. + +Why argparse (not ``fire``): the repo's other CLI surfaces use ``fire``, but +this command needs repeatable flags (``--image``/``--label``) and a +three-value ``--anisotropy-spacing`` flag, which argparse handles cleanly and +``fire`` does not. The ``trigger()`` entry point name is kept for consistency +with the existing ``pyproject.toml`` console script. + +Idempotency: the same shard spec always produces the same output paths and +content. ``--skip-existing`` skips a feature request whose output Parquet +already exists, so a re-run fills only missing shards/feature tables without +redoing finished ones. Writes are atomic (temp file + ``os.replace``) so a +crashed shard never leaves a partial file that ``--skip-existing`` would +mistake for a complete one. +""" + +from __future__ import annotations + +import argparse +import sys +from collections.abc import Sequence +from pathlib import Path + +from zedprofiler.featurization.colocalization import compute_colocalization +from zedprofiler.featurization.granularity import compute_granularity +from zedprofiler.featurization.intensity import compute_intensity +from zedprofiler.featurization.neighbors import compute_neighbors +from zedprofiler.featurization.texture import compute_texture +from zedprofiler.featurization.volumesizeshape import compute_volume_size_shape +from zedprofiler.identifiers import build_image_id +from zedprofiler.IO.feature_writing_utils import ( + FeatureMetadata, + format_morphology_feature_name, + save_features_as_parquet, +) +from zedprofiler.IO.loading_classes import ( + ImageSetLoader, + ObjectLoader, + TwoObjectLoader, + _image_loading, +) + +# CPU-backed featurizers; the ``cpu_or_gpu`` component of the output filename. +_CPU_OR_GPU = "cpu" + +# Minimum number of channels required for colocalization requests. +_MIN_CHANNELS_FOR_COLOCALIZATION = 2 + +# Feature types that consume a single channel + compartment via ObjectLoader. +_SINGLE_CHANNEL_TYPES = ( + "VolumeSizeShape", + "Intensity", + "Neighbors", + "Texture", + "Granularity", +) +# Feature types that require two channels via TwoObjectLoader. +_TWO_CHANNEL_TYPES = ("Colocalization",) +ALL_FEATURE_TYPES = (*_SINGLE_CHANNEL_TYPES, *_TWO_CHANNEL_TYPES) + +# Channel-agnostic features: their computation does not use the channel image, +# but like every ZedProfiler feature they are namespaced by a channel for +# warehouse organization, so a channel (for naming only) is still required. +_CHANNEL_AGNOSTIC_TYPES = ("VolumeSizeShape", "Neighbors") + + +def _parse_name_path(token: str) -> tuple[str, Path]: + """Parse a ``NAME=PATH`` flag value into a (name, path) pair.""" + if "=" not in token: + raise argparse.ArgumentTypeError( + f"Expected NAME=PATH, got {token!r}", + ) + name, raw_path = token.split("=", 1) + name = name.strip() + if not name: + raise argparse.ArgumentTypeError(f"NAME in {token!r} is empty") + return name, Path(raw_path) + + +def _parse_feature_spec(token: str) -> dict[str, object]: + """Parse a ``TYPE[,key=value,...]`` feature request into a dict. + + The first comma-separated token is the feature type; the rest are + ``key=value`` overrides for that feature's parameters. + """ + parts = [p.strip() for p in token.split(",") if p.strip()] + if not parts: + raise argparse.ArgumentTypeError(f"Empty feature spec: {token!r}") + feature_type = parts[0] + if feature_type not in ALL_FEATURE_TYPES: + valid = ", ".join(ALL_FEATURE_TYPES) + raise argparse.ArgumentTypeError( + f"Unknown feature type {feature_type!r}; valid: {valid}", + ) + request: dict[str, object] = {"type": feature_type} + for part in parts[1:]: + if "=" not in part: + raise argparse.ArgumentTypeError( + f"Expected key=value, got {part!r} in {token!r}", + ) + key, value = part.split("=", 1) + request[key.strip()] = value.strip() + return request + + +def _output_path( + out_dir: Path, + compartment: str, + channel: str, + feature_type: str, +) -> Path: + """The Parquet path a feature request will write to. + + Mirrors the naming inside ``save_features_as_parquet`` so the CLI can check + existence (for ``--skip-existing``) before running a featurizer. + """ + prefix = format_morphology_feature_name( + compartment, + channel, + feature_type, + _CPU_OR_GPU, + ) + return out_dir / f"{prefix}_features.parquet" + + +def _coerce_param(value: object, cast: type, key: str, spec: str) -> object: + """Cast a parsed string param, raising a friendly error on failure.""" + try: + return cast(value) # type: ignore[arg-type] + except (TypeError, ValueError) as exc: + raise argparse.ArgumentTypeError( + f"Parameter {key}={value!r} in {spec!r} is not a valid {cast.__name__}", + ) from exc + + +def _resolve_channel(request: dict[str, object], spec: str) -> str: + """Return the single channel for a single-channel feature request.""" + channel = request.get("channel") + if not isinstance(channel, str) or not channel: + raise argparse.ArgumentTypeError( + f"Feature spec {spec!r} requires a 'channel' key", + ) + return channel + + +def _resolve_compartment(request: dict[str, object], spec: str) -> str: + """Return the compartment for a feature request.""" + compartment = request.get("compartment") + if not isinstance(compartment, str) or not compartment: + raise argparse.ArgumentTypeError( + f"Feature spec {spec!r} requires a 'compartment' key", + ) + return compartment + + +def _run_single_channel( + image_set_loader: ImageSetLoader, + request: dict[str, object], +) -> tuple[str, str, object]: + """Run a single-channel featurizer; return (channel, feature_type, df).""" + feature_type = str(request["type"]) + spec = ",".join(f"{k}={v}" for k, v in request.items()) + compartment = _resolve_compartment(request, spec) + channel = _resolve_channel(request, spec) + object_loader = ObjectLoader( + image_set_loader=image_set_loader, + channel_name=channel, + compartment_name=compartment, + ) + if feature_type == "VolumeSizeShape": + df = compute_volume_size_shape( + image_set_loader=image_set_loader, + object_loader=object_loader, + ) + elif feature_type == "Intensity": + df = compute_intensity(object_loader) + elif feature_type == "Neighbors": + df = compute_neighbors( + object_loader, + distance_threshold=int( + request.get("distance_threshold", 10), + ), + anisotropy_factor=float( + request.get( + "anisotropy_factor", + image_set_loader.anisotropy_factor, + ), + ), + ) + elif feature_type == "Texture": + df = compute_texture( + object_loader, + distance=int(request.get("distance", 1)), + grayscale=int(request.get("grayscale", 256)), + ) + elif feature_type == "Granularity": + df = compute_granularity( + object_loader, + radius=int(request.get("radius", 10)), + granular_spectrum_length=int( + request.get("granular_spectrum_length", 16), + ), + subsample_size=float(request.get("subsample_size", 0.25)), + image_sample_size=float(request.get("image_sample_size", 0.25)), + ) + else: # pragma: no cover - exhaustive dispatch above + raise ValueError(f"Unhandled single-channel feature type: {feature_type}") + return channel, feature_type, df + + +def _run_colocalization( + image_set_loader: ImageSetLoader, + request: dict[str, object], +) -> tuple[str, str, object]: + """Run a two-channel colocalization request; return (channel, type, df).""" + spec = ",".join(f"{k}={v}" for k, v in request.items()) + compartment = _resolve_compartment(request, spec) + channel1 = request.get("channel1") + channel2 = request.get("channel2") + if not isinstance(channel1, str) or not isinstance(channel2, str): + raise argparse.ArgumentTypeError( + f"Colocalization spec {spec!r} requires 'channel1' and 'channel2' keys", + ) + two_object_loader = TwoObjectLoader( + image_set_loader=image_set_loader, + compartment=compartment, + channel1=channel1, + channel2=channel2, + ) + df = compute_colocalization( + two_object_loader, + thr=int(request.get("thr", 15)), + fast_costes=str(request.get("fast_costes", "Accurate")), + channel1=channel1, + channel2=channel2, + ) + return f"{channel1}-{channel2}", "Colocalization", df + + +def _request_output_identity( + request: dict[str, object], +) -> tuple[str, str, str]: + """Return (compartment, channel, feature_type) for a request's output path. + + Pure (no featurizer runs): lets ``run()`` compute the target Parquet path + for ``--skip-existing`` filtering *before* reading any image from disk, so + a re-run over finished shards skips image I/O entirely. + """ + feature_type = str(request["type"]) + if feature_type == "Colocalization": + channel = f"{request['channel1']}-{request['channel2']}" + else: + channel = str(request["channel"]) + return str(request["compartment"]), channel, feature_type + + +def _execute_request( + image_set_loader: ImageSetLoader, + request: dict[str, object], + out_dir: Path, + target: Path, +) -> Path: + """Run one feature request and atomically write its Parquet. + + The caller precomputes ``target`` (via ``_output_path``) and handles + ``--skip-existing`` filtering, so this always writes. + """ + feature_type = str(request["type"]) + if feature_type == "Colocalization": + channel, ran_type, df = _run_colocalization(image_set_loader, request) + else: + channel, ran_type, df = _run_single_channel(image_set_loader, request) + compartment = str(request["compartment"]) + out_dir.mkdir(parents=True, exist_ok=True) + save_features_as_parquet( + out_dir, + df, + FeatureMetadata( + compartment=compartment, + channel=channel, + feature_type=ran_type, + cpu_or_gpu=_CPU_OR_GPU, + ), + atomic=True, + ) + print(f"wrote: {target}", file=sys.stderr) + return target + + +def _colocalization_requests( + channels: list[str], + compartments: list[str], +) -> list[dict[str, object]]: + """One Colocalization request per ordered channel pair x compartment.""" + if len(channels) < _MIN_CHANNELS_FOR_COLOCALIZATION: + raise argparse.ArgumentTypeError( + "Colocalization requires at least two --image channels", + ) + return [ + { + "type": "Colocalization", + "channel1": channels[i], + "channel2": channels[j], + "compartment": compartment, + } + for i in range(len(channels)) + for j in range(i + 1, len(channels)) + for compartment in compartments + ] + + +def _channel_agnostic_requests( + feature_type: str, + channel: str, + compartments: list[str], +) -> list[dict[str, object]]: + """One request per compartment, using ``channel`` for naming only.""" + return [ + {"type": feature_type, "channel": channel, "compartment": compartment} + for compartment in compartments + ] + + +def _per_channel_requests( + feature_type: str, + channels: list[str], + compartments: list[str], +) -> list[dict[str, object]]: + """One request per channel x compartment (Intensity/Texture/Granularity).""" + return [ + {"type": feature_type, "channel": channel, "compartment": compartment} + for channel in channels + for compartment in compartments + ] + + +def _auto_requests( + channels: list[str], + compartments: list[str], + feature_types: list[str], +) -> list[dict[str, object]]: + """Generate requests over the channel x compartment cross-product. + + Used when ``--features`` is given without explicit ``--feature`` specs. + Channel-agnostic features (VolumeSizeShape, Neighbors) use the first + channel for naming only. Colocalization runs on each ordered channel pair + x compartment. + """ + if not channels: + raise argparse.ArgumentTypeError("No --image flags provided") + if not compartments: + raise argparse.ArgumentTypeError("No --label flags provided") + first_channel = channels[0] + requests: list[dict[str, object]] = [] + for feature_type in feature_types: + if feature_type == "Colocalization": + requests.extend(_colocalization_requests(channels, compartments)) + elif feature_type in _CHANNEL_AGNOSTIC_TYPES: + requests.extend( + _channel_agnostic_requests(feature_type, first_channel, compartments), + ) + else: # Intensity, Texture, Granularity + requests.extend(_per_channel_requests(feature_type, channels, compartments)) + return requests + + +def _resolve_requests( + channels: list[str], + compartments: list[str], + feature_specs: list[dict[str, object]], + features_filter: list[str] | None, +) -> list[dict[str, object]]: + """Determine the final list of feature requests to run. + + - If explicit ``--feature`` specs are given, use them (filtered by + ``--features`` if present), validating referenced channels/compartments. + - Else auto-generate requests from ``--features`` (or all applicable types). + """ + if feature_specs: + requests = feature_specs + if features_filter: + requests = [r for r in requests if str(r["type"]) in features_filter] + _validate_request_channels_compartments(requests, channels, compartments) + return requests + feature_types = features_filter if features_filter else list(_SINGLE_CHANNEL_TYPES) + if (features_filter is None) and len(channels) >= _MIN_CHANNELS_FOR_COLOCALIZATION: + feature_types = [*_SINGLE_CHANNEL_TYPES, "Colocalization"] + return _auto_requests(channels, compartments, feature_types) + + +def _validate_request_channels_compartments( + requests: list[dict[str, object]], + channels: list[str], + compartments: list[str], +) -> None: + """Ensure each explicit request references declared channels/compartments.""" + for request in requests: + feature_type = str(request["type"]) + compartment = request.get("compartment") + if compartment not in compartments: + raise argparse.ArgumentTypeError( + f"compartment {compartment!r} in {feature_type} request is not " + f"declared via --label (valid: {compartments})", + ) + if feature_type == "Colocalization": + for key in ("channel1", "channel2"): + ch = request.get(key) + if ch not in channels: + raise argparse.ArgumentTypeError( + f"{key} {ch!r} in {feature_type} request is not " + f"declared via --image (valid: {channels})", + ) + else: + ch = request.get("channel") + if ch not in channels: + raise argparse.ArgumentTypeError( + f"channel {ch!r} in {feature_type} request is not " + f"declared via --image (valid: {channels})", + ) + + +def _build_image_set_loader( + images: list[tuple[str, Path]], + labels: list[tuple[str, Path]], + anisotropy_spacing: tuple[float, float, float], + identifiers: tuple[str, str, str, str], +) -> tuple[ImageSetLoader, list[str], list[str]]: + """Read images/labels and build an ImageSetLoader with identifiers.""" + image_dict: dict[str, object] = {} + for name, path in images: + image_dict[name] = _image_loading(path) + compartment_names: list[str] = [] + for name, path in labels: + image_dict[name] = _image_loading(path) + compartment_names.append(name) + patient_tumor, plate, well, field = identifiers + image_set_name = build_image_id_from_identifiers(identifiers) + image_set_loader = ImageSetLoader.from_image_dict( + image_dict, + anisotropy_spacing=anisotropy_spacing, + image_set_name=image_set_name, + label_key_names=compartment_names, + patient_tumor=patient_tumor, + plate=plate, + well=well, + field=field, + ) + channels = [name for name, _ in images] + return image_set_loader, channels, compartment_names + + +def build_image_id_from_identifiers( + identifiers: tuple[str, str, str, str], +) -> str: + """Build the deterministic image set name from identifier fields.""" + patient_tumor, plate, well, field = identifiers + return build_image_id(patient_tumor, plate, well, field) + + +def run( # noqa: PLR0913, PLR0917 + images: list[tuple[str, Path]], + labels: list[tuple[str, Path]], + anisotropy_spacing: tuple[float, float, float], + identifiers: tuple[str, str, str, str], + out_dir: Path, + feature_specs: list[dict[str, object]] | None = None, + features_filter: list[str] | None = None, + skip_existing: bool = False, + force: bool = False, +) -> list[Path]: + """Run feature extraction for one shard and write Parquet outputs. + + Parameters + ---------- + images : list[tuple[str, Path]] + (channel name, path) pairs from ``--image``. + labels : list[tuple[str, Path]] + (compartment name, path) pairs from ``--label``. + anisotropy_spacing : tuple[float, float, float] + (z, y, x) spacing. + identifiers : tuple[str, str, str, str] + (patient_tumor, plate, well, field) imaging coordinates. + out_dir : Path + Shard output directory. + feature_specs : list[dict[str, object]] | None + Explicit ``--feature`` requests; None means auto-generate. + features_filter : list[str] | None + ``--features`` selector restricting which feature types run. + skip_existing : bool + Skip a request whose output Parquet already exists. + force : bool + Overwrite even when the output exists (still crash-safe/atomic). + + Returns + ------- + list[Path] + Paths written (or skipped) per request. + """ + channels = [name for name, _ in images] + compartments = [name for name, _ in labels] + requests = _resolve_requests( + channels, + compartments, + feature_specs or [], + features_filter, + ) + if not requests: + print("No feature requests to run", file=sys.stderr) + return [] + # Compute target paths and partition into skipped vs. pending *before* + # reading any image. When every requested output already exists and + # --skip-existing is set, image I/O is avoided entirely. + results: list[Path] = [] + pending: list[tuple[dict[str, object], Path]] = [] + for request in requests: + compartment, channel, feature_type = _request_output_identity(request) + target = _output_path(out_dir, compartment, channel, feature_type) + if target.exists() and skip_existing and not force: + print(f"skip-existing: {target}", file=sys.stderr) + results.append(target) + else: + pending.append((request, target)) + if not pending: + print( + "all requested outputs exist; nothing to do", + file=sys.stderr, + ) + return results + image_set_loader, _, _ = _build_image_set_loader( + images, + labels, + anisotropy_spacing, + identifiers, + ) + for request, target in pending: + results.append(_execute_request(image_set_loader, request, out_dir, target)) + return results + + +def _build_parser() -> argparse.ArgumentParser: + """Construct the top-level argument parser.""" + parser = argparse.ArgumentParser( + prog="ZedProfiler", + description="Per-shard 3D featurization for the NF1 profiling warehouse.", + ) + subparsers = parser.add_subparsers(dest="command", required=True) + + run_parser = subparsers.add_parser( + "run", + help="Run feature extraction for one well/FOV shard.", + description=( + "Run feature extraction for one well/FOV shard. Loads one image " + "set from explicit file paths, runs the selected featurizers, and " + "writes one Parquet per feature table to --out-dir." + ), + ) + run_parser.add_argument( + "--image", + action="append", + required=True, + metavar="NAME=PATH", + help="Channel image as NAME=PATH (repeatable; >=1 required).", + ) + run_parser.add_argument( + "--label", + action="append", + required=True, + metavar="NAME=PATH", + help="Compartment label mask as NAME=PATH (repeatable; >=1 required).", + ) + run_parser.add_argument( + "--anisotropy-spacing", + nargs=3, + required=True, + type=float, + metavar=("Z", "Y", "X"), + help="Z, Y, X voxel spacing.", + ) + run_parser.add_argument( + "--patient-tumor", + required=True, + help="Patient-tumor identifier (e.g. NF0014_T1).", + ) + run_parser.add_argument("--plate", required=True, help="Plate identifier.") + run_parser.add_argument("--well", required=True, help="Well identifier.") + run_parser.add_argument( + "--field", + required=True, + help="Field-of-view index or identifier.", + ) + run_parser.add_argument( + "--out-dir", + required=True, + help="Shard output directory (created if needed).", + ) + run_parser.add_argument( + "--features", + default=None, + help=( + "Comma-separated feature types to run (selector). With no " + "--feature flags, runs these types over the channel x compartment " + "cross-product. With --feature flags, restricts those requests by " + "type. Default: all applicable single-channel types (plus " + "Colocalization when >=2 channels)." + ), + ) + run_parser.add_argument( + "--feature", + action="append", + default=[], + metavar="TYPE[,key=value,...]", + help=( + "Explicit feature request as TYPE[,key=value,...] (repeatable). " + "Examples: 'Intensity,channel=DNA,compartment=Nuclei'; " + "'Colocalization,channel1=DNA1,channel2=DNA2,compartment=Nuclei," + "fast_costes=Faster'." + ), + ) + run_parser.add_argument( + "--skip-existing", + action="store_true", + help="Skip a feature request whose output Parquet already exists.", + ) + run_parser.add_argument( + "--force", + action="store_true", + help="Overwrite even when the output exists (writes are still atomic).", + ) + return parser + + +def main(argv: Sequence[str] | None = None) -> int: + """CLI entry point. Returns a process exit code.""" + parser = _build_parser() + args = parser.parse_args(argv) + if args.command != "run": + parser.error("a subcommand is required") + return 2 # pragma: no cover - parser.error exits + + images = [_parse_name_path(token) for token in args.image] + labels = [_parse_name_path(token) for token in args.label] + feature_specs = [_parse_feature_spec(token) for token in args.feature] + features_filter = ( + [f.strip() for f in args.features.split(",") if f.strip()] + if args.features + else None + ) + identifiers = ( + args.patient_tumor, + args.plate, + args.well, + args.field, + ) + run( + images=images, + labels=labels, + anisotropy_spacing=tuple(args.anisotropy_spacing), + identifiers=identifiers, + out_dir=Path(args.out_dir), + feature_specs=feature_specs, + features_filter=features_filter, + skip_existing=args.skip_existing, + force=args.force, + ) + return 0 + + +def trigger() -> None: + """Console-script entry point (matches the pyproject ``scripts`` target).""" + raise SystemExit(main()) + + +if __name__ == "__main__": # pragma: no cover + trigger() diff --git a/src/zedprofiler/featurization/colocalization.py b/src/zedprofiler/featurization/colocalization.py index 02f1bcb..5a20869 100644 --- a/src/zedprofiler/featurization/colocalization.py +++ b/src/zedprofiler/featurization/colocalization.py @@ -627,6 +627,7 @@ def compute_colocalization( # noqa: C901, PLR0912 if full_name not in ( "Metadata_Object_ObjectID", "Metadata_Experiment_ImageSet", + "Metadata_Imaging_ImageID", ): try: row[full_name] = numpy.float32(meas_val) @@ -635,11 +636,12 @@ def compute_colocalization( # noqa: C901, PLR0912 else: row[full_name] = meas_val - # ensure object_id and image_set are present and first + # ensure object_id and image identifiers are present and first row["Metadata_Object_ObjectID"] = object_id row["Metadata_Experiment_ImageSet"] = ( two_object_loader.image_set_loader.image_set_name ) + row["Metadata_Imaging_ImageID"] = two_object_loader.image_set_loader.image_id list_of_dfs.append(row) # Convert list of row-dicts into a dict-of-lists with stable ordering @@ -650,7 +652,11 @@ def compute_colocalization( # noqa: C901, PLR0912 other_keys: list[str] = [] for d in list_of_dfs: for k in d: - if k in ("Metadata_Object_ObjectID", "Metadata_Experiment_ImageSet"): + if k in ( + "Metadata_Object_ObjectID", + "Metadata_Experiment_ImageSet", + "Metadata_Imaging_ImageID", + ): continue if k not in other_keys: other_keys.append(k) @@ -658,6 +664,7 @@ def compute_colocalization( # noqa: C901, PLR0912 all_keys = [ "Metadata_Object_ObjectID", "Metadata_Experiment_ImageSet", + "Metadata_Imaging_ImageID", *other_keys, ] result: dict[str, list[object]] = { diff --git a/src/zedprofiler/featurization/granularity.py b/src/zedprofiler/featurization/granularity.py index 72c1f78..d479b29 100644 --- a/src/zedprofiler/featurization/granularity.py +++ b/src/zedprofiler/featurization/granularity.py @@ -459,6 +459,11 @@ def compute_granularity( # noqa: C901, PLR0912, PLR0913, PLR0915 }, inplace=True, ) + final_df.insert( + 0, + "Metadata_Imaging_ImageID", + object_loader.image_set_loader.image_id, + ) final_df.insert( 0, "Metadata_Experiment_ImageSet", diff --git a/src/zedprofiler/featurization/intensity.py b/src/zedprofiler/featurization/intensity.py index d153e2f..604faae 100644 --- a/src/zedprofiler/featurization/intensity.py +++ b/src/zedprofiler/featurization/intensity.py @@ -227,6 +227,11 @@ def compute_intensity( # noqa: PLR0915 inplace=True, ) + final_df.insert( + 0, + "Metadata_Imaging_ImageID", + object_loader.image_set_loader.image_id, + ) final_df.insert( 0, "Metadata_Experiment_ImageSet", diff --git a/src/zedprofiler/featurization/neighbors.py b/src/zedprofiler/featurization/neighbors.py index 9ab04ca..9c78059 100644 --- a/src/zedprofiler/featurization/neighbors.py +++ b/src/zedprofiler/featurization/neighbors.py @@ -199,6 +199,11 @@ def compute_neighbors( inplace=True, ) if not final_df.empty: + final_df.insert( + 0, + "Metadata_Imaging_ImageID", + object_loader.image_set_loader.image_id, + ) final_df.insert( 0, "Metadata_Experiment_ImageSet", diff --git a/src/zedprofiler/featurization/texture.py b/src/zedprofiler/featurization/texture.py index 71f83fd..30813b6 100644 --- a/src/zedprofiler/featurization/texture.py +++ b/src/zedprofiler/featurization/texture.py @@ -205,6 +205,11 @@ def compute_texture( # noqa: C901 }, inplace=True, ) + final_df.insert( + 0, + "Metadata_Imaging_ImageID", + object_loader.image_set_loader.image_id, + ) final_df.insert( 0, "Metadata_Experiment_ImageSet", diff --git a/src/zedprofiler/featurization/volumesizeshape.py b/src/zedprofiler/featurization/volumesizeshape.py index 123b5be..6e44b36 100644 --- a/src/zedprofiler/featurization/volumesizeshape.py +++ b/src/zedprofiler/featurization/volumesizeshape.py @@ -206,6 +206,11 @@ def measure_3D_volume_size_shape( inplace=True, ) + final_df.insert( + 0, + "Metadata_Imaging_ImageID", + object_loader.image_set_loader.image_id, + ) final_df.insert( 1, "Metadata_Experiment_ImageSet", diff --git a/src/zedprofiler/identifiers.py b/src/zedprofiler/identifiers.py new file mode 100644 index 0000000..0cf8179 --- /dev/null +++ b/src/zedprofiler/identifiers.py @@ -0,0 +1,55 @@ +"""Deterministic imaging identifiers for warehouse join keys. + +The NF1 bioimage profiling warehouse joins every image, object, feature +table, and annotation via stable identifiers (see the future processing +plan's identifier spec). The central one is ``Metadata_Imaging_ImageID``, +built deterministically from the four imaging coordinates: + + patient-tumor, plate, well, field + +Because a shard is dispatched per well/FOV, every feature table emitted by a +shard carries this single image id so downstream tables can rejoin without a +database service. + +This module is the single source of truth for the id format. Changing the +format here changes every shard's output ids at once. +""" + +from __future__ import annotations + +from beartype import beartype + + +@beartype +def build_image_id( + patient_tumor: str, + plate: str, + well: str, + field: int | str, +) -> str: + """Build a deterministic ``Metadata_Imaging_ImageID`` value. + + The id is a stable string assembled from the four imaging coordinates so + that the same well/FOV always produces the same id across runs, batches, + and reprocessing. Component order is fixed + (patient-tumor, plate, well, field) so ids sort and group naturally by + patient then plate then well then field. + + Parameters + ---------- + patient_tumor : str + Patient-tumor identifier (e.g. ``"NF0014_T1"``). + plate : str + Plate identifier (e.g. ``"PLATE01"``). + well : str + Well identifier (e.g. ``"A1"``). + field : int | str + Field-of-view index or identifier (e.g. ``1`` or ``"f1"``). + + Returns + ------- + str + The deterministic image id, e.g. ``"NF0014_T1_PLATE01_A1_field1"``. + + """ + return f"{patient_tumor}_{plate}_{well}_field{field}" diff --git a/tests/IO/test_loading_classes.py b/tests/IO/test_loading_classes.py index 725f5b0..ccda2d6 100644 --- a/tests/IO/test_loading_classes.py +++ b/tests/IO/test_loading_classes.py @@ -383,3 +383,123 @@ def test_two_object_loader_loads_images_and_ids(self) -> None: assert np.array_equal(obj.image2, image) assert np.array_equal(obj.label_image, label_image) assert obj.object_ids == [ONE_LABEL, TWO_LABEL] + + +class TestImageSetConfigImageId: + """Tests for ImageSetConfig.image_id identifier propagation.""" + + def test_image_id_is_none_when_any_field_missing(self) -> None: + """image_id stays None until all four identifier fields are set.""" + assert ImageSetConfig().image_id is None + assert ( + ImageSetConfig( + patient_tumor="NF0014_T1", + plate="PLATE01", + well="A1", + field=None, + ).image_id + is None + ) + + def test_image_id_built_when_all_fields_set(self) -> None: + """With all four fields set, image_id is the deterministic image id.""" + config = ImageSetConfig( + patient_tumor="NF0014_T1", + plate="PLATE01", + well="A1", + field=1, + ) + assert config.image_id == "NF0014_T1_PLATE01_A1_field1" + + +class TestFromImageDict: + """Tests for ImageSetLoader.from_image_dict (multi-channel shard loader).""" + + def _build_dict(self) -> dict[str, np.ndarray]: + label = np.array( + [ + [[ZERO_LABEL, ONE_LABEL], [TWO_LABEL, TWO_LABEL]], + [[ZERO_LABEL, ONE_LABEL], [TWO_LABEL, TWO_LABEL]], + ], + dtype=np.int32, + ) + return { + "DNA": np.ones((2, 2, 2), dtype=np.float32), + "AGP": np.full((2, 2, 2), 5.0, dtype=np.float32), + "Nuclei": label, + } + + def test_builds_working_multi_channel_loader(self) -> None: + """from_image_dict resolves compartments and channels from the dict.""" + loader = ImageSetLoader.from_image_dict( + self._build_dict(), + anisotropy_spacing=(2.0, 1.0, 1.0), + image_set_name="shard-01", + label_key_names=["Nuclei"], + ) + assert loader.image_set_name == "shard-01" + assert loader.anisotropy_factor == EXPECTED_ANISOTROPY + assert loader.compartments == ["Nuclei"] + assert sorted(loader.image_names) == ["AGP", "DNA"] + assert loader.unique_compartment_objects["Nuclei"] == [ + ONE_LABEL, + TWO_LABEL, + ] + # ObjectLoader built on the resulting loader resolves both channels. + dna_obj = ObjectLoader( + image_set_loader=loader, + channel_name="DNA", + compartment_name="Nuclei", + ) + assert dna_obj.object_ids == [ONE_LABEL, TWO_LABEL] + assert np.array_equal( + dna_obj.image, + np.ones((2, 2, 2), dtype=np.float32), + ) + + def test_image_id_populated_when_identifiers_provided(self) -> None: + """Passing all four identifier fields populates image_id deterministically.""" + loader = ImageSetLoader.from_image_dict( + self._build_dict(), + anisotropy_spacing=(1.0, 1.0, 1.0), + image_set_name="shard-01", + label_key_names=["Nuclei"], + patient_tumor="NF0014_T1", + plate="PLATE01", + well="A1", + field=1, + ) + assert loader.image_id == "NF0014_T1_PLATE01_A1_field1" + + def test_image_id_falls_back_to_image_set_name_without_identifiers( + self, + ) -> None: + """Without identifiers, image_id falls back to the image set name.""" + loader = ImageSetLoader.from_image_dict( + self._build_dict(), + anisotropy_spacing=(1.0, 1.0, 1.0), + image_set_name="legacy-set", + label_key_names=["Nuclei"], + ) + assert loader.image_id == "legacy-set" + + def test_two_object_loader_on_from_image_dict(self) -> None: + """A TwoObjectLoader resolves two channels from a from_image_dict loader.""" + loader = ImageSetLoader.from_image_dict( + self._build_dict(), + anisotropy_spacing=(1.0, 1.0, 1.0), + image_set_name="shard-02", + label_key_names=["Nuclei"], + patient_tumor="X", + plate="P", + well="A1", + field="1", + ) + two = TwoObjectLoader( + image_set_loader=loader, + compartment="Nuclei", + channel1="DNA", + channel2="AGP", + ) + assert two.object_ids == [ONE_LABEL, TWO_LABEL] + assert two.image_set_loader.image_id == "X_P_A1_field1" diff --git a/tests/featurization/test_colocalization.py b/tests/featurization/test_colocalization.py index 7b34252..94a4c1f 100644 --- a/tests/featurization/test_colocalization.py +++ b/tests/featurization/test_colocalization.py @@ -22,6 +22,8 @@ class ImageSetLoaderModel(BaseModel): model_config = ConfigDict(arbitrary_types_allowed=True) image_set_name: str = "coloc" + # mirrors ImageSetLoader.image_id (falls back to image_set_name) + image_id: str = "coloc" class TwoObjectLoaderModel(BaseModel): diff --git a/tests/featurization/test_granularity.py b/tests/featurization/test_granularity.py index c2f13b8..acb0769 100644 --- a/tests/featurization/test_granularity.py +++ b/tests/featurization/test_granularity.py @@ -20,6 +20,8 @@ class ImageSetLoaderModel(BaseModel): model_config = ConfigDict(arbitrary_types_allowed=True) image_set_name: str = "gran" + # mirrors ImageSetLoader.image_id (falls back to image_set_name) + image_id: str = "gran" class ObjectLoaderModel(BaseModel): @@ -94,7 +96,7 @@ class Dummy: image = img label_image = lab object_ids: ClassVar[list[int]] = [1] - image_set_loader = type("ISL", (), {"image_set_name": "s"})() + image_set_loader = type("ISL", (), {"image_set_name": "s", "image_id": "s"})() compartment = "Cell" channel = "Ch1" @@ -120,7 +122,7 @@ class Dummy: image = img label_image = lab object_ids: ClassVar[list[int]] = [1] - image_set_loader = type("ISL", (), {"image_set_name": "s"})() + image_set_loader = type("ISL", (), {"image_set_name": "s", "image_id": "s"})() compartment = "Cell" channel = "Ch1" @@ -150,7 +152,7 @@ class Dummy: image = img label_image = lab object_ids: ClassVar[list[int]] = [1] - image_set_loader = type("ISL", (), {"image_set_name": "s"})() + image_set_loader = type("ISL", (), {"image_set_name": "s", "image_id": "s"})() compartment = "Cell" channel = "Ch1" @@ -195,7 +197,7 @@ class Dummy: image = img label_image = lab object_ids: ClassVar[list[int]] = [1] - image_set_loader = type("ISL", (), {"image_set_name": "s"})() + image_set_loader = type("ISL", (), {"image_set_name": "s", "image_id": "s"})() compartment = "Cell" channel = "Ch1" @@ -222,7 +224,7 @@ class Dummy: image = img label_image = lab object_ids: ClassVar[list[int]] = [257, 514] - image_set_loader = type("ISL", (), {"image_set_name": "s"})() + image_set_loader = type("ISL", (), {"image_set_name": "s", "image_id": "s"})() compartment = "Cell" channel = "Ch1" diff --git a/tests/featurization/test_intensity.py b/tests/featurization/test_intensity.py index a2593ea..26e9310 100644 --- a/tests/featurization/test_intensity.py +++ b/tests/featurization/test_intensity.py @@ -12,6 +12,8 @@ class ImageSetLoaderModel(BaseModel): model_config = ConfigDict(arbitrary_types_allowed=True) image_set_name: str = "intensity" + # mirrors ImageSetLoader.image_id (falls back to image_set_name) + image_id: str = "intensity" class ObjectLoaderModel(BaseModel): diff --git a/tests/featurization/test_neighbors.py b/tests/featurization/test_neighbors.py index 398d12f..5167287 100644 --- a/tests/featurization/test_neighbors.py +++ b/tests/featurization/test_neighbors.py @@ -23,6 +23,8 @@ class ImageSetLoaderModel(BaseModel): model_config = ConfigDict(arbitrary_types_allowed=True) image_set_name: str = "neighbors" + # mirrors ImageSetLoader.image_id (falls back to image_set_name) + image_id: str = "neighbors" class ObjectLoaderModel(BaseModel): diff --git a/tests/featurization/test_neighbors_additional.py b/tests/featurization/test_neighbors_additional.py index 2da68d7..fa2c0c1 100644 --- a/tests/featurization/test_neighbors_additional.py +++ b/tests/featurization/test_neighbors_additional.py @@ -50,7 +50,7 @@ def test_compute_neighbors_distance_counts() -> None: class Dummy: label_image = lab object_ids = (1, 2, 3) - image_set_loader = type("ISL", (), {"image_set_name": "s"})() + image_set_loader = type("ISL", (), {"image_set_name": "s", "image_id": "s"})() compartment = "Cell" channel = "Ch1" diff --git a/tests/featurization/test_real_world_data.py b/tests/featurization/test_real_world_data.py index 8625e65..9839cef 100644 --- a/tests/featurization/test_real_world_data.py +++ b/tests/featurization/test_real_world_data.py @@ -279,21 +279,20 @@ def _load_colocalization_case( colocalization_case: RealColocalizationCase, ) -> TwoObjectLoader: label = tifffile.imread(colocalization_case.label_image_case.label_path) - object_ids = [int(x) for x in np.unique(label) if x != 0] - image_set_loader = ImageSetLoader.__new__(ImageSetLoader) - image_set_loader.image_set_name = colocalization_case.image_set_name - image_set_loader.image_set_dict = { - colocalization_case.first_channel: tifffile.imread( - colocalization_case.first_image_case.image_path, - ), - colocalization_case.second_channel: tifffile.imread( - colocalization_case.second_image_case.image_path, - ), - colocalization_case.compartment: label, - } - image_set_loader.unique_compartment_objects = { - colocalization_case.compartment: object_ids, - } + image_set_loader = ImageSetLoader.from_image_dict( + { + colocalization_case.first_channel: tifffile.imread( + colocalization_case.first_image_case.image_path, + ), + colocalization_case.second_channel: tifffile.imread( + colocalization_case.second_image_case.image_path, + ), + colocalization_case.compartment: label, + }, + anisotropy_spacing=(1.0, 1.0, 1.0), + image_set_name=colocalization_case.image_set_name, + label_key_names=[colocalization_case.compartment], + ) return TwoObjectLoader( image_set_loader=image_set_loader, diff --git a/tests/featurization/test_texture.py b/tests/featurization/test_texture.py index 659eb92..7f9c02e 100644 --- a/tests/featurization/test_texture.py +++ b/tests/featurization/test_texture.py @@ -18,6 +18,8 @@ class ImageSetLoaderModel(BaseModel): model_config = ConfigDict(arbitrary_types_allowed=True) image_set_name: str = "texture" + # mirrors ImageSetLoader.image_id (falls back to image_set_name) + image_id: str = "texture" class ObjectLoaderModel(BaseModel): diff --git a/tests/featurization/test_volumesizeshape.py b/tests/featurization/test_volumesizeshape.py index 02637f0..985959d 100644 --- a/tests/featurization/test_volumesizeshape.py +++ b/tests/featurization/test_volumesizeshape.py @@ -15,6 +15,8 @@ class ImageSetLoaderModel(BaseModel): model_config = ConfigDict(arbitrary_types_allowed=True) anisotropy_spacing: tuple[float, float, float] image_set_name: str = "testset" + # mirrors ImageSetLoader.image_id (falls back to image_set_name) + image_id: str = "testset" class ObjectLoaderModel(BaseModel): diff --git a/tests/test_cli.py b/tests/test_cli.py index 4f9aeed..92d20d6 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -1,8 +1,365 @@ -"""Package namespace tests replacing obsolete template CLI checks.""" +"""End-to-end and unit tests for the per-shard CLI (``zedprofiler.cli``). +The end-to-end tests exercise ``ZedProfiler run`` against the CellProfiler 3D +tutorial data (the same fixtures used by ``test_real_world_data.py``) and assert +the things that matter to the NF1 pipeline: one Parquet per requested feature +table, a deterministic ``Metadata_Imaging_ImageID``, ``--features`` selection, +``--skip-existing`` idempotency, and two-channel colocalization. + +The unit tests cover the request-selection and path-naming helpers directly so +the full default feature matrix (which includes slow Texture/Granularity +runs) does not have to run end-to-end here. +""" + +from __future__ import annotations + +import argparse +from pathlib import Path + +import pandas as pd +import pytest + +from zedprofiler.cli import ( + _auto_requests, + _output_path, + _parse_feature_spec, + _parse_name_path, + _resolve_requests, + main, +) from zedprofiler.featurization import texture +from zedprofiler.identifiers import build_image_id + +tifffile = pytest.importorskip("tifffile") + +TUTORIAL_ROOT = ( + Path(__file__).resolve().parent + / "data" + / "CP_tutorial_3D_noise_nuclei_segmentation" +) +IMAGE1 = TUTORIAL_ROOT / "input" / "nuclei1_out_c00_dr90_image.tif" +IMAGE2 = TUTORIAL_ROOT / "input" / "nuclei2_out_c90_dr90_image.tif" +LABEL1 = ( + TUTORIAL_ROOT + / "output" + / "masks" + / "nuclei1_out_c00_dr90_imageSegmentationMask.tiff" +) + +EXPECTED_OBJECT_COUNT = 5 +EXPECTED_CHANNEL_COUNT = 2 +PATIENT_TUMOR, PLATE, WELL, FIELD = "NF0014_T1", "PLATE01", "A1", "1" +EXPECTED_IMAGE_ID = build_image_id(PATIENT_TUMOR, PLATE, WELL, FIELD) def test_feature_namespace_import() -> None: """Lower-level feature namespace remains importable.""" assert texture.__name__ == "zedprofiler.featurization.texture" + + +def _run(argv: list[str]) -> int: + return main(argv) + + +# --------------------------------------------------------------------------- +# Unit tests: name/path + feature-spec parsing +# --------------------------------------------------------------------------- + + +def test_parse_name_path_splits_on_first_equals() -> None: + """NAME=PATH splits on the first '=' so paths may contain '='.""" + name, path = _parse_name_path("DNA=/tmp/a=b.tif") + assert name == "DNA" + assert path == Path("/tmp/a=b.tif") + + +def test_parse_name_path_rejects_missing_equals() -> None: + """A token without '=' is a user error, not a silent default.""" + with pytest.raises(argparse.ArgumentTypeError): + _parse_name_path("DNA.tif") + + +def test_parse_feature_spec_parses_type_and_overrides() -> None: + """TYPE is the first comma token; the rest are key=value overrides.""" + request = _parse_feature_spec("Intensity,channel=DNA,compartment=Nuclei") + assert request == { + "type": "Intensity", + "channel": "DNA", + "compartment": "Nuclei", + } + + +def test_parse_feature_spec_rejects_unknown_type() -> None: + """An unknown feature type raises a friendly argument error.""" + with pytest.raises(argparse.ArgumentTypeError, match="Unknown feature type"): + _parse_feature_spec("NotAFeature,channel=DNA") + + +def test_parse_feature_spec_rejects_bad_override() -> None: + """An override without '=' raises a friendly argument error.""" + with pytest.raises(argparse.ArgumentTypeError, match="key=value"): + _parse_feature_spec("Intensity,bogus") + + +def test_output_path_mirrors_save_features_as_parquet_naming() -> None: + """_output_path must match the path save_features_as_parquet writes to.""" + assert _output_path(Path("/out"), "Nuclei", "DNA", "Intensity") == Path( + "/out/Nuclei_DNA_Intensity_cpu_features.parquet" + ) + assert _output_path(Path("/out"), "Nuclei", "DNA1-DNA2", "Colocalization") == Path( + "/out/Nuclei_DNA1-DNA2_Colocalization_cpu_features.parquet" + ) + + +# --------------------------------------------------------------------------- +# Unit tests: request selection (no featurizers run) +# --------------------------------------------------------------------------- + + +def test_auto_requests_single_channel_no_colocalization() -> None: + """One channel x one compartment yields the 5 single-channel types only.""" + requests = _auto_requests( + ["DNA"], + ["Nuclei"], + ["VolumeSizeShape", "Intensity", "Neighbors", "Texture", "Granularity"], + ) + types = sorted(str(r["type"]) for r in requests) + assert types == sorted( + ["VolumeSizeShape", "Intensity", "Neighbors", "Texture", "Granularity"], + ) + # Channel-agnostic features use the first (only) channel for naming. + vol = next(r for r in requests if r["type"] == "VolumeSizeShape") + assert vol["channel"] == "DNA" + + +def test_auto_requests_two_channels_adds_colocalization_pairs() -> None: + """Two channels produce a single ordered colocalization pair x compartment.""" + requests = _auto_requests( + ["DNA1", "DNA2"], + ["Nuclei"], + ["Intensity", "Colocalization"], + ) + coloc = [r for r in requests if r["type"] == "Colocalization"] + assert len(coloc) == 1 + assert coloc[0]["channel1"] == "DNA1" + assert coloc[0]["channel2"] == "DNA2" + # Intensity runs per channel x compartment. + assert ( + sum(1 for r in requests if r["type"] == "Intensity") == EXPECTED_CHANNEL_COUNT + ) + + +def test_auto_requests_colocalization_requires_two_channels() -> None: + """Requesting Colocalization with one channel is a user error.""" + with pytest.raises(argparse.ArgumentTypeError, match="at least two"): + _auto_requests(["DNA"], ["Nuclei"], ["Colocalization"]) + + +def test_resolve_requests_default_two_channels_includes_colocalization() -> None: + """Default selection (no --features, no --feature) adds coloc for >=2 channels.""" + requests = _resolve_requests(["DNA1", "DNA2"], ["Nuclei"], [], None) + types = {str(r["type"]) for r in requests} + assert "Colocalization" in types + assert "Intensity" in types + + +def test_resolve_requests_features_filter_restricts_types() -> None: + """--features Intensity selects only Intensity from the cross-product.""" + requests = _resolve_requests( + ["DNA1", "DNA2"], + ["Nuclei"], + [], + ["Intensity"], + ) + assert all(r["type"] == "Intensity" for r in requests) + assert len(requests) == EXPECTED_CHANNEL_COUNT # 2 channels x 1 compartment + + +def test_resolve_requests_explicit_specs_validated_against_declared() -> None: + """An explicit --feature referencing an undeclared channel is rejected.""" + specs = [_parse_feature_spec("Intensity,channel=Ghost,compartment=Nuclei")] + with pytest.raises(argparse.ArgumentTypeError, match="declared via --image"): + _resolve_requests(["DNA"], ["Nuclei"], specs, None) + + +def test_resolve_requests_features_filter_applied_to_explicit_specs() -> None: + """--features restricts explicit --feature requests by type.""" + specs = [ + _parse_feature_spec("Intensity,channel=DNA,compartment=Nuclei"), + _parse_feature_spec("Texture,channel=DNA,compartment=Nuclei"), + ] + requests = _resolve_requests(["DNA"], ["Nuclei"], specs, ["Intensity"]) + assert len(requests) == 1 + assert requests[0]["type"] == "Intensity" + + +# --------------------------------------------------------------------------- +# End-to-end tests on the CellProfiler 3D tutorial data +# --------------------------------------------------------------------------- + +# The end-to-end tests read the CellProfiler 3D tutorial images/masks, which +# are added by a separate data commit and may be absent on some branches. Skip +# them when the data is not present so the CLI test module stays green +# everywhere; they run in full wherever the tutorial data is available. +requires_tutorial_data = pytest.mark.skipif( + not TUTORIAL_ROOT.exists(), + reason="CellProfiler 3D tutorial data not present on this branch", +) + + +def _base_run_args(out_dir: Path, *extra: str) -> list[str]: + return [ + "run", + f"--image=DNA={IMAGE1}", + f"--label=Nuclei={LABEL1}", + "--anisotropy-spacing", + "1.0", + "1.0", + "1.0", + f"--patient-tumor={PATIENT_TUMOR}", + f"--plate={PLATE}", + f"--well={WELL}", + f"--field={FIELD}", + f"--out-dir={out_dir}", + *extra, + ] + + +@requires_tutorial_data +def test_cli_run_intensity_writes_parquet_with_image_id( + tmp_path: Path, +) -> None: + """A single Intensity request writes one Parquet carrying the image id.""" + out_dir = tmp_path / "shard" + assert _run(_base_run_args(out_dir, "--features=Intensity")) == 0 + + parquet = _output_path(out_dir, "Nuclei", "DNA", "Intensity") + assert parquet.exists() + + df = pd.read_parquet(parquet) + assert len(df) == EXPECTED_OBJECT_COUNT + assert "Metadata_Imaging_ImageID" in df.columns + assert "Metadata_Experiment_ImageSet" in df.columns + assert set(df["Metadata_Imaging_ImageID"]) == {EXPECTED_IMAGE_ID} + # No leftover temp files from the atomic write. + assert not any(p.suffix == ".tmp" for p in out_dir.iterdir()) + + +@requires_tutorial_data +def test_cli_features_selector_restricts_outputs(tmp_path: Path) -> None: + """--features controls exactly which feature tables are written.""" + out_dir = tmp_path / "shard" + assert _run(_base_run_args(out_dir, "--features=VolumeSizeShape,Intensity")) == 0 + files = sorted(p.name for p in out_dir.glob("*.parquet")) + assert files == [ + "Nuclei_DNA_Intensity_cpu_features.parquet", + "Nuclei_DNA_VolumeSizeShape_cpu_features.parquet", + ] + + out_dir2 = tmp_path / "shard2" + assert _run(_base_run_args(out_dir2, "--features=Intensity")) == 0 + assert [p.name for p in out_dir2.glob("*.parquet")] == [ + "Nuclei_DNA_Intensity_cpu_features.parquet", + ] + + +@requires_tutorial_data +def test_cli_rerun_is_content_identical(tmp_path: Path) -> None: + """Re-running without --skip-existing reproduces identical feature content.""" + out_dir = tmp_path / "shard" + _run(_base_run_args(out_dir, "--features=Intensity")) + first = pd.read_parquet(_output_path(out_dir, "Nuclei", "DNA", "Intensity")) + + out_dir2 = tmp_path / "shard2" + _run(_base_run_args(out_dir2, "--features=Intensity")) + second = pd.read_parquet(_output_path(out_dir2, "Nuclei", "DNA", "Intensity")) + + pd.testing.assert_frame_equal(first, second) + + +@requires_tutorial_data +def test_cli_skip_existing_skips_recompute(tmp_path: Path) -> None: + """--skip-existing leaves finished outputs untouched and skips image I/O.""" + out_dir = tmp_path / "shard" + _run(_base_run_args(out_dir, "--features=Intensity")) + target = _output_path(out_dir, "Nuclei", "DNA", "Intensity") + first_bytes = target.read_bytes() + first_mtime_ns = target.stat().st_mtime_ns + + # Re-run with --skip-existing using *nonexistent* image paths: if the CLI + # tried to load images it would fail, proving skip happens before I/O. + assert ( + _run( + [ + "run", + "--image=DNA=/does/not/exist.tif", + "--label=Nuclei=/does/not/exist.tiff", + "--anisotropy-spacing", + "1.0", + "1.0", + "1.0", + f"--patient-tumor={PATIENT_TUMOR}", + f"--plate={PLATE}", + f"--well={WELL}", + f"--field={FIELD}", + f"--out-dir={out_dir}", + "--features=Intensity", + "--skip-existing", + ], + ) + == 0 + ) + assert target.read_bytes() == first_bytes + assert target.stat().st_mtime_ns == first_mtime_ns + + +@requires_tutorial_data +def test_cli_colocalization_two_channels(tmp_path: Path) -> None: + """An explicit colocalization request writes a DNA1-DNA2 Parquet.""" + out_dir = tmp_path / "shard" + assert ( + _run( + [ + "run", + f"--image=DNA1={IMAGE1}", + f"--image=DNA2={IMAGE2}", + f"--label=Nuclei={LABEL1}", + "--anisotropy-spacing", + "1.0", + "1.0", + "1.0", + f"--patient-tumor={PATIENT_TUMOR}", + f"--plate={PLATE}", + f"--well={WELL}", + f"--field={FIELD}", + f"--out-dir={out_dir}", + "--feature=Colocalization,channel1=DNA1,channel2=DNA2,compartment=Nuclei,fast_costes=Faster", + ], + ) + == 0 + ) + target = _output_path(out_dir, "Nuclei", "DNA1-DNA2", "Colocalization") + assert target.exists() + df = pd.read_parquet(target) + assert len(df) == EXPECTED_OBJECT_COUNT + assert set(df["Metadata_Imaging_ImageID"]) == {EXPECTED_IMAGE_ID} + assert any("Colocalization" in c for c in df.columns) + + +def test_cli_missing_required_arg_errors(tmp_path: Path) -> None: + """A run without --out-dir exits non-zero (argparse error).""" + argv = [ + "run", + f"--image=DNA={IMAGE1}", + f"--label=Nuclei={LABEL1}", + "--anisotropy-spacing", + "1.0", + "1.0", + "1.0", + f"--patient-tumor={PATIENT_TUMOR}", + f"--plate={PLATE}", + f"--well={WELL}", + f"--field={FIELD}", + ] + with pytest.raises(SystemExit): + _run(argv) diff --git a/tests/test_identifiers.py b/tests/test_identifiers.py new file mode 100644 index 0000000..363cefd --- /dev/null +++ b/tests/test_identifiers.py @@ -0,0 +1,44 @@ +"""Tests for the deterministic Metadata_Imaging_ImageID builder.""" + +from __future__ import annotations + +import pytest + +from zedprofiler.identifiers import build_image_id + + +def test_build_image_id_formats_all_fields() -> None: + """The image id joins patient-tumor, plate, well, and a fielded suffix.""" + assert ( + build_image_id("NF0014_T1", "PLATE01", "A1", 1) == "NF0014_T1_PLATE01_A1_field1" + ) + + +def test_build_image_id_is_deterministic() -> None: + """Repeated calls with the same inputs produce the same id.""" + args = ("NF0014_T1", "PLATE01", "A1", "2") + assert build_image_id(*args) == build_image_id(*args) + + +def test_build_image_id_field_accepts_string_or_int() -> None: + """Field may be supplied as either an int or its string form.""" + assert build_image_id("X", "P", "A1", 3) == "X_P_A1_field3" + assert build_image_id("X", "P", "A1", "3") == "X_P_A1_field3" + + +@pytest.mark.parametrize( + "patient_tumor,plate,well,field,expected", + [ + ("NF0014_T1", "PLATE01", "A1", "1", "NF0014_T1_PLATE01_A1_field1"), + ("NF0009_T2", "PLATE02", "B3", "7", "NF0009_T2_PLATE02_B3_field7"), + ], +) +def test_build_image_id_parametrized( + patient_tumor: str, + plate: str, + well: str, + field: str, + expected: str, +) -> None: + """Format holds across distinct imaging coordinates.""" + assert build_image_id(patient_tumor, plate, well, field) == expected From f4563607f6b2986b9a4d1b58f0fe9be486105190 Mon Sep 17 00:00:00 2001 From: d33bs Date: Thu, 6 Aug 2026 09:11:58 -0600 Subject: [PATCH 3/7] test coverage --- src/zedprofiler/cli.py | 2 +- tests/IO/test_loading_classes.py | 4 + tests/test_cli.py | 173 +++++++++++++++++++++++++++++++ 3 files changed, 178 insertions(+), 1 deletion(-) diff --git a/src/zedprofiler/cli.py b/src/zedprofiler/cli.py index 3b5c13a..1861b19 100644 --- a/src/zedprofiler/cli.py +++ b/src/zedprofiler/cli.py @@ -638,7 +638,7 @@ def main(argv: Sequence[str] | None = None) -> int: """CLI entry point. Returns a process exit code.""" parser = _build_parser() args = parser.parse_args(argv) - if args.command != "run": + if args.command != "run": # pragma: no cover - subparsers required=True parser.error("a subcommand is required") return 2 # pragma: no cover - parser.error exits diff --git a/tests/IO/test_loading_classes.py b/tests/IO/test_loading_classes.py index ccda2d6..258b4ee 100644 --- a/tests/IO/test_loading_classes.py +++ b/tests/IO/test_loading_classes.py @@ -63,6 +63,10 @@ def test_config_post_init_none_defaults(self) -> None: {"raw_image_key_name": "raw"}, "raw_image_key_name must be a list of strings or None", ), + ({"patient_tumor": 123}, "patient_tumor must be a string or None"), + ({"plate": 123}, "plate must be a string or None"), + ({"well": 123}, "well must be a string or None"), + ({"field": 1.5}, "field must be an int, str, or None"), ], ) def test_config_rejects_invalid_types( diff --git a/tests/test_cli.py b/tests/test_cli.py index 92d20d6..ece3ba4 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -16,19 +16,29 @@ import argparse from pathlib import Path +import numpy as np import pandas as pd import pytest +import zedprofiler.cli as cli_module from zedprofiler.cli import ( _auto_requests, + _coerce_param, _output_path, _parse_feature_spec, _parse_name_path, + _resolve_channel, + _resolve_compartment, _resolve_requests, + _run_colocalization, + _run_single_channel, + _validate_request_channels_compartments, main, + trigger, ) from zedprofiler.featurization import texture from zedprofiler.identifiers import build_image_id +from zedprofiler.IO.loading_classes import ImageSetLoader tifffile = pytest.importorskip("tifffile") @@ -48,6 +58,7 @@ EXPECTED_OBJECT_COUNT = 5 EXPECTED_CHANNEL_COUNT = 2 +EXPECTED_TINY_OBJECT_COUNT = 2 PATIENT_TUMOR, PLATE, WELL, FIELD = "NF0014_T1", "PLATE01", "A1", "1" EXPECTED_IMAGE_ID = build_image_id(PATIENT_TUMOR, PLATE, WELL, FIELD) @@ -193,6 +204,168 @@ def test_resolve_requests_features_filter_applied_to_explicit_specs() -> None: assert requests[0]["type"] == "Intensity" +# --------------------------------------------------------------------------- +# Unit tests: parser/validation error branches +# --------------------------------------------------------------------------- + + +def test_parse_name_path_rejects_empty_name() -> None: + """A NAME=PATH token with an empty name is a user error.""" + with pytest.raises(argparse.ArgumentTypeError, match="empty"): + _parse_name_path("=/tmp/a.tif") + + +@pytest.mark.parametrize("spec", ["", " ", " , , "]) +def test_parse_feature_spec_rejects_empty_spec(spec: str) -> None: + """An all-whitespace feature spec yields no tokens and is rejected.""" + with pytest.raises(argparse.ArgumentTypeError, match="Empty feature spec"): + _parse_feature_spec(spec) + + +def test_coerce_param_returns_cast_value_and_raises_friendly_error() -> None: + """_coerce_param casts on success and raises ArgumentTypeError on failure.""" + assert _coerce_param("5", int, "distance", "Texture,distance=5") == int("5") + with pytest.raises(argparse.ArgumentTypeError, match="not a valid int"): + _coerce_param("not-a-number", int, "distance", "Texture,distance=bad") + + +def test_resolve_channel_requires_nonempty_string() -> None: + """A single-channel request missing 'channel' is rejected.""" + with pytest.raises(argparse.ArgumentTypeError, match="requires a 'channel' key"): + _resolve_channel({"type": "Intensity", "compartment": "Nuclei"}, "Intensity") + + +def test_resolve_compartment_requires_nonempty_string() -> None: + """A request missing 'compartment' is rejected.""" + with pytest.raises( + argparse.ArgumentTypeError, match="requires a 'compartment' key" + ): + _resolve_compartment({"type": "Intensity", "channel": "DNA"}, "Intensity") + + +def test_auto_requests_rejects_no_channels() -> None: + """Auto-generation with zero --image channels is a user error.""" + with pytest.raises(argparse.ArgumentTypeError, match="No --image flags"): + _auto_requests([], ["Nuclei"], ["Intensity"]) + + +def test_auto_requests_rejects_no_compartments() -> None: + """Auto-generation with zero --label compartments is a user error.""" + with pytest.raises(argparse.ArgumentTypeError, match="No --label flags"): + _auto_requests(["DNA"], [], ["Intensity"]) + + +def test_validate_rejects_undeclared_compartment() -> None: + """An explicit request referencing an undeclared compartment is rejected.""" + with pytest.raises(argparse.ArgumentTypeError, match="not declared via --label"): + _validate_request_channels_compartments( + [{"type": "Intensity", "channel": "DNA", "compartment": "Ghost"}], + ["DNA"], + ["Nuclei"], + ) + + +def test_validate_rejects_undeclared_colocalization_channel() -> None: + """A colocalization request with an undeclared channel is rejected.""" + with pytest.raises(argparse.ArgumentTypeError, match="not declared via --image"): + _validate_request_channels_compartments( + [ + { + "type": "Colocalization", + "channel1": "Ghost", + "channel2": "DNA", + "compartment": "Nuclei", + }, + ], + ["DNA"], + ["Nuclei"], + ) + + +# --------------------------------------------------------------------------- +# Unit tests: featurizer dispatch on tiny in-memory data +# --------------------------------------------------------------------------- + + +def _tiny_image_set_loader() -> ImageSetLoader: + """A small multi-channel loader for exercising dispatch branches quickly.""" + rng = np.random.default_rng(0) + label = np.zeros((6, 6, 6), dtype=np.int32) + label[1:3, 1:3, 1:3] = 1 + label[4:6, 4:6, 4:6] = 2 + image = rng.integers(0, 200, size=(6, 6, 6)).astype(np.float32) + image2 = rng.integers(0, 200, size=(6, 6, 6)).astype(np.float32) + return ImageSetLoader.from_image_dict( + {"DNA": image, "AGP": image2, "Nuclei": label}, + anisotropy_spacing=(2.0, 1.0, 1.0), + image_set_name="tiny", + label_key_names=["Nuclei"], + ) + + +@pytest.mark.parametrize( + "feature_type", + ["Neighbors", "Texture", "Granularity", "VolumeSizeShape"], +) +def test_run_single_channel_dispatches_each_type(feature_type: str) -> None: + """Each single-channel dispatch branch runs and returns a framed result.""" + request = {"type": feature_type, "channel": "DNA", "compartment": "Nuclei"} + channel, ran_type, df = _run_single_channel(_tiny_image_set_loader(), dict(request)) + assert channel == "DNA" + assert ran_type == feature_type + assert len(df) == EXPECTED_TINY_OBJECT_COUNT # two objects in the tiny label mask + + +def test_run_colocalization_requires_channel_keys() -> None: + """A colocalization request missing channel1/channel2 is rejected.""" + with pytest.raises(argparse.ArgumentTypeError, match="requires 'channel1'"): + _run_colocalization( + _tiny_image_set_loader(), + {"type": "Colocalization", "compartment": "Nuclei"}, + ) + + +# --------------------------------------------------------------------------- +# Unit tests: run() / main() / trigger() control flow +# --------------------------------------------------------------------------- + + +def test_run_returns_empty_when_features_filter_drops_all_requests( + tmp_path: Path, +) -> None: + """A --features filter that excludes every --feature request writes nothing.""" + argv = [ + "run", + "--image=DNA=/does/not/exist.tif", + "--label=Nuclei=/does/not/exist.tiff", + "--anisotropy-spacing", + "1.0", + "1.0", + "1.0", + "--patient-tumor=NF0014_T1", + "--plate=PLATE01", + "--well=A1", + "--field=1", + f"--out-dir={tmp_path}", + "--feature=Intensity,channel=DNA,compartment=Nuclei", + "--features=Colocalization", + ] + # No requests survive the Colocalization filter, so no images are read and + # no Parquet is written; the command still exits 0. + assert _run(argv) == 0 + assert not list(tmp_path.glob("*.parquet")) + + +def test_trigger_raises_system_exit_with_main_exit_code( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """trigger() wraps main()'s exit code in a SystemExit.""" + monkeypatch.setattr(cli_module, "main", lambda argv=None: 0) + with pytest.raises(SystemExit) as exc: + trigger() + assert exc.value.code == 0 + + # --------------------------------------------------------------------------- # End-to-end tests on the CellProfiler 3D tutorial data # --------------------------------------------------------------------------- From 3ca83828f4f615af850b9f964d2e7d6482b47faa Mon Sep 17 00:00:00 2001 From: d33bs Date: Thu, 6 Aug 2026 09:39:00 -0600 Subject: [PATCH 4/7] fov labeling Co-Authored-By: Mike Lippincott <58147848+MikeLippincott@users.noreply.github.com> --- src/zedprofiler/IO/loading_classes.py | 16 ++++++++-------- src/zedprofiler/cli.py | 14 +++++++------- src/zedprofiler/identifiers.py | 14 +++++++------- tests/IO/test_loading_classes.py | 18 +++++++++--------- tests/test_cli.py | 14 +++++++------- tests/test_identifiers.py | 22 +++++++++++----------- 6 files changed, 49 insertions(+), 49 deletions(-) diff --git a/src/zedprofiler/IO/loading_classes.py b/src/zedprofiler/IO/loading_classes.py index 38e380c..d012860 100644 --- a/src/zedprofiler/IO/loading_classes.py +++ b/src/zedprofiler/IO/loading_classes.py @@ -51,7 +51,7 @@ class ImageSetConfig: patient_tumor: str | None = None plate: str | None = None well: str | None = None - field: int | str | None = None + field_of_view: int | str | None = None # validate the arg types def __post_init__(self) -> None: @@ -68,8 +68,8 @@ def __post_init__(self) -> None: raise TypeError("plate must be a string or None") if not isinstance(self.well, (str, type(None))): raise TypeError("well must be a string or None") - if not isinstance(self.field, (int, str, type(None))): - raise TypeError("field must be an int, str, or None") + if not isinstance(self.field_of_view, (int, str, type(None))): + raise TypeError("field_of_view must be an int, str, or None") if self.label_key_name is None: self.label_key_name = [] @@ -88,13 +88,13 @@ def image_id(self) -> str | None: self.patient_tumor is not None and self.plate is not None and self.well is not None - and self.field is not None + and self.field_of_view is not None ): return build_image_id( patient_tumor=self.patient_tumor, plate=self.plate, well=self.well, - field=self.field, + field_of_view=self.field_of_view, ) return None @@ -253,7 +253,7 @@ def from_image_dict( # noqa: PLR0913 patient_tumor: str | None = None, plate: str | None = None, well: str | None = None, - field: int | str | None = None, + field_of_view: int | str | None = None, ) -> ImageSetLoader: """Build an ImageSetLoader from an in-memory channel/label dict. @@ -277,7 +277,7 @@ def from_image_dict( # noqa: PLR0913 Keys in ``image_dict`` that are compartment labels (not channels). Used by ``get_compartments`` to distinguish compartments from raw channels. - patient_tumor, plate, well, field : optional + patient_tumor, plate, well, field_of_view : optional Imaging-coordinate identifier fields. When all four are set, the loader's ``image_id`` is the deterministic ``Metadata_Imaging_ImageID``; otherwise it falls back to @@ -310,7 +310,7 @@ def from_image_dict( # noqa: PLR0913 patient_tumor=patient_tumor, plate=plate, well=well, - field=field, + field_of_view=field_of_view, ) self.image_id = ( config.image_id if config.image_id is not None else config.image_set_name diff --git a/src/zedprofiler/cli.py b/src/zedprofiler/cli.py index 1861b19..e701c1e 100644 --- a/src/zedprofiler/cli.py +++ b/src/zedprofiler/cli.py @@ -438,7 +438,7 @@ def _build_image_set_loader( for name, path in labels: image_dict[name] = _image_loading(path) compartment_names.append(name) - patient_tumor, plate, well, field = identifiers + patient_tumor, plate, well, field_of_view = identifiers image_set_name = build_image_id_from_identifiers(identifiers) image_set_loader = ImageSetLoader.from_image_dict( image_dict, @@ -448,7 +448,7 @@ def _build_image_set_loader( patient_tumor=patient_tumor, plate=plate, well=well, - field=field, + field_of_view=field_of_view, ) channels = [name for name, _ in images] return image_set_loader, channels, compartment_names @@ -458,8 +458,8 @@ def build_image_id_from_identifiers( identifiers: tuple[str, str, str, str], ) -> str: """Build the deterministic image set name from identifier fields.""" - patient_tumor, plate, well, field = identifiers - return build_image_id(patient_tumor, plate, well, field) + patient_tumor, plate, well, field_of_view = identifiers + return build_image_id(patient_tumor, plate, well, field_of_view) def run( # noqa: PLR0913, PLR0917 @@ -484,7 +484,7 @@ def run( # noqa: PLR0913, PLR0917 anisotropy_spacing : tuple[float, float, float] (z, y, x) spacing. identifiers : tuple[str, str, str, str] - (patient_tumor, plate, well, field) imaging coordinates. + (patient_tumor, plate, well, field_of_view) imaging coordinates. out_dir : Path Shard output directory. feature_specs : list[dict[str, object]] | None @@ -589,7 +589,7 @@ def _build_parser() -> argparse.ArgumentParser: run_parser.add_argument("--plate", required=True, help="Plate identifier.") run_parser.add_argument("--well", required=True, help="Well identifier.") run_parser.add_argument( - "--field", + "--fov", required=True, help="Field-of-view index or identifier.", ) @@ -654,7 +654,7 @@ def main(argv: Sequence[str] | None = None) -> int: args.patient_tumor, args.plate, args.well, - args.field, + args.fov, ) run( images=images, diff --git a/src/zedprofiler/identifiers.py b/src/zedprofiler/identifiers.py index 0cf8179..aed6ae7 100644 --- a/src/zedprofiler/identifiers.py +++ b/src/zedprofiler/identifiers.py @@ -5,7 +5,7 @@ plan's identifier spec). The central one is ``Metadata_Imaging_ImageID``, built deterministically from the four imaging coordinates: - patient-tumor, plate, well, field + patient-tumor, plate, well, field of view Because a shard is dispatched per well/FOV, every feature table emitted by a shard carries this single image id so downstream tables can rejoin without a @@ -25,15 +25,15 @@ def build_image_id( patient_tumor: str, plate: str, well: str, - field: int | str, + field_of_view: int | str, ) -> str: """Build a deterministic ``Metadata_Imaging_ImageID`` value. The id is a stable string assembled from the four imaging coordinates so that the same well/FOV always produces the same id across runs, batches, and reprocessing. Component order is fixed - (patient-tumor, plate, well, field) so ids sort and group naturally by - patient then plate then well then field. + (patient-tumor, plate, well, field of view) so ids sort and group + naturally by patient then plate then well then field of view. Parameters ---------- @@ -43,13 +43,13 @@ def build_image_id( Plate identifier (e.g. ``"PLATE01"``). well : str Well identifier (e.g. ``"A1"``). - field : int | str + field_of_view : int | str Field-of-view index or identifier (e.g. ``1`` or ``"f1"``). Returns ------- str - The deterministic image id, e.g. ``"NF0014_T1_PLATE01_A1_field1"``. + The deterministic image id, e.g. ``"NF0014_T1_PLATE01_A1_fov1"``. """ - return f"{patient_tumor}_{plate}_{well}_field{field}" + return f"{patient_tumor}_{plate}_{well}_fov{field_of_view}" diff --git a/tests/IO/test_loading_classes.py b/tests/IO/test_loading_classes.py index 258b4ee..5ae6ec6 100644 --- a/tests/IO/test_loading_classes.py +++ b/tests/IO/test_loading_classes.py @@ -66,7 +66,7 @@ def test_config_post_init_none_defaults(self) -> None: ({"patient_tumor": 123}, "patient_tumor must be a string or None"), ({"plate": 123}, "plate must be a string or None"), ({"well": 123}, "well must be a string or None"), - ({"field": 1.5}, "field must be an int, str, or None"), + ({"field_of_view": 1.5}, "field_of_view must be an int, str, or None"), ], ) def test_config_rejects_invalid_types( @@ -74,7 +74,7 @@ def test_config_rejects_invalid_types( kwargs: dict[str, object], message: str, ) -> None: - """ImageSetConfig should validate field types during initialization.""" + """ImageSetConfig should validate identifier field types.""" with pytest.raises(TypeError, match=message): ImageSetConfig(**kwargs) @@ -400,7 +400,7 @@ def test_image_id_is_none_when_any_field_missing(self) -> None: patient_tumor="NF0014_T1", plate="PLATE01", well="A1", - field=None, + field_of_view=None, ).image_id is None ) @@ -411,9 +411,9 @@ def test_image_id_built_when_all_fields_set(self) -> None: patient_tumor="NF0014_T1", plate="PLATE01", well="A1", - field=1, + field_of_view=1, ) - assert config.image_id == "NF0014_T1_PLATE01_A1_field1" + assert config.image_id == "NF0014_T1_PLATE01_A1_fov1" class TestFromImageDict: @@ -471,9 +471,9 @@ def test_image_id_populated_when_identifiers_provided(self) -> None: patient_tumor="NF0014_T1", plate="PLATE01", well="A1", - field=1, + field_of_view=1, ) - assert loader.image_id == "NF0014_T1_PLATE01_A1_field1" + assert loader.image_id == "NF0014_T1_PLATE01_A1_fov1" def test_image_id_falls_back_to_image_set_name_without_identifiers( self, @@ -497,7 +497,7 @@ def test_two_object_loader_on_from_image_dict(self) -> None: patient_tumor="X", plate="P", well="A1", - field="1", + field_of_view="1", ) two = TwoObjectLoader( image_set_loader=loader, @@ -506,4 +506,4 @@ def test_two_object_loader_on_from_image_dict(self) -> None: channel2="AGP", ) assert two.object_ids == [ONE_LABEL, TWO_LABEL] - assert two.image_set_loader.image_id == "X_P_A1_field1" + assert two.image_set_loader.image_id == "X_P_A1_fov1" diff --git a/tests/test_cli.py b/tests/test_cli.py index ece3ba4..89a08fd 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -59,8 +59,8 @@ EXPECTED_OBJECT_COUNT = 5 EXPECTED_CHANNEL_COUNT = 2 EXPECTED_TINY_OBJECT_COUNT = 2 -PATIENT_TUMOR, PLATE, WELL, FIELD = "NF0014_T1", "PLATE01", "A1", "1" -EXPECTED_IMAGE_ID = build_image_id(PATIENT_TUMOR, PLATE, WELL, FIELD) +PATIENT_TUMOR, PLATE, WELL, FOV = "NF0014_T1", "PLATE01", "A1", "1" +EXPECTED_IMAGE_ID = build_image_id(PATIENT_TUMOR, PLATE, WELL, FOV) def test_feature_namespace_import() -> None: @@ -345,7 +345,7 @@ def test_run_returns_empty_when_features_filter_drops_all_requests( "--patient-tumor=NF0014_T1", "--plate=PLATE01", "--well=A1", - "--field=1", + "--fov=1", f"--out-dir={tmp_path}", "--feature=Intensity,channel=DNA,compartment=Nuclei", "--features=Colocalization", @@ -392,7 +392,7 @@ def _base_run_args(out_dir: Path, *extra: str) -> list[str]: f"--patient-tumor={PATIENT_TUMOR}", f"--plate={PLATE}", f"--well={WELL}", - f"--field={FIELD}", + f"--fov={FOV}", f"--out-dir={out_dir}", *extra, ] @@ -474,7 +474,7 @@ def test_cli_skip_existing_skips_recompute(tmp_path: Path) -> None: f"--patient-tumor={PATIENT_TUMOR}", f"--plate={PLATE}", f"--well={WELL}", - f"--field={FIELD}", + f"--fov={FOV}", f"--out-dir={out_dir}", "--features=Intensity", "--skip-existing", @@ -504,7 +504,7 @@ def test_cli_colocalization_two_channels(tmp_path: Path) -> None: f"--patient-tumor={PATIENT_TUMOR}", f"--plate={PLATE}", f"--well={WELL}", - f"--field={FIELD}", + f"--fov={FOV}", f"--out-dir={out_dir}", "--feature=Colocalization,channel1=DNA1,channel2=DNA2,compartment=Nuclei,fast_costes=Faster", ], @@ -532,7 +532,7 @@ def test_cli_missing_required_arg_errors(tmp_path: Path) -> None: f"--patient-tumor={PATIENT_TUMOR}", f"--plate={PLATE}", f"--well={WELL}", - f"--field={FIELD}", + f"--fov={FOV}", ] with pytest.raises(SystemExit): _run(argv) diff --git a/tests/test_identifiers.py b/tests/test_identifiers.py index 363cefd..60f0101 100644 --- a/tests/test_identifiers.py +++ b/tests/test_identifiers.py @@ -8,9 +8,9 @@ def test_build_image_id_formats_all_fields() -> None: - """The image id joins patient-tumor, plate, well, and a fielded suffix.""" + """The image id joins patient-tumor, plate, well, and a field-of-view suffix.""" assert ( - build_image_id("NF0014_T1", "PLATE01", "A1", 1) == "NF0014_T1_PLATE01_A1_field1" + build_image_id("NF0014_T1", "PLATE01", "A1", 1) == "NF0014_T1_PLATE01_A1_fov1" ) @@ -20,25 +20,25 @@ def test_build_image_id_is_deterministic() -> None: assert build_image_id(*args) == build_image_id(*args) -def test_build_image_id_field_accepts_string_or_int() -> None: - """Field may be supplied as either an int or its string form.""" - assert build_image_id("X", "P", "A1", 3) == "X_P_A1_field3" - assert build_image_id("X", "P", "A1", "3") == "X_P_A1_field3" +def test_build_image_id_field_of_view_accepts_string_or_int() -> None: + """The field of view may be supplied as either an int or its string form.""" + assert build_image_id("X", "P", "A1", 3) == "X_P_A1_fov3" + assert build_image_id("X", "P", "A1", "3") == "X_P_A1_fov3" @pytest.mark.parametrize( - "patient_tumor,plate,well,field,expected", + "patient_tumor,plate,well,field_of_view,expected", [ - ("NF0014_T1", "PLATE01", "A1", "1", "NF0014_T1_PLATE01_A1_field1"), - ("NF0009_T2", "PLATE02", "B3", "7", "NF0009_T2_PLATE02_B3_field7"), + ("NF0014_T1", "PLATE01", "A1", "1", "NF0014_T1_PLATE01_A1_fov1"), + ("NF0009_T2", "PLATE02", "B3", "7", "NF0009_T2_PLATE02_B3_fov7"), ], ) def test_build_image_id_parametrized( patient_tumor: str, plate: str, well: str, - field: str, + field_of_view: str, expected: str, ) -> None: """Format holds across distinct imaging coordinates.""" - assert build_image_id(patient_tumor, plate, well, field) == expected + assert build_image_id(patient_tumor, plate, well, field_of_view) == expected From e31238be37165c62d8adb50c9a379489a791e4c4 Mon Sep 17 00:00:00 2001 From: d33bs Date: Thu, 6 Aug 2026 09:48:23 -0600 Subject: [PATCH 5/7] add cli docs to readme Co-Authored-By: Mike Lippincott <58147848+MikeLippincott@users.noreply.github.com> --- README.md | 57 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) diff --git a/README.md b/README.md index e8e5ce2..1b4faa7 100644 --- a/README.md +++ b/README.md @@ -31,6 +31,63 @@ Accepted image formats (order matters): - Single channel: `(z, y, x)` +## Command-line interface + +`ZedProfiler run` extracts features for a single well/field-of-view (FOV) shard: it loads one image set from explicit file paths, runs a selected subset of featurizers, and writes one Parquet per feature table to an output directory. It is the command a workflow manager (for example Nextflow via SLURM `sbatch`) dispatches once per shard. + +After `uv sync` (or `pip install .`) the `ZedProfiler` console script is available; from a checkout you can also use `uv run ZedProfiler run ...`. Run `ZedProfiler run --help` for the authoritative, up-to-date list of arguments. + +```bash +ZedProfiler run \ + --image=DNA=/path/to/channel1.tif \ + --label=Nuclei=/path/to/nuclei_mask.tiff \ + --anisotropy-spacing 1.0 1.0 1.0 \ + --patient-tumor NF0014_T1 \ + --plate PLATE01 \ + --well A1 \ + --fov 1 \ + --out-dir ./shard_output \ + --features Intensity +``` + +### Arguments + +| Argument | Required | Description | +| -------------------------------- | --------------------- | --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | +| `--image=NAME=PATH` | yes (>=1, repeatable) | A channel image as `NAME=PATH`. Repeat for multi-channel shards. | +| `--label=NAME=PATH` | yes (>=1, repeatable) | A compartment label mask as `NAME=PATH`. Repeat for multiple compartments. | +| `--anisotropy-spacing Z Y X` | yes | Z, Y, X voxel spacing (three floats). | +| `--patient-tumor` | yes | Patient-tumor identifier (e.g. `NF0014_T1`). | +| `--plate` | yes | Plate identifier. | +| `--well` | yes | Well identifier (e.g. `A1`). | +| `--fov` | yes | Field-of-view index or identifier. | +| `--out-dir` | yes | Shard output directory (created if needed). | +| `--features` | no | Comma-separated feature types to run (selector). With no `--feature` flags, runs these types over the channel x compartment cross-product; with `--feature` flags, restricts those requests by type. Default: all single-channel types, plus `Colocalization` when >=2 channels are declared. | +| `--feature=TYPE[,key=value,...]` | no (repeatable) | An explicit feature request, e.g. `Intensity,channel=DNA,compartment=Nuclei` or `Colocalization,channel1=DNA1,channel2=DNA2,compartment=Nuclei,fast_costes=Faster`. | +| `--skip-existing` | no | Skip a feature request whose output Parquet already exists. | +| `--force` | no | Overwrite even when the output exists (writes are still atomic). | + +### Feature types + +`VolumeSizeShape`, `Intensity`, `Neighbors`, `Texture`, and `Granularity` are single-channel features run per channel x compartment. `Colocalization` is a two-channel feature run per ordered channel pair x compartment. + +### Outputs + +Each request writes `{compartment}_{channel}_{feature_type}_cpu_features.parquet` into `--out-dir`. Every table carries `Metadata_Imaging_ImageID` (deterministically built from the patient-tumor, plate, well, and FOV coordinates) and `Metadata_Experiment_ImageSet` so downstream tables can rejoin. Writes are atomic (temp file + replace), so a crashed shard never leaves a partial file that `--skip-existing` would mistake for a complete one. + +### Two-channel colocalization example + +```bash +ZedProfiler run \ + --image=DNA1=/path/to/channel1.tif \ + --image=DNA2=/path/to/channel2.tif \ + --label=Nuclei=/path/to/nuclei_mask.tiff \ + --anisotropy-spacing 1.0 1.0 1.0 \ + --patient-tumor NF0014_T1 --plate PLATE01 --well A1 --fov 1 \ + --out-dir ./shard_output \ + --feature=Colocalization,channel1=DNA1,channel2=DNA2,compartment=Nuclei,fast_costes=Faster +``` + ## Quality Gates We lint and format code with our pre-commit configuration. From 9e70fef543e3a21f0e5d1b8e251d69716b2bde5b Mon Sep 17 00:00:00 2001 From: d33bs Date: Thu, 6 Aug 2026 09:58:43 -0600 Subject: [PATCH 6/7] feature specs Co-Authored-By: Mike Lippincott <58147848+MikeLippincott@users.noreply.github.com> --- src/zedprofiler/cli.py | 48 +++++++++++++++++++++++++++++------------- 1 file changed, 33 insertions(+), 15 deletions(-) diff --git a/src/zedprofiler/cli.py b/src/zedprofiler/cli.py index e701c1e..2bf6c06 100644 --- a/src/zedprofiler/cli.py +++ b/src/zedprofiler/cli.py @@ -51,22 +51,40 @@ # Minimum number of channels required for colocalization requests. _MIN_CHANNELS_FOR_COLOCALIZATION = 2 -# Feature types that consume a single channel + compartment via ObjectLoader. -_SINGLE_CHANNEL_TYPES = ( - "VolumeSizeShape", - "Intensity", - "Neighbors", - "Texture", - "Granularity", +# Each feature type is declared once, with the two independent axes the CLI +# cares about: +# two_channel - True if the feature reads two channel images via +# TwoObjectLoader; False if it reads one via ObjectLoader. +# channel_agnostic - True if the computation ignores the channel pixel data +# (it only uses the compartment mask). Such features are +# still namespaced by a channel for warehouse organization, +# so a channel (for naming only) is still required. These +# are always a subset of the single-channel types: a +# two-channel feature by definition uses its channel +# images, so it cannot be channel-agnostic. +_FEATURE_SPECS: dict[str, dict[str, bool]] = { + "VolumeSizeShape": {"two_channel": False, "channel_agnostic": True}, + "Intensity": {"two_channel": False, "channel_agnostic": False}, + "Neighbors": {"two_channel": False, "channel_agnostic": True}, + "Texture": {"two_channel": False, "channel_agnostic": False}, + "Granularity": {"two_channel": False, "channel_agnostic": False}, + "Colocalization": {"two_channel": True, "channel_agnostic": False}, +} + +# Loader axis: how many channel images each feature reads. +_SINGLE_CHANNEL_TYPES = tuple( + name for name, spec in _FEATURE_SPECS.items() if not spec["two_channel"] +) +_TWO_CHANNEL_TYPES = tuple( + name for name, spec in _FEATURE_SPECS.items() if spec["two_channel"] +) +ALL_FEATURE_TYPES = tuple(_FEATURE_SPECS) + +# Computation axis: features that ignore the channel pixels (a subset of +# _SINGLE_CHANNEL_TYPES). +_CHANNEL_AGNOSTIC_TYPES = tuple( + name for name, spec in _FEATURE_SPECS.items() if spec["channel_agnostic"] ) -# Feature types that require two channels via TwoObjectLoader. -_TWO_CHANNEL_TYPES = ("Colocalization",) -ALL_FEATURE_TYPES = (*_SINGLE_CHANNEL_TYPES, *_TWO_CHANNEL_TYPES) - -# Channel-agnostic features: their computation does not use the channel image, -# but like every ZedProfiler feature they are namespaced by a channel for -# warehouse organization, so a channel (for naming only) is still required. -_CHANNEL_AGNOSTIC_TYPES = ("VolumeSizeShape", "Neighbors") def _parse_name_path(token: str) -> tuple[str, Path]: From a93da014748229f91077e63a13208b1c689049e4 Mon Sep 17 00:00:00 2001 From: d33bs Date: Thu, 6 Aug 2026 10:12:59 -0600 Subject: [PATCH 7/7] cli spec Co-Authored-By: Mike Lippincott <58147848+MikeLippincott@users.noreply.github.com> --- src/zedprofiler/cli.py | 66 ++++++++++++++++++++---------------------- 1 file changed, 32 insertions(+), 34 deletions(-) diff --git a/src/zedprofiler/cli.py b/src/zedprofiler/cli.py index 2bf6c06..0a38a1b 100644 --- a/src/zedprofiler/cli.py +++ b/src/zedprofiler/cli.py @@ -51,39 +51,32 @@ # Minimum number of channels required for colocalization requests. _MIN_CHANNELS_FOR_COLOCALIZATION = 2 -# Each feature type is declared once, with the two independent axes the CLI -# cares about: -# two_channel - True if the feature reads two channel images via -# TwoObjectLoader; False if it reads one via ObjectLoader. -# channel_agnostic - True if the computation ignores the channel pixel data -# (it only uses the compartment mask). Such features are -# still namespaced by a channel for warehouse organization, -# so a channel (for naming only) is still required. These -# are always a subset of the single-channel types: a -# two-channel feature by definition uses its channel -# images, so it cannot be channel-agnostic. -_FEATURE_SPECS: dict[str, dict[str, bool]] = { - "VolumeSizeShape": {"two_channel": False, "channel_agnostic": True}, - "Intensity": {"two_channel": False, "channel_agnostic": False}, - "Neighbors": {"two_channel": False, "channel_agnostic": True}, - "Texture": {"two_channel": False, "channel_agnostic": False}, - "Granularity": {"two_channel": False, "channel_agnostic": False}, - "Colocalization": {"two_channel": True, "channel_agnostic": False}, -} - -# Loader axis: how many channel images each feature reads. -_SINGLE_CHANNEL_TYPES = tuple( - name for name, spec in _FEATURE_SPECS.items() if not spec["two_channel"] +# Feature types partitioned into three mutually exclusive groups (their union +# is ALL_FEATURE_TYPES), so each type appears in exactly one list: +# - _SINGLE_CHANNEL_TYPES: read one channel + one compartment via ObjectLoader +# and actually use the channel pixel data. +# - _TWO_CHANNEL_TYPES: read two channels via TwoObjectLoader. +# - _CHANNEL_AGNOSTIC_TYPES: read one channel via ObjectLoader but ignore the +# channel pixel data (they only use the compartment mask). They are still +# namespaced by a channel for warehouse organization, so a channel (for +# naming only) is still required. +_SINGLE_CHANNEL_TYPES = ( + "Intensity", + "Texture", + "Granularity", ) -_TWO_CHANNEL_TYPES = tuple( - name for name, spec in _FEATURE_SPECS.items() if spec["two_channel"] -) -ALL_FEATURE_TYPES = tuple(_FEATURE_SPECS) - -# Computation axis: features that ignore the channel pixels (a subset of -# _SINGLE_CHANNEL_TYPES). -_CHANNEL_AGNOSTIC_TYPES = tuple( - name for name, spec in _FEATURE_SPECS.items() if spec["channel_agnostic"] +# Feature types that require two channels via TwoObjectLoader. +_TWO_CHANNEL_TYPES = ("Colocalization",) + +# Channel-agnostic features: their computation does not use the channel image, +# but like every ZedProfiler feature they are namespaced by a channel for +# warehouse organization, so a channel (for naming only) is still required. +_CHANNEL_AGNOSTIC_TYPES = ("VolumeSizeShape", "Neighbors") + +ALL_FEATURE_TYPES = ( + *_SINGLE_CHANNEL_TYPES, + *_TWO_CHANNEL_TYPES, + *_CHANNEL_AGNOSTIC_TYPES, ) @@ -405,9 +398,14 @@ def _resolve_requests( requests = [r for r in requests if str(r["type"]) in features_filter] _validate_request_channels_compartments(requests, channels, compartments) return requests - feature_types = features_filter if features_filter else list(_SINGLE_CHANNEL_TYPES) + # Default (no --features, no --feature): run every non-colocalization type + # (single-channel + channel-agnostic), plus Colocalization when there are + # enough channels. _SINGLE_CHANNEL_TYPES alone is not the full default + # because the channel-agnostic types live in their own group. + default_types = [*_SINGLE_CHANNEL_TYPES, *_CHANNEL_AGNOSTIC_TYPES] + feature_types = features_filter if features_filter else default_types if (features_filter is None) and len(channels) >= _MIN_CHANNELS_FOR_COLOCALIZATION: - feature_types = [*_SINGLE_CHANNEL_TYPES, "Colocalization"] + feature_types = [*default_types, "Colocalization"] return _auto_requests(channels, compartments, feature_types)