From 324a2b0cdff71973d6a6fe2ad979bf11893745ff Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 6 Sep 2026 07:10:57 +0200 Subject: [PATCH] integration: the aperture guard watches the aperture, and the width law can be read at any resolution Two changes to the same file, both no-ops on every run today, measured bit-identical across the corpus on every merged statistic. The starvation guard remembers the pre-adaptive settings so it has something to fall back to, and it remembered them only when r1 moved. r1 is a proxy for the aperture, not the aperture: a prototype that widened the background radii at high angle while leaving r1 alone took the neighbour-starved fraction of a dense low-symmetry crystal to eight times the bound this guard exists to enforce, and the guard could not fire, because nothing had been remembered. It now tests all three radii. Today they move together, so nothing changes; the hole was latent, and latent is where it should be closed. The width fit computes how much wider a crystal's spots get at high angle - the 1/d coefficient of r80(d) - and threw it away, returning only the value at one resolution. The law is now returned whole and held at the ends of the measured range rather than extrapolated, since two coefficients may rest on as few as two bands. Read at that one resolution it reproduces the old number bit for bit, which is why the coefficients are carried in double. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N --- rugnux/Rugnux.cpp | 14 ++++++++++++-- rugnux/SpotWidth.cpp | 36 +++++++++++++++++++++++++++--------- rugnux/SpotWidth.h | 25 ++++++++++++++++++++++--- 3 files changed, 61 insertions(+), 14 deletions(-) diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index c08a7617e..0912de2a4 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1076,9 +1076,19 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru const float r2 = r1 + 2.0f; const float r3 = std::sqrt(r2 * r2 + 133.0f); BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings(); - // Remembered only where the radius actually moved, so the guard in RunAllPasses has + // Remembered where the APERTURE moved, so the starvation guard in RunAllPasses has // something to give back and cannot announce a fallback that changes nothing. - if (r1 != bis.GetR1()) bragg_before_adaptive_ = bis; + // + // The test is on all three radii and not on r1 alone. It is the same test today - r2 and + // r3 are taken from r1 just above, and the adaptive rule does not run at all when the + // radii were given explicitly - so this changes no run. It is written this way because + // what the guard is about is the aperture, and r1 is only a proxy for it. A prototype + // that widened r2 and r3 at high angle while leaving r1 at the shipped 4 took the + // neighbour-starved fraction of a dense low-symmetry crystal from 2.4% to 9.0%, eight + // times the bound this guard exists to enforce, and the guard could not fire: nothing + // had set bragg_before_adaptive_, because r1 had not moved. + if (r1 != bis.GetR1() || r2 != bis.GetR2() || r3 != bis.GetR3()) + bragg_before_adaptive_ = bis; bis.R1(r1).R2(r2).R3(r3); experiment_.ImportBraggIntegrationSettings(bis); logger.Info("Spot width: r80 = {:.2f} px at {:.0f} A ({} spots from {} images) => " diff --git a/rugnux/SpotWidth.cpp b/rugnux/SpotWidth.cpp index c31d1ccde..9095f6dbf 100644 --- a/rugnux/SpotWidth.cpp +++ b/rugnux/SpotWidth.cpp @@ -272,7 +272,11 @@ void MeasureSpotFluxCurves(const ImagePreprocessorBuffer &image, int width, int } } -std::optional spot_width::R80AtReference(const std::vector &curves) { +float spot_width::R80Fit::At(double d_A) const { + return static_cast(std::clamp(c0 + c1 / d_A, lo, hi)); +} + +std::optional spot_width::FitR80(const std::vector &curves) { if (curves.size() < MIN_SPOTS_TOTAL) return std::nullopt; // One point per band: the median curve of the band, the radius it holds 80 % of its flux at, and @@ -302,7 +306,18 @@ std::optional spot_width::R80AtReference(const std::vector &cu static_cast(band_d.size())}); } if (points.empty()) return std::nullopt; - if (points.size() == 1) return static_cast(points[0].r80); + + // Never extrapolate outside what the bands actually measured. A single band measures no slope, so + // its own value is the whole law; the bounds then bracket it and At() returns it unchanged. + R80Fit fit; + double lo = std::numeric_limits::max(), hi = 0.0; + for (const auto &p : points) { lo = std::min(lo, p.r80); hi = std::max(hi, p.r80); } + fit.lo = 0.8 * lo; + fit.hi = 1.25 * hi; + if (points.size() == 1) { + fit.c0 = points[0].r80; + return fit; + } // The mosaic contribution to the detector footprint grows as 1/d, so r80 is linear in 1/d. double sw = 0.0, sx = 0.0, sxx = 0.0, sy = 0.0, sxy = 0.0; @@ -314,15 +329,18 @@ std::optional spot_width::R80AtReference(const std::vector &cu sxy += p.weight * p.inv_d * p.r80; } const double det = sw * sxx - sx * sx; - double value = sy / sw; + fit.c0 = sy / sw; if (std::abs(det) > 1e-12) { - const double c1 = (sw * sxy - sx * sy) / det; - value = (sy - c1 * sx) / sw + c1 / D_REF_A; + fit.c1 = (sw * sxy - sx * sy) / det; + fit.c0 = (sy - fit.c1 * sx) / sw; } - // Never extrapolate outside what the bands actually measured. - double lo = std::numeric_limits::max(), hi = 0.0; - for (const auto &p : points) { lo = std::min(lo, p.r80); hi = std::max(hi, p.r80); } - return static_cast(std::clamp(value, 0.8 * lo, 1.25 * hi)); + return fit; +} + +std::optional spot_width::R80AtReference(const std::vector &curves) { + const auto fit = FitR80(curves); + if (!fit) return std::nullopt; + return fit->At(D_REF_A); } float spot_width::R1ForWidth(float r80) { diff --git a/rugnux/SpotWidth.h b/rugnux/SpotWidth.h index e5b450c7e..00fdb3d60 100644 --- a/rugnux/SpotWidth.h +++ b/rugnux/SpotWidth.h @@ -52,9 +52,28 @@ struct FluxCurve { std::array c{}; }; -// r80 at D_REF_A from a pool of curves, or nothing where too few spots were measurable. Fitted as -// r80(d) = c0 + c1/d over resolution bands, weighted by the spots in each, so crystals whose strong -// spots sit at different resolutions are compared like with like. +// The fitted width law: r80(d) = c0 + c1/d over resolution bands, weighted by the spots in each, so +// crystals whose strong spots sit at different resolutions are compared like with like. The fit is +// over 1/d because the mosaic contribution to the detector footprint grows that way, so c1 is the +// crystal's own answer to how much wider its spots get at high angle - a quantity the fit has always +// computed and nothing has ever been able to read, because only the value at one resolution was +// returned. +struct R80Fit { + // Kept in double, and evaluated in double, because At() has to reproduce the single number the + // caller used to get bit for bit: rounding the coefficients to float first and adding them after + // moves the last bit, and r1 = round(2*r80) can turn on it. + double c0 = 0.0, c1 = 0.0; + // What the bands actually measured, widened the way this fit has always widened it. The law is + // HELD at these outside the measured range rather than extrapolated - it has two coefficients and + // as few as two points, so a slope read off a short lever arm must not be allowed to run away. + double lo = 0.0, hi = 0.0; + [[nodiscard]] float At(double d_A) const; +}; + +// The width law from a pool of curves, or nothing where too few spots were measurable. +[[nodiscard]] std::optional FitR80(const std::vector &curves); + +// r80 at D_REF_A, which is the law above read at one resolution. [[nodiscard]] std::optional R80AtReference(const std::vector &curves); // The signal radius the measured width asks for: r1 = clamp(round(2*r80), 4, 6). The floor is the