Files
Jungfraujoch/image_analysis/IceScore.h
T
leonarski_fandClaude Opus 5 62ca2b24f6 image analysis: two per-image detection scores say whether there is protein and whether there is ice
proteinScore and iceScore are saturating scores in [0,1] that answer presence, not quality: a superb
crystal and a barely-diffracting one both read near 1, and neither the spot count nor the resolution
enters either of them as a term. iceRingScore already reports a magnitude - a ratio, unbounded - which
is not a number that can be thresholded; these are.

The protein score counts distinct d SHELLS above 5 A rather than spots, so a parasitic ring in the
low-resolution band cannot accumulate evidence, and weights each spot against the frame's own median
so a scattering of the weakest detections cannot fill a shell either. The ice score carries two
channels and takes the stronger: a radial one over the azimuthal profile, which runs the hexagonal
and the CUBIC phase as separate hypotheses and decides between them at the end (flash-cooled loops
show cubic or stacking-disordered ice at least as often as hexagonal, the two share only three
lines, and dropping the cubic hypothesis costs about 5 pp on iced loops), and a spot one that reads
an excess on the ice radii against the same band slid to every ice-free offset, which is what
catches ice arriving as discrete crystallites and leaving the radial profile flat.

Both read d out of the geometry, so both move with a beam-centre error; the centre is not fitted
here, and the one they were computed with is written beside them as scoreBeamCenterX/Y so a later
rescoring can tell an algorithm disagreement from a geometry one.

Ported from validated prototypes and checked against them frame by frame on stored data: mean
absolute difference 2.7e-5 (protein), 1.3e-8 (ice radial) and 3.4e-4 (ice spots). On a 41-loop
battery the protein score reaches 98.4% of confirmed-protein frames and 0.00% of water frames, and
finds no cluster on any water or ice raster. Cost is 0.01 ms/frame for the protein score and
0.08-0.32 ms/frame for the ice score.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
2026-09-07 23:59:46 +02:00

62 lines
4.2 KiB
C++

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#pragma once
#include <vector>
#include "../common/AzimuthalIntegrationSettings.h"
#include "../common/SpotToSave.h"
// Is there CRYSTALLINE ICE on this image? A detection score in [0,1] that saturates: a loop buried in
// ice and one carrying a single detectable ring both come out near 1. Unlike ice_ring_score - which is
// a ratio, unbounded, and answers "how strong is the worst ring" - this answers only "is ice present",
// and it is the number to threshold.
//
// Ice reaches the frame two ways and they need different evidence, so two channels are computed and the
// stronger one wins. Neither is a subset of the other: fine polycrystalline ice makes smooth powder
// rings and leaves few extra spots, while ice in large crystallites makes discrete spots on the same
// radii and leaves the radial profile flat.
//
// Both read d from the geometry, so both move with a beam-centre error. The centre is not fitted here -
// that belongs to geometry refinement - and the centre the scores were computed with is written beside
// them in the file.
// Channel 1 - the radial profile. Two ice phases are carried as separate hypotheses and the decision is
// taken at the end (the larger score wins), because flash-cooled loops show cubic or stacking-disordered
// ice at least as often as hexagonal and the two phases share only three lines. Each band is read as a
// standardised excess over a running median, in units of the bin mean's own error, and is compared with
// the same statistic measured on every profile bin that belongs to no band of either phase - so a grainy
// profile raises its own null as much as its own band values. Two statistics are formed against that
// null, an amplitude and a band-count concordance, and the SMALLER is taken: a single elevated bin then
// fails, because real ice shows a whole pattern.
//
// profile / profile_std / profile_count are q_bins long, or q_bins x azimuthal bins, in which case the
// first two are averaged over azimuth and the third summed. The scale of an excess is the bin MEAN's own
// error, std / sqrt(count) - the profile is a mean of many pixels, so its plain standard deviation is
// the wrong yardstick by two orders of magnitude. Returns 0 when profile_std carries nothing usable: the
// FPGA azimuthal integration does not produce one (its profile can be recomputed on the CPU -
// ForceCPUinFPGAWorkflow - which does).
float IceScoreRadial(const std::vector<float> &profile, const std::vector<float> &profile_std,
const std::vector<uint64_t> &profile_count, int32_t q_bins,
const AzimuthalIntegrationSettings &settings);
// Channel 2 - the spot population. Evidence is an EXCESS of found spots on the hexagonal-ice radii over
// what this frame's own radial spot density predicts. The null is not the two flanks either side of each
// band (a ratio of two ~1-count numbers, which is what spot_count_ice_control is) but each band slid to
// every ice-free offset within +-0.45 A^-1, in +-delta pairs so that the fall-off of spot density with q
// cancels to first order, each count divided by the live detector area at that radius. That turns a
// 1-count control into an average over ~100 of them. The excess is then read as a quasi-Poisson upper
// tail AND as a ratio, and the smaller of the two is taken: the tail alone fires on a 10 % band
// enrichment when a frame has 800 spots, and the ratio alone fires on 2 spots out of 2.
//
// profile_count is the azimuthal integration's live pixel count per bin, q_bins long or q_bins x
// azimuthal bins; it is what makes the offsets comparable where the detector edge cuts a radius short.
float IceScoreSpots(const std::vector<SpotToSave> &spots, const std::vector<uint64_t> &profile_count,
int32_t q_bins, const AzimuthalIntegrationSettings &settings);
// The score itself: whichever channel sees more.
float IceScore(const std::vector<float> &profile, const std::vector<float> &profile_std,
const std::vector<uint64_t> &profile_count, int32_t q_bins,
const AzimuthalIntegrationSettings &settings, const std::vector<SpotToSave> &spots);