diff --git a/pyproject.toml b/pyproject.toml index be84cf3..1bd64dc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,7 +10,7 @@ readme = "README.md" requires-python = ">=3.11" dependencies = [ "pyyaml", - "jfjoch-client", + "jfjoch-client>=1.0.0rc166", "pydantic", "numpy", "scipy", diff --git a/src/aarecommon/math/jfjoch_gridscan_union.py b/src/aarecommon/math/jfjoch_gridscan_union.py index 1dfadc4..070a672 100644 --- a/src/aarecommon/math/jfjoch_gridscan_union.py +++ b/src/aarecommon/math/jfjoch_gridscan_union.py @@ -40,7 +40,8 @@ the picture if the database wants one. from __future__ import annotations import io -from typing import Any, Optional, Sequence +from collections.abc import Sequence +from typing import Any import numpy as np from pydantic import BaseModel, Field @@ -68,7 +69,7 @@ except ImportError: # pragma: no cover z: float = 0.0 -__all__ = ["analyse", "GridScanResult", "Centre", "Size", "Counts", "Thresholds"] +__all__ = ["Centre", "Counts", "GridScanResult", "Size", "Thresholds", "analyse"] #: DAQStatusModel, or anything else carrying the beam: its ``.geom`` is a #: SampleGeometryModel, and a ScanResultPayloadModel has ``beam_size_mm`` itself. @@ -98,10 +99,10 @@ class Centre(BaseModel): nx: float ny: float image_number: int = Field(description="Nearest collected image, for addressing") - x_um: Optional[float] = Field( + x_um: float | None = Field( None, description="Signed offset from the centre of cell (0,0), along the grid's own axes" ) - y_um: Optional[float] = None + y_um: float | None = None def offset_mm(self) -> Coordinate: """Offset in mm from the centre of cell (0,0), as an aarecommon Coordinate. @@ -123,17 +124,17 @@ class Size(BaseModel): narrower than the beam still reports a size instead of nothing. """ - x_um: Optional[float] = None - y_um: Optional[float] = None - z_um: Optional[float] = Field(None, description="From an orthogonal scan, if passed") - x_um_deconv: Optional[float] = None - y_um_deconv: Optional[float] = None - z_um_deconv: Optional[float] = None - beam_x_um: Optional[float] = Field(None, description="Beam FWHM used to deconvolve x") - beam_y_um: Optional[float] = Field(None, description="Beam FWHM used to deconvolve y") - area_um2: Optional[float] = None - equiv_diameter_um: Optional[float] = None - volume_pl: Optional[float] = Field(None, description="Ellipsoid; needs all three axes") + x_um: float | None = None + y_um: float | None = None + z_um: float | None = Field(None, description="From an orthogonal scan, if passed") + x_um_deconv: float | None = None + y_um_deconv: float | None = None + z_um_deconv: float | None = None + beam_x_um: float | None = Field(None, description="Beam FWHM used to deconvolve x") + beam_y_um: float | None = Field(None, description="Beam FWHM used to deconvolve y") + area_um2: float | None = None + equiv_diameter_um: float | None = None + volume_pl: float | None = Field(None, description="Ellipsoid; needs all three axes") volume_deconvolved: bool = Field( False, description="True if volume_pl used the beam-removed axes" ) @@ -168,44 +169,50 @@ class Counts(BaseModel): n_noise: int = Field(description="Spots present but neither protein nor ice") n_blank: int peak_protein_spots: int - mean_protein_in_contour: Optional[float] = None - best_resolution_A: Optional[float] = Field(None, description="Best res inside the contour") - median_bkg: Optional[float] = None - ice_in_contour: Optional[float] = Field( + mean_protein_in_contour: float | None = None + best_resolution_A: float | None = Field(None, description="Best res inside the contour") + median_bkg: float | None = None + ice_in_contour: float | None = Field( None, description="Fraction of frames inside the contour showing ice" ) n_fragments: int = Field( 0, description="Separate blobs at the 50 % level; >1 means more than one region" ) ice_fraction: float = Field(description="Fraction of frames showing ice") - unit_cell_agreement: Optional[float] = Field( + unit_cell_agreement: float | None = Field( None, description="Fraction of indexed frames in the contour sharing one cell" ) class GridScanResult(BaseModel): found: bool - reason: Optional[str] = Field(None, description="Why not, when found is False") - file_prefix: Optional[str] = None + reason: str | None = Field(None, description="Why not, when found is False") + file_prefix: str | None = None n_fast: int = 0 n_slow: int = 0 channel: str = Field("", description="Which spot channel built the map") - centre: Optional[Centre] = None - size: Optional[Size] = None - counts: Optional[Counts] = None - unit_cell: Optional[list[float]] = Field( + centre: Centre | None = None + size: Size | None = None + counts: Counts | None = None + unit_cell: list[float] | None = Field( None, description="Median cell of the contour, when check_unit_cell is on" ) warnings: list[str] = Field( default_factory=list, description="Reasons to distrust a found=True result" ) - jpeg: Optional[bytes] = Field( + contour_cells: list[tuple[int, int]] | None = Field( + None, + description="(nx, ny) of every collected cell inside the chosen 50 % " + "blob -- the decision's own footprint, for drawing it over the sample " + "image without recomputing the map", + ) + jpeg: bytes | None = Field( None, exclude=True, repr=False, description="The picture, in whichever single style was asked for", ) - jpeg_cells: Optional[bytes] = Field( + jpeg_cells: bytes | None = Field( None, exclude=True, repr=False, @@ -225,10 +232,10 @@ class Thresholds(BaseModel): "the bkg-only gate then rejects the true crystal", ) ice_spots: int = Field(5, description="spots_ice at or above this counts as ice") - ice_ring: Optional[float] = Field( + ice_ring: float | None = Field( None, description="Also call ice if scan_result 'ice' exceeds this" ) - beam_um: Optional[float] = Field( + beam_um: float | None = Field( None, description="Beam FWHM in um, used only when no DAQ status is passed. " "A single number is taken as a square beam", @@ -284,7 +291,7 @@ def _upsample(a: np.ndarray, f: int) -> np.ndarray: return np.asarray(im.resize((nx * fx, ny * fy), Image.BICUBIC), dtype=float) -def _norm(a: np.ndarray, floor: float = 0.0, ref: Optional[np.ndarray] = None) -> np.ndarray: +def _norm(a: np.ndarray, floor: float = 0.0, ref: np.ndarray | None = None) -> np.ndarray: """Normalise for colour against the channel's own range, but never below floor. ``ref`` supplies the range: the bounds come from the frames actually measured, @@ -372,7 +379,7 @@ _BEAM_PATHS = ( ) -def _beam_um(obj: Optional[BeamSource]) -> tuple[Optional[float], Optional[float]]: +def _beam_um(obj: BeamSource | None) -> tuple[float | None, float | None]: """Beam FWHM (x, y) in micrometres, dug out of a DAQ status or scan payload. The beam is a Coordinate, not one number: on a beamline with an 80 x 20 um @@ -406,7 +413,7 @@ def _beam_um(obj: Optional[BeamSource]) -> tuple[Optional[float], Optional[float return None, None -def _deconvolve(fwhm: Optional[float], beam: Optional[float]) -> Optional[float]: +def _deconvolve(fwhm: float | None, beam: float | None) -> float | None: """Quadrature removal of the beam. None when the beam swallows the feature.""" if fwhm is None or beam is None: return fwhm @@ -489,11 +496,11 @@ def _classify(F: dict, th: Thresholds) -> Counts: def analyse( scan_result: ScanResult, - grid_scan: Optional[GridScan] = None, + grid_scan: GridScan | None = None, *, - daq: Optional[BeamSource] = None, - thresholds: Optional[Thresholds] = None, - z_um: Optional[float] = None, + daq: BeamSource | None = None, + thresholds: Thresholds | None = None, + z_um: float | None = None, jpeg: bool = True, ) -> GridScanResult: """Analyse one grid scan. @@ -610,6 +617,9 @@ def analyse( if keep.shape == inside.shape: inside = inside & keep counts.n_fragments = n_frag + if cell.any(): + yy_c, xx_c = np.nonzero(cell) + out.contour_cells = [(int(x), int(y)) for y, x in zip(yy_c, xx_c)] if line: # One position wide: the 50 % crossings either side of the peak. The outer @@ -640,7 +650,7 @@ def analyse( w_cells = ((xx.max() - xx.min() + 1) / mx * nx_n, (yy.max() - yy.min() + 1) / my * ny_n) # the nearest image actually collected, for the DAQ to address - iy, ix = int(round(np.clip(gy, 0, ny_n - 1))), int(round(np.clip(gx, 0, nx_n - 1))) + iy, ix = round(np.clip(gy, 0, ny_n - 1)), round(np.clip(gx, 0, nx_n - 1)) num = int(F["number"][iy, ix]) if num < 0: # that cell was never collected yy2, xx2 = np.nonzero(F["number"] >= 0) @@ -742,7 +752,7 @@ def _edge(mask: np.ndarray, width: int = 1) -> np.ndarray: return mask & ~core -def _pictures(F, inside, centre, res: "GridScanResult", th: Thresholds, obj_level: float) -> None: +def _pictures(F, inside, centre, res: GridScanResult, th: Thresholds, obj_level: float) -> None: """Fill res.jpeg, and res.jpeg_cells too when both styles were asked for.""" want = ("smooth", "cells") if th.jpeg_style == "both" else (th.jpeg_style,) res.jpeg = _render(F, inside, centre, res, th, obj_level, want[0]) @@ -751,7 +761,7 @@ def _pictures(F, inside, centre, res: "GridScanResult", th: Thresholds, obj_leve def _render( - F, inside, centre, res: "GridScanResult", th: Thresholds, obj_level: float, style: str + F, inside, centre, res: GridScanResult, th: Thresholds, obj_level: float, style: str ) -> bytes: """Map, contours, crosshair and a one-line caption, as JPEG. diff --git a/src/aarecommon/models/gridscan_decision.py b/src/aarecommon/models/gridscan_decision.py new file mode 100644 index 0000000..b8cfeda --- /dev/null +++ b/src/aarecommon/models/gridscan_decision.py @@ -0,0 +1,36 @@ +"""The grid-scan centring decision, as taken at the beamline. + +The DAQ runs the analysis (``aarecommon.math.jfjoch_gridscan_union``) to pick +the collection position; this wraps its result verbatim with provenance and +rides on ``RasterPayloadModel.decision`` into AareDB, whose Results viewer +displays THE decision instead of recomputing one -- two implementations of the +same analysis drifting apart is how a UI ends up contradicting the beamline. + +The pictures are excluded from the result's serialization by the model itself; +the DAQ ships them through the ordinary image pipeline. +""" + +from typing import Any + +from pydantic import BaseModel, Field + +from aarecommon.math.jfjoch_gridscan_union import GridScanResult + +GRIDSCAN_DECISION_CONTRACT = 1 + + +class GridScanDecision(BaseModel): + contract: int = GRIDSCAN_DECISION_CONTRACT + algorithm: str = Field( + description="Which analysis took the decision, e.g. a " + "GridscanAnalysisMode value such as 'jfjoch_gridscan_union'" + ) + algorithm_version: str | None = Field( + None, description="aarecommon version (or commit) the DAQ ran" + ) + thresholds: dict[str, Any] | None = Field( + None, + description="Thresholds.model_dump() used at decision time, so a " + "stored decision states its own tuning", + ) + result: GridScanResult diff --git a/src/aarecommon/models/raster_grid.py b/src/aarecommon/models/raster_grid.py index ed9a1b0..db37e4a 100644 --- a/src/aarecommon/models/raster_grid.py +++ b/src/aarecommon/models/raster_grid.py @@ -6,6 +6,7 @@ from pydantic import AfterValidator, BaseModel, Field from aarecommon.math.coordinate import Coordinate, SmargonCoordinate, positive_coords from aarecommon.math.sample_geometry import SampleGeometryModel +from aarecommon.models.gridscan_decision import GridScanDecision class RasterGridRequest(BaseModel): @@ -92,3 +93,6 @@ class RasterPayloadModel(BaseModel): cell_size_pxl: Annotated[Coordinate, AfterValidator(positive_coords)] beam_mark_pxl: tuple[float, float] beam_size_mm: Annotated[Coordinate, AfterValidator(positive_coords)] + # The centring decision as taken at the beamline (None from DAQs that + # have not adopted it yet); AareDB stores and displays it verbatim. + decision: GridScanDecision | None = None diff --git a/tests/test_gridscan_decision.py b/tests/test_gridscan_decision.py new file mode 100644 index 0000000..cbd1cfa --- /dev/null +++ b/tests/test_gridscan_decision.py @@ -0,0 +1,68 @@ +"""GridScanDecision wire contract: wraps the analysis result verbatim.""" + +import json + +import numpy as np + +from aarecommon.math.jfjoch_gridscan_union import Thresholds, analyse +from aarecommon.models.gridscan_decision import GridScanDecision +from aarecommon.models.raster_grid import RasterPayloadModel + + +def _synthetic_scan(ny_n=20, nx_n=20): + """A gaussian protein blob on a flat loop, as plain dicts (jfjoch-free).""" + rng = np.random.default_rng(0) + yy, xx = np.mgrid[0:ny_n, 0:nx_n] + blob = 80 * np.exp(-(((xx - 7) ** 2 + (yy - 12) ** 2) / 8.0)) + images = [ + { + "number": int(y * nx_n + x), + "nx": int(x), + "ny": int(y), + "spots": int(blob[y, x] + 10 + rng.integers(0, 3)), + "spots_low_res": int(4 + rng.integers(0, 2)), + "spots_ice": 0, + "bkg": float(30 + 40 * np.exp(-(((x - 9) ** 2 + (y - 10) ** 2) / 90.0))), + "res": 2.0, + } + for y in range(ny_n) + for x in range(nx_n) + ] + return {"file_prefix": "decision-test", "images": images} + + +def test_decision_round_trips_without_pictures(request): + th = Thresholds.model_validate({}) + r = analyse(_synthetic_scan(), {"step_x_um": 10.0, "step_y_um": 10.0}, thresholds=th) + assert r.found + assert r.contour_cells, "the chosen 50 % blob must ship its cells" + assert all(0 <= x < r.n_fast and 0 <= y < r.n_slow for x, y in r.contour_cells) + # the centre lies within the contour's bounding box + assert r.centre is not None + xs = [x for x, _ in r.contour_cells] + ys = [y for _, y in r.contour_cells] + assert min(xs) - 1 <= r.centre.nx <= max(xs) + 1 + assert min(ys) - 1 <= r.centre.ny <= max(ys) + 1 + + d = GridScanDecision( + algorithm="jfjoch_gridscan_union", + algorithm_version="test", + thresholds=th.model_dump(), + result=r, + ) + wire = json.loads(d.model_dump_json()) + assert wire["contract"] == 1 + assert "jpeg" not in wire["result"], "pictures never ride the JSON" + assert wire["result"]["centre"]["image_number"] == r.centre.image_number + assert wire["thresholds"]["object_union_spots"] is True + back = GridScanDecision.model_validate(wire) + assert back.result.centre is not None + assert back.result.contour_cells is not None + assert back.result.centre.nx == r.centre.nx + assert [tuple(c) for c in back.result.contour_cells] == r.contour_cells + + +def test_payload_decision_is_optional(): + fields = RasterPayloadModel.model_fields + assert "decision" in fields + assert fields["decision"].default is None