Jungfraujoch/image_analysis/IceScore.h

// 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_ratio - 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.
//
// This channel takes no band width: it reads the single bin the geometry predicts, searched over a small
// tolerance for the beam-centre error (see CENTRE_SMEAR_Q in the .cpp). The one band width in this file
// is the spot channel's, and it is the pipeline's.
//
// 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.
// half_width_q is the band half-width, the same one the spot finder marks ice rings with
// (SpotFindingSettings::ice_ring_width_Q_recipA) - one width, measured over the corpus, not two.
float IceScoreSpots(const std::vector<float> &spot_d_A, const std::vector<uint64_t> &profile_count,
                    int32_t q_bins, const AzimuthalIntegrationSettings &settings, float half_width_q);

// The score itself: whichever channel sees more. See the .cpp for why the maximum and not something
// gentler, and for the lever to reach for if the rate ever has to come down.
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<float> &spot_d_A,
               float half_width_q);