diff --git a/broker/OpenAPIConvert.cpp b/broker/OpenAPIConvert.cpp index 99290eeb8..cf40188ca 100644 --- a/broker/OpenAPIConvert.cpp +++ b/broker/OpenAPIConvert.cpp @@ -899,7 +899,7 @@ PlotType ConvertPlotType(const std::optional& input) { "Plot type is compulsory paramater"); if (input == "bkg_estimate") return PlotType::BkgEstimate; if (input == "spindle_blind_fraction") return PlotType::SpindleBlindFraction; - if (input == "ice_ring_score") return PlotType::IceRingScore; + if (input == "ice_ring_score") return PlotType::IceRingRatio; if (input == "protein_score") return PlotType::ProteinScore; if (input == "ice_score") return PlotType::IceScore; if (input == "azint") return PlotType::AzInt; @@ -1152,8 +1152,8 @@ org::openapitools::server::model::Scan_result Convert(const ScanResult& input) { tmp.setSpots(i.spot_count.value()); if (i.spot_count_ice.has_value()) tmp.setSpotsIce(i.spot_count_ice.value()); - if (i.ice_ring_score.has_value()) - tmp.setIce(i.ice_ring_score.value()); + if (i.ice_ring_ratio.has_value()) + tmp.setIce(i.ice_ring_ratio.value()); if (i.protein_score.has_value()) tmp.setProteinScore(i.protein_score.value()); if (i.ice_score.has_value()) diff --git a/common/AzimuthalIntegrationProfile.cpp b/common/AzimuthalIntegrationProfile.cpp index 412eea018..65d9ba053 100644 --- a/common/AzimuthalIntegrationProfile.cpp +++ b/common/AzimuthalIntegrationProfile.cpp @@ -178,12 +178,12 @@ float AzimuthalIntegrationProfile::GetBkgEstimate(const AzimuthalIntegrationSett return GetMeanValueOfBins(min_bin, max_bin); } -float AzimuthalIntegrationProfile::GetIceRingScore(const AzimuthalIntegrationSettings &settings, +float AzimuthalIntegrationProfile::GetIceRingRatio(const AzimuthalIntegrationSettings &settings, float half_width_q) const { - return IceRingScore(GetResult1D(), q_bins, settings, half_width_q); + return IceRingRatio(GetResult1D(), q_bins, settings, half_width_q); } -float AzimuthalIntegrationProfile::IceRingScore(const std::vector &profile, int32_t q_bins, +float AzimuthalIntegrationProfile::IceRingRatio(const std::vector &profile, int32_t q_bins, const AzimuthalIntegrationSettings &settings, float half_width_q) { // Strongest hexagonal-ice ring's intensity relative to the background *under* it (1 = no ice). The diff --git a/common/AzimuthalIntegrationProfile.h b/common/AzimuthalIntegrationProfile.h index 7eb6fcd2d..2ab523237 100644 --- a/common/AzimuthalIntegrationProfile.h +++ b/common/AzimuthalIntegrationProfile.h @@ -46,11 +46,11 @@ public: float GetBkgEstimate(const AzimuthalIntegrationSettings& settings) const; // Single per-image ice indicator: the strongest hexagonal-ice ring's intensity relative to the smooth // radial background interpolated under it (1 = no ice, >1 = ice). See the .cpp for the background fit. - float GetIceRingScore(const AzimuthalIntegrationSettings& settings, float half_width_q) const; + float GetIceRingRatio(const AzimuthalIntegrationSettings& settings, float half_width_q) const; // The same score for a radial profile that is not this object's: the adaptive spot finder computes a // peak-excluded per-ring background in these same bins, which is a much cleaner input (see the .cpp). // profile is q_bins long, or q_bins x azimuthal bins, in which case it is averaged over azimuth first. - static float IceRingScore(const std::vector& profile, int32_t q_bins, + static float IceRingRatio(const std::vector& profile, int32_t q_bins, const AzimuthalIntegrationSettings& settings, float half_width_q); MultiLinePlot GetPlot(bool force_1d = false, PlotAzintUnit plot_unit = PlotAzintUnit::Q_recipA) const; AzimuthalIntegrationProfile& operator+=(const AzimuthalIntegrationProfile& profile); // Not thread safe diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index eccd8f4ff..514d2863d 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -114,7 +114,7 @@ struct DataMessage { // bands' own q width - so spot_count_ice_rings / spot_count_ice_control is 1 when spots are spread // evenly and > 1 when they pile up on the rings. Pooled over the run (a ratio of two per-image // counts this small is meaningless on its own image), it is the second ice channel: TEXTURED ice - // arrives as discrete spots and leaves the radial profile - and so ice_ring_score - flat. + // arrives as discrete spots and leaves the radial profile - and so ice_ring_ratio - flat. std::optional spot_count_ice_control; std::optional spot_count_low_res; @@ -127,7 +127,10 @@ struct DataMessage { std::vector az_int_profile_count; std::optional bkg_estimate; - std::optional ice_ring_score; // strongest ice ring over the smooth radial background (1 = none) + // Strongest ice ring over the smooth radial background. A RATIO, so 1 means no ice and it is + // unbounded above - the opposite convention to the two *_score members below, which is why the + // name says ratio. Still travels the wire and the file as "ice_ring_score"/iceRingScore. + std::optional ice_ring_ratio; // How much of a single sweep's blind cone this orientation makes unrecoverable: 0 = one sweep // about the spindle reaches everything the point group can give, 1 = a short lattice row lies on // the spindle and the whole cone is lost coherently - each shell losing its own 1 - cos theta(d), @@ -140,7 +143,7 @@ struct DataMessage { // 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 + // term. ice_score answers the question ice_ring_ratio 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. @@ -470,12 +473,12 @@ struct EndMessage { std::vector image_scale_factor; std::vector image_scale_cc; std::vector image_scale_mosaicity; - // Per-image ice strength. Note the name carries no v_ prefix, unlike v_bkg_estimate above - it - // was shipped that way and the CBOR key is part of the stream format. - std::vector ice_ring_score; + // Per-image ice ring ratio (1 = no ice). Note the name carries no v_ prefix, unlike + // v_bkg_estimate above - it was shipped that way and the CBOR key is part of the stream format. + std::vector ice_ring_ratio; // Run mean of the above, the single "how icy was this dataset" number (1 = no ice). The // bkg_estimate scalar's counterpart; written to /entry/MX/iceRingScoreMean. - std::optional ice_ring_score_mean; + std::optional ice_ring_ratio_mean; // Per-image sweep-quality code: 0 = the image falls in no flagged range, otherwise 1 + the // SweepQualityReason of the range it falls in (image_analysis/scale_merge/Merge.h). Written to diff --git a/common/JFJochReceiverPlots.cpp b/common/JFJochReceiverPlots.cpp index 2bdb0ec28..92c0294c9 100644 --- a/common/JFJochReceiverPlots.cpp +++ b/common/JFJochReceiverPlots.cpp @@ -62,7 +62,7 @@ void JFJochReceiverPlots::Setup(const DiffractionExperiment &experiment, const A } bkg_estimate.Clear(r); spindle_blind_fraction.Clear(r); - ice_ring_score.Clear(r); + ice_ring_ratio.Clear(r); protein_score.Clear(r); ice_score.Clear(r); spot_count.Clear(r); @@ -132,7 +132,7 @@ 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); spindle_blind_fraction.AddElement(msg.number, msg.spindle_blind_fraction); - ice_ring_score.AddElement(msg.number, msg.ice_ring_score); + ice_ring_ratio.AddElement(msg.number, msg.ice_ring_ratio); protein_score.AddElement(msg.number, msg.protein_score); ice_score.AddElement(msg.number, msg.ice_score); resolution_estimate.AddElement(msg.number, msg.resolution_estimate); @@ -286,8 +286,8 @@ MultiLinePlot JFJochReceiverPlots::GetPlots(const PlotRequest &request) { case PlotType::SpindleBlindFraction: ret = spindle_blind_fraction.GetMeanPlot(nbins, start, incr, request.fill_value); break; - case PlotType::IceRingScore: - ret = ice_ring_score.GetMeanPlot(nbins, start, incr, request.fill_value); + case PlotType::IceRingRatio: + ret = ice_ring_ratio.GetMeanPlot(nbins, start, incr, request.fill_value); break; case PlotType::ProteinScore: ret = protein_score.GetMeanPlot(nbins, start, incr, request.fill_value); @@ -520,16 +520,16 @@ std::optional JFJochReceiverPlots::GetIceRingSpotRatio() const { return ring / control; } -std::optional JFJochReceiverPlots::GetIceRingScore() const { - auto tmp = ice_ring_score.Mean(); +std::optional JFJochReceiverPlots::GetIceRingRatio() const { + auto tmp = ice_ring_ratio.Mean(); if (std::isfinite(tmp)) return tmp; else return {}; } -std::vector JFJochReceiverPlots::GetIceRingScoreArray() const { - return ice_ring_score.ExportArray(); +std::vector JFJochReceiverPlots::GetIceRingRatioArray() const { + return ice_ring_ratio.ExportArray(); } std::optional JFJochReceiverPlots::GetProteinScore() const { @@ -623,8 +623,8 @@ void JFJochReceiverPlots::GetPlotRaw(std::vector &v, PlotType type, const case PlotType::SpindleBlindFraction: v = spindle_blind_fraction.ExportArray(); break; - case PlotType::IceRingScore: - v = ice_ring_score.ExportArray(); + case PlotType::IceRingRatio: + v = ice_ring_ratio.ExportArray(); break; case PlotType::ProteinScore: v = protein_score.ExportArray(); diff --git a/common/JFJochReceiverPlots.h b/common/JFJochReceiverPlots.h index 05a6270eb..13384cf76 100644 --- a/common/JFJochReceiverPlots.h +++ b/common/JFJochReceiverPlots.h @@ -43,7 +43,7 @@ class JFJochReceiverPlots { StatusVector bkg_estimate; StatusVector spindle_blind_fraction; - StatusVector ice_ring_score; + StatusVector ice_ring_ratio; StatusVector protein_score; StatusVector ice_score; StatusVector spot_count; @@ -127,7 +127,7 @@ public: std::optional GetBkgEstimate() const; // Run mean of the per-image spindle severity, over the frames that had one; nothing when none did. std::optional GetSpindleBlindFraction() const; - std::optional GetIceRingScore() const; + std::optional GetIceRingRatio() const; // Pooled over the run: spots on the hexagonal rings over the same q width of ice-free control // flanks. 1 = spots spread evenly, > 1 = they pile up on the rings (textured ice). [[nodiscard]] std::optional GetIceRingSpotRatio() const; @@ -135,7 +135,7 @@ public: // not a mean - a few blank or badly-shot frames pull a mean a long way and leave the median where // the run actually is. [[nodiscard]] std::optional GetResolutionEstimate() const; - std::vector GetIceRingScoreArray() const; + std::vector GetIceRingRatioArray() 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; diff --git a/common/Plot.h b/common/Plot.h index 72a22277b..c38e01800 100644 --- a/common/Plot.h +++ b/common/Plot.h @@ -15,7 +15,7 @@ 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, IceRingRatio, ProteinScore, IceScore }; diff --git a/common/ScanResult.h b/common/ScanResult.h index d88a114ec..02ce9af4b 100644 --- a/common/ScanResult.h +++ b/common/ScanResult.h @@ -44,7 +44,7 @@ struct ScanResultElem { std::optional integrated_reflections; std::optional image_scale_factor; std::optional image_scale_cc; - std::optional ice_ring_score; + std::optional ice_ring_ratio; std::optional protein_score; std::optional ice_score; }; diff --git a/common/ScanResultGenerator.cpp b/common/ScanResultGenerator.cpp index 7c6ffc420..c510e2f1a 100644 --- a/common/ScanResultGenerator.cpp +++ b/common/ScanResultGenerator.cpp @@ -64,7 +64,7 @@ void ScanResultGenerator::Add(const DataMessage &message) { v[image_number].integrated_reflections = message.integrated_reflections; 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].ice_ring_ratio = message.ice_ring_ratio; v[image_number].protein_score = message.protein_score; v[image_number].ice_score = message.ice_score; if (message.lattice_type) @@ -119,7 +119,7 @@ void ScanResultGenerator::FillEndMessage(EndMessage &message) const { message.error_pixel_count.resize(n); message.image_scale_factor.resize(n, NAN); message.image_scale_cc.resize(n, NAN); - message.ice_ring_score.resize(n, NAN); + message.ice_ring_ratio.resize(n, NAN); message.v_protein_score.resize(n, NAN); message.v_ice_score.resize(n, NAN); message.integrated_reflections.resize(n); @@ -154,7 +154,7 @@ void ScanResultGenerator::FillEndMessage(EndMessage &message) const { message.error_pixel_count[number] = static_cast(value_or_zero(e.err_pixels)); 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.ice_ring_ratio[number] = e.ice_ring_ratio.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)); diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 6d66d1dcf..062beb5be 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -34,6 +34,25 @@ across them and are the ones the rest of the documentation cites. twinning and translational pseudo-symmetry, the resolution cutoff, diffraction anisotropy, and model-based validation. +## Naming: score, ratio, count + +Per-image quantities are named for the *kind* of number they are, because several of them describe +the same physics on opposite conventions and the name is the only thing that says which: + +| Suffix | Range | What 1 means | Examples | +|---|---|---|---| +| `*_score` | bounded, 0 to 1, saturating | **certainly present** | `protein_score`, `ice_score` | +| `*_ratio` | unbounded above | **none** | `ice_ring_ratio`, the ice-band spot ratio | +| `*_count` | integer | one occurrence | `spot_count`, `spot_count_ice_rings`, `spot_count_ice_control` | + +So a score is bounded and 1 is certainty; a ratio is unbounded and 1 is nothing. `ice_score` and +`ice_ring_ratio` both describe ice and their extremes are opposite ends of the scale. + +The C++ identifiers follow this. The **stored and wire names do not yet**: the ratio is still +written as `/entry/MX/iceRingScore`, sent as the CBOR key `ice_ring_score`, served as the +`ice_ring_score` plot type, and gated by `--ice-min-score`. Renaming those changes stored files, the +stream format, the REST API and a CLI flag, so it is a separate decision. + ## References The methods draw on, and in places reimplement, solutions from: diff --git a/frame_serialize/CBORStream2Deserializer.cpp b/frame_serialize/CBORStream2Deserializer.cpp index 96a9a989c..b37f54281 100644 --- a/frame_serialize/CBORStream2Deserializer.cpp +++ b/frame_serialize/CBORStream2Deserializer.cpp @@ -823,7 +823,7 @@ namespace { else if (key == "bkg_estimate") message.bkg_estimate = GetCBORFloat(value); else if (key == "ice_ring_score") - message.ice_ring_score = GetCBORFloat(value); + message.ice_ring_ratio = GetCBORFloat(value); else if (key == "protein_score") message.protein_score = GetCBORFloat(value); else if (key == "ice_score") @@ -1518,11 +1518,11 @@ namespace { else if (key == "v_spindle_blind_fraction") GetCBORFloatArray(value, message.v_spindle_blind_fraction); else if (key == "ice_ring_score") - GetCBORFloatArray(value, message.ice_ring_score); + GetCBORFloatArray(value, message.ice_ring_ratio); else if (key == "spot_count_ice_control") GetCBORFloatArray(value, message.spot_count_ice_control); else if (key == "ice_ring_score_mean") - message.ice_ring_score_mean = GetCBORFloat(value); + message.ice_ring_ratio_mean = GetCBORFloat(value); else if (key == "profile_radius") GetCBORFloatArray(value, message.profile_radius); else if (key == "mosaicity") diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index c27e24070..fb7c5f2ff 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -829,8 +829,8 @@ void CBORStream2Serializer::SerializeSequenceEnd(const EndMessage& message) { CBOR_ENC(mapEncoder, "protein_score", message.protein_score); CBOR_ENC(mapEncoder, "ice_score", message.ice_score); CBOR_ENC(mapEncoder, "v_spindle_blind_fraction", message.v_spindle_blind_fraction); - 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, "ice_ring_score", message.ice_ring_ratio); + CBOR_ENC(mapEncoder, "ice_ring_score_mean", message.ice_ring_ratio_mean); CBOR_ENC(mapEncoder, "spot_count_ice_control", message.spot_count_ice_control); CBOR_ENC(mapEncoder, "profile_radius", message.profile_radius); CBOR_ENC(mapEncoder, "mosaicity", message.mosaicity); @@ -926,7 +926,7 @@ void CBORStream2Serializer::SerializeImageInternal(CborEncoder &mapEncoder, cons CBOR_ENC(mapEncoder, "packets_expected", message.packets_expected); 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, "ice_ring_score", message.ice_ring_ratio); CBOR_ENC(mapEncoder, "protein_score", message.protein_score); CBOR_ENC(mapEncoder, "ice_score", message.ice_score); CBOR_ENC(mapEncoder, "spindle_blind_fraction", message.spindle_blind_fraction); diff --git a/image_analysis/IceScore.cpp b/image_analysis/IceScore.cpp index 3e4aca46d..6257e4568 100644 --- a/image_analysis/IceScore.cpp +++ b/image_analysis/IceScore.cpp @@ -26,18 +26,17 @@ namespace { 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; + // How far either side of a predicted band the radial channel looks for the bin that carries it. + // This is NOT a second band half-width competing with the spot channel's half_width_q: it is a + // tolerance for a mis-set beam centre, which smears and splits a ring, and it is set a priori from + // the measured centre error (+-0.6 % in d) rather than fitted. There is one band width in this + // file, and it is the pipeline's ice_ring_width_Q_recipA, passed in. + constexpr float CENTRE_SMEAR_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; @@ -144,20 +143,20 @@ namespace { 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) { + // Live pixels within +-w of q. + double AreaAt(const std::vector &count_q, float low_q, float dq, float q, float w) { 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)); + const int lo = static_cast(std::ceil((q - w - low_q) / dq - 0.5f)); + const int hi = static_cast(std::floor((q + w - 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); + // Spots within +-w of q, from the sorted spot q list. + double CountAt(const std::vector &sorted_q, float q, float w) { + const auto lo = std::lower_bound(sorted_q.begin(), sorted_q.end(), q - w); + const auto hi = std::upper_bound(sorted_q.begin(), sorted_q.end(), q + w); return static_cast(std::distance(lo, hi)); } } @@ -186,7 +185,7 @@ float IceScoreRadial(const std::vector &profile, const std::vector 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))); + const int bh = std::max(1, static_cast(std::lround(CENTRE_SMEAR_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. @@ -260,11 +259,11 @@ float IceScoreRadial(const std::vector &profile, const std::vector } float IceScoreSpots(const std::vector &spots, const std::vector &profile_count, - int32_t q_bins, const AzimuthalIntegrationSettings &settings) { + int32_t q_bins, const AzimuthalIntegrationSettings &settings, float half_width_q) { 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()) + if (nq < 2 || !(dq > 0.0f) || !(half_width_q > 0.0f) || spots.empty() || profile_count.empty()) return 0.0f; const std::vector count_q = FoldCountToQ(profile_count, nq); @@ -282,23 +281,23 @@ float IceScoreSpots(const std::vector &spots, const std::vector 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) { + for (float delta = 2.0f * half_width_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) + if (nearest < 2.0f * half_width_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); + const double area_hi = AreaAt(count_q, low_q, dq, qb + delta, half_width_q); + const double area_lo = AreaAt(count_q, low_q, dq, qb - delta, half_width_q); 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)); + control.push_back(0.5 * area_band * (CountAt(spot_q, qb + delta, half_width_q) / area_hi + + CountAt(spot_q, qb - delta, half_width_q) / area_lo)); } if (static_cast(control.size()) < SPOT_MIN_OFFSETS) continue; @@ -309,7 +308,7 @@ float IceScoreSpots(const std::vector &spots, const std::vector(control.size()); } @@ -327,7 +326,26 @@ float IceScoreSpots(const std::vector &spots, 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) { + const AzimuthalIntegrationSettings &settings, const std::vector &spots, + float half_width_q) { + // The larger of the two, deliberately, and this is a decided question rather than an open one. + // + // The radial channel cannot see ice that arrives as large crystallites: that ice diffracts as + // discrete spots and leaves the profile flat, so there is no radial evidence to read. Removing the + // spot channel would make the score blind to the regime the second channel exists for, by + // construction and not merely on this corpus. + // + // Measured over the corpus, the maximum takes the fraction of frames called icy in the protein + // classes labelled clean from 7% (radial alone) to 16%. That looks like a specificity cost and is + // not yet known to be one: the loops where the spot channel is the ONLY detector are precisely the + // ones whose ice status the labels do not settle, and several loops labelled clean - and three of + // four water loops - were independently found to carry a full hexagonal pattern. Until the corpus + // is relabelled from measurement rather than from directory names, the corpus cannot say whether + // those frames are false positives or correct detections nobody had recorded. + // + // If the rate does have to come down afterwards, the lever is to drop each hypothesis's worst band + // before scoring it (measured on the prototype: clean protein 20.5 -> 12.2%, ice 88.8 -> ~83%), not + // to remove a channel. return std::max(IceScoreRadial(profile, profile_std, profile_count, q_bins, settings), - IceScoreSpots(spots, profile_count, q_bins, settings)); + IceScoreSpots(spots, profile_count, q_bins, settings, half_width_q)); } diff --git a/image_analysis/IceScore.h b/image_analysis/IceScore.h index 31c7f53cd..a3cdc522a 100644 --- a/image_analysis/IceScore.h +++ b/image_analysis/IceScore.h @@ -31,6 +31,10 @@ // null, an amplitude and a band-count concordance, and the SMALLER is taken: a single elevated bin then // fails, because real ice shows a whole pattern. // +// This channel takes no band width: it reads the single bin the geometry predicts, searched over a small +// tolerance for the beam-centre error (see CENTRE_SMEAR_Q in the .cpp). The one band width in this file +// is the spot channel's, and it is the pipeline's. +// // profile / profile_std / profile_count are q_bins long, or q_bins x azimuthal bins, in which case the // first two are averaged over azimuth and the third summed. The scale of an excess is the bin MEAN's own // error, std / sqrt(count) - the profile is a mean of many pixels, so its plain standard deviation is @@ -52,10 +56,14 @@ float IceScoreRadial(const std::vector &profile, const std::vector // // profile_count is the azimuthal integration's live pixel count per bin, q_bins long or q_bins x // azimuthal bins; it is what makes the offsets comparable where the detector edge cuts a radius short. +// half_width_q is the band half-width, the same one the spot finder marks ice rings with +// (SpotFindingSettings::ice_ring_width_Q_recipA) - one width, measured over the corpus, not two. float IceScoreSpots(const std::vector &spots, const std::vector &profile_count, - int32_t q_bins, const AzimuthalIntegrationSettings &settings); + int32_t q_bins, const AzimuthalIntegrationSettings &settings, float half_width_q); -// The score itself: whichever channel sees more. +// The score itself: whichever channel sees more. See the .cpp for why the maximum and not something +// gentler, and for the lever to reach for if the rate ever has to come down. float IceScore(const std::vector &profile, const std::vector &profile_std, const std::vector &profile_count, int32_t q_bins, - const AzimuthalIntegrationSettings &settings, const std::vector &spots); + const AzimuthalIntegrationSettings &settings, const std::vector &spots, + float half_width_q); diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 8f21bb823..a0fefefd9 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -179,10 +179,10 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, const std::vector &ring_bkg = adaptiveSpotFinder->GetRingBackground(); if (ring_bkg.empty()) return; - output.ice_ring_score = AzimuthalIntegrationProfile::IceRingScore( + output.ice_ring_ratio = AzimuthalIntegrationProfile::IceRingRatio( ring_bkg, integration.GetQBinCount(), integration.Settings(), spot_finding_settings.ice_ring_width_Q_recipA); - bragg_engine->BackgroundRadial(*output.ice_ring_score + bragg_engine->BackgroundRadial(*output.ice_ring_ratio >= experiment.GetScalingSettings().GetIceMinScore()); }; @@ -289,7 +289,7 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, const std::vector &ring_bkg = adaptiveSpotFinder->GetRingBackground(); const bool have_ring_bkg = spot_finding_settings.enable && spot_finding_settings.adaptive_threshold && !ring_bkg.empty(); - output.ice_ring_score = AzimuthalIntegrationProfile::IceRingScore( + output.ice_ring_ratio = AzimuthalIntegrationProfile::IceRingRatio( have_ring_bkg ? ring_bkg : profile.GetResult1D(), integration.GetQBinCount(), integration.Settings(), spot_finding_settings.ice_ring_width_Q_recipA); @@ -298,7 +298,8 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, // 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); + integration.Settings(), output.spots, + spot_finding_settings.ice_ring_width_Q_recipA); } ImageSpotFinder &MXAnalysisWithoutFPGA::FixedThresholdFinder() { diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index 0f47ca16b..30653eeed 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -575,7 +575,7 @@ 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->spindle_blind_fraction = master_file->ReadOptVector("/entry/MX/spindleBlindFraction"); - dataset->ice_ring_score = master_file->ReadOptVector("/entry/MX/iceRingScore"); + dataset->ice_ring_ratio = 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"); @@ -726,7 +726,7 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen data_file, "/entry/MX/spindleBlindFraction", number_of_images, fimages); - ReadVector(dataset->ice_ring_score, + ReadVector(dataset->ice_ring_ratio, data_file, "/entry/MX/iceRingScore", number_of_images, fimages); @@ -1380,8 +1380,8 @@ void HDF5MetadataSource::FillPerImage(DataMessage &message, int64_t requested_im if (dataset->spindle_blind_fraction.size() > image_number && std::isfinite(dataset->spindle_blind_fraction[image_number])) message.spindle_blind_fraction = dataset->spindle_blind_fraction[image_number]; - if (dataset->ice_ring_score.size() > image_number) - message.ice_ring_score = dataset->ice_ring_score[image_number]; + if (dataset->ice_ring_ratio.size() > image_number) + message.ice_ring_ratio = dataset->ice_ring_ratio[image_number]; // Same convention as spindle_blind_fraction above: a frame that produced no score is stored as // NaN, and 0 is a real answer ("nothing detected here"), so a missing one has to come back absent // rather than as a NaN that a threshold would then compare against. diff --git a/reader/JFJochReaderDataset.h b/reader/JFJochReaderDataset.h index 1f1525a0f..8bf91072c 100644 --- a/reader/JFJochReaderDataset.h +++ b/reader/JFJochReaderDataset.h @@ -58,7 +58,7 @@ struct JFJochReaderDataset { std::vector bkg_estimate; // NaN entries are frames with no value (CANNOT SAY), not zeros. std::vector spindle_blind_fraction; - std::vector ice_ring_score; + std::vector ice_ring_ratio; // 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 diff --git a/receiver/JFJochReceiver.cpp b/receiver/JFJochReceiver.cpp index 18c16515b..bb00230c0 100644 --- a/receiver/JFJochReceiver.cpp +++ b/receiver/JFJochReceiver.cpp @@ -163,7 +163,7 @@ void JFJochReceiver::SendEndMessage() { message.bkg_estimate = plots.GetBkgEstimate(); message.spindle_blind_fraction = plots.GetSpindleBlindFraction(); - message.ice_ring_score_mean = plots.GetIceRingScore(); + message.ice_ring_ratio_mean = plots.GetIceRingRatio(); message.protein_score = plots.GetProteinScore(); message.ice_score = plots.GetIceScore(); message.indexing_rate = plots.GetIndexingRate(); diff --git a/receiver/JFJochReceiverFPGA.cpp b/receiver/JFJochReceiverFPGA.cpp index f5e4013d7..72fc23be8 100644 --- a/receiver/JFJochReceiverFPGA.cpp +++ b/receiver/JFJochReceiverFPGA.cpp @@ -435,7 +435,7 @@ void JFJochReceiverFPGA::FrameTransformationThread(uint32_t threadid) { if (force_cpu_azint) message.az_int_profile_std = az_int_profile_image.GetStd(); message.bkg_estimate = az_int_profile_image.GetBkgEstimate(experiment.GetAzimuthalIntegrationSettings()); - message.ice_ring_score = az_int_profile_image.GetIceRingScore( + message.ice_ring_ratio = az_int_profile_image.GetIceRingRatio( 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 @@ -443,7 +443,8 @@ void JFJochReceiverFPGA::FrameTransformationThread(uint32_t threadid) { 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); + experiment.GetAzimuthalIntegrationSettings(), message.spots, + spot_finding_settings.ice_ring_width_Q_recipA); scan_result.Add(message); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index c44d7e2c2..3fd6e078f 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -3440,13 +3440,14 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b msg.az_int_profile_count = profile.GetPixelCount(); msg.az_int_profile_std = profile.GetStd(); msg.bkg_estimate = profile.GetBkgEstimate(mapping.Settings()); - msg.ice_ring_score = profile.GetIceRingScore(mapping.Settings(), + msg.ice_ring_ratio = profile.GetIceRingRatio(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); + mapping.Settings(), msg.spots, + config_.spot_finding.ice_ring_width_Q_recipA); msg.run_number = experiment_.GetRunNumber(); msg.run_name = experiment_.GetRunName(); @@ -3651,8 +3652,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b end_msg.bkg_estimate = plots.GetBkgEstimate(); end_msg.spindle_blind_fraction = plots.GetSpindleBlindFraction(); result.spot_resolution_estimate_A = plots.GetResolutionEstimate(); - end_msg.ice_ring_score = plots.GetIceRingScoreArray(); - end_msg.ice_ring_score_mean = plots.GetIceRingScore(); + end_msg.ice_ring_ratio = plots.GetIceRingRatioArray(); + end_msg.ice_ring_ratio_mean = plots.GetIceRingRatio(); end_msg.v_protein_score = plots.GetProteinScoreArray(); end_msg.v_ice_score = plots.GetIceScoreArray(); end_msg.protein_score = plots.GetProteinScore(); @@ -3720,24 +3721,24 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // ring over the smooth radial background, measured on the peak-excluded profile so Bragg peaks // cannot pose as ice). The score is logged either way, so a run that gates off still reports why. // Ice is looked for on TWO channels, because it arrives in two forms. Fine-grained ice is a - // smooth powder ring and shows up in the radial profile (ice_ring_score). Ice in large + // smooth powder ring and shows up in the radial profile (ice_ring_ratio). Ice in large // crystallites arrives as discrete reflections instead, leaves the profile flat, and is only // visible as a pile-up of found SPOTS on the ring positions - measured against the same q width // of ice-free flanks beside them, so the fall-off of spot density with resolution cancels. On // the rotation battery the two channels barely overlap: the crystals with smooth ice read // 2.1-2.4 / ~1.0, the ones whose spots are mostly ice read ~1.1 / 3.8-14.6, and a clean crystal // reads 1.04 on both. Either channel alone is blind to half the cases. - const auto ice_score = plots.GetIceRingScore(); + const auto ice_ring_ratio = plots.GetIceRingRatio(); const auto ice_spot_ratio = plots.GetIceRingSpotRatio(); const float ice_min_score = experiment_.GetScalingSettings().GetIceMinScore(); const float ice_min_spot_ratio = experiment_.GetScalingSettings().GetIceMinSpotRatio(); - const bool ice_present = (!ice_score.has_value() || *ice_score >= ice_min_score) + const bool ice_present = (!ice_ring_ratio.has_value() || *ice_ring_ratio >= ice_min_score) || (ice_min_spot_ratio > 0.0f && ice_spot_ratio.has_value() && *ice_spot_ratio >= ice_min_spot_ratio); if (experiment_.IsDetectIceRings() && !ice_present) { logger.Info("Ice-ring handling: measured ice score {:.2f} < {:.2f} and spot ratio {:.2f} < " "{:.2f}, no ice detected - ice-ring handling skipped entirely", - *ice_score, ice_min_score, ice_spot_ratio.value_or(NAN), ice_min_spot_ratio); + *ice_ring_ratio, ice_min_score, ice_spot_ratio.value_or(NAN), ice_min_spot_ratio); } else if (experiment_.IsDetectIceRings()) { const float ice_width = config_.spot_finding.ice_ring_width_Q_recipA; size_t total = 0, flagged = 0; @@ -3752,7 +3753,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b logger.Info("Ice-ring handling: ice score {:.2f} (gate {:.2f}), spot ratio {:.2f} (gate " "{:.2f}); flagged {} of {} reflections on ice rings (half-width {:.3f} A^-1); " "excluded from scaling, kept for merging", - ice_score.value_or(NAN), ice_min_score, ice_spot_ratio.value_or(NAN), + ice_ring_ratio.value_or(NAN), ice_min_score, ice_spot_ratio.value_or(NAN), ice_min_spot_ratio, flagged, total, ice_width); } diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index c08b77e2e..64d2374e6 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1630,11 +1630,12 @@ static int RunRugnux(int argc, char **argv) { // otherwise --scale re-scales a dataset the writing run had scaled without those reflections, // and the per-image scales come out of a different fit than the ones in the file. const float ice_width = SpotFindingSettings().ice_ring_width_Q_recipA; - // ...and gated the same way, on the per-image ice score the writing run stored in the file, so + // ...and gated the same way, on the per-image ice ring ratio the writing run stored in the file, + // so // --scale reaches the same verdict on the same data as the pipeline that produced it. double ice_sum = 0.0; size_t ice_n = 0; - for (const float s : dataset->ice_ring_score) + for (const float s : dataset->ice_ring_ratio) if (std::isfinite(s)) { ice_sum += s; ++ice_n; diff --git a/tests/AzimuthalIntegrationTest.cpp b/tests/AzimuthalIntegrationTest.cpp index 98dc3e8de..4b701189e 100644 --- a/tests/AzimuthalIntegrationTest.cpp +++ b/tests/AzimuthalIntegrationTest.cpp @@ -547,10 +547,10 @@ TEST_CASE("AzimuthalIntegrationMapping_Threading_ConvertedGeometry_18Modules", " CheckAzimuthalIntegrationMappingThreadingExact(x); } -// The ice score on an explicit profile: flat means no ice (1.0), and a bump planted on a hexagonal +// The ice ring ratio on an explicit profile: flat means no ice (1.0), and a bump planted on a hexagonal // ring is reported at its own height over the background. Uses the static entry point, which is what // feeds the score the peak-excluded per-ring background instead of the plain profile. -TEST_CASE("AzimuthalIntegrationProfile_IceRingScore","[AzimuthalIntegration]") { +TEST_CASE("AzimuthalIntegrationProfile_IceRingRatio","[AzimuthalIntegration]") { AzimuthalIntegrationSettings settings; settings.QSpacing_recipA(0.01f).QRange_recipA(0.1f, 4.5f); @@ -563,21 +563,21 @@ TEST_CASE("AzimuthalIntegrationProfile_IceRingScore","[AzimuthalIntegration]") { }; std::vector profile(q_bins, 100.0f); - CHECK(AzimuthalIntegrationProfile::IceRingScore(profile, q_bins, settings, 0.03f) == Catch::Approx(1.0f)); + CHECK(AzimuthalIntegrationProfile::IceRingRatio(profile, q_bins, settings, 0.03f) == Catch::Approx(1.0f)); // A bump at the 2.249 A hexagonal ring, on a background the running median still reads as 100. profile[bin_of(2.249f)] = 250.0f; - CHECK(AzimuthalIntegrationProfile::IceRingScore(profile, q_bins, settings, 0.03f) == Catch::Approx(2.5f)); + CHECK(AzimuthalIntegrationProfile::IceRingRatio(profile, q_bins, settings, 0.03f) == Catch::Approx(2.5f)); // A bump of the same size well away from every ring is not ice and must not be reported. std::vector off_ring(q_bins, 100.0f); off_ring[bin_of(2.500f)] = 250.0f; - CHECK(AzimuthalIntegrationProfile::IceRingScore(off_ring, q_bins, settings, 0.03f) == Catch::Approx(1.0f)); + CHECK(AzimuthalIntegrationProfile::IceRingRatio(off_ring, q_bins, settings, 0.03f) == Catch::Approx(1.0f)); } // A profile given per (q, azimuth) bin is averaged over azimuth first, so a ring seen in every sector // scores the same as the equivalent 1-D profile. -TEST_CASE("AzimuthalIntegrationProfile_IceRingScore_Azimuthal","[AzimuthalIntegration]") { +TEST_CASE("AzimuthalIntegrationProfile_IceRingRatio_Azimuthal","[AzimuthalIntegration]") { AzimuthalIntegrationSettings settings; settings.QSpacing_recipA(0.01f).QRange_recipA(0.1f, 4.5f).AzimuthalBinCount(4); @@ -591,6 +591,6 @@ TEST_CASE("AzimuthalIntegrationProfile_IceRingScore_Azimuthal","[AzimuthalIntegr for (int az = 0; az < 4; az++) sectors[static_cast(az) * q_bins + ring] = 250.0f; - CHECK(AzimuthalIntegrationProfile::IceRingScore(sectors, q_bins, settings, 0.03f) - == Catch::Approx(AzimuthalIntegrationProfile::IceRingScore(flat, q_bins, settings, 0.03f))); + CHECK(AzimuthalIntegrationProfile::IceRingRatio(sectors, q_bins, settings, 0.03f) + == Catch::Approx(AzimuthalIntegrationProfile::IceRingRatio(flat, q_bins, settings, 0.03f))); } diff --git a/tests/DetectionScoreTest.cpp b/tests/DetectionScoreTest.cpp index 7f8b79d93..64830cf7a 100644 --- a/tests/DetectionScoreTest.cpp +++ b/tests/DetectionScoreTest.cpp @@ -9,6 +9,8 @@ namespace { constexpr float TWO_PI = 6.283185307f; + // The pipeline's own ice band half-width (SpotFindingSettings::ice_ring_width_Q_recipA). + constexpr float ICE_W = 0.03f; SpotToSave spot(float d_A, float intensity) { return SpotToSave{.intensity = intensity, .d_A = d_A}; @@ -158,23 +160,23 @@ TEST_CASE("IceScoreSpots_ExcessOnTheIceRadii") { 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); + CHECK(IceScoreSpots(even, count, q_bins, settings, ICE_W) < 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); + CHECK(IceScoreSpots(with_ice, count, q_bins, settings, ICE_W) > 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(two, count, q_bins, settings, ICE_W) < 0.5f); - CHECK(IceScoreSpots({}, count, q_bins, settings) == 0.0f); + CHECK(IceScoreSpots({}, count, q_bins, settings, ICE_W) == 0.0f); } TEST_CASE("IceScore_TakesTheStrongerChannel") { @@ -190,7 +192,7 @@ TEST_CASE("IceScore_TakesTheStrongerChannel") { if (b >= 0 && b < q_bins) powder[b] = 130.0f; } - CHECK(IceScore(powder, sigma, count, q_bins, settings, {}) > 0.5f); + CHECK(IceScore(powder, sigma, count, q_bins, settings, {}, ICE_W) > 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); @@ -201,5 +203,5 @@ TEST_CASE("IceScore_TakesTheStrongerChannel") { 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); + CHECK(IceScore(flat, sigma, count, q_bins, settings, textured, ICE_W) > 0.5f); } diff --git a/viewer/JFJochHttpReader.cpp b/viewer/JFJochHttpReader.cpp index b4c6e0f98..0a9db1022 100644 --- a/viewer/JFJochHttpReader.cpp +++ b/viewer/JFJochHttpReader.cpp @@ -268,7 +268,7 @@ std::shared_ptr JFJochHttpReader::UpdateDataset_i() { dataset->bkg_estimate = GetPlot_i("bkg_estimate"); dataset->spindle_blind_fraction = GetPlot_i("spindle_blind_fraction"); - dataset->ice_ring_score = GetPlot_i("ice_ring_score"); + dataset->ice_ring_ratio = GetPlot_i("ice_ring_score"); dataset->spot_count = GetPlot_i("spot_count"); dataset->spot_count_ice_rings = GetPlot_i("spot_count_ice"); dataset->spot_count_low_res = GetPlot_i("spot_count_low_res"); diff --git a/viewer/JFJochViewerDatasetInfo.cpp b/viewer/JFJochViewerDatasetInfo.cpp index 237836734..12abb776c 100644 --- a/viewer/JFJochViewerDatasetInfo.cpp +++ b/viewer/JFJochViewerDatasetInfo.cpp @@ -237,7 +237,7 @@ std::vector JFJochViewerDatasetInfo::ExtractMetric(const JFJochReaderData else if (val == 11) data = ds.image_scale_cc; else if (val == 12) data = ds.integrated_reflections; else if (val == 13) data = ds.indexing_lattice_count; - else if (val == 14) data = ds.ice_ring_score; + else if (val == 14) data = ds.ice_ring_ratio; else if (val == 15) data = ds.spindle_blind_fraction; else if (val >= 100) { const int roi_index = (val - 100) / 4; diff --git a/writer/HDF5DataFilePluginMX.cpp b/writer/HDF5DataFilePluginMX.cpp index fc66225eb..71696afe8 100644 --- a/writer/HDF5DataFilePluginMX.cpp +++ b/writer/HDF5DataFilePluginMX.cpp @@ -53,7 +53,7 @@ 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); + ice_ring_ratio.reserve(images_per_file); protein_score.reserve(images_per_file); ice_score.reserve(images_per_file); spindle_blind_fraction.reserve(images_per_file); @@ -102,8 +102,8 @@ void HDF5DataFilePluginMX::OpenFile(HDF5File &data_file, const DataMessage &msg, void HDF5DataFilePluginMX::Write(const DataMessage &msg, uint64_t image_number) { if (msg.bkg_estimate.has_value()) 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.ice_ring_ratio.has_value()) + ice_ring_ratio[image_number] = msg.ice_ring_ratio.value(); if (msg.protein_score.has_value()) protein_score[image_number] = msg.protein_score.value(); if (msg.ice_score.has_value()) @@ -259,8 +259,8 @@ void HDF5DataFilePluginMX::WriteFinal(HDF5File &data_file) { if (!bkg_estimate.empty()) 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 (!ice_ring_ratio.empty()) + data_file.SaveVector("/entry/MX/iceRingScore", ice_ring_ratio.vec()); if (!protein_score.empty()) data_file.SaveVector("/entry/MX/proteinScore", protein_score.vec()); if (!ice_score.empty()) diff --git a/writer/HDF5DataFilePluginMX.h b/writer/HDF5DataFilePluginMX.h index 61822e1f2..1463aa479 100644 --- a/writer/HDF5DataFilePluginMX.h +++ b/writer/HDF5DataFilePluginMX.h @@ -47,7 +47,7 @@ class HDF5DataFilePluginMX : public HDF5DataFilePlugin { // bkg_estimate AutoIncrVector bkg_estimate{NAN}; - AutoIncrVector ice_ring_score{NAN}; + AutoIncrVector ice_ring_ratio{NAN}; AutoIncrVector protein_score{NAN}; AutoIncrVector ice_score{NAN}; // NaN = the frame had no value (CANNOT SAY), which is not a zero - see SpindleTrigger. diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp index 3dcfb48bb..6fa941a91 100644 --- a/writer/HDF5NXmx.cpp +++ b/writer/HDF5NXmx.cpp @@ -1186,8 +1186,8 @@ void NXmx::Finalize(const EndMessage &end) { SaveScalar(*hdf5_file, "/entry/MX/spindleLostUniqueFraction", end.spindle_lost_unique_fraction.value()); } - if (end.ice_ring_score_mean) { - SaveScalar(*hdf5_file, "/entry/MX/iceRingScoreMean", end.ice_ring_score_mean.value()); + if (end.ice_ring_ratio_mean) { + SaveScalar(*hdf5_file, "/entry/MX/iceRingScoreMean", end.ice_ring_ratio_mean.value()); } if (end.protein_score) { SaveScalar(*hdf5_file, "/entry/MX/proteinScoreMean", end.protein_score.value()); @@ -1282,7 +1282,7 @@ void NXmx::EndResultVectors(const EndMessage &end) { 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/spindleBlindFraction", end.v_spindle_blind_fraction); - SaveVectorIfMissing(*hdf5_file, "/entry/MX/iceRingScore", end.ice_ring_score); + SaveVectorIfMissing(*hdf5_file, "/entry/MX/iceRingScore", end.ice_ring_ratio); SaveVectorIfMissing(*hdf5_file, "/entry/MX/profileRadius", end.profile_radius, "Angstrom^-1"); SaveVectorIfMissing(*hdf5_file, "/entry/MX/mosaicity", end.mosaicity, "deg"); SaveVectorIfMissing(*hdf5_file, "/entry/MX/bFactor", end.bFactor, "Angstrom^2");