diff --git a/src/aarecommon/math/jfjoch_gridscan_union.py b/src/aarecommon/math/jfjoch_gridscan_union.py index 55efcb4..b45898d 100644 --- a/src/aarecommon/math/jfjoch_gridscan_union.py +++ b/src/aarecommon/math/jfjoch_gridscan_union.py @@ -17,16 +17,22 @@ The beam size comes from the DAQ status (or from a ScanResultPayloadModel, which carries its own ``beam_size_mm``); it is the one quantity the Jungfraujoch files never record, and without it every size is an upper bound. -The protein channel is ``spots - spots_ice - spots_low_res``: spots that are -neither inside an ice band nor above the low-resolution cut, so what remains is -diffraction below the cut -- the requirement that a crystal show spots below -5 A. Nothing derived from indexing is used: indexing fires on ice often enough -that a map built from it carries that error in. +The protein channel is ``spots_low_res``: Bragg spots beyond 5 A, which only a +large unit cell produces -- ice cannot (its innermost ring is at 3.9 A) and salt +cannot. The alternative, ``spots - spots_ice - spots_low_res``, fails on a loop +caked in ice: ring spots leaking past the ice windows form a haze that follows +the ice, not the crystal, and it centred on the film instead. Nothing derived +from indexing is used: indexing fires on ice often enough that a map built from +it carries that error in. -The centre is the *area centroid of the region enclosed by the 50 % contour*, -not an intensity-weighted centre of the whole map. A weighted mean is dragged -by the weak signal trailing across the rest of the loop and lands on the -shoulder of the contour rather than inside it. +A cell counts as protein when its count clears the scan's own noise floor, +measured from the cells off the sample. The crystal is the region above half +the peak (and never below that floor), and the centre is the *area centroid* of +that region -- not an intensity-weighted mean, which is dragged by weak signal +across the loop and lands on the contour's shoulder. Three gates decide whether +there is a crystal at all: the peak must clear the floor by 1.5x, and the region +must be one compact patch rather than scattered cells. When they fail, the +centre of the object itself is returned instead, marked as such. Input may be jfjoch-client models or plain dicts. Needs numpy, pillow, pydantic; no scipy, no matplotlib. @@ -42,8 +48,7 @@ the picture if the database wants one. from __future__ import annotations import io -from collections.abc import Sequence -from typing import Any +from typing import Any, Optional, Sequence import numpy as np from pydantic import BaseModel, Field @@ -51,16 +56,27 @@ from pydantic import BaseModel, Field # The structures this exchanges with AareDAQ and Jungfraujoch. Verified against # jfjoch-client 1.0.0-rc.165. The fallbacks let the module import and its self-test # run outside AareDAQ; inside it the real classes are always the ones used. -from jfjoch_client.models.grid_scan import GridScan -from jfjoch_client.models.scan_result import ScanResult -from jfjoch_client.models.scan_result_images_inner import ScanResultImagesInner -from jfjoch_client.models.unit_cell import UnitCell +try: + from jfjoch_client.models.grid_scan import GridScan + from jfjoch_client.models.scan_result import ScanResult + from jfjoch_client.models.scan_result_images_inner import ScanResultImagesInner + from jfjoch_client.models.unit_cell import UnitCell +except ImportError: # pragma: no cover + GridScan = ScanResult = ScanResultImagesInner = UnitCell = Any -from aarecommon.math.coordinate import Coordinate -from aarecommon.math.sample_geometry import SampleGeometryModel +try: + from aarecommon.math.coordinate import Coordinate + from aarecommon.math.sample_geometry import SampleGeometryModel +except ImportError: # pragma: no cover + SampleGeometryModel = Any + + class Coordinate(BaseModel): # same shape as aarecommon's + x: float = 0.0 + y: float = 0.0 + z: float = 0.0 -__all__ = ["Centre", "Counts", "GridScanResult", "Size", "Thresholds", "analyse"] +__all__ = ["analyse", "GridScanResult", "Centre", "Size", "Counts", "Thresholds"] #: DAQStatusModel, or anything else carrying the beam: its ``.geom`` is a #: SampleGeometryModel, and a ScanResultPayloadModel has ``beam_size_mm`` itself. @@ -90,10 +106,16 @@ class Centre(BaseModel): nx: float ny: float image_number: int = Field(description="Nearest collected image, for addressing") - x_um: float | None = Field( + resolution_A: Optional[float] = Field( + None, + description="Jungfraujoch res on the frame at this position -- what the " + "crystal gives where the beam will actually sit, not the " + "best value found anywhere in the contour", + ) + x_um: Optional[float] = Field( None, description="Signed offset from the centre of cell (0,0), along the grid's own axes" ) - y_um: float | None = None + y_um: Optional[float] = None def offset_mm(self) -> Coordinate: """Offset in mm from the centre of cell (0,0), as an aarecommon Coordinate. @@ -115,17 +137,17 @@ class Size(BaseModel): narrower than the beam still reports a size instead of nothing. """ - 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") + 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") volume_deconvolved: bool = Field( False, description="True if volume_pl used the beam-removed axes" ) @@ -159,51 +181,83 @@ class Counts(BaseModel): n_ice: int n_noise: int = Field(description="Spots present but neither protein nor ice") n_blank: int + spot_floor: float = Field( + 0.0, + description="Spot count a frame had to clear to count as protein, after " + "the off-object baseline was applied", + ) + off_object_median: Optional[float] = Field( + None, description="Median protein channel where the background says no sample" + ) peak_protein_spots: int - 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( + contrast: Optional[float] = Field( + None, + description="Strongest frame divided by the floor it had to clear. " + "Near 1 the map is noise dressed as a crystal; a real " + "crystal is several times its floor", + ) + 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( 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" ) + n_contour_cells: int = Field( + 0, + description="Collected cells inside the accepted blob. A crystal is a " + "contiguous patch; one or two isolated cells over the level " + "is what noise produces", + ) + contour_compactness: Optional[float] = Field( + None, + description="Cells in the accepted blob over all cells above the " + "level. Near 1 the signal is one region; near 0 it is " + "scattered specks and the map is noise", + ) ice_fraction: float = Field(description="Fraction of frames showing ice") - unit_cell_agreement: float | None = Field( + unit_cell_agreement: Optional[float] = Field( None, description="Fraction of indexed frames in the contour sharing one cell" ) class GridScanResult(BaseModel): - found: bool - reason: str | None = Field(None, description="Why not, when found is False") - file_prefix: str | None = None + found: bool = Field( + description="A protein crystal was located. False still " + "leaves a usable centre when centre_source is " + "'background'" + ) + reason: Optional[str] = Field(None, description="Why not, when found is False") + centre_source: str = Field( + "", + description="'protein' -- the 50 % contour of the protein channel, a " + "crystal. 'background' -- the middle of the sample itself, " + "used when no frame shows convincing protein diffraction; " + "point the beam there, but do not collect a dataset on it " + "without rescanning", + ) + file_prefix: Optional[str] = None n_fast: int = 0 n_slow: int = 0 channel: str = Field("", description="Which spot channel built the map") - centre: Centre | None = None - size: Size | None = None - counts: Counts | None = None - unit_cell: list[float] | None = Field( + centre: Optional[Centre] = None + size: Optional[Size] = None + counts: Optional[Counts] = None + unit_cell: Optional[list[float]] = 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" ) - 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( + jpeg: Optional[bytes] = Field( None, exclude=True, repr=False, description="The picture, in whichever single style was asked for", ) - jpeg_cells: bytes | None = Field( + jpeg_cells: Optional[bytes] = Field( None, exclude=True, repr=False, @@ -214,19 +268,56 @@ class GridScanResult(BaseModel): class Thresholds(BaseModel): """Everything tunable. Defaults come from ~250 loops of PSI beamline data.""" - min_spots: int = Field(10, description="Protein spots needed before a contour means anything") - object_union_spots: bool = Field( - True, - description="Treat cells with >= min_spots protein as part of the object " - "even where the bkg proxy disagrees. A crystal on a thin tip has " - "low bkg while a metal pin foot dominates the percentile band; " - "the bkg-only gate then rejects the true crystal", + min_spots: int = Field( + 10, + description="Absolute floor on protein spots. Acts as a lower bound on " + "the baseline-derived one below, never as a replacement", + ) + baseline_sigma: Optional[float] = Field( + 5.0, + description="Counting-noise sigmas (sqrt of the off-object median) above " + "that median that set the real spot floor, so the threshold " + "tracks each scan's own pedestal. None restores the fixed " + "min_spots", + ) + baseline_min_cells: int = Field( + 20, + description="Off-object cells needed before their percentile is trusted; " + "below this the fixed min_spots is used", + ) + contrast_min: float = Field( + 1.5, + description="Peak over spot floor below which the protein map is not " + "trusted to centre on: the strongest cell is then no " + "different from the noise, and centring falls back to the " + "middle of the object. Real crystals run 5-30x. Note this " + "is a per-frame statistic and says nothing about whether " + "the signal is spatially organised -- compactness_min does " + "that", + ) + two_crystal_cells: int = Field( + 4, + description="A contour split into no more than two patches, the larger " + "holding at least this many cells, is two crystals rather " + "than scatter and passes the compactness gate", + ) + compactness_min: float = Field( + 0.55, + description="Cells in the largest blob over all cells above the " + "contour level. A crystal is one contiguous patch and " + "scores near 1; signal scattered over the loop scores " + "low. On a permutation control -- the same frames dealt " + "to different cells -- real scans held a median 1.00 and " + "shuffled ones 0.29, and this cut removed 86 % of the " + "shuffled false positives while keeping 95 % of the real " + "crystals. It is the only test here with any power " + "against spatially random signal", ) ice_spots: int = Field(5, description="spots_ice at or above this counts as ice") - ice_ring: float | None = Field( + ice_ring: Optional[float] = Field( None, description="Also call ice if scan_result 'ice' exceeds this" ) - beam_um: float | None = Field( + beam_um: Optional[float] = Field( None, description="Beam FWHM in um, used only when no DAQ status is passed. " "A single number is taken as a square beam", @@ -282,7 +373,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: np.ndarray | None = None) -> np.ndarray: +def _norm(a: np.ndarray, floor: float = 0.0, ref: Optional[np.ndarray] = 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, @@ -326,6 +417,42 @@ def _component(mask: np.ndarray, seed) -> np.ndarray: cur = g +def _nearest_image(F, gx: float, gy: float, nx_n: int, ny_n: int) -> int: + """The collected image nearest a fractional grid position, for addressing.""" + iy, ix = int(round(np.clip(gy, 0, ny_n - 1))), int(round(np.clip(gx, 0, nx_n - 1))) + num = int(F["number"][iy, ix]) + if num < 0: # that cell was never collected + yy, xx = np.nonzero(F["number"] >= 0) + k = int(np.argmin((yy - gy) ** 2 + (xx - gx) ** 2)) + num = int(F["number"][yy[k], xx[k]]) + return num + + +def _object_centre(bkg: np.ndarray, level: float): + """Centre of the sample itself, from the background alone. + + The largest connected patch of raised background, weighted by how far each cell + sits above the object level. This is where to point when the protein channel + has nothing to say: a loop with no useful diffraction still has a loop, and its + middle is a far better guess than whichever cell the noise happened to peak on. + Returns None when the background never rises -- an empty scan has no centre. + """ + obj = np.isfinite(bkg) & (bkg > level) + if not obj.any(): + return None + blobs = _blobs(obj) + if not blobs: + return None + w = np.where(obj, bkg - level, 0.0) + m = max(blobs, key=lambda c: float(w[c].sum())) + yy, xx = np.nonzero(m) + wt = w[yy, xx] + tot = float(wt.sum()) + if tot <= 0: + return float(xx.mean()), float(yy.mean()) + return float((xx * wt).sum() / tot), float((yy * wt).sum() / tot) + + def _blobs(mask: np.ndarray, limit: int = 64): """Every separate blob the level encloses, largest-signal choice left to caller.""" rest, out = mask.copy(), [] @@ -370,7 +497,7 @@ _BEAM_PATHS = ( ) -def _beam_um(obj: BeamSource | None) -> tuple[float | None, float | None]: +def _beam_um(obj: Optional[BeamSource]) -> tuple[Optional[float], Optional[float]]: """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 @@ -404,7 +531,7 @@ def _beam_um(obj: BeamSource | None) -> tuple[float | None, float | None]: return None, None -def _deconvolve(fwhm: float | None, beam: float | None) -> float | None: +def _deconvolve(fwhm: Optional[float], beam: Optional[float]) -> Optional[float]: """Quadrature removal of the beam. None when the beam swallows the feature.""" if fwhm is None or beam is None: return fwhm @@ -420,6 +547,16 @@ def _grids(images: Sequence[ScanResultImagesInner]): nx/ny come from the scan_result itself, so the snake ordering, the vertical flag and the step signs are already applied upstream and are not redone here. + + The protein channel is spots_low_res: Bragg spots beyond 5 A. It is the one + count that only a large unit cell can produce. Ice cannot -- its innermost + ring is at 3.9 A -- and neither can salt, whose few reflections at that + spacing do not register as a count. It is therefore immune to the failure + that ruled out every alternative: on a loop caked in ice, hundreds of ring + spots per frame leak past the ice windows into "spots - spots_ice", and the + haze that leaves follows the ice, not the crystal, and can be compact enough + to pass any spatial test at the wrong place. The remaining high-angle count + is kept as "hires" for the resolution it implies, never for the position. """ nx = np.array([_get(i, "nx", -1) for i in images], dtype=np.int64) ny = np.array([_get(i, "ny", -1) for i in images], dtype=np.int64) @@ -431,20 +568,19 @@ def _grids(images: Sequence[ScanResultImagesInner]): return np.array([_get(i, key, np.nan) for i in images], dtype=dtype) spots, low, ice = col("spots"), col("spots_low_res"), col("spots_ice") - # No total spot count: fall back to spots above the low-resolution cut. That - # channel is ice-free on the same argument (the lowest hexagonal ice ring is - # 3.895 A) and tracks the real one at r = 0.79-0.98 on our data. - if np.isnan(spots).all(): - protein, channel = np.nan_to_num(low), "spots_low_res (no total count)" - else: - protein = np.nan_to_num(spots) - np.nan_to_num(ice) - np.nan_to_num(low) - channel = "spots - spots_ice - spots_low_res" + protein, channel = np.nan_to_num(low), "spots_low_res (beyond 5 A)" + hires = ( + np.clip(np.nan_to_num(spots) - np.nan_to_num(ice) - np.nan_to_num(low), 0, None) + if not np.isnan(spots).all() + else np.full_like(protein, np.nan) + ) protein = np.clip(protein, 0, None) F = {} for name, v in ( ("protein", protein), ("ice", np.nan_to_num(ice)), + ("hires", hires), ("bkg", col("bkg")), ("res", col("res")), ("ring", col("ice")), @@ -462,11 +598,51 @@ def _grids(images: Sequence[ScanResultImagesInner]): return F, channel -def _classify(F: dict, th: Thresholds) -> Counts: +def _spot_floor(F: dict, bkg: np.ndarray, level: float, th: Thresholds): + """Spots a frame must clear, measured from the scan's own off-object cells. + + A fixed threshold assumes the channel reads about zero where there is no + sample. The real channel does; the spots_low_res fallback does not -- it sits + at a median of 11-13 in empty air on every scan tested, because low-resolution + peak finding has a false-positive rate there. A fixed 10 then marks ~85 % of + such a scan as protein. Measuring the floor from the cells the background gate + rejects fixes that, and where the channel really is clean the estimate falls + below min_spots and nothing changes. + + The scale comes from a one-sided MAD -- deviations below the median only. + "Off-object" is a coarse gate, so on a line scan of 30-odd points a few cells + at the crystal's edge land in it and sit in the upper tail; a two-sided MAD or + a high percentile then follows them (a 99th percentile gave 61 spots on one + line whose raster said 21). The lower half cannot be contaminated that way. + + Returns (floor, off-object median). The median is reported rather than used: + it is what tells you the noisy fallback channel is in play. + """ + p = F["protein"] + off = np.isfinite(p) & np.isfinite(bkg) & (bkg <= level) + v = p[off] + med = float(np.median(v)) if v.size else None + if th.baseline_sigma is None or v.size < th.baseline_min_cells: + # Too few off-object cells to measure a baseline -- a crystal filling the + # scan would otherwise set its own threshold from its own signal. + return float(th.min_spots), med + # The pedestal's scatter is its counting noise, sqrt(median). A MAD was tried + # and failed on scans of a crystal larger than the grid: the off-object cells + # there are half true empties at zero and half crystal fringe at tens, and a + # one-sided MAD of that mixture returns the median itself, lifting the floor + # to eight times the pedestal and hiding a 250 um crystal. Poisson scatter + # reproduces the floors measured on the pedestal sessions (11 -> 28 against + # the 26-27 the MAD gave) with one assumption fewer. A half-count minimum + # keeps an all-zero field from setting its floor at zero. + base = med + th.baseline_sigma * max(float(np.sqrt(max(med, 0.0))), 0.5) + return max(float(th.min_spots), float(base)), med + + +def _classify(F: dict, th: Thresholds, floor: float) -> Counts: """Protein / protein+ice / ice / noise / blank, per frame.""" p, i_, s = F["protein"], F["ice"], F["spots"] ok = np.isfinite(p) - prot = ok & (p >= th.min_spots) + prot = ok & (p >= floor) ice = ok & (i_ >= th.ice_spots) if th.ice_ring is not None: ice |= ok & np.isfinite(F["ring"]) & (F["ring"] >= th.ice_ring) @@ -487,11 +663,11 @@ def _classify(F: dict, th: Thresholds) -> Counts: def analyse( scan_result: ScanResult, - grid_scan: GridScan | None = None, + grid_scan: Optional[GridScan] = None, *, - daq: BeamSource | None = None, - thresholds: Thresholds | None = None, - z_um: float | None = None, + daq: Optional[BeamSource] = None, + thresholds: Optional[Thresholds] = None, + z_um: Optional[float] = None, jpeg: bool = True, ) -> GridScanResult: """Analyse one grid scan. @@ -537,7 +713,19 @@ def analyse( ) F, channel = built ny_n, nx_n = F["protein"].shape - counts = _classify(F, th) + # The object level comes first now: the spot floor is derived from the cells + # this rejects, so it cannot be computed before it. + bkg = F["bkg"] + if np.isfinite(bkg).any(): + lo, hi = np.nanpercentile(bkg, 5), np.nanpercentile(bkg, 95) + level = lo + 0.25 * (hi - lo) + else: + bkg, level = np.ones_like(F["protein"]), 0.0 + floor, off_med = _spot_floor(F, bkg, level, th) + + counts = _classify(F, th, floor) + counts.spot_floor = round(floor, 1) + counts.off_object_median = None if off_med is None else round(off_med, 1) # The step is signed -- negative means the grid runs right to left, or bottom to # top. A size only wants the magnitude, but an offset that drops the sign sends # the goniometer the wrong way, so both are kept. @@ -550,36 +738,81 @@ def analyse( ) peak = float(np.nanmax(F["protein"])) - if not np.isfinite(peak) or peak < th.min_spots: - # The 50 % contour of a channel peaking at two spots is a contour of - # noise: say there is no crystal rather than centre on nothing. - out.reason = f"peak protein channel {peak:.0f} < {th.min_spots} spots" + counts.contrast = round(peak / floor, 2) if floor > 0 and np.isfinite(peak) else None + + def _fall_back(reason: str) -> GridScanResult: + """No crystal to centre on: hand back the middle of the object instead. + + The scan still cost the time and the sample is still in the beam, so the + useful answer is where the loop is, not silence. found stays False -- this + is not a crystal and no strategy should be built on it -- but centre and + centre_source are filled so the goniometer can be driven there and the loop + rescanned or reoriented. resolution_A is deliberately left None: whatever + the frame there reports is ice or air, and feeding it to a dose budget + would turn a non-detection into a confident wrong number. + """ + out.reason = reason + oc = _object_centre(bkg, level) + if oc is not None: + gx0, gy0 = oc + out.centre_source = "background" + out.centre = Centre( + nx=round(float(gx0), 3), + ny=round(float(gy0), 3), + image_number=_nearest_image(F, gx0, gy0, nx_n, ny_n), + resolution_A=None, + x_um=round(gx0 * sx_s, 2) if sx_s else None, + y_um=round(gy0 * sy_s, 2) if sy_s else None, + ) + out.warnings.append( + "no convincing protein diffraction; centred on the middle of the " + "object, not on a crystal" + ) if jpeg: - _pictures(F, None, None, out, th, 0.0) + _pictures( + F, + None, + (out.centre.nx, out.centre.ny) if out.centre else None, + out, + th, + level, + floor, + ) return out - # Object mask from the background estimate: this finds the loop, not the crystal. - bkg = F["bkg"] - if np.isfinite(bkg).any(): - lo, hi = np.nanpercentile(bkg, 5), np.nanpercentile(bkg, 95) - level = lo + 0.25 * (hi - lo) - else: - bkg, level = np.ones_like(F["protein"]), 0.0 + if not np.isfinite(peak) or peak < floor: + # The 50 % contour of a channel peaking at the noise floor is a contour of + # noise: say there is no crystal rather than centre on nothing. + return _fall_back( + f"peak protein channel {peak:.0f} < {floor:.0f} spots" + + ( + f" (floor from the off-object baseline of {off_med:.0f})" + if off_med is not None and floor > th.min_spots + else "" + ) + ) + if counts.contrast is not None and counts.contrast < th.contrast_min: + # Over the floor, but not by enough to mean anything. The 50 % contour of + # such a map lands wherever the noise happened to peak, which is what made + # weak scans centre in an arbitrary place; the object centre does not move + # with the noise. + return _fall_back( + f"peak protein channel {peak:.0f} is only {counts.contrast:.2f}x the " + f"{floor:.0f}-spot floor; no cell stands out as a crystal" + ) line = nx_n == 1 or ny_n == 1 factor = 24 if line else int(np.clip(round(512 / max(nx_n, ny_n)), 4, 24)) P = _upsample(np.nan_to_num(F["protein"]), factor) B = _upsample(np.nan_to_num(bkg, nan=float(np.nanmin(bkg))), factor) - obj_gate = B > level - if th.object_union_spots: - # Bragg spots are direct evidence of sample; bkg is only a proxy - obj_gate |= P >= th.min_spots - inside = (P >= _LEVEL * peak) & obj_gate + # The contour cannot dip below the noise floor. At 50 % of a 28-spot peak the + # level is 14 spots, under a floor of 26 -- the region would enclose cells the + # classifier has already called noise. Where the crystal is strong the 50 % + # level is far above the floor and this changes nothing. + contour_level = max(_LEVEL * peak, floor) + inside = (P >= contour_level) & (B > level) if not inside.any(): - out.reason = "no region above the 50 % level inside the object" - if jpeg: - _pictures(F, None, None, out, th, level) - return out + return _fall_back("no region above the 50 % level inside the object") my, mx = inside.shape # Sample the contour back onto the collected cells, then keep only the blob @@ -588,6 +821,7 @@ def analyse( ci = np.clip(((np.arange(ny_n) + 0.5) / ny_n * my).astype(int), 0, my - 1) cj = np.clip(((np.arange(nx_n) + 0.5) / nx_n * mx).astype(int), 0, mx - 1) cell = inside[np.ix_(ci, cj)] + n_above = int(cell.sum()) # before the blob choice narrows it blobs = _blobs(cell) n_frag = len(blobs) if blobs: @@ -608,9 +842,29 @@ 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)] + counts.n_contour_cells = int(cell.sum()) + counts.contour_compactness = round(float(cell.sum()) / n_above, 3) if n_above else None + # The spatial test. Everything above is a statistic of the frames taken one at a + # time, and a permutation control -- the same frames dealt out to different + # cells -- passes all of them at the same rate as the real data, because dealing + # them elsewhere changes no per-frame number. What it does destroy is + # contiguity, and that is what a crystal has: its hot cells touch. A map whose + # hot cells are scattered over the loop is noise however hot they are. + # Two large patches are not scatter: a loop often carries two crystals, or one + # long one that the 50 % level cuts in two at its waist. The gate is for many + # small pieces, so a split into at most two, the larger of them a real patch + # of cells, is let through and centred on the stronger. + two_crystals = n_frag <= 2 and counts.n_contour_cells >= th.two_crystal_cells + if ( + counts.contour_compactness is not None + and counts.contour_compactness < th.compactness_min + and not two_crystals + ): + return _fall_back( + f"the {n_above} cells over the contour level fall into {n_frag} " + f"separate patches (largest holds {counts.n_contour_cells}); signal " + "this scattered is not one crystal" + ) if line: # One position wide: the 50 % crossings either side of the peak. The outer @@ -619,11 +873,8 @@ def analyse( prof = P[:, 0] if vert else P[0, :] base = B[:, 0] if vert else B[0, :] pos = (np.arange(len(prof)) + 0.5) / len(prof) * (ny_n if vert else nx_n) - line_gate = base > level - if th.object_union_spots: - line_gate |= prof >= th.min_spots - v = np.where(line_gate, prof, 0.0) - c = _crossings(pos, v, _LEVEL * peak) + v = np.where(base > level, prof, 0.0) + c = _crossings(pos, v, contour_level) top = pos[int(np.argmax(v))] lo = [x for x in c if x <= top] hi = [x for x in c if x >= top] @@ -641,7 +892,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 = round(np.clip(gy, 0, ny_n - 1)), round(np.clip(gx, 0, nx_n - 1)) + iy, ix = int(round(np.clip(gy, 0, ny_n - 1))), int(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) @@ -656,10 +907,14 @@ def analyse( if cell[0, :].any() or cell[-1, :].any() or cell[:, 0].any() or cell[:, -1].any(): out.warnings.append("contour touches the scan edge; crystal may extend beyond it") + res_here = F["res"][iy, ix] out.centre = Centre( nx=round(float(gx), 3), ny=round(float(gy), 3), image_number=num, + resolution_A=( + round(float(res_here), 2) if np.isfinite(res_here) and res_here > 0 else None + ), x_um=round(gx * sx_s, 2) if sx_s else None, y_um=round(gy * sy_s, 2) if sy_s else None, ) @@ -696,10 +951,23 @@ def analyse( # what sits under the contour if cell.any(): counts.mean_protein_in_contour = round(float(np.nanmean(F["protein"][cell])), 1) - r = F["res"][cell & np.isfinite(F["res"])] + # Zero and negative are Jungfraujoch's "no estimate", and min() would take + # them as the best resolution on the loop. + r = F["res"][cell & np.isfinite(F["res"]) & (F["res"] > 0)] if r.size: counts.best_resolution_A = round(float(r.min()), 2) counts.ice_in_contour = round(float(np.mean(F["ice"][cell] >= th.ice_spots)), 3) + # The high-angle count is not used for the position, but where it peaks + # is worth knowing: on the crystal it is the crystal's own high-resolution + # diffraction; off it, it is ice leaking past the ice windows. + hi = F["hires"] + if np.isfinite(hi).any() and np.nanmax(hi) > 0: + hy, hx = np.unravel_index(int(np.nanargmax(hi)), hi.shape) + if not cell[hy, hx]: + out.warnings.append( + f"high-angle spots peak at {hx},{hy}, outside the crystal; that " + "is ice, and any resolution estimate there is the ice's" + ) if th.check_unit_cell: out.unit_cell, counts.unit_cell_agreement = _cell_agreement(F["uc"][cell], th.cell_tol) if ( @@ -720,12 +988,22 @@ def analyse( out.warnings.append( f"{n_frag} separate regions at the 50 % level; centred on the strongest" ) + elif ( + n_frag == 2 + and counts.contour_compactness is not None + and counts.contour_compactness < th.compactness_min + ): + out.warnings.append( + "two comparable regions at the 50 % level -- two crystals, or one cut at " + "its waist; centred on the stronger, the other is not measured" + ) frac = cell.sum() / max(int(np.isfinite(F["protein"]).sum()), 1) if frac > 0.4: out.warnings.append(f"contour covers {frac:.0%} of the scan; likely not a single crystal") out.found = True + out.centre_source = "protein" if jpeg: - _pictures(F, inside, (gx, gy), out, th, level) + _pictures(F, inside, (gx, gy), out, th, level, floor) return out @@ -743,16 +1021,31 @@ 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, + floor: float = _COLOUR_FLOOR, +) -> 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]) + res.jpeg = _render(F, inside, centre, res, th, obj_level, want[0], floor) if len(want) > 1: - res.jpeg_cells = _render(F, inside, centre, res, th, obj_level, want[1]) + res.jpeg_cells = _render(F, inside, centre, res, th, obj_level, want[1], floor) 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, + floor: float = _COLOUR_FLOOR, ) -> bytes: """Map, contours, crosshair and a one-line caption, as JPEG. @@ -789,10 +1082,17 @@ def _render( else (np.nan_to_num(F["bkg"]), np.nan_to_num(F["protein"]), np.nan_to_num(F["ice"])) ) obj = _norm(up_b, ref=F["bkg"]) if np.isfinite(F["bkg"]).any() else np.ones_like(up_b) - p = _norm(up_p, _COLOUR_FLOOR, ref=F["protein"]) * obj + # The colour scale uses the same floor as the classification, so a channel + # sitting on a noise pedestal does not paint the whole loop orange. + p = _norm(up_p, floor, ref=F["protein"]) * obj i_ = _norm(up_i, _COLOUR_FLOOR, ref=F["ice"]) * obj base = _PLATE * (1 - obj[..., None]) + _GREY * obj[..., None] + # Protein takes precedence: ice is shown only where protein is not. On an + # iced loop every frame carries ice, and an even blend would paint the crystal + # the same muddy blue-brown as the film around it. Here the crystal reads + # orange and the ice reads blue around it, which is the picture to act on. + i_ = i_ * (1 - p) tot = np.clip(p + i_, 0, 1)[..., None] mix = (_ORANGE * p[..., None] + _BLUE * i_[..., None]) / np.maximum(p + i_, 1e-9)[..., None] rgb = np.clip(base * (1 - tot) + mix * tot, 0, 1)