diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 6b4c101f7..ea03b1ad7 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -237,7 +237,8 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, // plain profile carries the Bragg peaks and cannot support an absolute threshold. // // On stills the correction also takes its CURVE from that background, and auto applies it to - // every image that has one. A serial frame has a few dozen predicted reflections, so their own + // every image that has one (rugnux's pre-scan turns it off for a run whose empty apertures say + // the radial model does not hold there). A serial frame has a few dozen predicted reflections, so their own // ring pixels - where the curve comes from otherwise - sample the radial background sparsely, // while the finder's rings hold every pixel of the frame. And the correction is the measured // curvature itself: where the background has none it is zero, so there is nothing for a gate diff --git a/image_analysis/bragg_integration/BackgroundAperture.cpp b/image_analysis/bragg_integration/BackgroundAperture.cpp index a5b696756..14e89dc97 100644 --- a/image_analysis/bragg_integration/BackgroundAperture.cpp +++ b/image_analysis/bragg_integration/BackgroundAperture.cpp @@ -6,6 +6,8 @@ #include #include +#include "../../common/JFJochMath.h" + std::vector background_aperture::CandidateR3(float r2, float max_r3) { std::vector out; for (float r3 = r2 + 2.0f; r3 <= max_r3 + 1e-3f; r3 += 1.0f) @@ -16,10 +18,12 @@ std::vector background_aperture::CandidateR3(float r2, float max_r3) { void background_aperture::MeasureFrame(ImagePreprocessorBuffer &image, int width, int height, const std::vector &spots, float r1, float r2, const std::vector &r3, float clip_nsigma, float gain, - Tally &tally) { + const RingBackground &rings, const DiffractionGeometry &geometry, + const DetectorModules &modules, Tally &tally) { if (r3.empty()) return; tally.sum_d2.resize(r3.size(), 0.0); + tally.sum_d2_radial.resize(r3.size(), 0.0); const float r_max = r3.back(); const int half = static_cast(std::ceil(r_max)); const int side = 2 * half + 1; @@ -48,45 +52,57 @@ void background_aperture::MeasureFrame(ImagePreprocessorBuffer &image, int width return; const int32_t *tiles = image.GatherTiles(cx, cy, width, half); - std::vector ring_sum(r3.size()); - std::vector ring_n(r3.size()); - std::vector ring_px; - std::vector ring_d2; + struct Pixel { + int32_t v; + float d2; // squared distance from the aperture centre + double b; // the ring background at this pixel's resolution + }; + std::vector disk, ring; for (size_t a = 0; a < cx.size(); ++a) { const int32_t *t = tiles + a * static_cast(side) * side; - double disk_sum = 0.0; - int disk_n = 0; + const int64_t panel = modules.ModuleOfImagePixel(cx[a], cy[a]); bool readable = true; - ring_px.clear(); - ring_d2.clear(); + disk.clear(); + ring.clear(); for (int dy = -half; dy <= half && readable; ++dy) for (int dx = -half; dx <= half; ++dx) { const float d2 = static_cast(dx * dx + dy * dy); - if (d2 >= r_max * r_max) + if (d2 >= r_max * r_max || (d2 >= r1 * r1 && d2 < r2 * r2)) continue; const int32_t v = t[(dy + half) * side + dx + half]; if (v == INT32_MIN || v == INT32_MAX) { readable = false; break; } - if (d2 < r1 * r1) { - disk_sum += v; - ++disk_n; - } else if (d2 >= r2 * r2) { - ring_px.push_back(v); - ring_d2.push_back(d2); + double b = 0.0; + if (!rings.empty()) { + const auto [px, py] = modules.ImageToPlane(panel, cx[a] + dx, cy[a] + dy); + b = rings.At(2.0 * PI / geometry.PxlToRes(static_cast(px), static_cast(py))); } + (d2 < r1 * r1 ? disk : ring).push_back({v, d2, b}); } - if (!readable || disk_n == 0) + if (!readable || disk.empty()) continue; - // Every candidate ring as the integrator reads it: its mean, then one high-side clip pass. + double disk_v = 0.0, disk_b = 0.0; + for (const auto &p : disk) { + disk_v += p.v; + disk_b += p.b; + } + disk_v /= static_cast(disk.size()); + disk_b /= static_cast(disk.size()); + + // Every candidate ring as the integrator reads it: the mean, then one high-side clip pass. The + // radial correction is taken over the ring's geometry, as the integrator's kernel is. + std::vector d_plain(r3.size()), d_radial(r3.size()); bool measured = true; for (size_t k = 0; k < r3.size() && measured; ++k) { - double sum = 0.0; + const float r3_sq = r3[k] * r3[k]; + double sum = 0.0, sum_b = 0.0; int n = 0; - for (size_t p = 0; p < ring_px.size(); ++p) - if (ring_d2[p] < r3[k] * r3[k]) { - sum += ring_px[p]; + for (const auto &p : ring) + if (p.d2 < r3_sq) { + sum += p.v; + sum_b += p.b; ++n; } if (n == 0) { @@ -95,28 +111,55 @@ void background_aperture::MeasureFrame(ImagePreprocessorBuffer &image, int width } const double mean = sum / n; const double thr = mean + clip_nsigma * std::sqrt(gain * std::max(mean, 1.0)); - ring_sum[k] = 0.0; - ring_n[k] = 0; - for (size_t p = 0; p < ring_px.size(); ++p) - if (ring_d2[p] < r3[k] * r3[k] && (clip_nsigma <= 0.0f || ring_px[p] <= thr)) { - ring_sum[k] += ring_px[p]; - ++ring_n[k]; + double kept = 0.0; + int n_kept = 0; + for (const auto &p : ring) + if (p.d2 < r3_sq && (clip_nsigma <= 0.0f || p.v <= thr)) { + kept += p.v; + ++n_kept; } - measured = ring_n[k] > 0; + measured = n_kept > 0; + if (measured) { + d_plain[k] = disk_v - kept / n_kept; + d_radial[k] = d_plain[k] - (disk_b - sum_b / n); + } } if (!measured) continue; - const double disk = disk_sum / disk_n; for (size_t k = 0; k < r3.size(); ++k) { - const double d = disk - ring_sum[k] / ring_n[k]; - tally.sum_d2[k] += d * d; + tally.sum_d2[k] += d_plain[k] * d_plain[k]; + if (!rings.empty()) + tally.sum_d2_radial[k] += d_radial[k] * d_radial[k]; } ++tally.apertures; + if (!rings.empty()) + ++tally.apertures_radial; } } -int background_aperture::Best(const Tally &tally) { - if (tally.apertures == 0 || tally.sum_d2.empty()) +void background_aperture::Tally::Add(const Tally &other) { + sum_d2.resize(std::max(sum_d2.size(), other.sum_d2.size()), 0.0); + sum_d2_radial.resize(std::max(sum_d2_radial.size(), other.sum_d2_radial.size()), 0.0); + for (size_t k = 0; k < other.sum_d2.size(); ++k) + sum_d2[k] += other.sum_d2[k]; + for (size_t k = 0; k < other.sum_d2_radial.size(); ++k) + sum_d2_radial[k] += other.sum_d2_radial[k]; + apertures += other.apertures; + apertures_radial += other.apertures_radial; +} + +std::vector background_aperture::Tally::MeanSquare(bool radial) const { + const size_t n = radial ? apertures_radial : apertures; + std::vector out; + if (n == 0) + return out; + for (const double v : radial ? sum_d2_radial : sum_d2) + out.push_back(v / static_cast(n)); + return out; +} + +int background_aperture::Best(const std::vector &mean_square) { + if (mean_square.empty()) return -1; - return static_cast(std::min_element(tally.sum_d2.begin(), tally.sum_d2.end()) - tally.sum_d2.begin()); + return static_cast(std::min_element(mean_square.begin(), mean_square.end()) - mean_square.begin()); } diff --git a/image_analysis/bragg_integration/BackgroundAperture.h b/image_analysis/bragg_integration/BackgroundAperture.h index 7227301d6..e927fae29 100644 --- a/image_analysis/bragg_integration/BackgroundAperture.h +++ b/image_analysis/bragg_integration/BackgroundAperture.h @@ -7,8 +7,11 @@ #include #include +#include "../../common/DetectorModules.h" +#include "../../common/DiffractionGeometry.h" #include "../../common/DiffractionSpot.h" #include "../image_preprocessing/ImagePreprocessorBuffer.h" +#include "RingBackground.h" // How far out the background ring should reach, measured on apertures with nothing in them. // @@ -27,12 +30,24 @@ // a reflection: its noise plus its bias. The disk's own noise is the same for every candidate, so it // does not move the choice. The ring is clipped as the integrator clips it (high side, at // mean + nsigma * sqrt(gain * mean)). +// +// The same apertures also say whether the radial background-curvature correction helps: each difference +// is taken a second time with the correction applied - the frame's own ring background (RingBackground) +// averaged over the ring's pixels minus over the disk's, which is what the integrator's kernel takes off a +// reflection. Where the background really is a function of the radius the corrected ring predicts the +// disk better; where it is not - a strongly polarised beam, panels of different response - it predicts it +// worse, and the run measures which. The ring is chosen on the plain differences first, so its reach does +// not rest on the radial model; the correction is then kept where it improves that ring. namespace background_aperture { -// The squared disk-minus-ring differences, summed per candidate r3, over every aperture measured. +// The squared disk-minus-ring differences, summed per candidate r3 over every aperture measured: plain, +// and with the radial correction (counted separately: a frame without a ring background adds none). struct Tally { - std::vector sum_d2; - size_t apertures = 0; + std::vector sum_d2, sum_d2_radial; + size_t apertures = 0, apertures_radial = 0; + void Add(const Tally &other); + // The mean square difference of each candidate, plain or corrected; empty where nothing was measured. + [[nodiscard]] std::vector MeanSquare(bool radial) const; }; // The candidate outer radii: from r2 + 2 up to max_r3, one pixel apart. @@ -40,13 +55,16 @@ struct Tally { // The apertures of one frame, appended to `tally`. Centres lie on a grid over the frame, and an aperture // is used only where every pixel out to the largest candidate is readable and no found spot is within -// that radius plus SPOT_CLEARANCE_PX of its centre. +// that radius plus SPOT_CLEARANCE_PX of its centre. `rings` is this frame's ring background, or empty; +// `geometry` and `modules` place the pixels on the plane to read it at their resolution. void MeasureFrame(ImagePreprocessorBuffer &image, int width, int height, const std::vector &spots, float r1, float r2, - const std::vector &r3, float clip_nsigma, float gain, Tally &tally); + const std::vector &r3, float clip_nsigma, float gain, + const RingBackground &rings, const DiffractionGeometry &geometry, + const DetectorModules &modules, Tally &tally); -// The candidate with the smallest mean square difference, or -1 where no aperture was measured. -[[nodiscard]] int Best(const Tally &tally); +// The candidate with the smallest mean square difference, or -1 where nothing was measured. +[[nodiscard]] int Best(const std::vector &mean_square); constexpr float SPOT_CLEARANCE_PX = 4.0f; // At most this many apertures a frame: a 16 Mpx detector would otherwise lay twenty thousand. diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp index fec16cea1..23647ad9c 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp @@ -13,6 +13,7 @@ #include "../../common/JFJochMath.h" // PI (M_PI is not standard, and MSVC does not define it) #include "../SensorAbsorption.h" +#include "RingBackground.h" namespace { @@ -221,45 +222,16 @@ void BraggIntegrationEngine::BuildRadialKernel(float grow) { void BraggIntegrationEngine::RadialBackgroundFromRings(const std::vector &ring_bkg, int32_t q_bins, float low_q, float dq) { radial_curve.clear(); - if (q_bins <= 0 || !(dq > 0.0f)) + const RingBackground rings(ring_bkg, q_bins, low_q, dq); + if (rings.empty()) return; - // The rings averaged over their azimuthal bins, where the profile has any. - std::vector sum(q_bins, 0.0); - std::vector cnt(q_bins, 0); - for (size_t i = 0; i < ring_bkg.size(); ++i) - if (std::isfinite(ring_bkg[i])) { - sum[i % q_bins] += ring_bkg[i]; - ++cnt[i % q_bins]; - } - std::vector q, b; - for (int i = 0; i < q_bins; ++i) - if (cnt[i] > 0) { - q.push_back(low_q + (i + 0.5) * dq); - b.push_back(sum[i] / cnt[i]); - } - if (q.size() < 2) - return; - // The same bins as the engines' own curve: one per pixel of distance from the beam, out to the far - // corner. Each is read off the rings by linear interpolation in q; beyond the first and the last - // ring the end value is held, which has no curvature to correct. + // corner. const int n = static_cast(std::ceil(r_max)) + 2; const double two_k = 4.0 * PI / geom.GetWavelength_A(); radial_curve.resize(n); - size_t j = 0; - for (int r = 0; r < n; ++r) { - const double qr = two_k * std::sin(0.5 * std::atan(r / F_px)); - while (j + 1 < q.size() && q[j + 1] < qr) - ++j; - double v; - if (qr <= q.front()) - v = b.front(); - else if (j + 1 >= q.size()) - v = b.back(); - else - v = b[j] + (b[j + 1] - b[j]) * (qr - q[j]) / (q[j + 1] - q[j]); - radial_curve[r] = static_cast(v); - } + for (int r = 0; r < n; ++r) + radial_curve[r] = static_cast(rings.At(two_k * std::sin(0.5 * std::atan(r / F_px)))); } std::vector BraggIntegrationEngine::Panels(const std::vector &predicted, diff --git a/image_analysis/bragg_integration/CMakeLists.txt b/image_analysis/bragg_integration/CMakeLists.txt index 3a632ffd7..d90545394 100644 --- a/image_analysis/bragg_integration/CMakeLists.txt +++ b/image_analysis/bragg_integration/CMakeLists.txt @@ -20,6 +20,8 @@ ADD_LIBRARY(JFJochBraggIntegration STATIC BackgroundGain.h BackgroundAperture.cpp BackgroundAperture.h + RingBackground.cpp + RingBackground.h ) TARGET_LINK_LIBRARIES(JFJochBraggIntegration JFJochImagePreprocessing JFJochCommon) diff --git a/image_analysis/bragg_integration/RingBackground.cpp b/image_analysis/bragg_integration/RingBackground.cpp new file mode 100644 index 000000000..746ba8e21 --- /dev/null +++ b/image_analysis/bragg_integration/RingBackground.cpp @@ -0,0 +1,33 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "RingBackground.h" + +#include +#include + +RingBackground::RingBackground(const std::vector &ring_bkg, int32_t q_bins, float low_q, float dq) { + if (q_bins <= 0 || !(dq > 0.0f)) + return; + std::vector sum(q_bins, 0.0); + std::vector cnt(q_bins, 0); + for (size_t i = 0; i < ring_bkg.size(); ++i) + if (std::isfinite(ring_bkg[i])) { + sum[i % q_bins] += ring_bkg[i]; + ++cnt[i % q_bins]; + } + for (int i = 0; i < q_bins; ++i) + if (cnt[i] > 0) { + q.push_back(low_q + (i + 0.5) * dq); + b.push_back(sum[i] / cnt[i]); + } +} + +double RingBackground::At(double q_recipA) const { + if (q_recipA <= q.front()) + return b.front(); + if (q_recipA >= q.back()) + return b.back(); + const size_t j = std::upper_bound(q.begin(), q.end(), q_recipA) - q.begin() - 1; + return b[j] + (b[j + 1] - b[j]) * (q_recipA - q[j]) / (q[j + 1] - q[j]); +} diff --git a/image_analysis/bragg_integration/RingBackground.h b/image_analysis/bragg_integration/RingBackground.h new file mode 100644 index 000000000..58b2bb206 --- /dev/null +++ b/image_analysis/bragg_integration/RingBackground.h @@ -0,0 +1,20 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +// A frame's background as a function of q, from a peak-excluded background per resolution ring of the +// whole detector (the adaptive spot finder's GetRingBackground): ring i holds low_q + (i + 0.5) * dq, +// entries past the first q_bins are further azimuthal bins and are averaged in, and a ring with no +// value (NAN) is bridged by its neighbours. Between rings the value is interpolated linearly; beyond +// the first and the last ring the end value is held, which has no curvature. +struct RingBackground { + std::vector q, b; + + RingBackground(const std::vector &ring_bkg, int32_t q_bins, float low_q, float dq); + [[nodiscard]] bool empty() const { return q.size() < 2; } + [[nodiscard]] double At(double q_recipA) const; +}; diff --git a/rugnux/RugnuxPreScan.cpp b/rugnux/RugnuxPreScan.cpp index 1eca770ad..060cb54c3 100644 --- a/rugnux/RugnuxPreScan.cpp +++ b/rugnux/RugnuxPreScan.cpp @@ -8,6 +8,7 @@ #include "../image_analysis/bragg_integration/SpotWidth.h" #include "../image_analysis/bragg_integration/BackgroundGain.h" #include "../image_analysis/bragg_integration/BackgroundAperture.h" +#include "../image_analysis/bragg_integration/RingBackground.h" #include "../image_analysis/bragg_integration/SpotFootprint.h" #include "../image_analysis/bragg_integration/BraggStencil.h" // BRAGG_FOOTPRINT_NSIGMA #include "../image_analysis/hot_pixels/HotPixels.h" @@ -425,11 +426,21 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru MeasureSpotsOnTiles(*w.preprocessed, static_cast(prescan_x.GetXPixelsNumConv()), static_cast(prescan_x.GetYPixelsNumConv()), prescan_x.GetDiffractionGeometry(), prescan_x.GetModules(), spots, curves, footprint); - if (for_width && !aperture_r3.empty()) + if (for_width && !aperture_r3.empty()) { + // The frame's ring background, as the integration of this frame will read it - from the + // adaptive finder, which has it whole-frame and peak-excluded. + const std::vector no_rings; + const RingBackground rings(config_.spot_finding.adaptive_threshold ? w.spot_finder->GetRingBackground() + : no_rings, + prescan_mapping->GetQBinCount(), + prescan_mapping->Settings().GetLowQ_recipA(), + prescan_mapping->Settings().GetQSpacing_recipA()); background_aperture::MeasureFrame(*w.preprocessed, static_cast(prescan_x.GetXPixelsNumConv()), static_cast(prescan_x.GetYPixelsNumConv()), spots, aperture_bis.GetR1(), aperture_bis.GetR2(), aperture_r3, - aperture_bis.GetBackgroundClipNSigma(), aperture_gain, apertures); + aperture_bis.GetBackgroundClipNSigma(), aperture_gain, rings, + prescan_x.GetDiffractionGeometry(), prescan_x.GetModules(), apertures); + } }; std::vector prescan_spot_q; @@ -730,13 +741,9 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru if (!shape_only[i]) prescan_spot_q.insert(prescan_spot_q.end(), spot_q_of[i].begin(), spot_q_of[i].end()); // In sample order, over the frames the width was taken on. - aperture_tally.sum_d2.assign(aperture_r3.size(), 0.0); for (size_t i = 0; i < ordinals.size(); i++) - if (position[i] < width_images && !apertures_of[i].sum_d2.empty()) { - for (size_t k = 0; k < aperture_r3.size(); k++) - aperture_tally.sum_d2[k] += apertures_of[i].sum_d2[k]; - aperture_tally.apertures += apertures_of[i].apertures; - } + if (position[i] < width_images) + aperture_tally.Add(apertures_of[i]); // In sample order, so the table does not depend on how the workers interleaved. for (const auto &f : footprint_of) footprint_spots.insert(footprint_spots.end(), f.begin(), f.end()); @@ -779,26 +786,44 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // Stills keep their signal disk and take the ring's outer radius from the empty apertures // instead (BackgroundAperture.h): the candidate whose ring predicted the background under the // disk best. Where the noise wins that is the widest, which is the ring the run started from. + // The same apertures decide the radial curvature correction, where it was left to auto: on + // where the corrected ring predicts the disk better. if (want_width && stills) { spot_width_measured_ = true; - const int best = background_aperture::Best(aperture_tally); + BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings(); + const auto radial_setting = bis.GetBackgroundRadialCorrection(); + // The ring is chosen on the plain differences, so its reach does not rest on the radial + // model; the correction is then kept, where it was left to auto, if it makes that ring predict + // the disk better. Read only where every frame gave a ring background. + const std::vector plain = aperture_tally.MeanSquare(false); + const std::vector radial = aperture_tally.apertures_radial == aperture_tally.apertures + ? aperture_tally.MeanSquare(true) : std::vector(); + const int best = background_aperture::Best(plain); + const bool use_radial = radial_setting.has_value() + ? *radial_setting + : best >= 0 && !radial.empty() && radial[best] < plain[best]; if (best < 0) { logger.Info("Background ring: no empty aperture measured on {} images, keeping r3={:.1f}", - width_images, experiment_.GetBraggIntegrationSettings().GetR3()); + width_images, bis.GetR3()); } else { - std::string table; - for (size_t k = 0; k < aperture_r3.size(); k++) - table += fmt::format(" {:.0f}:{:.4g}", aperture_r3[k], - aperture_tally.sum_d2[k] / static_cast(aperture_tally.apertures)); - BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings(); + const auto table = [&](const std::vector &m) { + std::string t; + for (size_t k = 0; k < m.size(); k++) + t += fmt::format(" {:.0f}:{:.4g}", aperture_r3[k], m[k]); + return t; + }; bis.R3(aperture_r3[best]); + if (!radial_setting) + bis.BackgroundRadialCorrection(use_radial ? std::nullopt : std::optional(false)); experiment_.ImportBraggIntegrationSettings(bis); logger.Info("Background ring: on {} empty apertures from {} images the r{:.0f}..r3 ring predicts " - "the r{:.0f} disk with a mean square error of (r3:value){} => integration radii " - "r1={:.1f} r2={:.1f} r3={:.1f}", aperture_tally.apertures, width_images, - bis.GetR2(), bis.GetR1(), table, bis.GetR1(), bis.GetR2(), bis.GetR3()); - Note(fmt::format("background ring {:.0f}..{:.0f} px, measured on empty apertures", - bis.GetR2(), bis.GetR3())); + "the r{:.0f} disk with a mean square error of (r3:value){}, and with the radial " + "curvature correction{} => integration radii r1={:.1f} r2={:.1f} r3={:.1f}, radial " + "correction {}", aperture_tally.apertures, width_images, bis.GetR2(), bis.GetR1(), + table(plain), radial.empty() ? std::string(" (no ring background)") : table(radial), + bis.GetR1(), bis.GetR2(), bis.GetR3(), use_radial ? "on" : "off"); + Note(fmt::format("background ring {:.0f}..{:.0f} px, radial correction {}, measured on empty " + "apertures", bis.GetR2(), bis.GetR3(), use_radial ? "on" : "off")); } } if (want_width && !stills) { diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 3256a54d7..95cf46310 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -254,7 +254,7 @@ void print_usage() { std::cout << " --integration-high-resolution High resolution limit for prediction/integration. If omitted (or 0), integration extends as far as the detector reaches" << std::endl; std::cout << " --max-hkl Predict reflections with |h|,|k|,|l| <= n. Default: derived per crystal from the refined cell (ceil(longest axis / d_min) + 1), which is the exact bound - set it only to override that" << std::endl; std::cout << " --background-clip High-side clip of the background ring at mean + n*sqrt(mean) (default 4; 0 = off). This is the default background estimator - it rejects neighbour cores and zingers without the symmetric trim's Poisson skew bias. Ignored by --integrator boxsum" << std::endl; - std::cout << " --background-radial[=on|off|auto] Correct the background ring for the CURVATURE of the radial background (default auto). The signal disk and the background ring are concentric, so a background linear in position cancels between them and only curvature survives - which on a smooth ice ring reaches +26 counts on a single reflection. =auto applies it per image where that image's ice score shows a smooth powder ring, which is where a radius-only background model holds; on ice made of discrete crystallite spots there is no such ring and the correction makes the bias worse. On stills =auto applies it to every image, with the curve taken from the spot finder's peak-excluded background of the whole frame: a smooth solvent or carrier hump curves the background as much as ice does. Costs one short dot product per reflection and no extra pixel reads" << std::endl; + std::cout << " --background-radial[=on|off|auto] Correct the background ring for the CURVATURE of the radial background (default auto). The signal disk and the background ring are concentric, so a background linear in position cancels between them and only curvature survives - which on a smooth ice ring reaches +26 counts on a single reflection. =auto applies it per image where that image's ice score shows a smooth powder ring, which is where a radius-only background model holds; on ice made of discrete crystallite spots there is no such ring and the correction makes the bias worse. On stills the curve is taken from the spot finder's peak-excluded background of the whole frame, and =auto applies it to every image of a run where it makes the background ring predict the background of empty apertures better (measured in the pre-scan): a smooth solvent or carrier hump curves the background as much as ice does. Costs one short dot product per reflection and no extra pixel reads" << std::endl; std::cout << " --background-trim Use the old symmetric trimmed mean for the background ring instead of the clip (0<=f<0.5; 0.10 was the former default). Switches --background-clip off. A symmetric trim is biased low on Poisson data and adds ~5 counts to every partial, so this is for back compatibility only; 0 = plain ring mean" << std::endl; std::cout << " --overlap What to do where two predicted reflections share signal pixels: off|reject|exclude (default exclude). A pixel inside two signal disks belongs to the nearer centre; before this, nothing kept a neighbour's flux out of a reflection's own disk, so on a dense pattern a crowded reflection read high. exclude drops the shared PIXELS from the profile fit, which renormalises itself, so the reflection is kept; reject instead drops the whole reflection when less than --overlap-minpk of its expected profile is cleanly its own (what XDS calls MINPK). --integrator boxsum has no profile to renormalise with, so exclude does nothing there and only reject acts" << std::endl; std::cout << " --overlap-minpk Least fraction of a reflection's expected profile that must be usable for the reflection to be kept (default 0.75, XDS MINPK). Governs both ways part of a profile is lost: the fraction that must be READABLE - not masked, untrusted, in a detector gap or overloaded - in every profile mode, and, under --overlap reject, the fraction that must be cleanly the reflection's own. --integrator boxsum has no profile to renormalise with, so there a signal disk with any unreadable pixel is still discarded outright and the reject fraction is by disk AREA, which cuts harder" << std::endl; @@ -3429,7 +3429,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { radial.has_value() ? (*radial ? "on" : "off") : rotation_indexing ? "auto (per image, where the ice score sees a smooth ring)" - : "auto (every image, from the radial background of the whole frame)"); + : "auto (decided in the pre-scan on empty apertures; the curve from the whole frame)"); } // Detection defaults differ by workflow (RugnuxDefaultSpotFinding); each flag always wins. diff --git a/tests/BackgroundApertureTest.cpp b/tests/BackgroundApertureTest.cpp index 3121546ed..bedf72ee5 100644 --- a/tests/BackgroundApertureTest.cpp +++ b/tests/BackgroundApertureTest.cpp @@ -26,15 +26,25 @@ ImagePreprocessorBuffer Frame(double level, double ripple, uint32_t seed) { return img; } +const RingBackground no_rings({}, 0, 0.0f, 0.0f); + +DiffractionGeometry Geometry() { + DiffractionGeometry g; + g.BeamX_pxl(300.0f).BeamY_pxl(300.0f).DetectorDistance_mm(60.0f).PixelSize_mm(0.075f).Wavelength_A(1.0f); + return g; +} + int BestR3(double level, double ripple) { const auto r3 = background_aperture::CandidateR3(8.0f, 14.0f); background_aperture::Tally tally; for (uint32_t seed = 1; seed <= 4; seed++) { auto img = Frame(level, ripple, seed); - background_aperture::MeasureFrame(img, W, H, {}, 6.0f, 8.0f, r3, 4.0f, 1.0f, tally); + background_aperture::MeasureFrame(img, W, H, {}, 6.0f, 8.0f, r3, 4.0f, 1.0f, no_rings, Geometry(), + DetectorModules(), tally); } REQUIRE(tally.apertures > 100); - return static_cast(r3[background_aperture::Best(tally)]); + CHECK(tally.apertures_radial == 0); + return static_cast(r3[background_aperture::Best(tally.MeanSquare(false))]); } } // namespace @@ -62,11 +72,47 @@ TEST_CASE("BackgroundAperture_SkipsUnreadablePixels", "[Integration]") { const auto r3 = background_aperture::CandidateR3(8.0f, 14.0f); auto img = Frame(1.0, 0.0, 7); background_aperture::Tally all; - background_aperture::MeasureFrame(img, W, H, {}, 6.0f, 8.0f, r3, 4.0f, 1.0f, all); + background_aperture::MeasureFrame(img, W, H, {}, 6.0f, 8.0f, r3, 4.0f, 1.0f, no_rings, Geometry(), + DetectorModules(), all); for (size_t i = 0; i < static_cast(W) * H; i += 97) img[i] = INT32_MIN; background_aperture::Tally masked; - background_aperture::MeasureFrame(img, W, H, {}, 6.0f, 8.0f, r3, 4.0f, 1.0f, masked); + background_aperture::MeasureFrame(img, W, H, {}, 6.0f, 8.0f, r3, 4.0f, 1.0f, no_rings, Geometry(), + DetectorModules(), masked); CHECK(masked.apertures == 0); CHECK(all.apertures > 0); } + +// A bright background that is a curved function of the radius alone: the radial correction, read off +// the frame's ring background, makes the ring predict the disk better than it does without it. +TEST_CASE("BackgroundAperture_RadialCorrectionHelpsOnARadialBackground", "[Integration]") { + const DiffractionGeometry g = Geometry(); + const auto background = [](double r) { return 2000.0 + 0.2 * (r - 150.0) * (r - 150.0); }; + std::mt19937 rng(3); + ImagePreprocessorBuffer img(static_cast(W) * H); + for (int y = 0; y < H; y++) + for (int x = 0; x < W; x++) { + std::poisson_distribution p(background(std::hypot(x - 300.0, y - 300.0))); + img[static_cast(y) * W + x] = p(rng); + } + const float low_q = 0.05f, dq = 0.005f; + std::vector profile(1000); + for (size_t i = 0; i < profile.size(); i++) { + const double q = low_q + (i + 0.5) * dq; + const double r = 60.0 / 0.075 * std::tan(2.0 * std::asin(q * 1.0 / (4.0 * PI))); + profile[i] = static_cast(background(r)); + } + const RingBackground rings(profile, static_cast(profile.size()), low_q, dq); + const auto r3 = background_aperture::CandidateR3(8.0f, 14.0f); + background_aperture::Tally tally; + // Unclipped: a gradient this steep would put the high side of the ring above a photon clip. + background_aperture::MeasureFrame(img, W, H, {}, 6.0f, 8.0f, r3, 0.0f, 1.0f, rings, g, DetectorModules(), tally); + REQUIRE(tally.apertures_radial == tally.apertures); + REQUIRE(tally.apertures > 100); + const auto plain = tally.MeanSquare(false), radial = tally.MeanSquare(true); + // The pre-scan's question: at the ring the plain differences choose, does the correction help. + const int best = background_aperture::Best(plain); + CHECK(radial[best] < 0.8 * plain[best]); + // With the curvature taken out, the noise wins again and a wide ring is best. + CHECK(r3[background_aperture::Best(radial)] >= 12.0f); +}