From 5a80d2df53bfea59733feee4d88c739677afdb04 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 31 Aug 2026 16:18:27 +0200 Subject: [PATCH] calibration: fix four ways the powder fit quietly loses its input None of these changes the answer on a well-separated cubic standard - the LaB6 distance series is bit-identical by both methods - but each one is a case where input is dropped or mis-assigned without saying so. The circumcentre vote grid was a fixed 4000x4000 box, and the caller never passed anything else. That allocated 128 MB whatever the detector, and on a detector larger than 4000 px in either direction it put the beam centre outside the grid, so every vote was discarded and the guess failed with "Beam center not found". Span the spots' own bounding box instead: a powder ring encloses its centre, so that is where the answer has to be. uint32 votes while there - the most any bin can take is C(500,3). Spots were assigned to the FIRST calibrant ring within a fixed 0.1 1/A, not the nearest. Silver behenate's orders sit 0.108 1/A apart and hexagonal ice has three rings inside 0.06, so for those two standards the window reaches the neighbour and every point lands on the lower-q ring of the pair, biasing the distance. Take the nearest ring, and clamp the window to half the gap to the neighbour - which is what the profile path already did inline, now shared as RingMatchWindow and covered by a test that checks it actually narrows on the crowded standards and not on LaB6. A profile bin no pixel fell in is NaN. SectorPeakQ dropped such a sector by accident, through NaN comparisons falling false; check the four background bins and return explicitly. Ice-ring handling is switched off in calibration mode. Flagged spots are sorted last by the spot budget and so discarded first, which for --calibrant ice throws away exactly what is being calibrated on. The two per-ring std::cout lines in GuessGeometry are gone: a library has no business writing to a terminal, and constructing a Logger to keep them would emit a version banner from inside a fit. What matched belongs in the result struct, which the quality gating still to come needs anyway. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01NfuDvf5ipV3Hi8TiCUKD27 --- .../geom_refinement/AssignSpotsToRings.cpp | 88 ++++++++++++------- .../geom_refinement/AssignSpotsToRings.h | 20 ++++- .../geom_refinement/RingsFromProfile.cpp | 16 ++-- rugnux/rugnux_cli.cpp | 6 ++ tests/CalibrationTest.cpp | 40 +++++++++ 5 files changed, 132 insertions(+), 38 deletions(-) diff --git a/image_analysis/geom_refinement/AssignSpotsToRings.cpp b/image_analysis/geom_refinement/AssignSpotsToRings.cpp index d6152385c..2bc6a0c87 100644 --- a/image_analysis/geom_refinement/AssignSpotsToRings.cpp +++ b/image_analysis/geom_refinement/AssignSpotsToRings.cpp @@ -4,23 +4,42 @@ #include "../../common/JFJochMath.h" #include "AssignSpotsToRings.h" +#include +#include #include #include #include #include -#include #include "../../common/CrystalLattice.h" -FindCircleCenterResult FindCircleCenter(const std::vector &v, int64_t width, int64_t height, int64_t max_spots) { - if ((width <= 0) || (height <= 0) || (max_spots <= 0)) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Invalid image size"); - - std::vector vote(width * height, 0); +FindCircleCenterResult FindCircleCenter(const std::vector &v, int64_t max_spots) { + if (max_spots <= 0) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Invalid spot limit"); + if (v.size() < 3) + return {.total_votes = 0, .votes_for_beam_center = 0, .x = 0.0f, .y = 0.0f}; // Limit to only first 250 spots, given the algorithm is N^3 auto task_size = std::min(v.size(), max_spots); + // The vote grid, one bin per pixel over the spots' bounding box (see the header). + int64_t x0 = std::numeric_limits::max(), x1 = std::numeric_limits::min(); + int64_t y0 = x0, y1 = x1; + for (size_t i = 0; i < task_size; i++) { + x0 = std::min(x0, std::floor(v[i].x)); + x1 = std::max(x1, std::ceil(v[i].x)); + y0 = std::min(y0, std::floor(v[i].y)); + y1 = std::max(y1, std::ceil(v[i].y)); + } + const int64_t width = x1 - x0 + 1; + const int64_t height = y1 - y0 + 1; + if ((width <= 0) || (height <= 0)) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Invalid image size"); + + // uint32_t, not int64_t: the most votes any bin can take is C(max_spots,3), which for the 500-spot + // limit is 2.1e7 - three orders below what a uint32 holds - and it halves the allocation. + std::vector vote(width * height, 0); + for (int i = 0; i < task_size; i++) { for (int j = i+1; j < task_size; j++) { for (int k = j+1; k < task_size; k++) { @@ -36,8 +55,8 @@ FindCircleCenterResult FindCircleCenter(const std::vector &v, int64_ float bx = x1_sq * (v[j].y - v[k].y) + x2_sq * (v[k].y - v[i].y) + x3_sq * (v[i].y - v[j].y); float by = x1_sq * (v[k].x - v[j].x) + x2_sq * (v[i].x - v[k].x) + x3_sq * (v[j].x - v[i].x); - int64_t cx = std::lround(bx / (2.0f * a)); - int64_t cy = std::lround(by / (2.0f * a)); + const int64_t cx = std::lround(bx / (2.0f * a)) - x0; + const int64_t cy = std::lround(by / (2.0f * a)) - y0; if ((cx >= 0) && (cx < width) && (cy >= 0) && (cy < height)) vote[cx + cy * width]++; @@ -63,8 +82,8 @@ FindCircleCenterResult FindCircleCenter(const std::vector &v, int64_ return { .total_votes = total_votes, .votes_for_beam_center = max_votes, - .x = static_cast(cx), - .y = static_cast(cy) + .x = static_cast(cx + x0), + .y = static_cast(cy + y0) }; } @@ -250,27 +269,22 @@ void GuessGeometry(DiffractionGeometry &geom, const std::vector &v, int ring_idx = 0; int cluster_idx = 0; + // Walk the observed clusters and the calibrant's rings together, in ascending q. Which of them + // paired up used to be printed to stdout from here; it is diagnostics a caller wants back as data, + // not something a library should write to a terminal. while (cluster_idx < cluster_annot.size() && ring_idx < ring_q.size()) { - float obs_q = 2 * PI / geom.PxlToRes(cluster_annot[cluster_idx].R_obs); - if (std::fabs(ring_q[ring_idx] - obs_q) < 0.1) { - std::cout << "Found ring " << ring_idx - << " with q_ring = " << ring_q[ring_idx] - << " and q_obs = " << obs_q - << " diff = " << std::fabs(ring_q[ring_idx] - obs_q) << std::endl; + const float obs_q = 2 * PI / geom.PxlToRes(cluster_annot[cluster_idx].R_obs); + if (std::fabs(ring_q[ring_idx] - obs_q) + < RingMatchWindow(ring_q, ring_idx, RING_MATCH_Q_RECIPA)) { for (const auto &spot: cluster_annot[cluster_idx].spots) optimizer_input.push_back({v[spot].x, v[spot].y, ring_q[ring_idx]}); ring_idx++; cluster_idx++; + } else if (ring_q[ring_idx] < obs_q) { + ring_idx++; } else { - std::cout << "Cannot match " << ring_idx - << " with q_ring = " << ring_q[ring_idx] - << " and q_obs = " << obs_q - << " diff = " << std::fabs(ring_q[ring_idx] - obs_q) << std::endl; - if (ring_q[ring_idx] < obs_q) - ring_idx++; - else - cluster_idx++; + cluster_idx++; } } @@ -278,19 +292,33 @@ void GuessGeometry(DiffractionGeometry &geom, const std::vector &v, geom = optimizer.Run(optimizer_input); } +float RingMatchWindow(const std::vector &ring_q, size_t i, float max_window) { + float window = max_window; + if (i > 0) + window = std::min(window, 0.5f * (ring_q[i] - ring_q[i - 1])); + if (i + 1 < ring_q.size()) + window = std::min(window, 0.5f * (ring_q[i + 1] - ring_q[i])); + return window; +} + std::vector AssignSpotsToRings(const DiffractionGeometry &geom, const std::vector &v, const std::vector &ring_q) { std::vector optimizer_input; for (const auto& s: v) { - float q_obs = 2 * PI / geom.PxlToRes(s.x, s.y); - for (const auto &q : ring_q) { - if (std::fabs(q - q_obs) < 0.1) { - optimizer_input.push_back({s.x, s.y, q}); - break; - } + const float q_obs = 2 * PI / geom.PxlToRes(s.x, s.y); + // The NEAREST ring, then the window - not the first ring within a fixed window. Where two rings + // are closer together than that window, taking the first match assigns both of them to the + // lower-q one and biases the distance it fits. + size_t nearest = 0; + for (size_t i = 1; i < ring_q.size(); i++) { + if (std::fabs(ring_q[i] - q_obs) < std::fabs(ring_q[nearest] - q_obs)) + nearest = i; } + if (!ring_q.empty() + && std::fabs(ring_q[nearest] - q_obs) < RingMatchWindow(ring_q, nearest, RING_MATCH_Q_RECIPA)) + optimizer_input.push_back({s.x, s.y, ring_q[nearest]}); } return optimizer_input; } diff --git a/image_analysis/geom_refinement/AssignSpotsToRings.h b/image_analysis/geom_refinement/AssignSpotsToRings.h index 4a806107e..69cd79ef0 100644 --- a/image_analysis/geom_refinement/AssignSpotsToRings.h +++ b/image_analysis/geom_refinement/AssignSpotsToRings.h @@ -21,9 +21,12 @@ struct FindCircleCenterResult { float y; }; +// Beam centre from the circumcentres of spot triples, voted into a pixel grid. The grid spans the +// spots' own bounding box: a powder ring encloses the beam centre, so that is where the answer has to +// be, and it is the only bound the spots themselves justify. It used to be a fixed 4000x4000 box, +// which both allocated 128 MB whatever the detector and put the centre outside the grid - so every +// vote was discarded - on any detector larger than that in either direction. FindCircleCenterResult FindCircleCenter(const std::vector &v, - int64_t width = 4000, - int64_t height = 4000, int64_t max_spots = 500); std::vector> ClusterSpotsIntoRings(const std::vector& r, float eps = 1.5, int minPts = 5); @@ -56,8 +59,19 @@ void GuessGeometry(DiffractionGeometry &geom, const std::vector &v, void OptimizeGeometry(DiffractionGeometry &geom, const std::vector &v, const std::vector &ring_q, bool refine_tilt = true, RingFitUncertainty *unc = nullptr); +// Widest a ring's match window is allowed to be, in q. The clamp below narrows it where the rings +// crowd; this is the value for a well-separated ring. +constexpr float RING_MATCH_Q_RECIPA = 0.1f; + +// Half-width of the window ring i may be matched in, clamped so it can never reach the neighbouring +// ring. Silver behenate's orders sit 0.108 1/A apart and hexagonal ice has three rings within 0.06, so +// a fixed window takes in the neighbour's flank - and a search that stops at its first match then puts +// every point on the lower-q ring of the pair. +float RingMatchWindow(const std::vector &ring_q, size_t i, float max_window); + // Each spot paired with the calibrant ring nearest its observed q, as the points RingOptimizer fits. -// Spots more than 0.1 1/A from every ring are dropped rather than forced onto the closest one. +// Spots further from every ring than RingMatchWindow allows are dropped rather than forced onto the +// closest one. std::vector AssignSpotsToRings(const DiffractionGeometry &geom, const std::vector &v, const std::vector &ring_q); diff --git a/image_analysis/geom_refinement/RingsFromProfile.cpp b/image_analysis/geom_refinement/RingsFromProfile.cpp index b33f26d1c..9402c357c 100644 --- a/image_analysis/geom_refinement/RingsFromProfile.cpp +++ b/image_analysis/geom_refinement/RingsFromProfile.cpp @@ -5,6 +5,7 @@ #include #include "RingsFromProfile.h" +#include "AssignSpotsToRings.h" // RingMatchWindow #include "../../common/JFJochMath.h" namespace { @@ -22,6 +23,15 @@ float SectorPeakQ(const std::vector &profile, int32_t q_bins, int phi_bin const auto value = [&](int i) { return profile[row + static_cast(i)]; }; const auto q_of = [&](int i) { return low_q + (static_cast(i) + 0.5f) * q_spacing; }; + // A bin no pixel fell in is NaN, not zero (AzimuthalIntegrationProfile::GetResult), and the four + // background bins are where a module gap or the beam stop shows up first. Say so rather than + // relying on NaN comparisons to fail the peak test further down: a sector whose background cannot + // be measured has no measurable peak either. + for (int i : {lo_bin, lo_bin + 1, hi_bin - 1, hi_bin}) { + if (!std::isfinite(value(i))) + return NAN; + } + 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) { @@ -99,11 +109,7 @@ std::vector RingsFromAzimuthalProfile(const std::vector 0) - window = std::min(window, 0.5f * (q_ring - calibrant_ring_q[i - 1])); - if (i + 1 < calibrant_ring_q.size()) - window = std::min(window, 0.5f * (calibrant_ring_q[i + 1] - q_ring)); + const float window = RingMatchWindow(calibrant_ring_q, i, q_window_recipA); if (!(q_ring - window > low_q) || !(q_ring + window < high_q)) continue; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 3e3c661be..dd3465172 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1988,6 +1988,12 @@ static int RunRugnux(int argc, char **argv) { config.fit_spindle = fit_spindle; config.write_process_h5 = false; // the .poni below is the output of this mode + // Ice-ring handling has no place in a calibration run and does active harm. Flagged spots are + // sorted LAST by the spot budget (FilterSpotsByCount), so they are the first thrown away - which + // for --calibrant ice discards exactly the spots being calibrated on, and for every other + // standard discards whichever of its rings happen to fall in the fixed hexagonal bands. + experiment.DetectIceRings(false); + // Spot finding for --calibration spots. Indexing is off: a calibration wants the spot positions // and nothing else, and the calibrant is a powder with no lattice to index. config.spot_finding.enable = true; diff --git a/tests/CalibrationTest.cpp b/tests/CalibrationTest.cpp index 2226b7a37..1f6bda6c4 100644 --- a/tests/CalibrationTest.cpp +++ b/tests/CalibrationTest.cpp @@ -10,6 +10,7 @@ #include "../common/Definitions.h" #include "../common/JFJochMath.h" +#include "../image_analysis/geom_refinement/AssignSpotsToRings.h" #include "../image_analysis/geom_refinement/Calibrants.h" #include "../image_analysis/geom_refinement/PowderCalibration.h" @@ -128,3 +129,42 @@ TEST_CASE("Calibration_PoniFileAxisConvention", "[DetGeomCalib]") { + std::to_string(x.GetXPixelsNumConv()) + "]"; CHECK(keys["Detector_config"].find(shape) != std::string::npos); } + +// The match window may never reach the neighbouring ring, for any calibrant. Where two rings are closer +// together than twice the nominal window, a fixed window takes in the neighbour's flank - which the +// rings path reads as this ring's background, and which the spots path (before it took the NEAREST ring) +// resolved by assigning both to the lower-q one. +TEST_CASE("Calibration_RingMatchWindowNeverReachesTheNeighbour", "[DetGeomCalib]") { + for (const auto &c : Calibrants()) { + const auto q = CalibrantRings(c.name); + REQUIRE(q.size() > 1); + for (size_t i = 0; i < q.size(); ++i) { + const float w = RingMatchWindow(q, i, RING_MATCH_Q_RECIPA); + CHECK(w <= RING_MATCH_Q_RECIPA); + CHECK(w > 0.0f); + if (i > 0) + CHECK(q[i] - w >= 0.5f * (q[i] + q[i - 1]) - 1e-6f); + if (i + 1 < q.size()) + CHECK(q[i] + w <= 0.5f * (q[i] + q[i + 1]) + 1e-6f); + } + } +} + +// ...and the clamp is not a no-op. Silver behenate's orders sit about 0.108 1/A apart and hexagonal ice +// has rings inside 0.06, so both are narrowed below the nominal window - while LaB6, whose rings are +// well separated at low q, keeps it. Without a standard that actually crowds, the test above would pass +// on a clamp that never fired. +TEST_CASE("Calibration_CrowdedStandardsNarrowTheWindow", "[DetGeomCalib]") { + auto narrowed = [](const std::string &name) { + const auto q = CalibrantRings(name); + size_t n = 0; + for (size_t i = 0; i < q.size(); ++i) + if (RingMatchWindow(q, i, RING_MATCH_Q_RECIPA) < RING_MATCH_Q_RECIPA) ++n; + return n; + }; + CHECK(narrowed("agbh") > 0); + CHECK(narrowed("ice") > 0); + // The innermost LaB6 rings are more than 0.2 1/A apart, so nothing narrows them. + const auto lab6 = CalibrantRings("lab6"); + CHECK(RingMatchWindow(lab6, 0, RING_MATCH_Q_RECIPA) == Catch::Approx(RING_MATCH_Q_RECIPA)); +}