mx/AareCommon
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 29 Wiki Activity

fix: update gridscan analysis script #37

Merged
perl_d merged 1 commits from fix/update_gridscan_analysis into main 2026-09-16 13:10:47 +02:00
Conversation 0 Commits 1 Files Changed 1
+415 -115
1 changed files with 415 additions and 115 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/aarecommon/math/jfjoch_gridscan_union.py
+415 -115
View File
@@ -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)