Powder calibration: read the rings off an azimuthal profile, not off a spot list

The ring calibration already here (AssignSpotsToRings + RingOptimizer, driven from
the viewer's powder panel) is given a SPOT LIST from a single image. A powder ring
is not a set of spots - it is a smooth arc - so a spot finder samples it wherever
its threshold happens to bite, and one image carries only the counts that image
collected. An azimuthally-binned profile summed over a run measures the same ring
directly, at every azimuth, with the whole run behind it.

RingsFromAzimuthalProfile turns such a profile into the (x, y, q_expected) triples
RingOptimizer already consumes, so nothing downstream changes: for each calibrant
ring and each azimuthal sector it fits the radial peak against a locally
interpolated background, and maps the measured (q, phi) back through the current
geometry to the pixel it came from.

What this is for is the BEAM CENTRE. A powder ring is a conic centred on the beam,
so a wrong centre makes its apparent radius oscillate once per turn and a detector
tilt twice - and neither depends on the calibrant's d-spacings or on the detector
distance. That matters, because the beam centre is otherwise the weakest parameter
we have: fitted from Bragg spots it is gauge-coupled to the crystal orientation,
which is why PostRefine has to restrain it toward the header and commit only a
sub-1 % move, and why XtalOptimizer carries a soft prior noting the beam is "only
LaB6-monitored to ~a few px". A ring does not know about the crystal.

Two things the peak fit is careful about, both of which would otherwise show up as
a spurious cos(phi) - i.e. as a beam-centre shift:

 - the sector's CENTRE is used, not its lower edge. GetBin() floors phi into the
   sector, so a bin stands for [j, j+1), and taking its edge rotates every ring
   point by half a sector.
 - a peak has to stand clear of the scatter of the background either side of it,
   or a sector with no ring in it contributes its largest noise excursion as
   though it were a measurement.

Refuses a single-azimuthal-bin profile outright: that is a plain radial profile,
the ring has been averaged over every direction, and there is nothing left to say
where its centre is.

Tested by round trip against the forward model, as the existing calibration tests
are: synthesise the profile the azimuthal integration would build with the rings
where a shifted geometry puts them but every pixel binned with the unshifted one,
then extract and fit. A 6.0 / -4.0 px beam offset is recovered as 6.13 / -4.03
from 192 ring points. Only the beam centre is exercised here; the tilt path is
covered by the existing DetGeomCalibTest round trips.

This is the extraction only - nothing calls it yet, and the run-scoped accumulator
it is meant to read (JFJochReceiverPlots::az_int_profile, already summed over a run
and written to /entry/azint/dataset) is still integrated with one azimuthal bin by
default.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

image_analysis/geom_refinement/CMakeLists.txt

@@ -4,6 +4,8 @@ ADD_LIBRARY(JFJochGeomRefinement STATIC
         RingOptimizer.h
         AssignSpotsToRings.cpp
         AssignSpotsToRings.h
+        RingsFromProfile.cpp
+        RingsFromProfile.h
         XtalOptimizer.cpp
         XtalOptimizer.h
         XtalResidual.h

image_analysis/geom_refinement/RingsFromProfile.cpp

@@ -0,0 +1,125 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <algorithm>
+#include <cmath>
+
+#include "RingsFromProfile.h"
+#include "AssignSpotsToRings.h"   // CalculateXtalRings
+#include "../../common/JFJochMath.h"
+
+namespace {
+
+// Peak position of one ring in one azimuthal sector, in q, or NaN if there is no peak worth using.
+//
+// The window is narrow and centred on where the ring is expected, so the background under it is close
+// to a straight line: take it from the two bins at each end and interpolate. The position itself is the
+// intensity-weighted centroid of everything above half the peak height, which is insensitive to the
+// exact half-maximum crossing and needs no line-shape assumption - a powder ring is not Gaussian, it is
+// the instrumental profile convolved with whatever strain and size broadening the standard has.
+float SectorPeakQ(const std::vector<float> &profile, int32_t q_bins, int phi_bin,
+                  int lo_bin, int hi_bin, float low_q, float q_spacing, float min_peak_over_noise) {
+    const size_t row = static_cast<size_t>(phi_bin) * static_cast<size_t>(q_bins);
+    const auto value = [&](int i) { return profile[row + static_cast<size_t>(i)]; };
+    const auto q_of = [&](int i) { return low_q + (static_cast<float>(i) + 0.5f) * q_spacing; };
+
+    const float bkg_lo = 0.5f * (value(lo_bin) + value(lo_bin + 1));
+    const float bkg_hi = 0.5f * (value(hi_bin) + value(hi_bin - 1));
+    const auto bkg_at = [&](int i) {
+        const float t = static_cast<float>(i - lo_bin) / static_cast<float>(hi_bin - lo_bin);
+        return bkg_lo + t * (bkg_hi - bkg_lo);
+    };
+
+    int peak = -1;
+    float peak_height = 0.0f;
+    for (int i = lo_bin + 2; i <= hi_bin - 2; ++i) {
+        const float h = value(i) - bkg_at(i);
+        if (h > peak_height) { peak_height = h; peak = i; }
+    }
+    if (peak < 0)
+        return NAN;
+
+    // Scatter of the background shoulders, as the noise this peak has to stand clear of. A sector with
+    // no ring in it has a "peak" that is just the largest background fluctuation, and this is what
+    // rejects it - the alternative, an absolute intensity cut, would need a value per detector and beam.
+    float s = 0.0f;
+    int n = 0;
+    for (int i : {lo_bin, lo_bin + 1, hi_bin - 1, hi_bin}) {
+        const float r = value(i) - bkg_at(i);
+        s += r * r;
+        ++n;
+    }
+    const float noise = std::sqrt(s / static_cast<float>(n));
+    if (!(peak_height > min_peak_over_noise * noise))
+        return NAN;
+
+    const float half = 0.5f * peak_height;
+    double sum_wq = 0.0, sum_w = 0.0;
+    for (int i = peak; i >= lo_bin && value(i) - bkg_at(i) >= half; --i) {
+        const double w = value(i) - bkg_at(i);
+        sum_wq += w * q_of(i);
+        sum_w += w;
+    }
+    for (int i = peak + 1; i <= hi_bin && value(i) - bkg_at(i) >= half; ++i) {
+        const double w = value(i) - bkg_at(i);
+        sum_wq += w * q_of(i);
+        sum_w += w;
+    }
+    if (!(sum_w > 0.0))
+        return NAN;
+    return static_cast<float>(sum_wq / sum_w);
+}
+
+}   // namespace
+
+std::vector<RingOptimizerInput> RingsFromAzimuthalProfile(const std::vector<float> &profile,
+                                                          const AzimuthalIntegrationMapping &mapping,
+                                                          const DiffractionGeometry &geom,
+                                                          const UnitCell &calibrant,
+                                                          float q_window_recipA,
+                                                          float min_peak_over_noise) {
+    std::vector<RingOptimizerInput> out;
+
+    const int32_t q_bins = mapping.GetQBinCount();
+    const int32_t azim_bins = mapping.GetAzimuthalBinCount();
+    // One azimuthal bin is a plain radial profile: the ring is averaged over every direction at once, so
+    // nothing remains to say where its centre is. This needs the run to have been integrated with
+    // azimuthal bins (jfjoch_broker azim_int_settings.azimuthal_bins, rugnux --azim-phi-bins).
+    if (azim_bins < 4 || q_bins < 8
+        || profile.size() != static_cast<size_t>(q_bins) * static_cast<size_t>(azim_bins))
+        return out;
+
+    const auto &settings = mapping.Settings();
+    const float low_q = settings.GetLowQ_recipA();
+    const float q_spacing = settings.GetQSpacing_recipA();
+    const float high_q = low_q + static_cast<float>(q_bins) * q_spacing;
+    const int window_bins = std::max(3, static_cast<int>(std::lround(q_window_recipA / q_spacing)));
+
+    for (const float q_ring : CalculateXtalRings(calibrant)) {
+        if (!(q_ring - q_window_recipA > low_q) || !(q_ring + q_window_recipA < high_q))
+            continue;
+        const int centre_bin = static_cast<int>((q_ring - low_q) / q_spacing);
+        const int lo_bin = std::max(0, centre_bin - window_bins);
+        const int hi_bin = std::min(q_bins - 1, centre_bin + window_bins);
+        if (hi_bin - lo_bin < 6)
+            continue;
+
+        for (int phi_bin = 0; phi_bin < azim_bins; ++phi_bin) {
+            const float q_obs = SectorPeakQ(profile, q_bins, phi_bin, lo_bin, hi_bin,
+                                            low_q, q_spacing, min_peak_over_noise);
+            if (!std::isfinite(q_obs))
+                continue;
+
+            // The sector's CENTRE, not its lower edge: GetBin() floors phi into the sector, so a bin
+            // stands for [j, j+1) and taking its edge would rotate every ring point by half a sector -
+            // which is exactly the cos(phi) signal the beam centre is read from.
+            const float phi_rad = static_cast<float>((static_cast<double>(phi_bin) + 0.5)
+                                                     * 2.0 * PI / static_cast<double>(azim_bins));
+            const auto [x, y] = geom.ResPhiToPxl(static_cast<float>(2.0 * PI) / q_obs, phi_rad);
+            if (!std::isfinite(x) || !std::isfinite(y))
+                continue;
+            out.push_back({x, y, q_ring});
+        }
+    }
+    return out;
+}

image_analysis/geom_refinement/RingsFromProfile.h

@@ -0,0 +1,37 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <vector>
+
+#include "../../common/AzimuthalIntegrationMapping.h"
+#include "../../common/DiffractionGeometry.h"
+#include "../../common/UnitCell.h"
+#include "RingOptimizer.h"
+
+// Turn an accumulated (q x azimuth) powder profile into ring points for RingOptimizer.
+//
+// The calibration this feeds already exists (AssignSpotsToRings + RingOptimizer); what it has always
+// been given is a SPOT LIST from a single image. A powder ring is not a set of spots - it is a smooth
+// arc - so a spot finder samples it wherever its threshold happens to bite, and one image carries only
+// as much of the ring as that image's counting statistics allow. An azimuthally-binned profile summed
+// over a run measures the same ring directly, at every azimuth, with the whole run's counts behind it.
+//
+// Where the ring falls is what carries the geometry. A powder ring is a conic centred on the beam, so
+// if the beam centre is wrong its apparent radius oscillates once per turn (a cos(phi) term), and if
+// the detector is tilted, twice (cos(2 phi)). Neither depends on the calibrant's d-spacings, which is
+// why the beam centre is the one thing a powder pattern determines without assuming anything about the
+// standard - the distance, by contrast, is only as good as the lattice constant it is measured against.
+//
+// profile is the mean intensity per bin (AzimuthalIntegrationProfile::GetResult()): q_bins x azimuthal
+// bins, indexed bin = q_bin + phi_bin * q_bins. geom supplies the CURRENT geometry, used only to turn a
+// measured (q, phi) back into the pixel it came from - RingOptimizer then refines that geometry so the
+// q it predicts at that pixel matches the calibrant's. Rings outside the profile's q range, and sectors
+// where no peak stands clear of the local background, are skipped rather than guessed at.
+std::vector<RingOptimizerInput> RingsFromAzimuthalProfile(const std::vector<float> &profile,
+                                                          const AzimuthalIntegrationMapping &mapping,
+                                                          const DiffractionGeometry &geom,
+                                                          const UnitCell &calibrant,
+                                                          float q_window_recipA = 0.06f,
+                                                          float min_peak_over_noise = 3.0f);

tests/CMakeLists.txt

@@ -59,6 +59,7 @@ ADD_EXECUTABLE(jfjoch_test
         JFJochReceiverPlotsTest.cpp
         GoniometerAxisTest.cpp
         DetGeomCalibTest.cpp
+        RingsFromProfileTest.cpp
         XtalOptimizerTest.cpp
         CrystalLatticeTest.cpp
         FPGAPTPTest.cpp

tests/RingsFromProfileTest.cpp

@@ -0,0 +1,123 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <catch2/catch_all.hpp>
+
+#include "../image_analysis/geom_refinement/RingsFromProfile.h"
+#include "../image_analysis/geom_refinement/AssignSpotsToRings.h"
+#include "../common/Definitions.h"
+#include "../common/JFJochMath.h"
+
+namespace {
+
+constexpr UnitCell LAB6{LAB6_CELL_A, LAB6_CELL_A, LAB6_CELL_A, 90.0f, 90.0f, 90.0f};
+
+// A (q x azimuth) powder profile as the azimuthal integration would build it: the rings sit where
+// geom_true puts them, but every pixel is binned with geom_assumed - which is the whole point, since a
+// wrong assumed geometry is what makes a ring's apparent q wander with azimuth.
+std::vector<float> SynthesiseProfile(const AzimuthalIntegrationMapping &mapping,
+                                     const DiffractionGeometry &geom_assumed,
+                                     const DiffractionGeometry &geom_true) {
+    const auto &settings = mapping.Settings();
+    const int32_t q_bins = mapping.GetQBinCount();
+    const int32_t azim_bins = mapping.GetAzimuthalBinCount();
+    std::vector<float> profile(static_cast<size_t>(q_bins) * azim_bins, 100.0f);   // flat background
+
+    for (const float q_ring : CalculateXtalRings(LAB6)) {
+        const float d = static_cast<float>(2.0 * PI) / q_ring;
+        if (d <= geom_true.GetWavelength_A() / 2.0f)
+            continue;
+        for (int t = 0; t < 3600; ++t) {
+            const float phi_true = static_cast<float>(2.0 * PI * t / 3600.0);
+            const auto [px, py] = geom_true.ResPhiToPxl(d, phi_true);
+            if (!std::isfinite(px) || !std::isfinite(py))
+                continue;
+            const float q_obs = geom_assumed.PxlToQ(px, py);
+            float phi_deg = geom_assumed.Phi_rad(px, py) * 180.0f / static_cast<float>(PI);
+            if (phi_deg < 0.0f)
+                phi_deg += 360.0f;
+            const uint16_t bin = settings.GetBin(q_obs, phi_deg);
+            if (bin == UINT16_MAX)
+                continue;
+            // Lay a narrow peak over the neighbouring q bins of this azimuthal row.
+            const int q_bin = bin % q_bins, phi_bin = bin / q_bins;
+            for (int k = -3; k <= 3; ++k) {
+                const int b = q_bin + k;
+                if (b < 0 || b >= q_bins)
+                    continue;
+                profile[static_cast<size_t>(phi_bin) * q_bins + b] +=
+                        2000.0f * std::exp(-0.5f * static_cast<float>(k * k) / (1.2f * 1.2f));
+            }
+        }
+    }
+    return profile;
+}
+
+}   // namespace
+
+// The measurement this is for: a powder ring is a conic centred on the beam, so a wrong beam centre
+// makes its apparent radius oscillate once per turn. Recovering the centre from that needs neither the
+// calibrant's lattice constant nor the detector distance - only that the ring be round.
+TEST_CASE("RingsFromProfile_RecoversBeamCenter", "[DetGeomCalib]") {
+    DiffractionExperiment x(DetJF4M());
+    x.QSpacingForAzimInt_recipA(0.004).QRangeForAzimInt_recipA(0.5, 4.0);
+    auto azint = x.GetAzimuthalIntegrationSettings();
+    azint.AzimuthalBinCount(32);
+    x.ImportAzimuthalIntegrationSettings(azint);
+
+    PixelMask pixel_mask(x);
+    AzimuthalIntegrationMapping mapping(x, pixel_mask);
+
+    const DiffractionGeometry geom_assumed = x.GetDiffractionGeometry();
+    DiffractionGeometry geom_true = geom_assumed;
+    geom_true.BeamX_pxl(geom_assumed.GetBeamX_pxl() + 6.0f)
+             .BeamY_pxl(geom_assumed.GetBeamY_pxl() - 4.0f);
+
+    const auto profile = SynthesiseProfile(mapping, geom_assumed, geom_true);
+    const auto rings = RingsFromAzimuthalProfile(profile, mapping, geom_assumed, LAB6);
+
+    // Several rings, sampled all the way round: without azimuthal coverage there is no centre to find.
+    REQUIRE(rings.size() > 64);
+
+    RingOptimizer optimizer(geom_assumed);
+    const auto fitted = optimizer.Run(rings);
+
+    CHECK(fitted.GetBeamX_pxl() == Catch::Approx(geom_true.GetBeamX_pxl()).margin(0.5));
+    CHECK(fitted.GetBeamY_pxl() == Catch::Approx(geom_true.GetBeamY_pxl()).margin(0.5));
+    // The starting point was wrong by 6 and 4 pixels, so a fit that did nothing would fail the above -
+    // but check explicitly that it moved toward the truth rather than merely landing near it.
+    CHECK(std::abs(fitted.GetBeamX_pxl() - geom_true.GetBeamX_pxl())
+          < std::abs(geom_assumed.GetBeamX_pxl() - geom_true.GetBeamX_pxl()));
+}
+
+// One azimuthal bin is a plain radial profile: the ring has been averaged over every direction, so
+// nothing is left to say where its centre is. Refuse rather than return points that cannot constrain it.
+TEST_CASE("RingsFromProfile_NeedsAzimuthalBins", "[DetGeomCalib]") {
+    DiffractionExperiment x(DetJF4M());
+    x.QSpacingForAzimInt_recipA(0.004).QRangeForAzimInt_recipA(0.5, 4.0);
+
+    PixelMask pixel_mask(x);
+    AzimuthalIntegrationMapping mapping(x, pixel_mask);
+    REQUIRE(mapping.GetAzimuthalBinCount() == 1);
+
+    const std::vector<float> profile(static_cast<size_t>(mapping.GetQBinCount()), 1000.0f);
+    CHECK(RingsFromAzimuthalProfile(profile, mapping, x.GetDiffractionGeometry(), LAB6).empty());
+}
+
+// A profile with no rings in it must yield no ring points: the peak has to stand clear of the scatter
+// of the background either side of it, or every azimuthal sector would contribute its largest noise
+// excursion as though it were a measurement.
+TEST_CASE("RingsFromProfile_FlatProfileGivesNothing", "[DetGeomCalib]") {
+    DiffractionExperiment x(DetJF4M());
+    x.QSpacingForAzimInt_recipA(0.004).QRangeForAzimInt_recipA(0.5, 4.0);
+    auto azint = x.GetAzimuthalIntegrationSettings();
+    azint.AzimuthalBinCount(32);
+    x.ImportAzimuthalIntegrationSettings(azint);
+
+    PixelMask pixel_mask(x);
+    AzimuthalIntegrationMapping mapping(x, pixel_mask);
+
+    const std::vector<float> profile(
+            static_cast<size_t>(mapping.GetQBinCount()) * mapping.GetAzimuthalBinCount(), 100.0f);
+    CHECK(RingsFromAzimuthalProfile(profile, mapping, x.GetDiffractionGeometry(), LAB6).empty());
+}