diff --git a/broker/OpenAPIConvert.cpp b/broker/OpenAPIConvert.cpp index e44b9fca2..85ef537f8 100644 --- a/broker/OpenAPIConvert.cpp +++ b/broker/OpenAPIConvert.cpp @@ -899,6 +899,8 @@ PlotType ConvertPlotType(const std::optional& 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()) diff --git a/broker/jfjoch_api.yaml b/broker/jfjoch_api.yaml index 7b3fe7587..26051b7e5 100644 --- a/broker/jfjoch_api.yaml +++ b/broker/jfjoch_api.yaml @@ -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 diff --git a/common/Definitions.h b/common/Definitions.h index 140d48836..00eda7822 100644 --- a/common/Definitions.h +++ b/common/Definitions.h @@ -83,6 +83,15 @@ constexpr std::array 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 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. diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index ef2cd0601..cd0de3f66 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -129,6 +129,16 @@ struct DataMessage { std::optional bkg_estimate; std::optional 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 protein_score; + std::optional ice_score; + std::optional indexing_result; std::optional indexing_lattice; std::vector indexing_extra_lattices; @@ -369,6 +379,10 @@ struct EndMessage { std::optional efficiency; std::optional indexing_rate; std::optional bkg_estimate; + // Run means of the two per-image detection scores; written to /entry/MX/proteinScoreMean and + // /entry/MX/iceScoreMean. + std::optional protein_score; + std::optional ice_score; std::optional end_date; @@ -419,6 +433,8 @@ struct EndMessage { std::vector image_indexed; std::vector indexed_lattice_count; std::vector v_bkg_estimate; + std::vector v_protein_score; + std::vector v_ice_score; std::vector profile_radius; std::vector mosaicity; std::vector bFactor; diff --git a/common/JFJochReceiverPlots.cpp b/common/JFJochReceiverPlots.cpp index 239606614..e972a0ad7 100644 --- a/common/JFJochReceiverPlots.cpp +++ b/common/JFJochReceiverPlots.cpp @@ -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 JFJochReceiverPlots::GetIceRingScoreArray() const { return ice_ring_score.ExportArray(); } +std::optional JFJochReceiverPlots::GetProteinScore() const { + auto tmp = protein_score.Mean(); + if (std::isfinite(tmp)) + return tmp; + else + return {}; +} + +std::optional JFJochReceiverPlots::GetIceScore() const { + auto tmp = ice_score.Mean(); + if (std::isfinite(tmp)) + return tmp; + else + return {}; +} + +std::vector JFJochReceiverPlots::GetProteinScoreArray() const { + return protein_score.ExportArray(); +} + +std::vector 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 &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; diff --git a/common/JFJochReceiverPlots.h b/common/JFJochReceiverPlots.h index a7ebdddf3..642071e7d 100644 --- a/common/JFJochReceiverPlots.h +++ b/common/JFJochReceiverPlots.h @@ -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 GetResolutionEstimate() const; std::vector GetIceRingScoreArray() const; + // Run means of the two per-image detection scores, and the per-image arrays behind them. + [[nodiscard]] std::optional GetProteinScore() const; + [[nodiscard]] std::optional GetIceScore() const; + [[nodiscard]] std::vector GetProteinScoreArray() const; + [[nodiscard]] std::vector GetIceScoreArray() const; std::vector GetAzIntProfile() const; // The run-summed profile object itself, rather than the array GetAzIntProfile() flattens it to. diff --git a/common/Plot.h b/common/Plot.h index b629c1b44..af2321d9e 100644 --- a/common/Plot.h +++ b/common/Plot.h @@ -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 { diff --git a/common/ScanResult.h b/common/ScanResult.h index 92a409bdc..c695b714b 100644 --- a/common/ScanResult.h +++ b/common/ScanResult.h @@ -44,6 +44,8 @@ struct ScanResultElem { std::optional image_scale_factor; std::optional image_scale_cc; std::optional ice_ring_score; + std::optional protein_score; + std::optional ice_score; }; struct ScanResult { diff --git a/common/ScanResultGenerator.cpp b/common/ScanResultGenerator.cpp index 1510da683..b49854083 100644 --- a/common/ScanResultGenerator.cpp +++ b/common/ScanResultGenerator.cpp @@ -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(value_or_zero(e.integrated_reflections)); message.niggli_class[number] = static_cast(value_or_zero(e.niggli_class)); message.pixel_sum[number] = value_or_zero(e.pixel_sum); diff --git a/docs/CBOR.md b/docs/CBOR.md index 095c256a8..5ea0201ac 100644 --- a/docs/CBOR.md +++ b/docs/CBOR.md @@ -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\] | | diff --git a/docs/CPU_DATA_ANALYSIS_IMAGE.md b/docs/CPU_DATA_ANALYSIS_IMAGE.md index d3e9aef3e..2df1851ad 100644 --- a/docs/CPU_DATA_ANALYSIS_IMAGE.md +++ b/docs/CPU_DATA_ANALYSIS_IMAGE.md @@ -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. diff --git a/docs/HDF5.md b/docs/HDF5.md index 1bd2c6ebb..c8faa43b3 100644 --- a/docs/HDF5.md +++ b/docs/HDF5.md @@ -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 | diff --git a/frame_serialize/CBORStream2Deserializer.cpp b/frame_serialize/CBORStream2Deserializer.cpp index 658e5a1cb..7bee5b144 100644 --- a/frame_serialize/CBORStream2Deserializer.cpp +++ b/frame_serialize/CBORStream2Deserializer.cpp @@ -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") diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index 7a17da3f9..2de5d4c7a 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -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); diff --git a/frontend/src/components/DataProcessingPlot.tsx b/frontend/src/components/DataProcessingPlot.tsx index 8843dc1b2..1643d1497 100644 --- a/frontend/src/components/DataProcessingPlot.tsx +++ b/frontend/src/components/DataProcessingPlot.tsx @@ -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: diff --git a/frontend/src/components/DataProcessingPlots.tsx b/frontend/src/components/DataProcessingPlots.tsx index 7e950811c..378864a21 100644 --- a/frontend/src/components/DataProcessingPlots.tsx +++ b/frontend/src/components/DataProcessingPlots.tsx @@ -51,6 +51,8 @@ function DataProcessingPlots({type: initialType, height}: MyProps) { Spot count indexed Spot count ice ring Ice ring score + Protein detection score + Ice detection score Azimuthal integration profile Azimuthal integration profile (1D) Background estimate diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index 2bf22bd3a..017616aac 100644 --- a/image_analysis/CMakeLists.txt +++ b/image_analysis/CMakeLists.txt @@ -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 diff --git a/image_analysis/IceScore.cpp b/image_analysis/IceScore.cpp new file mode 100644 index 000000000..3e4aca46d --- /dev/null +++ b/image_analysis/IceScore.cpp @@ -0,0 +1,333 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include + +#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 FoldToQ(const std::vector &in, int nq) { + std::vector out(nq, 0.0f); + std::vector n(nq, 0); + for (size_t i = 0; i < in.size() && nq > 0; i++) { + const int q = static_cast(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(n[q]) : NAN; + return out; + } + + // The same fold for the live pixel count, which adds rather than averages. + std::vector FoldCountToQ(const std::vector &in, int nq) { + std::vector out(nq, 0.0); + for (size_t i = 0; i < in.size() && nq > 0; i++) + out[i % nq] += static_cast(in[i]); + return out; + } + + // Running median, NaN where the window holds nothing finite. + std::vector RunningMedian(const std::vector &v, int half) { + std::vector out(v.size(), NAN); + std::vector window; + for (int i = 0; i < static_cast(v.size()); i++) { + window.clear(); + const int lo = std::max(0, i - half); + const int hi = std::min(static_cast(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 &c, const std::vector &band_q, + const std::vector &band_weight, int min_bands, + float low_q, float dq, int bh, + float mu0, float v0, float p0) { + const int nq = static_cast(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(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(nbands); + const float s_cnt = (static_cast(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(j) * log_mu - std::lgamma(static_cast(j) + 1.0)); + sum += term; + if (static_cast(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 &count_q, float low_q, float dq, float q) { + double area = 0.0; + const int lo = static_cast(std::ceil((q - SPOT_BAND_HALF_Q - low_q) / dq - 0.5f)); + const int hi = static_cast(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(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 &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(std::distance(lo, hi)); + } +} + +float IceScoreRadial(const std::vector &profile, const std::vector &profile_std, + const std::vector &profile_count, int32_t q_bins, + const AzimuthalIntegrationSettings &settings) { + const int nq = std::max(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 prof = FoldToQ(profile, nq); + std::vector sigma = FoldToQ(profile_std, nq); + const std::vector 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(std::sqrt(count[i])); + else + sigma[i] = NAN; + } + + const std::vector background = RunningMedian(prof, BACKGROUND_HALF_WINDOW); + const std::vector sigma_smooth = RunningMedian(sigma, SIGMA_HALF_WINDOW); + const int bh = std::max(1, static_cast(std::lround(BAND_HALF_Q / dq))); + auto q_of = [&](int i) { return low_q + (static_cast(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 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 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 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(bh) * dq + NULL_EXCLUSION_Q) + null_c.push_back(c[i]); + } + if (static_cast(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(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 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 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 &spots, const std::vector &profile_count, + int32_t q_bins, const AzimuthalIntegrationSettings &settings) { + const int nq = std::max(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 count_q = FoldCountToQ(profile_count, nq); + + std::vector 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 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 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(control.size()) < SPOT_MIN_OFFSETS) + continue; + double mean = 0.0; + for (const double v: control) + mean += v; + mean /= static_cast(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(control.size()); + } + + mu = std::max(mu, static_cast(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(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(std::max(0.0, std::min(s_sig, s_size))); +} + +float IceScore(const std::vector &profile, const std::vector &profile_std, + const std::vector &profile_count, int32_t q_bins, + const AzimuthalIntegrationSettings &settings, const std::vector &spots) { + return std::max(IceScoreRadial(profile, profile_std, profile_count, q_bins, settings), + IceScoreSpots(spots, profile_count, q_bins, settings)); +} diff --git a/image_analysis/IceScore.h b/image_analysis/IceScore.h new file mode 100644 index 000000000..31c7f53cd --- /dev/null +++ b/image_analysis/IceScore.h @@ -0,0 +1,61 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#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 &profile, const std::vector &profile_std, + const std::vector &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 &spots, const std::vector &profile_count, + int32_t q_bins, const AzimuthalIntegrationSettings &settings); + +// The score itself: whichever channel sees more. +float IceScore(const std::vector &profile, const std::vector &profile_std, + const std::vector &profile_count, int32_t q_bins, + const AzimuthalIntegrationSettings &settings, const std::vector &spots); diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 57849fb15..8f21bb823 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -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() { diff --git a/image_analysis/spot_finding/SpotUtils.cpp b/image_analysis/spot_finding/SpotUtils.cpp index 433eb6610..23b78ee77 100644 --- a/image_analysis/spot_finding/SpotUtils.cpp +++ b/image_analysis/spot_finding/SpotUtils.cpp @@ -1,6 +1,8 @@ // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only +#include + #include "../../common/JFJochMath.h" #include "SpotUtils.h" #include "../../common/ResolutionShells.h" @@ -167,6 +169,53 @@ std::optional GetResolution(const std::vector &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 &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 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 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(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 &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. diff --git a/image_analysis/spot_finding/SpotUtils.h b/image_analysis/spot_finding/SpotUtils.h index 1eaa5f03e..5b367cd07 100644 --- a/image_analysis/spot_finding/SpotUtils.h +++ b/image_analysis/spot_finding/SpotUtils.h @@ -43,6 +43,25 @@ void FilterSpuriousHighResolutionSpots(std::vector &spots, float thr // Returns nothing when the image has too few spots to have a fall-off at all. std::optional GetResolution(const std::vector &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 &spots); + void SpotAnalyze(const DiffractionExperiment &experiment, const SpotFindingSettings &settings, const std::vector &spots, diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index 8b12c1fab..7a4d465c3 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -558,6 +558,10 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen dataset->indexing_result = master_file->ReadOptVector("/entry/MX/imageIndexed"); dataset->bkg_estimate = master_file->ReadOptVector("/entry/MX/bkgEstimate"); dataset->ice_ring_score = master_file->ReadOptVector("/entry/MX/iceRingScore"); + dataset->protein_score = master_file->ReadOptVector("/entry/MX/proteinScore"); + dataset->ice_score = master_file->ReadOptVector("/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("/entry/MX/resolutionEstimate"); dataset->profile_radius = master_file->ReadOptVector("/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) diff --git a/reader/JFJochReaderDataset.h b/reader/JFJochReaderDataset.h index ab6d36e8a..e33da758d 100644 --- a/reader/JFJochReaderDataset.h +++ b/reader/JFJochReaderDataset.h @@ -57,6 +57,14 @@ struct JFJochReaderDataset { std::vector indexing_lattice_count; std::vector bkg_estimate; std::vector 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 protein_score; + std::vector ice_score; + std::optional score_beam_center_x; + std::optional score_beam_center_y; std::vector resolution_estimate; std::vector efficiency; std::vector profile_radius; diff --git a/receiver/JFJochReceiver.cpp b/receiver/JFJochReceiver.cpp index 7c1a21799..8c64778ea 100644 --- a/receiver/JFJochReceiver.cpp +++ b/receiver/JFJochReceiver.cpp @@ -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(); diff --git a/receiver/JFJochReceiverFPGA.cpp b/receiver/JFJochReceiverFPGA.cpp index 23288259e..f5e4013d7 100644 --- a/receiver/JFJochReceiverFPGA.cpp +++ b/receiver/JFJochReceiverFPGA.cpp @@ -6,6 +6,7 @@ #include #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); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index a75631ee9..9fb564f54 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -2,6 +2,7 @@ // SPDX-License-Identifier: GPL-3.0-only #include +#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; diff --git a/tests/CBORTest.cpp b/tests/CBORTest.cpp index 4c36ea360..6425bead0 100644 --- a/tests/CBORTest.cpp +++ b/tests/CBORTest.cpp @@ -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); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index a4af137b6..7fcfdb5f2 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -78,6 +78,7 @@ ADD_EXECUTABLE(jfjoch_test SpotExtractorGPUParityTest.cpp CalcBraggPredictionTest.cpp SpotUtilsTest.cpp + DetectionScoreTest.cpp LatticeSearchTest.cpp TimeTest.cpp RotationIndexerTest.cpp diff --git a/tests/DetectionScoreTest.cpp b/tests/DetectionScoreTest.cpp new file mode 100644 index 000000000..7f8b79d93 --- /dev/null +++ b/tests/DetectionScoreTest.cpp @@ -0,0 +1,205 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#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(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 spots; + for (int i = 0; i < 200; i++) + spots.push_back(spot(1.0f + 0.015f * static_cast(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 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 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 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 many; + for (int i = 0; i < 400; i++) + many.push_back(spot(5.1f + 0.08f * static_cast(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 flat(q_bins, 100.0f); + // Per-pixel standard deviation and pixel count: the score uses std / sqrt(count) = 2 photons. + const std::vector sigma(q_bins, 20.0f); + const std::vector 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 sigma(q_bins, 20.0f); + const std::vector 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 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 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 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 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 flat(q_bins, 100.0f); + std::vector sigma(q_bins, 20.0f); + std::vector 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 sectors(static_cast(q_bins) * 4); + std::vector sectors_sigma(static_cast(q_bins) * 4); + std::vector sectors_count(static_cast(q_bins) * 4); + for (int az = 0; az < 4; az++) + for (int q = 0; q < q_bins; q++) { + sectors[static_cast(az) * q_bins + q] = flat[q]; + sectors_sigma[static_cast(az) * q_bins + q] = sigma[q]; + sectors_count[static_cast(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 count(q_bins, 10000); + + // 400 spots spread evenly in q: whatever lands on an ice radius is what the control predicts. + std::vector 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(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 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 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 count(q_bins, 10000); + const std::vector sigma(q_bins, 200.0f); + + // Powder ice, no spots at all: the radial channel carries it on its own. + std::vector 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 flat(q_bins, 100.0f); + std::vector textured; + for (int i = 0; i < 400; i++) + textured.push_back(spot(TWO_PI / (1.3f + 2.9f * static_cast(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); +} diff --git a/writer/HDF5DataFilePluginMX.cpp b/writer/HDF5DataFilePluginMX.cpp index 808a2006a..4a3f1823a 100644 --- a/writer/HDF5DataFilePluginMX.cpp +++ b/writer/HDF5DataFilePluginMX.cpp @@ -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()) diff --git a/writer/HDF5DataFilePluginMX.h b/writer/HDF5DataFilePluginMX.h index 68a568edb..5f3bbc043 100644 --- a/writer/HDF5DataFilePluginMX.h +++ b/writer/HDF5DataFilePluginMX.h @@ -48,6 +48,8 @@ class HDF5DataFilePluginMX : public HDF5DataFilePlugin { // bkg_estimate AutoIncrVector bkg_estimate{NAN}; AutoIncrVector ice_ring_score{NAN}; + AutoIncrVector protein_score{NAN}; + AutoIncrVector ice_score{NAN}; // resolution_estimation AutoIncrVector resolution_estimate{NAN}; diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp index d1eecac63..829e5e001 100644 --- a/writer/HDF5NXmx.cpp +++ b/writer/HDF5NXmx.cpp @@ -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");