From 49dd01331b85f8649f9964a065e6d5210e5ec791 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 31 Aug 2026 17:58:26 +0200 Subject: [PATCH] calibration: let the spots vote on the geometry in ring mode too The rings and the spots fail in different regimes, which is the whole reason to carry both. The header is right on a well-configured instrument and is the thing a calibration is run to check. The profile's own estimates are exact while the error stays small and stop meaning anything beyond that - past about half a ring spacing each ring reaches the azimuthally averaged profile as two horns rather than one peak, and the distance search reads a list of horns as a list of rings. The circle through the spots reads nothing from the header at all: measured on a 110 mm LaB6 exposure, --calibration spots returns the same geometry from a header 400 px out in the centre AND eight times out in distance. So ring mode now finds spots as well - half a second - and offers what they make of the geometry as one more starting hypothesis, fitted like the others with the residual left to choose. It is added whole rather than as a centre alone: GuessGeometry votes for the circle centre, clusters the radii into rings and takes the distance from the innermost, and those two belong together. Taking only its centre would not have helped, because the distance candidates come from a profile averaged about the header's centre, and a ring smeared over hundreds of pixels cannot be un-smeared by reading its bins differently. Measured on the 110 mm exposure, whose true PONI is 765.90 at 110.03 mm. A header 40 px out in the centre now lands within 0.6 px (it landed 41 px away before). The cases with BOTH wrong, which failed before this and equally before the beam centre work, now come out: 250 mm / 780 px gives 110.041 mm and 766.10 px against 42.3 mm and 782.7; 250 mm / 867 px gives 110.065 and 765.48. All five datasets are unchanged from their correct headers and the distance still recovers from any header between 25 and 1200 mm. Past about a hundred pixels nothing rescues ring mode, and the reason is the profile rather than the seeding: binned about a centre that far out it shows each ring smeared across its own sectors, so even started from the exact answer the fit leaves rms 4.8 px and drifts. Re-binning would fix it and would need the images read a second time; --calibration spots, which never touches the profile, already covers it. That regime is now visible rather than silent. The spots' beam centre is printed beside the fitted one as a cross-check - two methods sharing no assumption, so a reader can see at a glance whether they agree. It costs nothing, the spots having been found already, and it separates cleanly: 0.3-0.4 px on the good runs against 201.4 px on the exposure whose profile could not be fitted at all. Reported as a fact and not gated on, since at long distance both methods weaken together and the honest thing is to show the number (5.2 px at 300 mm, 12.6 at 500). Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01NfuDvf5ipV3Hi8TiCUKD27 --- .../geom_refinement/PowderCalibration.cpp | 74 +++++++++++++++---- .../geom_refinement/PowderCalibration.h | 22 +++++- rugnux/Rugnux.cpp | 15 ++-- rugnux/rugnux_cli.cpp | 12 +++ 4 files changed, 101 insertions(+), 22 deletions(-) diff --git a/image_analysis/geom_refinement/PowderCalibration.cpp b/image_analysis/geom_refinement/PowderCalibration.cpp index adfbb5192..a9f514ac5 100644 --- a/image_analysis/geom_refinement/PowderCalibration.cpp +++ b/image_analysis/geom_refinement/PowderCalibration.cpp @@ -10,7 +10,7 @@ #include "PowderCalibration.h" #include "../../common/GitInfo.h" #include "../../common/JFJochMath.h" -#include "AssignSpotsToRings.h" +#include "AssignSpotsToRings.h" // FindCircleCenter #include "RingOptimizer.h" #include "RingsFromProfile.h" @@ -73,7 +73,8 @@ CalibrationResult CalibrateFromProfile(const std::vector &profile, const AzimuthalIntegrationMapping &mapping, const DiffractionGeometry &geom, const std::vector &calibrant_ring_q, - bool refine_tilt) { + bool refine_tilt, + const std::vector &spots) { // Where to start. The ring search below is local - each ring is looked for inside a window a few // pixels of radius wide - so a header distance more than a percent or two out puts every ring // outside its own window, and what the fit then converges on is noise. Ask the rings what the @@ -90,15 +91,47 @@ CalibrationResult CalibrateFromProfile(const std::vector &profile, // whatever sectors are left, which is how a 20 px error used to end 31 px wrong. The rings answer // this without a calibrant and without a distance: a powder ring is a conic centred on the beam, so // a wrong centre makes EVERY ring's radius oscillate once per turn by the same amount. - DiffractionGeometry seed_geometry = geom; - auto centre_offset = BeamCentreOffsetFromProfile(profile, mapping, geom, observed); - // Under a pixel it is not a different hypothesis, it is the same one - and fitting it as well would - // double the work for two answers that cannot be told apart. - if (centre_offset && std::hypot(centre_offset->first, centre_offset->second) < 1.0f) - centre_offset.reset(); - if (centre_offset) - seed_geometry.BeamX_pxl(geom.GetBeamX_pxl() + centre_offset->first) - .BeamY_pxl(geom.GetBeamY_pxl() + centre_offset->second); + // Three places a beam centre can come from, and they fail in different regimes - which is the whole + // reason to carry all of them. The header is right on a well-configured instrument and is what a + // calibration is run to check. The rings' own wobble is exact while the error stays under about half + // a ring spacing, and stops meaning anything beyond that, where each ring reaches the azimuthally + // averaged profile as two horns rather than one peak. The circle through the spots reads nothing + // from the header whatsoever, so it holds where both of the others have given up - measured, it + // returns the same answer from a header 400 px out in the centre and eight times out in distance. + // None of them is trusted; each is fitted and the residual chooses. + std::vector centres = {geom}; + const auto add_centre = [&](float x, float y) { + // Within a pixel of one already on the list it is not another hypothesis, it is the same one, + // and fitting it again would only double the work for two answers nothing can tell apart. + for (const auto &c : centres) + if (std::hypot(x - c.GetBeamX_pxl(), y - c.GetBeamY_pxl()) < 1.0f) + return; + DiffractionGeometry candidate = geom; + centres.push_back(candidate.BeamX_pxl(x).BeamY_pxl(y)); + }; + + if (const auto offset = BeamCentreOffsetFromProfile(profile, mapping, geom, observed)) + add_centre(geom.GetBeamX_pxl() + offset->first, geom.GetBeamY_pxl() + offset->second); + + // What the spots make of it, as a COMPLETE geometry rather than only a centre. Taking the circle + // centre alone is not enough: the distance candidates are read off the azimuthally averaged profile, + // which is averaged about the header's centre, and once that is badly wrong the profile no longer + // shows rings at all - a ring smeared over hundreds of pixels cannot be un-smeared by reading its + // bins differently. The spots have no such problem, because a spot's position is a fact about the + // image; GuessGeometry votes for the circle centre, clusters the radii into rings and takes the + // distance from the innermost one, which is exactly the path that makes --calibration spots immune + // to the header. Let it produce one whole starting geometry of its own. + std::optional from_spots; + if (spots.size() >= 3) { + try { + DiffractionGeometry guess = geom; + GuessGeometry(guess, spots, calibrant_ring_q, refine_tilt); + from_spots = guess; + } catch (const std::exception &) { + // No circle, or no ring clusters - the spots simply have nothing to say here. The profile's + // own hypotheses stand on their own, so there is nothing to report and nothing to stop. + } + } // One starting distance, fitted to convergence. A seed only has to land in the fit's basin, not on // the answer: it is measured from blended peaks in the azimuthally-averaged profile and is good to @@ -195,14 +228,17 @@ CalibrationResult CalibrateFromProfile(const std::vector &profile, struct Start { DiffractionGeometry geometry; bool tracked; Provenance provenance; }; std::vector starts; - for (int centre = 0; centre < (centre_offset ? 2 : 1); ++centre) { - const DiffractionGeometry &base = centre == 0 ? geom : seed_geometry; + // The spots' geometry is a start in its own right, not one to cross with the profile's distances - + // its centre and its distance are measured together and belong together. + if (from_spots) + starts.push_back({*from_spots, true, {from_spots->GetDetectorDistance_mm(), 0.0}}); + for (size_t centre = 0; centre < centres.size(); ++centre) { for (size_t i = 0; i <= candidates.size(); ++i) { const bool seeded_distance = i < candidates.size(); - DiffractionGeometry start = base; + DiffractionGeometry start = centres[centre]; if (seeded_distance) start.DetectorDistance_mm(candidates[i].distance_mm); - starts.push_back({start, seeded_distance || centre == 1, + starts.push_back({start, seeded_distance || centre > 0, {seeded_distance ? candidates[i].distance_mm : 0.0f, seeded_distance ? candidates[i].score : 0.0}}); } @@ -275,6 +311,14 @@ CalibrationResult CalibrateFromProfile(const std::vector &profile, } auto result = Summarize(best->geometry, best->points, best->uncertainty); + if (from_spots) { + result.spots_available = true; + result.spots_beam_x_pxl = from_spots->GetBeamX_pxl(); + result.spots_beam_y_pxl = from_spots->GetBeamY_pxl(); + result.spots_disagreement_pxl = + std::hypot(from_spots->GetBeamX_pxl() - best->geometry.GetBeamX_pxl(), + from_spots->GetBeamY_pxl() - best->geometry.GetBeamY_pxl()); + } result.tilt_refined = tilt_refined; result.tilt_significance = significance; result.seed_distance_mm = best_provenance.seed_mm; diff --git a/image_analysis/geom_refinement/PowderCalibration.h b/image_analysis/geom_refinement/PowderCalibration.h index c4d4cf392..16769e5d9 100644 --- a/image_analysis/geom_refinement/PowderCalibration.h +++ b/image_analysis/geom_refinement/PowderCalibration.h @@ -45,6 +45,17 @@ struct CalibrationResult { // the tilt was never a free parameter. bool tilt_refined = false; float tilt_significance = 0.0f; + // How far the beam centre the SPOTS vote for is from the one the rings were fitted to, in pixels, + // and where the spots put it. An independent cross-check that costs nothing, because the spots have + // already been found: the circle through them reads nothing from the header, so it holds in exactly + // the regime the profile does not - a profile binned about a badly wrong centre shows each ring + // smeared across its sectors, and no seeding recovers that. A large disagreement is therefore not a + // close call to arbitrate but a statement that this profile could not have been fitted, whatever + // the residual says. Zero when no spots were available. + float spots_beam_x_pxl = 0.0f; + float spots_beam_y_pxl = 0.0f; + float spots_disagreement_pxl = 0.0f; + bool spots_available = false; }; // How many of its own sigmas the fitted tilt is away from zero - the number the gate above reads. @@ -69,11 +80,20 @@ constexpr float TILT_MIN_SIGNIFICANCE = 3.0f; // given a geometry measured with the tilt pinned than one measured tilted and then flattened. // Fit the geometry to the rings found in a run-summed azimuthal profile (CalibrationMethod::Rings). +// +// spots, where they are given, add one more beam-centre hypothesis: the centre of the circle through +// them, voted for by every triple (FindCircleCenter). It reads nothing from the header at all - only +// where the spots landed - so it is the one estimate that survives a header the profile itself cannot +// correct from. The rings' own estimate stops working once the beam centre is out by more than about +// half a ring spacing, because past that the azimuthally averaged profile shows each ring as the two +// horns of a sinusoid rather than as a ring; the circle through the spots has no such limit. Measured +// on a 110 mm LaB6 exposure it returns the same geometry from a header 400 px and 8x in distance wrong. CalibrationResult CalibrateFromProfile(const std::vector &profile, const AzimuthalIntegrationMapping &mapping, const DiffractionGeometry &geom, const std::vector &calibrant_ring_q, - bool refine_tilt = true); + bool refine_tilt = true, + const std::vector &spots = {}); // Fit the geometry to a pooled spot list (CalibrationMethod::Spots): the beam centre from scratch off // the Hough circle centre, then the same ring fit. diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 90a6e41bb..98be31ba6 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1611,12 +1611,14 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } const bool full = (config_.mode == ProcessMode::FullAnalysis); - // Powder calibration by spots finds spots and nothing else, so it runs the same per-image engine as - // the full analysis (spot finding lives inside MXAnalysisWithoutFPGA) with indexing switched off in - // its settings; by rings it only needs the azimuthal profile, which is the azint worker exactly. + // Powder calibration runs the same per-image engine as the full analysis (spot finding lives inside + // MXAnalysisWithoutFPGA) with indexing switched off in its settings. By rings it needs only the + // azimuthal profile and the azint worker would do - but the spots are worth their half second even + // then, because the circle through them places the beam centre without reference to the header at + // all, which is the one hypothesis that survives a header the profile cannot correct. const bool calibration = (config_.mode == ProcessMode::Calibration); const bool calibration_spots = calibration && config_.calibration_method == CalibrationMethod::Spots; - const bool per_image_analysis = full || calibration_spots; + const bool per_image_analysis = full || calibration; const bool write_files = write_output && !config_.output_prefix.empty(); // Output/runtime invariants. Algorithm settings (indexing, scaling, integration, polarization, @@ -2812,7 +2814,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b msg.run_number = experiment_.GetRunNumber(); msg.run_name = experiment_.GetRunName(); - if (calibration_spots) { + if (calibration) { std::lock_guard lock(calibration_spot_m); for (size_t i = 0; i < msg.spots.size() && i < calibration_spots_per_image; ++i) calibration_spot_list.push_back(msg.spots[i]); @@ -2926,7 +2928,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.calibration = CalibrateFromProfile(summed->GetResult(), mapping, experiment_.GetDiffractionGeometry(), config_.calibrant_ring_q, - config_.calibration_refine_tilt); + config_.calibration_refine_tilt, + calibration_spot_list); } } diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index f9d271808..711dc1393 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -2092,6 +2092,18 @@ static int RunRugnux(int argc, char **argv) { const auto [beam_x, beam_y] = g.GetDirectBeam_pxl(); std::cout << fmt::format("Direct beam: {:.3f}, {:.3f} px", beam_x, beam_y) << std::endl; + // Where the spots put the beam, independently of everything above. The circle through them + // reads nothing from the header, so it holds in exactly the regime the summed profile does not: + // a profile binned about a badly wrong centre shows each ring smeared across its sectors, and no + // amount of seeding recovers that - the fit then converges, prettily, on the wrong geometry. + // Printed as a fact rather than gated on, because the reader can see at a glance whether two + // methods that share no assumption agree. + if (cal.spots_available) + std::cout << fmt::format("Spot cross-check: the circle through the spots puts the beam at " + "{:.1f}, {:.1f} px - {:.1f} px from the fitted centre", + cal.spots_beam_x_pxl, cal.spots_beam_y_pxl, + cal.spots_disagreement_pxl) << std::endl; + // What the fit says about itself. The scatter line above is about the MEASUREMENTS; this is // about the PARAMETERS, and the two disagree exactly where the calibration is worth doubting: // a fit with few rings can sit tightly on the points it has while leaving the tilt free to