diff --git a/broker/gen/model/Scan_result_images_inner.h b/broker/gen/model/Scan_result_images_inner.h index 9677d5e90..2a9beb4ed 100644 --- a/broker/gen/model/Scan_result_images_inner.h +++ b/broker/gen/model/Scan_result_images_inner.h @@ -125,7 +125,7 @@ public: bool spotsIceIsSet() const; void unsetSpots_ice(); /// - /// Strongest hexagonal-ice ring band/shoulder intensity ratio (1 = no ice) + /// Strongest hexagonal-ice ring intensity over the smooth radial background (1 = no ice) /// float getIceRingScore() const; void setIceRingScore(float const value); diff --git a/broker/jfjoch_api.yaml b/broker/jfjoch_api.yaml index c3f05238b..5d6bdafb6 100644 --- a/broker/jfjoch_api.yaml +++ b/broker/jfjoch_api.yaml @@ -1582,7 +1582,7 @@ components: ice_ring_score: type: number format: float - description: Strongest hexagonal-ice ring band/shoulder intensity ratio (1 = no ice) + description: Strongest hexagonal-ice ring intensity over the smooth radial background (1 = no ice) index: type: integer format: int64 diff --git a/common/AzimuthalIntegrationProfile.cpp b/common/AzimuthalIntegrationProfile.cpp index 91530dd6c..d4794e38a 100644 --- a/common/AzimuthalIntegrationProfile.cpp +++ b/common/AzimuthalIntegrationProfile.cpp @@ -176,23 +176,74 @@ float AzimuthalIntegrationProfile::GetBkgEstimate(const AzimuthalIntegrationSett float AzimuthalIntegrationProfile::GetIceRingScore(const AzimuthalIntegrationSettings &settings, float half_width_q) const { - // For each hexagonal-ice powder ring, the mean profile intensity in the ring band (+/- half_width in - // q = 2*pi/d) over a baseline from the two shoulders just outside it. Report the strongest ring's - // ratio (1 = no ice); rings whose band+shoulders fall off the measured q-range are skipped. + // Strongest hexagonal-ice ring's intensity relative to the background *under* it (1 = no ice). The + // background is a smooth whole-profile estimate: a running median of the NON-ice bins, interpolated to + // each ring position - not a couple of adjacent shoulder bins (the azint binning is coarser than the + // ring width, so a local shoulder is only ~1 bin and a narrow ratio is noisy and can double-count the + // ring's own edge). Clean profiles then sit at ~1 at every ring; ice makes the ring bin stand out. constexpr float two_pi = 6.283185307f; + const std::vector prof = GetResult1D(); + const int nq = static_cast(prof.size()); const float low_q = settings.GetLowQ_recipA(); - const float high_q = settings.GetHighQ_recipA(); + const float dq = settings.GetQSpacing_recipA(); + if (nq < 12 || !(dq > 0.0f)) + return 1.0f; + + auto q_of = [&](int i) { return low_q + (static_cast(i) + 0.5f) * dq; }; + auto on_ice = [&](float q) { + for (const float d : ICE_RING_RES_A) + if (std::fabs(q - two_pi / d) < 1.5f * half_width_q) + return true; + return false; + }; + + // Non-ice, finite, positive bins (ascending q) carry the background. + std::vector base; + for (int i = 0; i < nq; ++i) + if (std::isfinite(prof[i]) && prof[i] > 0.0f && !on_ice(q_of(i))) + base.push_back(i); + if (base.size() < 8) + return 1.0f; + + // Running median over the base bins => a smooth background robust to the ice peaks. + constexpr int K = 4; + std::vector base_bg(base.size()); + std::vector window; + for (int j = 0; j < static_cast(base.size()); ++j) { + const int lo = std::max(0, j - K); + const int hi = std::min(static_cast(base.size()), j + K + 1); + window.clear(); + for (int m = lo; m < hi; ++m) + window.push_back(prof[base[m]]); + std::sort(window.begin(), window.end()); + base_bg[j] = window[window.size() / 2]; + } + float score = 1.0f; - for (const float ice_d : ICE_RING_RES_A) { - const float q = two_pi / ice_d; - if (q - 2 * half_width_q < low_q || q + 2 * half_width_q > high_q) + for (const float d : ICE_RING_RES_A) { + const float qr = two_pi / d; + const int b = static_cast(std::lround((qr - low_q) / dq - 0.5f)); + if (b < 0 || b >= nq || !std::isfinite(prof[b]) || prof[b] <= 0.0f) continue; - const float ring = GetMeanValueOfBins(settings.QToBin(q - half_width_q), settings.QToBin(q + half_width_q)); - const float lo = GetMeanValueOfBins(settings.QToBin(q - 2 * half_width_q), settings.QToBin(q - half_width_q)); - const float hi = GetMeanValueOfBins(settings.QToBin(q + half_width_q), settings.QToBin(q + 2 * half_width_q)); - const float baseline = 0.5f * (lo + hi); - if (std::isfinite(baseline) && baseline > 0.0f && std::isfinite(ring)) - score = std::max(score, ring / baseline); + // Linear-interpolate the smooth background to the ring position. + float bg; + if (qr <= q_of(base.front())) + bg = base_bg.front(); + else if (qr >= q_of(base.back())) + bg = base_bg.back(); + else { + bg = NAN; + for (int j = 0; j + 1 < static_cast(base.size()); ++j) { + const float qa = q_of(base[j]), qb = q_of(base[j + 1]); + if (qa <= qr && qr <= qb) { + const float t = (qb > qa) ? (qr - qa) / (qb - qa) : 0.0f; + bg = base_bg[j] + t * (base_bg[j + 1] - base_bg[j]); + break; + } + } + } + if (std::isfinite(bg) && bg > 0.0f) + score = std::max(score, prof[b] / bg); } return score; } diff --git a/common/AzimuthalIntegrationProfile.h b/common/AzimuthalIntegrationProfile.h index 31e52df61..f7d998f37 100644 --- a/common/AzimuthalIntegrationProfile.h +++ b/common/AzimuthalIntegrationProfile.h @@ -44,8 +44,8 @@ public: float GetMeanValueOfBins(uint16_t min_bin, uint16_t max_bin) const; float GetBkgEstimate(const AzimuthalIntegrationSettings& settings) const; - // Single per-image ice indicator: the strongest hexagonal-ice ring's band/shoulder intensity ratio - // (1 = no ice, >1 = ice above the local background). Max over the rings in range; 1 if none. + // 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; 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 d8ced99d4..f186a65bb 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -119,7 +119,7 @@ struct DataMessage { std::vector az_int_profile_count; std::optional bkg_estimate; - std::optional ice_ring_score; // strongest ice-ring band/shoulder ratio from the azint profile (1 = none) + std::optional ice_ring_score; // strongest ice ring over the smooth radial background (1 = none) std::optional indexing_result; std::optional indexing_lattice; diff --git a/docs/CBOR.md b/docs/CBOR.md index 79bacf647..5f3fd2c55 100644 --- a/docs/CBOR.md +++ b/docs/CBOR.md @@ -215,7 +215,7 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | packets_expected | uint64 | Number of packets expected per image (in units of 2 kB) | | | | 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 band/shoulder intensity ratio from the azint profile (1 = no ice) | | | +| ice_ring_score | float | Strongest hexagonal-ice ring intensity over the smooth radial background (1 = no ice) | | | | 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 | | image_scale_factor | float | Scaling result: Image scale factor (g) | | X | @@ -306,7 +306,7 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | spot_count_indexed | Array(int32) | Per-image number of spots fitting indexing solution | | | 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 band/shoulder intensity ratio (1 = no ice) | | +| ice_ring_score | Array(float) | Per-image strongest ice-ring intensity over the smooth radial background (1 = no ice) | | | 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.md b/docs/CPU_DATA_ANALYSIS.md index df352a123..f7617b168 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -180,7 +180,7 @@ Special cases: Spot finding can be restricted to a resolution range $[d_\mathrm{high}, d_\mathrm{low}]$ by masking pixels outside the range. Optionally, pixels in identified ice-ring regions can be tagged so that subsequent indexing/refinement may include or exclude them (see §4 and §6). -A single per-image **ice-ring score** is derived from the azimuthally-integrated radial profile: for each hexagonal-ice powder ring (positions $d$ from Moreau *et al.*, Acta Cryst D77, 2021), the mean profile intensity in the ring band ($\pm$ ice-ring half-width in $q$) is divided by a baseline interpolated from the two shoulders just outside it, and the strongest ring's ratio is reported (1 = no ice, $>1$ = ice above background). It is stored per image (`ice_ring_score`, HDF5 `/entry/MX/iceRingScore`) as a monitoring quantity. Note this is distinct from the merge-time ice masking, which is data-driven from the per-ring merged CC1/2 rather than this background ratio. +A single per-image **ice-ring score** is derived from the azimuthally-integrated radial profile: for each hexagonal-ice powder ring (positions $d$ from Moreau *et al.*, Acta Cryst D77, 2021), the profile intensity at the ring is divided by a smooth background estimated from the *whole* profile — a running median of the non-ice bins, interpolated under each ring — and the strongest ring's ratio is reported (1 = no ice, $>1$ = ice above background). A whole-profile background is used rather than a couple of adjacent shoulder bins because the radial binning is typically coarser than the ring width (a local shoulder is only ~1 bin), so a narrow ratio is noisy and can double-count the ring's own edge. (A significance/z-score was considered but is uninformative here: with many photons any real ice ring is highly significant, so the discriminating quantity is the ice *magnitude*, i.e. this ratio.) It is stored per image (`ice_ring_score`, HDF5 `/entry/MX/iceRingScore`) as a monitoring quantity, distinct from the merge-time ice masking, which is data-driven from the per-ring merged CC1/2. A further optional safeguard removes isolated high-resolution “spur” spots by detecting large gaps in $1/d$ (or $q$) space and discarding spots beyond the gap. This is intended for macromolecular diffraction where edge-of-detector backgrounds can be extremely low. diff --git a/docs/HDF5.md b/docs/HDF5.md index c98f49159..4a9e01dd5 100644 --- a/docs/HDF5.md +++ b/docs/HDF5.md @@ -251,7 +251,7 @@ In legacy/VDS mode these live in the data files and are linked/virtual-stacked i | `resolutionEstimate` | Å | diffraction resolution estimate | | `integratedReflections` | | number of integrated reflections | | `bkgEstimate` | photons | mean background in the 3–5 Å resolution band | -| `iceRingScore` | ratio | strongest hexagonal-ice ring band/shoulder intensity ratio (1 = no ice) | +| `iceRingScore` | ratio | strongest hexagonal-ice ring intensity over the smooth radial background (1 = no ice) | | `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 | diff --git a/frontend/src/client/types.gen.ts b/frontend/src/client/types.gen.ts index 54f29f163..f5db882a4 100644 --- a/frontend/src/client/types.gen.ts +++ b/frontend/src/client/types.gen.ts @@ -1011,7 +1011,7 @@ export type scan_result = { */ spots_ice?: number; /** - * Strongest hexagonal-ice ring band/shoulder intensity ratio (1 = no ice) + * Strongest hexagonal-ice ring intensity over the smooth radial background (1 = no ice) */ ice_ring_score?: number; /**