Merge branch 'worktree-agent-aa55c81220288558b' into worktree-agent-ab3a3e173261a95d8
This commit is contained in:
@@ -899,6 +899,8 @@ PlotType ConvertPlotType(const std::optional<std::string>& input) {
|
||||
"Plot type is compulsory paramater");
|
||||
if (input == "bkg_estimate") return PlotType::BkgEstimate;
|
||||
if (input == "ice_ring_score") return PlotType::IceRingScore;
|
||||
if (input == "protein_score") return PlotType::ProteinScore;
|
||||
if (input == "ice_score") return PlotType::IceScore;
|
||||
if (input == "azint") return PlotType::AzInt;
|
||||
if (input == "azint_1d") return PlotType::AzInt1D;
|
||||
if (input == "spot_count") return PlotType::SpotCount;
|
||||
@@ -1149,6 +1151,10 @@ org::openapitools::server::model::Scan_result Convert(const ScanResult& input) {
|
||||
tmp.setSpotsIce(i.spot_count_ice.value());
|
||||
if (i.ice_ring_score.has_value())
|
||||
tmp.setIce(i.ice_ring_score.value());
|
||||
if (i.protein_score.has_value())
|
||||
tmp.setProteinScore(i.protein_score.value());
|
||||
if (i.ice_score.has_value())
|
||||
tmp.setIceScore(i.ice_score.value());
|
||||
if (i.spot_count_low_res.has_value())
|
||||
tmp.setSpotsLowRes(i.spot_count_low_res.value());
|
||||
if (i.spot_count_indexed.has_value())
|
||||
|
||||
@@ -121,6 +121,8 @@ components:
|
||||
- image_scale_cc
|
||||
- compression_ratio
|
||||
- ice_ring_score
|
||||
- protein_score
|
||||
- ice_score
|
||||
roi:
|
||||
in: query
|
||||
name: roi
|
||||
@@ -1661,6 +1663,14 @@ components:
|
||||
type: number
|
||||
format: float
|
||||
description: Strongest hexagonal-ice ring intensity over the smooth radial background (1 = no ice)
|
||||
protein_score:
|
||||
type: number
|
||||
format: float
|
||||
description: Protein diffraction detection score (0 = none, 1 = certain); saturating, not a quality measure
|
||||
ice_score:
|
||||
type: number
|
||||
format: float
|
||||
description: Crystalline ice detection score (0 = none, 1 = certain); saturating, not a quality measure
|
||||
index:
|
||||
type: integer
|
||||
format: int64
|
||||
|
||||
@@ -83,6 +83,15 @@ constexpr std::array<float, 19> ICE_RING_RES_A = {3.895, 3.661, 3.438, 2.667, 2.
|
||||
1.472, 1.443, 1.371, 1.366, 1.298, 1.261, 1.224,
|
||||
1.170};
|
||||
|
||||
// Cubic-ice (Ic) ring positions, d = a / sqrt(h^2+k^2+l^2) with a = 6.358 A and the diamond-lattice
|
||||
// reflection conditions (hkl all odd, or all even with h+k+l = 4n). Flash-cooled loops show cubic or
|
||||
// stacking-disordered ice at least as often as hexagonal, and the two phases share only three lines,
|
||||
// so a detector that only looks for ICE_RING_RES_A above misses it. Used by the ice detection score
|
||||
// as a second, independent hypothesis - not by the spot ice-ring flag, which stays hexagonal because
|
||||
// masking on a phase that is usually absent would throw away good reflections.
|
||||
constexpr std::array<float, 9> ICE_RING_CUBIC_RES_A = {3.670, 2.248, 1.917, 1.835, 1.590,
|
||||
1.459, 1.298, 1.224, 1.124};
|
||||
|
||||
// True when resolution d (Angstrom) sits within half_width of a hexagonal-ice powder ring, in the
|
||||
// q = 2*pi/d units the spot-finder uses (ice_ring_width_Q_recipA). Used to drop ice-contaminated
|
||||
// reflections from scaling/merging when ice-ring handling is enabled.
|
||||
|
||||
@@ -129,6 +129,16 @@ struct DataMessage {
|
||||
std::optional<float> bkg_estimate;
|
||||
std::optional<float> ice_ring_score; // strongest ice ring over the smooth radial background (1 = none)
|
||||
|
||||
// Two DETECTION scores in [0,1], saturating: is there protein diffraction on this image, and is
|
||||
// there crystalline ice. Not quality measures - resolution, B factor and the indexing rate say how
|
||||
// good the diffraction is, and neither the spot count nor the resolution enters either score as a
|
||||
// term. ice_score answers the question ice_ring_score above only reports on: it combines the radial
|
||||
// powder evidence with a spot excess on the ice radii, and it saturates, so it can be thresholded.
|
||||
// Both are computed from the geometry's beam centre; see image_analysis/IceScore.h and
|
||||
// image_analysis/spot_finding/SpotUtils.h.
|
||||
std::optional<float> protein_score;
|
||||
std::optional<float> ice_score;
|
||||
|
||||
std::optional<bool> indexing_result;
|
||||
std::optional<CrystalLattice> indexing_lattice;
|
||||
std::vector<CrystalLattice> indexing_extra_lattices;
|
||||
@@ -369,6 +379,10 @@ struct EndMessage {
|
||||
std::optional<float> efficiency;
|
||||
std::optional<float> indexing_rate;
|
||||
std::optional<float> bkg_estimate;
|
||||
// Run means of the two per-image detection scores; written to /entry/MX/proteinScoreMean and
|
||||
// /entry/MX/iceScoreMean.
|
||||
std::optional<float> protein_score;
|
||||
std::optional<float> ice_score;
|
||||
|
||||
std::optional<std::string> end_date;
|
||||
|
||||
@@ -419,6 +433,8 @@ struct EndMessage {
|
||||
std::vector<uint8_t> image_indexed;
|
||||
std::vector<int32_t> indexed_lattice_count;
|
||||
std::vector<float> v_bkg_estimate;
|
||||
std::vector<float> v_protein_score;
|
||||
std::vector<float> v_ice_score;
|
||||
std::vector<float> profile_radius;
|
||||
std::vector<float> mosaicity;
|
||||
std::vector<float> bFactor;
|
||||
|
||||
@@ -62,6 +62,8 @@ void JFJochReceiverPlots::Setup(const DiffractionExperiment &experiment, const A
|
||||
}
|
||||
bkg_estimate.Clear(r);
|
||||
ice_ring_score.Clear(r);
|
||||
protein_score.Clear(r);
|
||||
ice_score.Clear(r);
|
||||
spot_count.Clear(r);
|
||||
spot_count_low_res.Clear(r);
|
||||
spot_count_indexed.Clear(r);
|
||||
@@ -129,6 +131,8 @@ void JFJochReceiverPlots::Setup(const DiffractionExperiment &experiment, const A
|
||||
void JFJochReceiverPlots::Add(const DataMessage &msg, const AzimuthalIntegrationProfile &profile) {
|
||||
bkg_estimate.AddElement(msg.number, msg.bkg_estimate);
|
||||
ice_ring_score.AddElement(msg.number, msg.ice_ring_score);
|
||||
protein_score.AddElement(msg.number, msg.protein_score);
|
||||
ice_score.AddElement(msg.number, msg.ice_score);
|
||||
resolution_estimate.AddElement(msg.number, msg.resolution_estimate);
|
||||
spot_count.AddElement(msg.number, msg.spot_count);
|
||||
spot_count_low_res.AddElement(msg.number, msg.spot_count_low_res);
|
||||
@@ -280,6 +284,12 @@ MultiLinePlot JFJochReceiverPlots::GetPlots(const PlotRequest &request) {
|
||||
case PlotType::IceRingScore:
|
||||
ret = ice_ring_score.GetMeanPlot(nbins, start, incr, request.fill_value);
|
||||
break;
|
||||
case PlotType::ProteinScore:
|
||||
ret = protein_score.GetMeanPlot(nbins, start, incr, request.fill_value);
|
||||
break;
|
||||
case PlotType::IceScore:
|
||||
ret = ice_score.GetMeanPlot(nbins, start, incr, request.fill_value);
|
||||
break;
|
||||
case PlotType::ResolutionEstimate:
|
||||
ret = resolution_estimate.GetMeanPlot(nbins, start, incr, request.fill_value);
|
||||
break;
|
||||
@@ -509,6 +519,30 @@ std::vector<float> JFJochReceiverPlots::GetIceRingScoreArray() const {
|
||||
return ice_ring_score.ExportArray();
|
||||
}
|
||||
|
||||
std::optional<float> JFJochReceiverPlots::GetProteinScore() const {
|
||||
auto tmp = protein_score.Mean();
|
||||
if (std::isfinite(tmp))
|
||||
return tmp;
|
||||
else
|
||||
return {};
|
||||
}
|
||||
|
||||
std::optional<float> JFJochReceiverPlots::GetIceScore() const {
|
||||
auto tmp = ice_score.Mean();
|
||||
if (std::isfinite(tmp))
|
||||
return tmp;
|
||||
else
|
||||
return {};
|
||||
}
|
||||
|
||||
std::vector<float> JFJochReceiverPlots::GetProteinScoreArray() const {
|
||||
return protein_score.ExportArray();
|
||||
}
|
||||
|
||||
std::vector<float> JFJochReceiverPlots::GetIceScoreArray() const {
|
||||
return ice_score.ExportArray();
|
||||
}
|
||||
|
||||
MeanProcessingTime JFJochReceiverPlots::GetMeanProcessingTime() const {
|
||||
MeanProcessingTime ret{};
|
||||
ret.compression = compression_time.Mean();
|
||||
@@ -576,6 +610,12 @@ void JFJochReceiverPlots::GetPlotRaw(std::vector<float> &v, PlotType type, const
|
||||
case PlotType::IceRingScore:
|
||||
v = ice_ring_score.ExportArray();
|
||||
break;
|
||||
case PlotType::ProteinScore:
|
||||
v = protein_score.ExportArray();
|
||||
break;
|
||||
case PlotType::IceScore:
|
||||
v = ice_score.ExportArray();
|
||||
break;
|
||||
case PlotType::ResolutionEstimate:
|
||||
v = resolution_estimate.ExportArray();
|
||||
break;
|
||||
|
||||
@@ -43,6 +43,8 @@ class JFJochReceiverPlots {
|
||||
|
||||
StatusVector bkg_estimate;
|
||||
StatusVector ice_ring_score;
|
||||
StatusVector protein_score;
|
||||
StatusVector ice_score;
|
||||
StatusVector spot_count;
|
||||
StatusVector spot_count_low_res;
|
||||
StatusVector spot_count_indexed;
|
||||
@@ -131,6 +133,11 @@ public:
|
||||
// the run actually is.
|
||||
[[nodiscard]] std::optional<float> GetResolutionEstimate() const;
|
||||
std::vector<float> GetIceRingScoreArray() const;
|
||||
// Run means of the two per-image detection scores, and the per-image arrays behind them.
|
||||
[[nodiscard]] std::optional<float> GetProteinScore() const;
|
||||
[[nodiscard]] std::optional<float> GetIceScore() const;
|
||||
[[nodiscard]] std::vector<float> GetProteinScoreArray() const;
|
||||
[[nodiscard]] std::vector<float> GetIceScoreArray() const;
|
||||
|
||||
std::vector<float> GetAzIntProfile() const;
|
||||
// The run-summed profile object itself, rather than the array GetAzIntProfile() flattens it to.
|
||||
|
||||
+2
-1
@@ -15,7 +15,8 @@ enum class PlotType {
|
||||
ROISum, ROIMean, ROIMaxCount, ROIPixels, ROIWeightedX, ROIWeightedY, PacketsReceived, MaxValue,
|
||||
ResolutionEstimate, ProfileRadius, Mosaicity, BFactor, PixelSum, StrongPixels,
|
||||
RefinementBeamX, RefinementBeamY, ImageProcessingTime, IntegratedReflections,
|
||||
ImageScaleFactor, ImageScaleCC, CompressionRatio, IndexingLatticeCount, IceRingScore
|
||||
ImageScaleFactor, ImageScaleCC, CompressionRatio, IndexingLatticeCount, IceRingScore,
|
||||
ProteinScore, IceScore
|
||||
};
|
||||
|
||||
enum class PlotAzintUnit {
|
||||
|
||||
@@ -64,6 +64,8 @@ void ScanResultGenerator::Add(const DataMessage &message) {
|
||||
v[image_number].image_scale_factor = message.image_scale_factor;
|
||||
v[image_number].image_scale_cc = message.image_scale_cc;
|
||||
v[image_number].ice_ring_score = message.ice_ring_score;
|
||||
v[image_number].protein_score = message.protein_score;
|
||||
v[image_number].ice_score = message.ice_score;
|
||||
if (message.lattice_type)
|
||||
v[image_number].niggli_class = message.lattice_type->niggli_class;
|
||||
}
|
||||
@@ -110,6 +112,8 @@ void ScanResultGenerator::FillEndMessage(EndMessage &message) const {
|
||||
message.image_scale_factor.resize(n);
|
||||
message.image_scale_cc.resize(n);
|
||||
message.ice_ring_score.resize(n);
|
||||
message.v_protein_score.resize(n);
|
||||
message.v_ice_score.resize(n);
|
||||
message.integrated_reflections.resize(n);
|
||||
message.niggli_class.resize(n);
|
||||
message.pixel_sum.resize(n);
|
||||
@@ -142,6 +146,8 @@ void ScanResultGenerator::FillEndMessage(EndMessage &message) const {
|
||||
message.image_scale_factor[number] = e.image_scale_factor.value_or(NAN);
|
||||
message.image_scale_cc[number] = e.image_scale_cc.value_or(NAN);
|
||||
message.ice_ring_score[number] = e.ice_ring_score.value_or(NAN);
|
||||
message.v_protein_score[number] = e.protein_score.value_or(NAN);
|
||||
message.v_ice_score[number] = e.ice_score.value_or(NAN);
|
||||
message.integrated_reflections[number] = static_cast<int32_t>(value_or_zero(e.integrated_reflections));
|
||||
message.niggli_class[number] = static_cast<uint8_t>(value_or_zero(e.niggli_class));
|
||||
message.pixel_sum[number] = value_or_zero(e.pixel_sum);
|
||||
|
||||
@@ -225,6 +225,8 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s
|
||||
| packets_received | uint64 | Number of packets received per image (in units of 2 kB) | | |
|
||||
| bkg_estimate | float | Mean value for pixels in resolution range from 3.0 to 5.0 A \[photons\] | | |
|
||||
| ice_ring_score | float | Strongest hexagonal-ice ring intensity over the smooth radial background (1 = no ice) | | |
|
||||
| protein_score | float | Protein diffraction detection score, 0 to 1, saturating (0 = none, 1 = certain) | | |
|
||||
| ice_score | float | Crystalline ice detection score, 0 to 1, saturating (0 = none, 1 = certain) | | |
|
||||
| spot_count_ice_control | float | Spots in the ice-free flanks beside the hexagonal rings, rescaled to the ring bands' own q width (control for spot_count_ice_rings) | | |
|
||||
| beam_corr_x | float | Beam center correction X applied during processing \[pixel\] | | X |
|
||||
| beam_corr_y | float | Beam center correction Y applied during processing \[pixel\] | | X |
|
||||
@@ -317,8 +319,12 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s
|
||||
| image_indexed | Array(uint8) | Per-image indexing result; 0 = not indexed, nonzero = indexed | |
|
||||
| v_bkg_estimate | Array(float) | Per-image background estimate | |
|
||||
| ice_ring_score | Array(float) | Per-image strongest ice-ring intensity over the smooth radial background (1 = no ice) | |
|
||||
| v_protein_score | Array(float) | Per-image protein diffraction detection score, 0 to 1 | |
|
||||
| v_ice_score | Array(float) | Per-image crystalline ice detection score, 0 to 1 | |
|
||||
| spot_count_ice_control | Array(float) | Per-image spot count in the ice-free flanks beside the hexagonal rings, rescaled to the ring bands' q width | |
|
||||
| ice_ring_score_mean | float | Mean ice-ring score for the whole run (1 = no ice) | |
|
||||
| protein_score | float | Mean protein detection score for the whole run | |
|
||||
| ice_score | float | Mean ice detection score for the whole run | |
|
||||
| profile_radius | Array(float) | Per-image profile radius \[Angstrom^-1\] | |
|
||||
| mosaicity | Array(float) | Per-image mosaicity \[degree\] | |
|
||||
| bFactor | Array(float) | Per-image estimated B-factor \[Angstrom^2\] | |
|
||||
|
||||
@@ -416,3 +416,68 @@ is kept; the frame is then integrated once at that min-pix. The fraction factor
|
||||
A per-image **resolution estimate** is read off the finished spot list. It predicts how far the *merged* data will reach, not how far the furthest spot on this image lies. Each non-ice spot is weighted by $\sqrt{I}$ — the intensity is a summed photon count, so $\sqrt{I}$ is its Poisson significance — the $1/d^2$ is found beyond which a fraction $f=0.30$ of that weight lies, and the estimate is that resolution taken $2.25\times$ further in $1/d$. It is deliberately **not** limited to what the detector records: the quantile sits in the middle of the fall-off, well inside the recorded range, so it goes on measuring the crystal where the detector stops before the diffraction does, and on such a run it reads finer than the detector corner. The dataset value is the median over images.
|
||||
|
||||
Both constants carry a mechanism. A quantile from the middle of the distribution measures the *shape* of the fall-off, which is the crystal's own $\exp(-B/2d^{2})$, where the extreme end of it measures where detection stops — a threshold that moves with the exposure and with how many reflections the unit cell puts on a frame. And merging averages many observations of each reflection, so intensities go on being measurable a fixed factor in $1/d$ past the point at which one image's spot finder still detects them; that factor is the $2.25$. Both are calibrated on rotation data against the resolution at which per-shell CC1/2 falls through 0.30, and the estimate is good to about 0.2 Å there. It is a prediction and not a measurement of what a run achieved: nothing downstream is cut on it, and it is reported alone (rugnux `SPOT_RESOLUTION_ESTIMATE`, and per image in the stream, the plots and HDF5).
|
||||
|
||||
### 3.7 Detection scores: is there protein here, is there ice here
|
||||
|
||||
`iceRingScore` above is a *magnitude* — a ratio, unbounded, answering "how strong is the worst ring".
|
||||
Two further per-image scalars answer a different question, the one a grid scan actually asks:
|
||||
`proteinScore` and `iceScore`, both in $[0,1]$ and both **saturating**, so a superb crystal and a
|
||||
barely-diffracting one score the same. They are **detection** scores, not quality measures — the
|
||||
resolution estimate (§3.6), the $B$ factor and the indexing rate already say how good the
|
||||
diffraction is — and neither the spot count nor the resolution enters either of them as a term.
|
||||
|
||||
**Protein.** The evidence is a resolution band nothing else can reach. No ice or salt phase a cryo
|
||||
loop can carry diffracts beyond about 4 Å (the largest hexagonal-ice spacing is 3.895 Å, the largest
|
||||
elemental-metal one about 2.9 Å), while every protein cell has an axis over 20 Å, so a spot at
|
||||
$d > 5$ Å is near-proof of a protein-scale repeat. What spoils that argument is hardware rather than
|
||||
physics — parasitic scatter, beam-stop haloes and detector artefacts also put narrow rings in the
|
||||
low-resolution band — so the evidence counts distinct $d$ **shells** rather than spots. Spots in
|
||||
$5 < d < 40$ Å are grouped into shells 2 % wide in $\ln(1/d)$, each spot weighted
|
||||
$\min(2I/(I+I_\mathrm{med}), 1)$ against the frame's own median spot intensity, and a shell
|
||||
carrying total weight $w$ contributes $1 - e^{-w}$: one shell is worth at most 1, so no single ring
|
||||
can accumulate a score, and a scattering of the weakest detections cannot fill a shell either. The
|
||||
evidence $E$ is the sum over shells and the score is $E/(E + k)$ with $k = 3.70$, the only fitted
|
||||
number, calibrated by leaving out one negative loop at a time and taking the smallest $k$ that holds
|
||||
the pooled false-positive rate over the remaining ones below $5\times10^{-4}$.
|
||||
|
||||
**Ice** reaches the frame two ways, and they need different evidence, so two channels are computed
|
||||
and the stronger one wins. The *radial* channel reads the plain azimuthal profile — not the
|
||||
peak-excluded background the `iceRingScore` uses, because it needs the profile's own standard
|
||||
deviation, which that background does not carry. Each band is read as an excess over a running
|
||||
median (half-window 6 bins, which rejects a 3–5 bin powder ring but follows the ~40-bin vitreous
|
||||
halo, so the halo cannot score), in units of the bin mean's own error $\sigma/\sqrt{n}$ smoothed
|
||||
over ±20 bins and floored at 1 % of the background, clipped to $[0, 6]$ after subtracting a floor of
|
||||
2 so that only a real $>2\sigma$ excess counts at all. The best bin within ±0.012 Å⁻¹ is taken,
|
||||
which absorbs the radial smear a mis-set beam centre produces. The same statistic is measured at
|
||||
every profile bin belonging to no band, giving the frame its own null — a grainy profile therefore
|
||||
raises its own null as much as its own band values — and two standardised statistics are formed
|
||||
against it, an amplitude and a band-count concordance; the **smaller** is taken, so a single elevated
|
||||
bin fails and only a whole pattern scores. **Two ice phases are carried as separate hypotheses and
|
||||
the decision is taken at the end** (the larger score wins): flash-cooled loops show cubic or
|
||||
stacking-disordered ice at least as often as hexagonal, the two phases share only three lines, and
|
||||
dropping the cubic hypothesis costs about 5 pp of detection on iced loops. Cubic ice is
|
||||
$Fd\bar{3}m$ with $a = 6.358$ Å and the diamond reflection conditions. Finally
|
||||
$I = S^2/(S^2 + 3^2)$, so the conventional $3\sigma$ detection maps to 0.5.
|
||||
|
||||
The *spot* channel catches ice that arrives as large crystallites, which diffracts as discrete spots
|
||||
and leaves the radial profile flat. Its evidence is an excess of found spots on the hexagonal radii
|
||||
over what the frame's own radial spot density predicts. The null is not the two flanks beside each
|
||||
band — that is `spot_count_ice_control` above, a ratio of two ~1-count numbers — but each band slid
|
||||
to every ice-free offset within ±0.45 Å⁻¹, in $\pm\delta$ pairs so the fall-off of spot density with
|
||||
$q$ cancels to first order, each count divided by the live detector area at that radius (taken from
|
||||
the azimuthal integration's own per-bin pixel count, so a radius the detector edge cuts short is not
|
||||
mistaken for one with no spots). That turns a 1-count control into an average over ~100. The excess
|
||||
is read both as a quasi-Poisson upper tail — the controls' own scatter setting the overdispersion —
|
||||
and as a ratio, and the **smaller** 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.
|
||||
|
||||
Both scores read $d$ out of the geometry, so both move with a beam-centre error; the centre is
|
||||
**not** fitted here, because that belongs to geometry refinement. The centre they were computed with
|
||||
is written beside them as `/entry/MX/scoreBeamCenterX`/`Y`, since the refined centre may later
|
||||
overwrite `/entry/instrument/detector/beam_center_x` and a rescoring would then have no way to tell
|
||||
an algorithm disagreement from a geometry one.
|
||||
|
||||
The radial ice channel needs a profile standard deviation, which the FPGA azimuthal integration does
|
||||
not produce; on that path it abstains and `iceScore` is the spot channel alone, unless the profile is
|
||||
recomputed on the CPU (`force_cpu_in_fpga_workflow`). In `--mode azint` no spots are looked for, so
|
||||
the radial channel is all there is.
|
||||
|
||||
@@ -346,6 +346,8 @@ In legacy/VDS mode these live in the data files and are linked/virtual-stacked i
|
||||
| `integratedReflections` | | number of integrated reflections |
|
||||
| `bkgEstimate` | photons | mean background in the 3–5 Å resolution band |
|
||||
| `iceRingScore` | ratio | strongest hexagonal-ice ring intensity over the smooth radial background (1 = no ice) |
|
||||
| `proteinScore` | | protein diffraction detection score, 0 to 1, saturating (0 = none, 1 = certain) |
|
||||
| `iceScore` | | crystalline ice detection score, 0 to 1, saturating (0 = none, 1 = certain) |
|
||||
| `beam_corr_x`, `beam_corr_y` | pixel | beam-center correction applied during processing |
|
||||
| `imageScaleFactor` | | on-the-fly per-image scale factor *g* |
|
||||
| `imageScaleCC` | | on-the-fly scaling correlation coefficient |
|
||||
@@ -368,6 +370,9 @@ variants.
|
||||
| `imageIndexedMean` | | mean indexing rate over the run |
|
||||
| `bkgEstimateMean` | photons | mean background over the run |
|
||||
| `iceRingScoreMean` | ratio | mean `iceRingScore` over the run — the single "how icy was this dataset" number (1 = no ice) |
|
||||
| `proteinScoreMean` | | mean `proteinScore` over the run |
|
||||
| `iceScoreMean` | | mean `iceScore` over the run |
|
||||
| `scoreBeamCenterX`, `scoreBeamCenterY` | pixel | the beam centre `proteinScore` and `iceScore` were computed with. Both read *d* out of the geometry, and `/entry/instrument/detector/beam_center_x`/`_y` carries the **refined** centre where refinement ran, so without this pair a rescoring could not tell an algorithm disagreement from a geometry one |
|
||||
| `indexedLatticeCount` | | per-image lattice count summary (master). *Note: data files use `indexingLatticeCount`; readers accept either.* |
|
||||
| `reindexMatrix` | | change of basis from the setting the per-image data are in to the setting of `/entry/sample/unit_cell` (`[9]`, `int32`, flattened 3×3, row major) — see below |
|
||||
|
||||
|
||||
@@ -820,6 +820,10 @@ namespace {
|
||||
message.bkg_estimate = GetCBORFloat(value);
|
||||
else if (key == "ice_ring_score")
|
||||
message.ice_ring_score = GetCBORFloat(value);
|
||||
else if (key == "protein_score")
|
||||
message.protein_score = GetCBORFloat(value);
|
||||
else if (key == "ice_score")
|
||||
message.ice_score = GetCBORFloat(value);
|
||||
else if (key == "adu_histogram")
|
||||
GetCBORUInt64Array(value, message.adu_histogram);
|
||||
else if (key == "beam_corr_x")
|
||||
@@ -1493,6 +1497,14 @@ namespace {
|
||||
GetCBORUInt8Array(value, message.image_indexed);
|
||||
else if (key == "v_bkg_estimate")
|
||||
GetCBORFloatArray(value, message.v_bkg_estimate);
|
||||
else if (key == "v_protein_score")
|
||||
GetCBORFloatArray(value, message.v_protein_score);
|
||||
else if (key == "v_ice_score")
|
||||
GetCBORFloatArray(value, message.v_ice_score);
|
||||
else if (key == "protein_score")
|
||||
message.protein_score = GetCBORFloat(value);
|
||||
else if (key == "ice_score")
|
||||
message.ice_score = GetCBORFloat(value);
|
||||
else if (key == "ice_ring_score")
|
||||
GetCBORFloatArray(value, message.ice_ring_score);
|
||||
else if (key == "spot_count_ice_control")
|
||||
|
||||
@@ -820,6 +820,10 @@ void CBORStream2Serializer::SerializeSequenceEnd(const EndMessage& message) {
|
||||
CBOR_ENC(mapEncoder, "spot_count_indexed", message.spot_count_indexed);
|
||||
CBOR_ENC(mapEncoder, "image_indexed", message.image_indexed);
|
||||
CBOR_ENC(mapEncoder, "v_bkg_estimate", message.v_bkg_estimate);
|
||||
CBOR_ENC(mapEncoder, "v_protein_score", message.v_protein_score);
|
||||
CBOR_ENC(mapEncoder, "v_ice_score", message.v_ice_score);
|
||||
CBOR_ENC(mapEncoder, "protein_score", message.protein_score);
|
||||
CBOR_ENC(mapEncoder, "ice_score", message.ice_score);
|
||||
CBOR_ENC(mapEncoder, "ice_ring_score", message.ice_ring_score);
|
||||
CBOR_ENC(mapEncoder, "ice_ring_score_mean", message.ice_ring_score_mean);
|
||||
CBOR_ENC(mapEncoder, "spot_count_ice_control", message.spot_count_ice_control);
|
||||
@@ -918,6 +922,8 @@ void CBORStream2Serializer::SerializeImageInternal(CborEncoder &mapEncoder, cons
|
||||
CBOR_ENC(mapEncoder, "packets_received", message.packets_received);
|
||||
CBOR_ENC(mapEncoder, "bkg_estimate", message.bkg_estimate);
|
||||
CBOR_ENC(mapEncoder, "ice_ring_score", message.ice_ring_score);
|
||||
CBOR_ENC(mapEncoder, "protein_score", message.protein_score);
|
||||
CBOR_ENC(mapEncoder, "ice_score", message.ice_score);
|
||||
CBOR_ENC(mapEncoder, "adu_histogram", message.adu_histogram);
|
||||
CBOR_ENC(mapEncoder, "roi_integrals", message.roi);
|
||||
CBOR_ENC(mapEncoder, "beam_corr_x", message.beam_corr_x);
|
||||
|
||||
@@ -54,6 +54,9 @@ function AxisTypeY(plot: plot_type) : string | ReactNode {
|
||||
return "Count";
|
||||
case plot_type.ICE_RING_SCORE:
|
||||
return "Ratio";
|
||||
case plot_type.PROTEIN_SCORE:
|
||||
case plot_type.ICE_SCORE:
|
||||
return "Detection score";
|
||||
case plot_type.AZINT:
|
||||
case plot_type.AZINT_1D:
|
||||
case plot_type.BKG_ESTIMATE:
|
||||
|
||||
@@ -51,6 +51,8 @@ function DataProcessingPlots({type: initialType, height}: MyProps) {
|
||||
<MenuItem value={plot_type.SPOT_COUNT_INDEXED}>Spot count indexed</MenuItem>
|
||||
<MenuItem value={plot_type.SPOT_COUNT_ICE}>Spot count ice ring</MenuItem>
|
||||
<MenuItem value={plot_type.ICE_RING_SCORE}>Ice ring score</MenuItem>
|
||||
<MenuItem value={plot_type.PROTEIN_SCORE}>Protein detection score</MenuItem>
|
||||
<MenuItem value={plot_type.ICE_SCORE}>Ice detection score</MenuItem>
|
||||
<MenuItem value={plot_type.AZINT}>Azimuthal integration profile</MenuItem>
|
||||
<MenuItem value={plot_type.AZINT_1D}>Azimuthal integration profile (1D)</MenuItem>
|
||||
<MenuItem value={plot_type.BKG_ESTIMATE}>Background estimate</MenuItem>
|
||||
|
||||
@@ -59,6 +59,8 @@ ADD_LIBRARY(JFJochImageAnalysis STATIC
|
||||
MXAnalysisAfterFPGA.cpp
|
||||
IndexAndRefine.cpp
|
||||
IndexAndRefine.h
|
||||
IceScore.cpp
|
||||
IceScore.h
|
||||
dark_mask_analysis/DarkMaskAnalysis.cpp
|
||||
dark_mask_analysis/DarkMaskAnalysis.h
|
||||
beam_stop/ShadowFinder.cpp
|
||||
|
||||
@@ -0,0 +1,333 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
|
||||
#include "IceScore.h"
|
||||
#include "../common/Definitions.h"
|
||||
|
||||
namespace {
|
||||
constexpr float TWO_PI = 6.283185307f;
|
||||
|
||||
// Below this q the profile is direct beam and beam-stop shadow, not diffraction. The lowest ice
|
||||
// line of either phase sits at 1.61 A^-1, so nothing is lost.
|
||||
constexpr float Q_MIN_EVAL = 1.20f;
|
||||
// Running-median half-window for the background under the bands. A median over 13 bins rejects a
|
||||
// 3-5 bin powder ring but follows the ~40-bin-wide vitreous-water halo, so the halo never reaches
|
||||
// the residual and cannot score.
|
||||
constexpr int BACKGROUND_HALF_WINDOW = 6;
|
||||
// Half-window over which the bin mean's own error is smoothed, so that a ring cannot inflate its
|
||||
// own denominator.
|
||||
constexpr int SIGMA_HALF_WINDOW = 20;
|
||||
// A band must be a real excess to contribute at all; a uniform low-level positive bias spread over
|
||||
// many bands then contributes nothing. Capped so that one saturated band cannot carry a hypothesis.
|
||||
constexpr float Z_FLOOR = 2.0f;
|
||||
constexpr float Z_CAP = 6.0f;
|
||||
// The conventional 3-sigma detection maps to a score of 0.5.
|
||||
constexpr float S_HALF = 3.0f;
|
||||
// Half-width of a band, in q. A mis-set beam centre smears and splits a ring, so the best bin
|
||||
// within this distance is taken rather than the one the geometry predicts.
|
||||
constexpr float BAND_HALF_Q = 0.012f;
|
||||
// Bins nearer than this to a band of EITHER phase are not part of the null.
|
||||
constexpr float NULL_EXCLUSION_Q = 0.025f;
|
||||
// Fewer null bins than this and the null is not measured.
|
||||
constexpr int NULL_MIN_BINS = 20;
|
||||
|
||||
// Band half-width for the spot channel, in q. Wider than the spot finder's own ice marking width
|
||||
// (ice_ring_width_Q_recipA, 0.03): measured, 0.03 costs 5 pp of specificity on protein frames for
|
||||
// no gain on ice, because a wider band collects more of the crystal's own reflections.
|
||||
constexpr float SPOT_BAND_HALF_Q = 0.02f;
|
||||
// How far a band may be slid for its control, and in what steps.
|
||||
constexpr float SPOT_OFFSET_MAX_Q = 0.45f;
|
||||
constexpr float SPOT_OFFSET_STEP_Q = 0.0025f;
|
||||
// A control radius the detector barely covers is not a control.
|
||||
constexpr float SPOT_MIN_AREA_FRACTION = 0.05f;
|
||||
constexpr int SPOT_MIN_OFFSETS = 5;
|
||||
// A control count below this is not measurable; do not divide by it.
|
||||
constexpr float SPOT_CONTROL_FLOOR = 0.5f;
|
||||
// -log10 p at which the significance term saturates, and the band/control ratio at which the size
|
||||
// term does.
|
||||
constexpr float SPOT_SIGNIFICANCE_SAT = 4.0f;
|
||||
constexpr float SPOT_RATIO_SAT = 3.0f;
|
||||
|
||||
// Average a q_bins x azimuthal-bins array over azimuth. A profile that is already 1-D passes
|
||||
// through. Bins with nothing finite in them come out NaN.
|
||||
std::vector<float> FoldToQ(const std::vector<float> &in, int nq) {
|
||||
std::vector<float> out(nq, 0.0f);
|
||||
std::vector<int> n(nq, 0);
|
||||
for (size_t i = 0; i < in.size() && nq > 0; i++) {
|
||||
const int q = static_cast<int>(i % nq);
|
||||
if (std::isfinite(in[i])) {
|
||||
out[q] += in[i];
|
||||
n[q]++;
|
||||
}
|
||||
}
|
||||
for (int q = 0; q < nq; q++)
|
||||
out[q] = n[q] ? out[q] / static_cast<float>(n[q]) : NAN;
|
||||
return out;
|
||||
}
|
||||
|
||||
// The same fold for the live pixel count, which adds rather than averages.
|
||||
std::vector<double> FoldCountToQ(const std::vector<uint64_t> &in, int nq) {
|
||||
std::vector<double> out(nq, 0.0);
|
||||
for (size_t i = 0; i < in.size() && nq > 0; i++)
|
||||
out[i % nq] += static_cast<double>(in[i]);
|
||||
return out;
|
||||
}
|
||||
|
||||
// Running median, NaN where the window holds nothing finite.
|
||||
std::vector<float> RunningMedian(const std::vector<float> &v, int half) {
|
||||
std::vector<float> out(v.size(), NAN);
|
||||
std::vector<float> window;
|
||||
for (int i = 0; i < static_cast<int>(v.size()); i++) {
|
||||
window.clear();
|
||||
const int lo = std::max(0, i - half);
|
||||
const int hi = std::min(static_cast<int>(v.size()), i + half + 1);
|
||||
for (int j = lo; j < hi; j++)
|
||||
if (std::isfinite(v[j]))
|
||||
window.push_back(v[j]);
|
||||
if (window.empty())
|
||||
continue;
|
||||
std::ranges::nth_element(window, window.begin() + window.size() / 2);
|
||||
out[i] = window[window.size() / 2];
|
||||
}
|
||||
return out;
|
||||
}
|
||||
|
||||
// One ice phase read as a hypothesis over the per-bin evidence c: the standardised amplitude of its
|
||||
// bands and the standardised count of bands that show anything, against the frame's own null.
|
||||
float HypothesisScore(const std::vector<float> &c, const std::vector<float> &band_q,
|
||||
const std::vector<float> &band_weight, int min_bands,
|
||||
float low_q, float dq, int bh,
|
||||
float mu0, float v0, float p0) {
|
||||
const int nq = static_cast<int>(c.size());
|
||||
float w1 = 0.0f, w2 = 0.0f, evidence = 0.0f;
|
||||
int nbands = 0, nhit = 0;
|
||||
for (size_t i = 0; i < band_q.size(); i++) {
|
||||
if (band_q[i] < Q_MIN_EVAL)
|
||||
continue;
|
||||
const int b = static_cast<int>(std::lround((band_q[i] - low_q) / dq - 0.5f));
|
||||
if (b < bh || b >= nq - bh || !std::isfinite(c[b]))
|
||||
continue;
|
||||
nbands++;
|
||||
w1 += band_weight[i];
|
||||
w2 += band_weight[i] * band_weight[i];
|
||||
evidence += c[b] * band_weight[i];
|
||||
if (c[b] > 0.0f)
|
||||
nhit++;
|
||||
}
|
||||
if (nbands < min_bands)
|
||||
return 0.0f;
|
||||
const float s_amp = (evidence - mu0 * w1) / std::sqrt(std::max(v0 * w2, 1e-6f));
|
||||
const float n = static_cast<float>(nbands);
|
||||
const float s_cnt = (static_cast<float>(nhit) - n * p0)
|
||||
/ std::sqrt(std::max(n * p0 * (1.0f - p0), 0.0f) + 0.25f);
|
||||
return std::max(std::min(s_amp, s_cnt), 0.0f);
|
||||
}
|
||||
|
||||
// P(X >= k) for X ~ Poisson(mu), summed upward from k. The lower tail would lose a small
|
||||
// probability to cancellation, and it is the small probabilities this is wanted for.
|
||||
double PoissonUpperTail(int64_t k, double mu) {
|
||||
if (k <= 0)
|
||||
return 1.0;
|
||||
if (!(mu > 0.0))
|
||||
return 0.0;
|
||||
const double log_mu = std::log(mu);
|
||||
double sum = 0.0;
|
||||
for (int64_t j = k; j < k + 100000; j++) {
|
||||
const double term = std::exp(-mu + static_cast<double>(j) * log_mu - std::lgamma(static_cast<double>(j) + 1.0));
|
||||
sum += term;
|
||||
if (static_cast<double>(j) > mu && term < 1e-18 * sum)
|
||||
break;
|
||||
}
|
||||
return std::min(sum, 1.0);
|
||||
}
|
||||
|
||||
// Live pixels within +-SPOT_BAND_HALF_Q of q.
|
||||
double AreaAt(const std::vector<double> &count_q, float low_q, float dq, float q) {
|
||||
double area = 0.0;
|
||||
const int lo = static_cast<int>(std::ceil((q - SPOT_BAND_HALF_Q - low_q) / dq - 0.5f));
|
||||
const int hi = static_cast<int>(std::floor((q + SPOT_BAND_HALF_Q - low_q) / dq - 0.5f));
|
||||
for (int i = std::max(lo, 0); i <= std::min(hi, static_cast<int>(count_q.size()) - 1); i++)
|
||||
area += count_q[i];
|
||||
return area;
|
||||
}
|
||||
|
||||
// Spots within +-SPOT_BAND_HALF_Q of q, from the sorted spot q list.
|
||||
double CountAt(const std::vector<float> &sorted_q, float q) {
|
||||
const auto lo = std::lower_bound(sorted_q.begin(), sorted_q.end(), q - SPOT_BAND_HALF_Q);
|
||||
const auto hi = std::upper_bound(sorted_q.begin(), sorted_q.end(), q + SPOT_BAND_HALF_Q);
|
||||
return static_cast<double>(std::distance(lo, hi));
|
||||
}
|
||||
}
|
||||
|
||||
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) {
|
||||
const int nq = std::max<int>(q_bins, 0);
|
||||
const float low_q = settings.GetLowQ_recipA();
|
||||
const float dq = settings.GetQSpacing_recipA();
|
||||
if (nq < 2 * SIGMA_HALF_WINDOW || !(dq > 0.0f) || profile_std.empty())
|
||||
return 0.0f;
|
||||
|
||||
std::vector<float> prof = FoldToQ(profile, nq);
|
||||
std::vector<float> sigma = FoldToQ(profile_std, nq);
|
||||
const std::vector<double> count = FoldCountToQ(profile_count, nq);
|
||||
for (int i = 0; i < nq; i++) {
|
||||
if (!(prof[i] > 0.0f))
|
||||
prof[i] = NAN;
|
||||
// The error of the bin mean, not of one pixel.
|
||||
if (sigma[i] > 0.0f && count[i] > 0.0)
|
||||
sigma[i] /= static_cast<float>(std::sqrt(count[i]));
|
||||
else
|
||||
sigma[i] = NAN;
|
||||
}
|
||||
|
||||
const std::vector<float> background = RunningMedian(prof, BACKGROUND_HALF_WINDOW);
|
||||
const std::vector<float> sigma_smooth = RunningMedian(sigma, SIGMA_HALF_WINDOW);
|
||||
const int bh = std::max(1, static_cast<int>(std::lround(BAND_HALF_Q / dq)));
|
||||
auto q_of = [&](int i) { return low_q + (static_cast<float>(i) + 0.5f) * dq; };
|
||||
|
||||
// Per-bin excess in units of the bin mean's own error, floored so that only a real excess counts.
|
||||
// The sigma floor at 1 % of the background keeps a bin whose error is reported as tiny - or not at
|
||||
// all - from turning noise into a detection.
|
||||
std::vector<float> z(nq, NAN);
|
||||
for (int i = 0; i < nq; i++) {
|
||||
if (q_of(i) < Q_MIN_EVAL || !std::isfinite(prof[i]) || !std::isfinite(background[i]))
|
||||
continue;
|
||||
const float floor_sigma = 0.01f * background[i];
|
||||
const float sg = std::max(std::isfinite(sigma_smooth[i]) ? sigma_smooth[i] : 0.0f, floor_sigma);
|
||||
if (sg > 0.0f)
|
||||
z[i] = (prof[i] - background[i]) / sg;
|
||||
}
|
||||
std::vector<float> c(nq, NAN);
|
||||
for (int i = 0; i < nq; i++) {
|
||||
float best = NAN;
|
||||
for (int o = -bh; o <= bh; o++) {
|
||||
const float v = z[std::clamp(i + o, 0, nq - 1)];
|
||||
if (std::isfinite(v) && (!std::isfinite(best) || v > best))
|
||||
best = v;
|
||||
}
|
||||
if (std::isfinite(best))
|
||||
c[i] = std::clamp(best - Z_FLOOR, 0.0f, Z_CAP);
|
||||
}
|
||||
|
||||
// The null: the same statistic at every bin belonging to no band of either phase.
|
||||
std::vector<float> null_c;
|
||||
for (int i = 0; i < nq; i++) {
|
||||
if (!std::isfinite(c[i]) || q_of(i) < Q_MIN_EVAL)
|
||||
continue;
|
||||
float nearest = 1e9f;
|
||||
for (const float d: ICE_RING_RES_A)
|
||||
nearest = std::min(nearest, std::fabs(q_of(i) - TWO_PI / d));
|
||||
for (const float d: ICE_RING_CUBIC_RES_A)
|
||||
nearest = std::min(nearest, std::fabs(q_of(i) - TWO_PI / d));
|
||||
if (nearest > static_cast<float>(bh) * dq + NULL_EXCLUSION_Q)
|
||||
null_c.push_back(c[i]);
|
||||
}
|
||||
if (static_cast<int>(null_c.size()) < NULL_MIN_BINS)
|
||||
return 0.0f;
|
||||
|
||||
float mu0 = 0.0f, p0 = 0.0f;
|
||||
for (const float v: null_c) {
|
||||
mu0 += v;
|
||||
p0 += (v > 0.0f) ? 1.0f : 0.0f;
|
||||
}
|
||||
const float kn = static_cast<float>(null_c.size());
|
||||
mu0 /= kn;
|
||||
p0 /= kn;
|
||||
float v0 = 0.0f;
|
||||
for (const float v: null_c)
|
||||
v0 += (v - mu0) * (v - mu0);
|
||||
v0 /= (kn - 1.0f);
|
||||
|
||||
// Hexagonal ice: the primary triplet above 3 A carries twice the weight of the rest.
|
||||
std::vector<float> hex_q, hex_w;
|
||||
for (const float d: ICE_RING_RES_A) {
|
||||
hex_q.push_back(TWO_PI / d);
|
||||
hex_w.push_back(d > 3.0f ? 2.0f : 1.0f);
|
||||
}
|
||||
std::vector<float> cubic_q, cubic_w;
|
||||
for (const float d: ICE_RING_CUBIC_RES_A) {
|
||||
cubic_q.push_back(TWO_PI / d);
|
||||
cubic_w.push_back(d > 3.0f ? 2.0f : 1.0f);
|
||||
}
|
||||
|
||||
const float s = std::max(HypothesisScore(c, hex_q, hex_w, 4, low_q, dq, bh, mu0, v0, p0),
|
||||
HypothesisScore(c, cubic_q, cubic_w, 3, low_q, dq, bh, mu0, v0, p0));
|
||||
return s * s / (s * s + S_HALF * S_HALF);
|
||||
}
|
||||
|
||||
float IceScoreSpots(const std::vector<SpotToSave> &spots, const std::vector<uint64_t> &profile_count,
|
||||
int32_t q_bins, const AzimuthalIntegrationSettings &settings) {
|
||||
const int nq = std::max<int>(q_bins, 0);
|
||||
const float low_q = settings.GetLowQ_recipA();
|
||||
const float dq = settings.GetQSpacing_recipA();
|
||||
if (nq < 2 || !(dq > 0.0f) || spots.empty() || profile_count.empty())
|
||||
return 0.0f;
|
||||
|
||||
const std::vector<double> count_q = FoldCountToQ(profile_count, nq);
|
||||
|
||||
std::vector<float> spot_q;
|
||||
spot_q.reserve(spots.size());
|
||||
for (const auto &s: spots)
|
||||
if (s.d_A > 0.0f)
|
||||
spot_q.push_back(TWO_PI / s.d_A);
|
||||
std::ranges::sort(spot_q);
|
||||
|
||||
std::vector<float> band_q;
|
||||
for (const float d: ICE_RING_RES_A)
|
||||
band_q.push_back(TWO_PI / d);
|
||||
|
||||
double k0 = 0.0, mu = 0.0, var = 0.0;
|
||||
for (const float qb: band_q) {
|
||||
const double area_band = AreaAt(count_q, low_q, dq, qb);
|
||||
if (!(area_band > 0.0))
|
||||
continue;
|
||||
std::vector<double> control;
|
||||
for (float delta = 2.0f * SPOT_BAND_HALF_Q; delta <= SPOT_OFFSET_MAX_Q; delta += SPOT_OFFSET_STEP_Q) {
|
||||
// An offset that lands on another ice band is not a control.
|
||||
float nearest = 1e9f;
|
||||
for (const float qo: band_q)
|
||||
nearest = std::min(nearest, std::min(std::fabs(qb + delta - qo), std::fabs(qb - delta - qo)));
|
||||
if (nearest < 2.0f * SPOT_BAND_HALF_Q)
|
||||
continue;
|
||||
const double area_hi = AreaAt(count_q, low_q, dq, qb + delta);
|
||||
const double area_lo = AreaAt(count_q, low_q, dq, qb - delta);
|
||||
if (area_hi < SPOT_MIN_AREA_FRACTION * area_band || area_lo < SPOT_MIN_AREA_FRACTION * area_band)
|
||||
continue;
|
||||
control.push_back(0.5 * area_band * (CountAt(spot_q, qb + delta) / area_hi
|
||||
+ CountAt(spot_q, qb - delta) / area_lo));
|
||||
}
|
||||
if (static_cast<int>(control.size()) < SPOT_MIN_OFFSETS)
|
||||
continue;
|
||||
double mean = 0.0;
|
||||
for (const double v: control)
|
||||
mean += v;
|
||||
mean /= static_cast<double>(control.size());
|
||||
double m2 = 0.0;
|
||||
for (const double v: control)
|
||||
m2 += (v - mean) * (v - mean);
|
||||
k0 += CountAt(spot_q, qb);
|
||||
mu += mean;
|
||||
var += m2 / static_cast<double>(control.size());
|
||||
}
|
||||
|
||||
mu = std::max(mu, static_cast<double>(SPOT_CONTROL_FLOOR));
|
||||
if (k0 <= mu)
|
||||
return 0.0f;
|
||||
// Quasi-Poisson: the controls' own scatter says how much wider than Poisson the count really is.
|
||||
const double phi = std::max(1.0, var / mu);
|
||||
const double p = PoissonUpperTail(static_cast<int64_t>(std::floor((k0 - 1.0) / phi)) + 1, mu / phi);
|
||||
const double s_sig = std::min(1.0, -std::log10(std::max(p, 1e-300)) / SPOT_SIGNIFICANCE_SAT);
|
||||
const double s_size = std::min(1.0, std::log(k0 / mu) / std::log(SPOT_RATIO_SAT));
|
||||
return static_cast<float>(std::max(0.0, std::min(s_sig, s_size)));
|
||||
}
|
||||
|
||||
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) {
|
||||
return std::max(IceScoreRadial(profile, profile_std, profile_count, q_bins, settings),
|
||||
IceScoreSpots(spots, profile_count, q_bins, settings));
|
||||
}
|
||||
@@ -0,0 +1,61 @@
|
||||
// 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);
|
||||
@@ -10,6 +10,7 @@
|
||||
#include "../compression/JFJochDecompress.h"
|
||||
|
||||
#include "spot_finding/SpotUtils.h"
|
||||
#include "IceScore.h"
|
||||
#include "bragg_prediction/BraggPredictionFactory.h"
|
||||
#include "image_preprocessing/ImagePreprocessorCPU.h"
|
||||
|
||||
@@ -291,6 +292,13 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output,
|
||||
output.ice_ring_score = AzimuthalIntegrationProfile::IceRingScore(
|
||||
have_ring_bkg ? ring_bkg : profile.GetResult1D(), integration.GetQBinCount(),
|
||||
integration.Settings(), spot_finding_settings.ice_ring_width_Q_recipA);
|
||||
|
||||
// The ice DETECTION score, unlike the ratio above, reads the plain profile: it measures each band
|
||||
// against the profile's own standard deviation, which the peak-excluded ring background does not
|
||||
// carry. Its second channel reads the spots, which SpotAnalyze has already put in the message.
|
||||
output.ice_score = IceScore(output.az_int_profile, output.az_int_profile_std,
|
||||
output.az_int_profile_count, integration.GetQBinCount(),
|
||||
integration.Settings(), output.spots);
|
||||
}
|
||||
|
||||
ImageSpotFinder &MXAnalysisWithoutFPGA::FixedThresholdFinder() {
|
||||
|
||||
@@ -1,6 +1,8 @@
|
||||
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include <map>
|
||||
|
||||
#include "../../common/JFJochMath.h"
|
||||
#include "SpotUtils.h"
|
||||
#include "../../common/ResolutionShells.h"
|
||||
@@ -167,6 +169,53 @@ std::optional<float> GetResolution(const std::vector<SpotToSave> &spots) {
|
||||
return 1.0f / (SPOT_RESOLUTION_MERGE_REACH * std::sqrt(one_over_d2));
|
||||
}
|
||||
|
||||
namespace {
|
||||
// The band nothing but a protein-scale repeat reaches. Its floor sits above every ice and salt
|
||||
// spacing; its ceiling is where a "spot" stops being a lattice reflection and starts being
|
||||
// beam-stop halo. Measured: dropping the ceiling costs a factor of two in false positives.
|
||||
constexpr float PROTEIN_BAND_LOW_A = 5.0f;
|
||||
constexpr float PROTEIN_BAND_HIGH_A = 40.0f;
|
||||
// Shell width for the saturation, in ln(1/d) - 2% in d. Spots are grouped into shells this wide and
|
||||
// each shell contributes at most 1, so a single narrow parasitic ring cannot accumulate evidence.
|
||||
constexpr float PROTEIN_SHELL_BIN = 0.02f;
|
||||
// Evidence at which the score reaches 0.5. The only fitted number in the score; calibrated by
|
||||
// leaving out one negative loop at a time and taking the smallest value that holds the pooled false
|
||||
// positive rate over the rest below 5e-4.
|
||||
constexpr float PROTEIN_SATURATION = 3.70f;
|
||||
}
|
||||
|
||||
float ProteinScore(const std::vector<SpotToSave> &spots) {
|
||||
// Intensity reference: the frame's own median spot. The weight below is 1 for a spot at or above it
|
||||
// and falls off smoothly beneath, so what counts is how a candidate stands against the rest of this
|
||||
// image - not an absolute photon count, which would make the score follow the exposure.
|
||||
std::vector<float> intensity;
|
||||
intensity.reserve(spots.size());
|
||||
for (const auto &s: spots)
|
||||
intensity.push_back(std::max(s.intensity, 0.0f));
|
||||
float i_ref = 0.0f;
|
||||
if (!intensity.empty()) {
|
||||
const size_t mid = intensity.size() / 2;
|
||||
std::ranges::nth_element(intensity, intensity.begin() + mid);
|
||||
i_ref = intensity[mid];
|
||||
}
|
||||
|
||||
// Sum the weights shell by shell, then saturate each shell separately.
|
||||
std::map<int, float> shell_weight;
|
||||
for (const auto &s: spots) {
|
||||
if (!(s.d_A > PROTEIN_BAND_LOW_A) || !(s.d_A < PROTEIN_BAND_HIGH_A))
|
||||
continue;
|
||||
const float i = std::max(s.intensity, 0.0f);
|
||||
const float w = (i_ref > 0.0f) ? std::min(2.0f * i / (i + i_ref), 1.0f) : 1.0f;
|
||||
shell_weight[static_cast<int>(std::log(1.0f / s.d_A) / PROTEIN_SHELL_BIN)] += w;
|
||||
}
|
||||
|
||||
float evidence = 0.0f;
|
||||
for (const auto &[shell, w]: shell_weight)
|
||||
evidence += 1.0f - std::exp(-w);
|
||||
|
||||
return evidence / (evidence + PROTEIN_SATURATION);
|
||||
}
|
||||
|
||||
void GenerateSpotPlot(DataMessage &msg, const std::vector<SpotToSave> &spots, float d_min_A) {
|
||||
const int nshells = 20;
|
||||
// The geometry gives no usable high-resolution corner (no distance or no wavelength), so there is
|
||||
@@ -233,6 +282,7 @@ void SpotAnalyze(const DiffractionExperiment &experiment,
|
||||
spot_d_min.value_or(0.0f) > 0 ? *spot_d_min : experiment.GetDetectorMaxResolution_A());
|
||||
|
||||
output.resolution_estimate = GetResolution(spots_out);
|
||||
output.protein_score = ProteinScore(spots_out);
|
||||
|
||||
// One decision drives both: if indexing is to use the ice-band spots, the spot budget must not
|
||||
// throw them away before it gets the chance.
|
||||
|
||||
@@ -43,6 +43,25 @@ void FilterSpuriousHighResolutionSpots(std::vector<SpotToSave> &spots, float thr
|
||||
// Returns nothing when the image has too few spots to have a fall-off at all.
|
||||
std::optional<float> GetResolution(const std::vector<SpotToSave> &spots);
|
||||
|
||||
// Is there PROTEIN diffraction on this image? A detection score in [0,1] that saturates: a superb
|
||||
// crystal and a barely-diffracting one both come out near 1. It is not a quality measure - resolution,
|
||||
// B factor and the indexing rate already say how good the diffraction is - and neither the spot count
|
||||
// nor the resolution enters it as a term.
|
||||
//
|
||||
// The physics is a resolution band nothing else can reach. No ice or salt phase a cryo loop can carry
|
||||
// diffracts beyond ~4 A (the largest hexagonal-ice spacing is 3.895 A, the largest elemental-metal one
|
||||
// ~2.9 A), while every protein cell has an axis over 20 A. A spot at d > 5 A is therefore near-proof of
|
||||
// a protein-scale repeat. What spoils that argument is hardware, not physics: parasitic scatter,
|
||||
// beam-stop haloes and detector artefacts also put narrow rings in the low-resolution band. So the
|
||||
// evidence counts distinct d SHELLS rather than spots - one shell is worth at most 1, so no single ring
|
||||
// can accumulate a score - and weights each spot by its intensity against the frame's own median, so a
|
||||
// scattering of the weakest detections cannot fill a shell either.
|
||||
//
|
||||
// Uses d_A as the geometry gave it. The band floor at 5 A is a resolution, so the score moves with a
|
||||
// beam-centre error; the centre is not fitted here, deliberately - that belongs to geometry refinement,
|
||||
// and the beam centre the scores were computed with is recorded beside them in the written file.
|
||||
float ProteinScore(const std::vector<SpotToSave> &spots);
|
||||
|
||||
void SpotAnalyze(const DiffractionExperiment &experiment,
|
||||
const SpotFindingSettings &settings,
|
||||
const std::vector<DiffractionSpot> &spots,
|
||||
|
||||
@@ -558,6 +558,10 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen
|
||||
dataset->indexing_result = master_file->ReadOptVector<float>("/entry/MX/imageIndexed");
|
||||
dataset->bkg_estimate = master_file->ReadOptVector<float>("/entry/MX/bkgEstimate");
|
||||
dataset->ice_ring_score = master_file->ReadOptVector<float>("/entry/MX/iceRingScore");
|
||||
dataset->protein_score = master_file->ReadOptVector<float>("/entry/MX/proteinScore");
|
||||
dataset->ice_score = master_file->ReadOptVector<float>("/entry/MX/iceScore");
|
||||
dataset->score_beam_center_x = master_file->GetOptFloat("/entry/MX/scoreBeamCenterX");
|
||||
dataset->score_beam_center_y = master_file->GetOptFloat("/entry/MX/scoreBeamCenterY");
|
||||
dataset->resolution_estimate = master_file->ReadOptVector<float>("/entry/MX/resolutionEstimate");
|
||||
dataset->profile_radius = master_file->ReadOptVector<float>("/entry/MX/profileRadius");
|
||||
// Master files write indexedLatticeCount; data files / the per-file MX
|
||||
@@ -704,6 +708,14 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen
|
||||
data_file, "/entry/MX/iceRingScore",
|
||||
number_of_images, fimages);
|
||||
|
||||
ReadVector(dataset->protein_score,
|
||||
data_file, "/entry/MX/proteinScore",
|
||||
number_of_images, fimages);
|
||||
|
||||
ReadVector(dataset->ice_score,
|
||||
data_file, "/entry/MX/iceScore",
|
||||
number_of_images, fimages);
|
||||
|
||||
ReadVector(dataset->profile_radius,
|
||||
data_file, "/entry/MX/profileRadius",
|
||||
number_of_images, fimages);
|
||||
@@ -1329,6 +1341,10 @@ void HDF5MetadataSource::FillPerImage(DataMessage &message, int64_t requested_im
|
||||
message.bkg_estimate = dataset->bkg_estimate[image_number];
|
||||
if (dataset->ice_ring_score.size() > image_number)
|
||||
message.ice_ring_score = dataset->ice_ring_score[image_number];
|
||||
if (dataset->protein_score.size() > image_number)
|
||||
message.protein_score = dataset->protein_score[image_number];
|
||||
if (dataset->ice_score.size() > image_number)
|
||||
message.ice_score = dataset->ice_score[image_number];
|
||||
if (dataset->efficiency.size() > image_number)
|
||||
message.image_collection_efficiency = dataset->efficiency[image_number];
|
||||
if (dataset->profile_radius.size() > image_number)
|
||||
|
||||
@@ -57,6 +57,14 @@ struct JFJochReaderDataset {
|
||||
std::vector<float> indexing_lattice_count;
|
||||
std::vector<float> bkg_estimate;
|
||||
std::vector<float> ice_ring_score;
|
||||
// The two per-image detection scores, in [0,1] and saturating: is there protein diffraction on
|
||||
// this image, and is there crystalline ice. From /entry/MX/proteinScore and /entry/MX/iceScore;
|
||||
// empty when the file predates them. score_beam_center_x/y is the beam centre they were computed
|
||||
// with, which is NOT necessarily the refined one the geometry above carries.
|
||||
std::vector<float> protein_score;
|
||||
std::vector<float> ice_score;
|
||||
std::optional<float> score_beam_center_x;
|
||||
std::optional<float> score_beam_center_y;
|
||||
std::vector<float> resolution_estimate;
|
||||
std::vector<float> efficiency;
|
||||
std::vector<float> profile_radius;
|
||||
|
||||
@@ -163,6 +163,8 @@ void JFJochReceiver::SendEndMessage() {
|
||||
|
||||
message.bkg_estimate = plots.GetBkgEstimate();
|
||||
message.ice_ring_score_mean = plots.GetIceRingScore();
|
||||
message.protein_score = plots.GetProteinScore();
|
||||
message.ice_score = plots.GetIceScore();
|
||||
message.indexing_rate = plots.GetIndexingRate();
|
||||
|
||||
message.az_int_result["dataset"] = plots.GetAzIntProfile();
|
||||
|
||||
@@ -6,6 +6,7 @@
|
||||
#include <thread>
|
||||
|
||||
#include "ImageMetadata.h"
|
||||
#include "../image_analysis/IceScore.h"
|
||||
|
||||
JFJochReceiverFPGA::JFJochReceiverFPGA(const DiffractionExperiment &in_experiment,
|
||||
const PixelMask &in_pixel_mask,
|
||||
@@ -436,6 +437,13 @@ void JFJochReceiverFPGA::FrameTransformationThread(uint32_t threadid) {
|
||||
message.bkg_estimate = az_int_profile_image.GetBkgEstimate(experiment.GetAzimuthalIntegrationSettings());
|
||||
message.ice_ring_score = az_int_profile_image.GetIceRingScore(
|
||||
experiment.GetAzimuthalIntegrationSettings(), spot_finding_settings.ice_ring_width_Q_recipA);
|
||||
// The radial channel of the ice score needs the profile's standard deviation, which the
|
||||
// FPGA azimuthal integration does not produce - it abstains there and only the spot
|
||||
// channel contributes. ForceCPUinFPGAWorkflow gives it the standard deviation back.
|
||||
message.ice_score = IceScore(message.az_int_profile, message.az_int_profile_std,
|
||||
message.az_int_profile_count,
|
||||
experiment.GetAzimuthalIntegrationSettings().GetQBinCount(),
|
||||
experiment.GetAzimuthalIntegrationSettings(), message.spots);
|
||||
|
||||
scan_result.Add(message);
|
||||
|
||||
|
||||
@@ -2,6 +2,7 @@
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include <cstdlib>
|
||||
#include "../image_analysis/IceScore.h"
|
||||
#include "Rugnux.h"
|
||||
#include "ModelValidation.h"
|
||||
#include "WriteModel.h"
|
||||
@@ -3291,6 +3292,11 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
|
||||
msg.bkg_estimate = profile.GetBkgEstimate(mapping.Settings());
|
||||
msg.ice_ring_score = profile.GetIceRingScore(mapping.Settings(),
|
||||
config_.spot_finding.ice_ring_width_Q_recipA);
|
||||
// Azimuthal integration only: no spots were looked for, so the radial channel is all the
|
||||
// ice score has.
|
||||
msg.ice_score = IceScore(msg.az_int_profile, msg.az_int_profile_std,
|
||||
msg.az_int_profile_count, mapping.GetQBinCount(),
|
||||
mapping.Settings(), msg.spots);
|
||||
msg.run_number = experiment_.GetRunNumber();
|
||||
msg.run_name = experiment_.GetRunName();
|
||||
|
||||
@@ -3496,6 +3502,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
|
||||
result.spot_resolution_estimate_A = plots.GetResolutionEstimate();
|
||||
end_msg.ice_ring_score = plots.GetIceRingScoreArray();
|
||||
end_msg.ice_ring_score_mean = plots.GetIceRingScore();
|
||||
end_msg.v_protein_score = plots.GetProteinScoreArray();
|
||||
end_msg.v_ice_score = plots.GetIceScoreArray();
|
||||
end_msg.protein_score = plots.GetProteinScore();
|
||||
end_msg.ice_score = plots.GetIceScore();
|
||||
end_msg.az_int_result["dataset"] = plots.GetAzIntProfile();
|
||||
end_msg.indexing_rate = result.indexing_rate;
|
||||
|
||||
|
||||
@@ -507,6 +507,10 @@ TEST_CASE("CBORSerialize_End", "[CBOR]") {
|
||||
},
|
||||
.rotation_lattice = CrystalLattice(40, 50, 60, 90, 90, 90)
|
||||
};
|
||||
message.protein_score = 0.75f;
|
||||
message.ice_score = 0.25f;
|
||||
message.v_protein_score = {0.5f, 0.625f};
|
||||
message.v_ice_score = {0.125f, 0.875f};
|
||||
|
||||
REQUIRE_NOTHROW(serializer.SerializeSequenceEnd(message));
|
||||
|
||||
@@ -526,6 +530,10 @@ TEST_CASE("CBORSerialize_End", "[CBOR]") {
|
||||
REQUIRE(output_message.run_number == message.run_number);
|
||||
REQUIRE(output_message.run_name == message.run_name);
|
||||
REQUIRE(output_message.az_int_result.empty());
|
||||
REQUIRE(output_message.protein_score == message.protein_score);
|
||||
REQUIRE(output_message.ice_score == message.ice_score);
|
||||
REQUIRE(output_message.v_protein_score == message.v_protein_score);
|
||||
REQUIRE(output_message.v_ice_score == message.v_ice_score);
|
||||
|
||||
REQUIRE(output_message.rotation_lattice_type.has_value());
|
||||
CHECK(output_message.rotation_lattice_type->centering == 'R');
|
||||
@@ -663,6 +671,8 @@ TEST_CASE("CBORSerialize_Image", "[CBOR]") {
|
||||
.spots = spots,
|
||||
.spot_count_ice_rings = 157,
|
||||
.bkg_estimate = 12.345f,
|
||||
.protein_score = 0.8125f,
|
||||
.ice_score = 0.375f,
|
||||
.indexing_result = true,
|
||||
.indexing_unit_cell = UnitCell{.a = 123, .b = 145, .c = 67.5, .alpha = 90, .beta = 120, .gamma = 134},
|
||||
.adu_histogram = {3, 4, 5, 8},
|
||||
@@ -716,6 +726,8 @@ TEST_CASE("CBORSerialize_Image", "[CBOR]") {
|
||||
REQUIRE(image_array.error_pixel_count == message.error_pixel_count);
|
||||
REQUIRE(image_array.strong_pixel_count == message.strong_pixel_count);
|
||||
REQUIRE(image_array.bkg_estimate == message.bkg_estimate);
|
||||
REQUIRE(image_array.protein_score == message.protein_score);
|
||||
REQUIRE(image_array.ice_score == message.ice_score);
|
||||
REQUIRE(image_array.image_collection_efficiency == message.image_collection_efficiency);
|
||||
REQUIRE(image_array.user_data == message.user_data);
|
||||
REQUIRE(image_array.original_number == message.original_number);
|
||||
|
||||
@@ -79,6 +79,7 @@ ADD_EXECUTABLE(jfjoch_test
|
||||
SpotExtractorGPUParityTest.cpp
|
||||
CalcBraggPredictionTest.cpp
|
||||
SpotUtilsTest.cpp
|
||||
DetectionScoreTest.cpp
|
||||
LatticeSearchTest.cpp
|
||||
TimeTest.cpp
|
||||
RotationIndexerTest.cpp
|
||||
|
||||
@@ -0,0 +1,205 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include <catch2/catch_all.hpp>
|
||||
|
||||
#include "../common/Definitions.h"
|
||||
#include "../image_analysis/IceScore.h"
|
||||
#include "../image_analysis/spot_finding/SpotUtils.h"
|
||||
|
||||
namespace {
|
||||
constexpr float TWO_PI = 6.283185307f;
|
||||
|
||||
SpotToSave spot(float d_A, float intensity) {
|
||||
return SpotToSave{.intensity = intensity, .d_A = d_A};
|
||||
}
|
||||
|
||||
AzimuthalIntegrationSettings ice_settings() {
|
||||
AzimuthalIntegrationSettings settings;
|
||||
settings.QSpacing_recipA(0.006f).QRange_recipA(0.1f, 4.5f);
|
||||
return settings;
|
||||
}
|
||||
|
||||
int bin_of(const AzimuthalIntegrationSettings &settings, float d_A) {
|
||||
return static_cast<int>(std::lround((TWO_PI / d_A - settings.GetLowQ_recipA())
|
||||
/ settings.GetQSpacing_recipA() - 0.5f));
|
||||
}
|
||||
}
|
||||
|
||||
TEST_CASE("ProteinScore_EmptyAndHighResolutionOnly") {
|
||||
CHECK(ProteinScore({}) == 0.0f);
|
||||
|
||||
// Nothing beyond 5 A: a salt or ice powder pattern, however strong, is not protein.
|
||||
std::vector<SpotToSave> spots;
|
||||
for (int i = 0; i < 200; i++)
|
||||
spots.push_back(spot(1.0f + 0.015f * static_cast<float>(i), 5000.0f));
|
||||
CHECK(ProteinScore(spots) == 0.0f);
|
||||
}
|
||||
|
||||
TEST_CASE("ProteinScore_SaturatesAndCountsShellsNotSpots") {
|
||||
// Six well-separated shells above 5 A, one strong spot each: enough evidence to be sure.
|
||||
std::vector<SpotToSave> shells;
|
||||
for (const float d: {5.5f, 6.5f, 8.0f, 10.0f, 14.0f, 20.0f})
|
||||
shells.push_back(spot(d, 1000.0f));
|
||||
const float six = ProteinScore(shells);
|
||||
CHECK(six > 0.5f);
|
||||
|
||||
// Ten times as many spots, all in the SAME shell: a parasitic ring, not a lattice. One shell is
|
||||
// worth at most 1, so it must score far below the six shells above.
|
||||
std::vector<SpotToSave> one_ring;
|
||||
for (int i = 0; i < 60; i++)
|
||||
one_ring.push_back(spot(5.5f, 1000.0f));
|
||||
CHECK(ProteinScore(one_ring) < 0.3f);
|
||||
CHECK(ProteinScore(one_ring) < six);
|
||||
|
||||
// Saturation: making every spot a hundred times stronger does not raise the score, because the
|
||||
// weight is measured against the frame's own median spot.
|
||||
std::vector<SpotToSave> strong;
|
||||
for (const float d: {5.5f, 6.5f, 8.0f, 10.0f, 14.0f, 20.0f})
|
||||
strong.push_back(spot(d, 100000.0f));
|
||||
CHECK(ProteinScore(strong) == Catch::Approx(six));
|
||||
|
||||
// And it stays inside [0, 1] however much evidence there is.
|
||||
std::vector<SpotToSave> many;
|
||||
for (int i = 0; i < 400; i++)
|
||||
many.push_back(spot(5.1f + 0.08f * static_cast<float>(i % 300), 1000.0f));
|
||||
CHECK(ProteinScore(many) > 0.9f);
|
||||
CHECK(ProteinScore(many) < 1.0f);
|
||||
}
|
||||
|
||||
TEST_CASE("IceScoreRadial_FlatProfileIsNotIce") {
|
||||
const auto settings = ice_settings();
|
||||
const int q_bins = settings.GetQBinCount();
|
||||
const std::vector<float> flat(q_bins, 100.0f);
|
||||
// Per-pixel standard deviation and pixel count: the score uses std / sqrt(count) = 2 photons.
|
||||
const std::vector<float> sigma(q_bins, 20.0f);
|
||||
const std::vector<uint64_t> count(q_bins, 100);
|
||||
CHECK(IceScoreRadial(flat, sigma, count, q_bins, settings) == 0.0f);
|
||||
|
||||
// No standard deviation, no radial channel: the FPGA azimuthal integration does not produce one.
|
||||
CHECK(IceScoreRadial(flat, {}, count, q_bins, settings) == 0.0f);
|
||||
}
|
||||
|
||||
TEST_CASE("IceScoreRadial_HexagonalAndCubicPatterns") {
|
||||
const auto settings = ice_settings();
|
||||
const int q_bins = settings.GetQBinCount();
|
||||
const std::vector<float> sigma(q_bins, 20.0f);
|
||||
const std::vector<uint64_t> count(q_bins, 100);
|
||||
|
||||
// A single strong band is not ice - real ice shows a whole pattern, and the band-count
|
||||
// concordance test is what refuses one bin.
|
||||
std::vector<float> one_band(q_bins, 100.0f);
|
||||
one_band[bin_of(settings, ICE_RING_RES_A[0])] = 400.0f;
|
||||
CHECK(IceScoreRadial(one_band, sigma, count, q_bins, settings) < 0.5f);
|
||||
|
||||
// The whole hexagonal pattern is.
|
||||
std::vector<float> hexagonal(q_bins, 100.0f);
|
||||
for (const float d: ICE_RING_RES_A) {
|
||||
const int b = bin_of(settings, d);
|
||||
if (b >= 0 && b < q_bins)
|
||||
hexagonal[b] = 130.0f;
|
||||
}
|
||||
CHECK(IceScoreRadial(hexagonal, sigma, count, q_bins, settings) > 0.5f);
|
||||
|
||||
// So is the cubic one, which shares only three lines with it - the phase that a hexagonal-only
|
||||
// detector misses entirely.
|
||||
std::vector<float> cubic(q_bins, 100.0f);
|
||||
for (const float d: ICE_RING_CUBIC_RES_A) {
|
||||
const int b = bin_of(settings, d);
|
||||
if (b >= 0 && b < q_bins)
|
||||
cubic[b] = 130.0f;
|
||||
}
|
||||
CHECK(IceScoreRadial(cubic, sigma, count, q_bins, settings) > 0.5f);
|
||||
|
||||
// The same excess spread over bins belonging to no phase is not ice.
|
||||
std::vector<float> off_band(q_bins, 100.0f);
|
||||
for (int i = 100; i < q_bins - 100; i += 37)
|
||||
off_band[i] = 130.0f;
|
||||
CHECK(IceScoreRadial(off_band, sigma, count, q_bins, settings) < 0.5f);
|
||||
}
|
||||
|
||||
TEST_CASE("IceScoreRadial_AzimuthalProfileFoldsToTheSameAnswer") {
|
||||
AzimuthalIntegrationSettings settings;
|
||||
settings.QSpacing_recipA(0.006f).QRange_recipA(0.1f, 4.5f).AzimuthalBinCount(4);
|
||||
const int q_bins = settings.GetQBinCount();
|
||||
|
||||
std::vector<float> flat(q_bins, 100.0f);
|
||||
std::vector<float> sigma(q_bins, 20.0f);
|
||||
std::vector<uint64_t> count(q_bins, 400);
|
||||
for (const float d: ICE_RING_RES_A) {
|
||||
const int b = bin_of(settings, d);
|
||||
if (b >= 0 && b < q_bins)
|
||||
flat[b] = 130.0f;
|
||||
}
|
||||
std::vector<float> sectors(static_cast<size_t>(q_bins) * 4);
|
||||
std::vector<float> sectors_sigma(static_cast<size_t>(q_bins) * 4);
|
||||
std::vector<uint64_t> sectors_count(static_cast<size_t>(q_bins) * 4);
|
||||
for (int az = 0; az < 4; az++)
|
||||
for (int q = 0; q < q_bins; q++) {
|
||||
sectors[static_cast<size_t>(az) * q_bins + q] = flat[q];
|
||||
sectors_sigma[static_cast<size_t>(az) * q_bins + q] = sigma[q];
|
||||
sectors_count[static_cast<size_t>(az) * q_bins + q] = count[q] / 4;
|
||||
}
|
||||
|
||||
CHECK(IceScoreRadial(sectors, sectors_sigma, sectors_count, q_bins, settings)
|
||||
== Catch::Approx(IceScoreRadial(flat, sigma, count, q_bins, settings)));
|
||||
}
|
||||
|
||||
TEST_CASE("IceScoreSpots_ExcessOnTheIceRadii") {
|
||||
const auto settings = ice_settings();
|
||||
const int q_bins = settings.GetQBinCount();
|
||||
// A detector that covers every radius equally, so the control offsets are directly comparable.
|
||||
const std::vector<uint64_t> count(q_bins, 10000);
|
||||
|
||||
// 400 spots spread evenly in q: whatever lands on an ice radius is what the control predicts.
|
||||
std::vector<SpotToSave> even;
|
||||
const float q_lo = 1.3f, q_hi = 4.2f;
|
||||
for (int i = 0; i < 400; i++) {
|
||||
const float q = q_lo + (q_hi - q_lo) * static_cast<float>(i) / 399.0f;
|
||||
even.push_back(spot(TWO_PI / q, 1000.0f));
|
||||
}
|
||||
CHECK(IceScoreSpots(even, count, q_bins, settings) < 0.5f);
|
||||
|
||||
// The same frame with 10 extra spots planted on each hexagonal radius.
|
||||
std::vector<SpotToSave> with_ice = even;
|
||||
for (const float d: ICE_RING_RES_A)
|
||||
for (int i = 0; i < 10; i++)
|
||||
with_ice.push_back(spot(d, 1000.0f));
|
||||
CHECK(IceScoreSpots(with_ice, count, q_bins, settings) > 0.5f);
|
||||
|
||||
// Two spots that both happen to sit on a ring are not ice: the ratio term refuses them even
|
||||
// though the Poisson tail alone would not.
|
||||
std::vector<SpotToSave> two;
|
||||
two.push_back(spot(ICE_RING_RES_A[0], 1000.0f));
|
||||
two.push_back(spot(ICE_RING_RES_A[1], 1000.0f));
|
||||
CHECK(IceScoreSpots(two, count, q_bins, settings) < 0.5f);
|
||||
|
||||
CHECK(IceScoreSpots({}, count, q_bins, settings) == 0.0f);
|
||||
}
|
||||
|
||||
TEST_CASE("IceScore_TakesTheStrongerChannel") {
|
||||
const auto settings = ice_settings();
|
||||
const int q_bins = settings.GetQBinCount();
|
||||
const std::vector<uint64_t> count(q_bins, 10000);
|
||||
const std::vector<float> sigma(q_bins, 200.0f);
|
||||
|
||||
// Powder ice, no spots at all: the radial channel carries it on its own.
|
||||
std::vector<float> powder(q_bins, 100.0f);
|
||||
for (const float d: ICE_RING_RES_A) {
|
||||
const int b = bin_of(settings, d);
|
||||
if (b >= 0 && b < q_bins)
|
||||
powder[b] = 130.0f;
|
||||
}
|
||||
CHECK(IceScore(powder, sigma, count, q_bins, settings, {}) > 0.5f);
|
||||
|
||||
// Ice as discrete crystallites: the profile is flat and only the spot channel sees it.
|
||||
const std::vector<float> flat(q_bins, 100.0f);
|
||||
std::vector<SpotToSave> textured;
|
||||
for (int i = 0; i < 400; i++)
|
||||
textured.push_back(spot(TWO_PI / (1.3f + 2.9f * static_cast<float>(i) / 399.0f), 1000.0f));
|
||||
for (const float d: ICE_RING_RES_A)
|
||||
for (int i = 0; i < 10; i++)
|
||||
textured.push_back(spot(d, 1000.0f));
|
||||
CHECK(IceScoreRadial(flat, sigma, count, q_bins, settings) == 0.0f);
|
||||
CHECK(IceScore(flat, sigma, count, q_bins, settings, textured) > 0.5f);
|
||||
}
|
||||
@@ -54,6 +54,8 @@ HDF5DataFilePluginMX::HDF5DataFilePluginMX(const StartMessage &msg)
|
||||
void HDF5DataFilePluginMX::OpenFile(HDF5File &data_file, const DataMessage &msg, size_t images_per_file) {
|
||||
bkg_estimate.reserve(images_per_file);
|
||||
ice_ring_score.reserve(images_per_file);
|
||||
protein_score.reserve(images_per_file);
|
||||
ice_score.reserve(images_per_file);
|
||||
|
||||
if (max_spots == 0)
|
||||
return;
|
||||
@@ -101,6 +103,10 @@ void HDF5DataFilePluginMX::Write(const DataMessage &msg, uint64_t image_number)
|
||||
bkg_estimate[image_number] = msg.bkg_estimate.value();
|
||||
if (msg.ice_ring_score.has_value())
|
||||
ice_ring_score[image_number] = msg.ice_ring_score.value();
|
||||
if (msg.protein_score.has_value())
|
||||
protein_score[image_number] = msg.protein_score.value();
|
||||
if (msg.ice_score.has_value())
|
||||
ice_score[image_number] = msg.ice_score.value();
|
||||
|
||||
if (max_spots == 0)
|
||||
return;
|
||||
@@ -252,6 +258,10 @@ void HDF5DataFilePluginMX::WriteFinal(HDF5File &data_file) {
|
||||
data_file.SaveVector("/entry/MX/bkgEstimate", bkg_estimate.vec());
|
||||
if (!ice_ring_score.empty())
|
||||
data_file.SaveVector("/entry/MX/iceRingScore", ice_ring_score.vec());
|
||||
if (!protein_score.empty())
|
||||
data_file.SaveVector("/entry/MX/proteinScore", protein_score.vec());
|
||||
if (!ice_score.empty())
|
||||
data_file.SaveVector("/entry/MX/iceScore", ice_score.vec());
|
||||
if (!profile_radius.empty())
|
||||
data_file.SaveVector("/entry/MX/profileRadius", profile_radius.vec())->Units("Angstrom^-1");
|
||||
if (!mosaicity_deg.empty())
|
||||
|
||||
@@ -48,6 +48,8 @@ class HDF5DataFilePluginMX : public HDF5DataFilePlugin {
|
||||
// bkg_estimate
|
||||
AutoIncrVector<float> bkg_estimate{NAN};
|
||||
AutoIncrVector<float> ice_ring_score{NAN};
|
||||
AutoIncrVector<float> protein_score{NAN};
|
||||
AutoIncrVector<float> ice_score{NAN};
|
||||
|
||||
// resolution_estimation
|
||||
AutoIncrVector<float> resolution_estimate{NAN};
|
||||
|
||||
@@ -1182,6 +1182,20 @@ void NXmx::Finalize(const EndMessage &end) {
|
||||
if (end.ice_ring_score_mean) {
|
||||
SaveScalar(*hdf5_file, "/entry/MX/iceRingScoreMean", end.ice_ring_score_mean.value());
|
||||
}
|
||||
if (end.protein_score) {
|
||||
SaveScalar(*hdf5_file, "/entry/MX/proteinScoreMean", end.protein_score.value());
|
||||
}
|
||||
if (end.ice_score) {
|
||||
SaveScalar(*hdf5_file, "/entry/MX/iceScoreMean", end.ice_score.value());
|
||||
}
|
||||
// The beam centre the two detection scores were computed with. Both read d out of the geometry,
|
||||
// and beam_center_x/y in /entry/instrument/detector above is the REFINED centre when refinement
|
||||
// ran - so without this pair a later rescoring could not tell a disagreement about the algorithm
|
||||
// from a disagreement about the geometry.
|
||||
if (!end.v_protein_score.empty() || !end.v_ice_score.empty()) {
|
||||
SaveScalar(*hdf5_file, "/entry/MX/scoreBeamCenterX", start_message.beam_center_x)->Units("pixel");
|
||||
SaveScalar(*hdf5_file, "/entry/MX/scoreBeamCenterY", start_message.beam_center_y)->Units("pixel");
|
||||
}
|
||||
|
||||
hdf5_file->Close();
|
||||
hdf5_file.reset();
|
||||
@@ -1258,6 +1272,8 @@ void NXmx::EndResultVectors(const EndMessage &end) {
|
||||
SaveVectorIfMissing(*hdf5_file, "/entry/MX/imageIndexed", end.image_indexed);
|
||||
SaveVectorIfMissing(*hdf5_file, "/entry/MX/indexedLatticeCount", end.indexed_lattice_count);
|
||||
SaveVectorIfMissing(*hdf5_file, "/entry/MX/bkgEstimate", end.v_bkg_estimate);
|
||||
SaveVectorIfMissing(*hdf5_file, "/entry/MX/proteinScore", end.v_protein_score);
|
||||
SaveVectorIfMissing(*hdf5_file, "/entry/MX/iceScore", end.v_ice_score);
|
||||
SaveVectorIfMissing(*hdf5_file, "/entry/MX/iceRingScore", end.ice_ring_score);
|
||||
SaveVectorIfMissing(*hdf5_file, "/entry/MX/profileRadius", end.profile_radius, "Angstrom^-1");
|
||||
SaveVectorIfMissing(*hdf5_file, "/entry/MX/mosaicity", end.mosaicity, "deg");
|
||||
|
||||
Reference in New Issue
Block a user