From fba5435c38854ee3a67b0171c0858eeb03ad98fb Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 31 Aug 2026 16:11:18 +0200 Subject: [PATCH] calibration: report what the ring fit knows about its own answer The ring fit reported the SCATTER of its measurements (rms, and the beam-centre standard error that follows from it) but nothing about how well the fit pinned each parameter. Those two part company exactly where a calibration is worth doubting: as the rings run out, the tilt and the beam centre stop being separable - both displace a ring's radius as cos(phi) and only the way that amplitude scales with radius tells them apart - so the fit can sit tightly on the few points it has while being free to spend tens of pixels of beam centre on a tilt the data do not support. Take the covariance of the converged problem from Ceres and report it. Measured on a LaB6 distance series, the fitted tilt is 50 sigma at 110 mm and 0.1 sigma at 500 mm, where only two rings reach the detector; at 500 mm the fit quotes its own beam centre to +-180 px and its tilt to +-2.9 deg on a 0.35 deg value, and the correlation between them is 1.000. Nothing acts on this yet - it is printed so the next change can gate on it. Ceres returns the bare (J'J)^-1 of an unweighted problem, so it is scaled by chi2 per degree of freedom; that leaves the sigmas in pixels, mm and radians whatever unit the residual is stated in. Cost is one 5x5 SVD per run, below the noise of the surrounding I/O. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01NfuDvf5ipV3Hi8TiCUKD27 --- .../geom_refinement/AssignSpotsToRings.cpp | 4 +- .../geom_refinement/AssignSpotsToRings.h | 2 +- .../geom_refinement/PowderCalibration.cpp | 13 ++- .../geom_refinement/PowderCalibration.h | 6 ++ .../geom_refinement/RingOptimizer.cpp | 90 ++++++++++++++++++- .../geom_refinement/RingOptimizer.h | 36 +++++++- rugnux/rugnux_cli.cpp | 22 +++++ 7 files changed, 160 insertions(+), 13 deletions(-) diff --git a/image_analysis/geom_refinement/AssignSpotsToRings.cpp b/image_analysis/geom_refinement/AssignSpotsToRings.cpp index a99e669a7..d6152385c 100644 --- a/image_analysis/geom_refinement/AssignSpotsToRings.cpp +++ b/image_analysis/geom_refinement/AssignSpotsToRings.cpp @@ -296,7 +296,7 @@ std::vector AssignSpotsToRings(const DiffractionGeometry &ge } void OptimizeGeometry(DiffractionGeometry &geom, const std::vector &v, const std::vector &ring_q, - bool refine_tilt) { + bool refine_tilt, RingFitUncertainty *unc) { RingOptimizer optimizer(geom, refine_tilt); - geom = optimizer.Run(AssignSpotsToRings(geom, v, ring_q)); + geom = optimizer.Run(AssignSpotsToRings(geom, v, ring_q), unc); } diff --git a/image_analysis/geom_refinement/AssignSpotsToRings.h b/image_analysis/geom_refinement/AssignSpotsToRings.h index 0dc5dd6dd..4a806107e 100644 --- a/image_analysis/geom_refinement/AssignSpotsToRings.h +++ b/image_analysis/geom_refinement/AssignSpotsToRings.h @@ -54,7 +54,7 @@ std::vector GuessInitialGeometry(DiffractionGeometry &geom, const void GuessGeometry(DiffractionGeometry &geom, const std::vector &v, const std::vector &ring_q, bool refine_tilt = true); void OptimizeGeometry(DiffractionGeometry &geom, const std::vector &v, const std::vector &ring_q, - bool refine_tilt = true); + bool refine_tilt = true, RingFitUncertainty *unc = nullptr); // 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. diff --git a/image_analysis/geom_refinement/PowderCalibration.cpp b/image_analysis/geom_refinement/PowderCalibration.cpp index fb2873330..67523d77c 100644 --- a/image_analysis/geom_refinement/PowderCalibration.cpp +++ b/image_analysis/geom_refinement/PowderCalibration.cpp @@ -24,9 +24,11 @@ namespace { // points of scatter s leaves the textbook var = 2 s^2 / n on each of dx and dy. That is the number that // separates a beam centre that was measured from one that was merely reported. CalibrationResult Summarize(const DiffractionGeometry &fitted, - const std::vector &points) { + const std::vector &points, + const RingFitUncertainty &unc) { CalibrationResult result; result.geometry = fitted; + result.uncertainty = unc; const float cx = fitted.GetBeamX_pxl(); const float cy = fitted.GetBeamY_pxl(); @@ -62,7 +64,9 @@ CalibrationResult CalibrateFromProfile(const std::vector &profile, if (points.empty()) throw JFJochException(JFJochExceptionCategory::CalibrationError, "No powder ring found in the summed azimuthal profile"); - return Summarize(RingOptimizer(geom, refine_tilt).Run(points), points); + RingFitUncertainty unc; + const auto fitted = RingOptimizer(geom, refine_tilt).Run(points, &unc); + return Summarize(fitted, points, unc); } CalibrationResult CalibrateFromSpots(const std::vector &spots, @@ -73,9 +77,10 @@ CalibrationResult CalibrateFromSpots(const std::vector &spots, // From scratch (Hough circle centre + ring clustering), then refined: the guess pins the centre to a // whole pixel and only sees the spots its clustering kept, so the refine re-matches every spot at // that geometry. + RingFitUncertainty unc; GuessGeometry(fitted, spots, calibrant_ring_q, refine_tilt); - OptimizeGeometry(fitted, spots, calibrant_ring_q, refine_tilt); - return Summarize(fitted, AssignSpotsToRings(fitted, spots, calibrant_ring_q)); + OptimizeGeometry(fitted, spots, calibrant_ring_q, refine_tilt, &unc); + return Summarize(fitted, AssignSpotsToRings(fitted, spots, calibrant_ring_q), unc); } void WritePoniFile(const std::string &path, const DiffractionExperiment &experiment, diff --git a/image_analysis/geom_refinement/PowderCalibration.h b/image_analysis/geom_refinement/PowderCalibration.h index 55401b778..f612e23cb 100644 --- a/image_analysis/geom_refinement/PowderCalibration.h +++ b/image_analysis/geom_refinement/PowderCalibration.h @@ -10,6 +10,7 @@ #include "../../common/DiffractionExperiment.h" #include "../../common/DiffractionGeometry.h" #include "../../common/SpotToSave.h" +#include "RingOptimizer.h" // RingFitUncertainty // How the powder rings the detector geometry is fitted to are measured (rugnux --calibration). enum class CalibrationMethod { @@ -25,6 +26,11 @@ struct CalibrationResult { // says so here, and that is the only warning a user gets. double rms_radial_pxl = 0.0; double beam_sigma_pxl = 0.0; + // What the fit itself says about how well each parameter is determined, and how badly the tilt is + // correlated with the beam centre. rms/beam_sigma above describe the SCATTER of the measurements; + // this describes the FIT, and the two part company exactly where it matters - a two-ring tilt can + // leave a small rms while being free to move tens of pixels of beam centre with it. + RingFitUncertainty uncertainty; }; // Both fits take the detector tilt as a free parameter unless refine_tilt is false, which holds diff --git a/image_analysis/geom_refinement/RingOptimizer.cpp b/image_analysis/geom_refinement/RingOptimizer.cpp index 733cddc31..097a8e7d5 100644 --- a/image_analysis/geom_refinement/RingOptimizer.cpp +++ b/image_analysis/geom_refinement/RingOptimizer.cpp @@ -1,6 +1,7 @@ // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only +#include #include #include "../../common/DetectorOrientation.h" @@ -72,7 +73,81 @@ struct RingResidual { RingOptimizer::RingOptimizer(const DiffractionGeometry& geom, bool refine_tilt) : reference(geom), refine_tilt(refine_tilt) {} -DiffractionGeometry RingOptimizer::Run(const std::vector &input) { +namespace { + +// Covariance of the converged fit, scaled to the residual scatter this fit actually left. Ceres hands +// back the bare (J^T J)^-1 of an unweighted problem, which carries the shape of the correlations but +// not their size; multiplying by chi2 per degree of freedom is what turns it into a sigma. The +// residual is in q^2, so the scale is in q^2 too - and it cancels, because (J^T J)^-1 is in +// parameter^2 per residual^2. The sigmas therefore come out in pixels, mm and radians whatever the +// residual is measured in, which is why this does not need to wait on the residual being restated. +void ComputeUncertainty(ceres::Problem &problem, const ceres::Solver::Summary &summary, + const std::vector &free_blocks, + double *const centre_x, double *const centre_y, double *const distance, + double *const rot1, double *const rot2, + RingFitUncertainty &unc) { + const int p = static_cast(free_blocks.size()); + const int n = static_cast(summary.num_residuals); + unc.free_parameters = p; + if (n <= p) + return; + + // DENSE_SVD rather than the SPARSE_QR default: five parameters is not a sparse problem, and SVD is + // the one that survives a near-degenerate normal matrix long enough to report how degenerate it is. + ceres::Covariance::Options options; + options.algorithm_type = ceres::DENSE_SVD; + options.num_threads = 1; + ceres::Covariance covariance(options); + + std::vector> pairs; + for (const double *a : free_blocks) + for (const double *b : free_blocks) + pairs.emplace_back(a, b); + // A failure here is not an error to report upwards: it means the problem is rank-deficient at the + // solution, i.e. some direction in parameter space costs the fit nothing at all. valid = false + // carries that, and it is the strongest statement this struct can make. + if (!covariance.Compute(pairs, &problem)) + return; + + std::vector cov(static_cast(p) * p); + if (!covariance.GetCovarianceMatrix(free_blocks, cov.data())) + return; + + // 2 * final_cost is sum of squared residuals - Ceres' cost is half of it. + unc.chi2_per_dof = 2.0 * summary.final_cost / static_cast(n - p); + for (auto &v : cov) + v *= unc.chi2_per_dof; + + const auto index_of = [&](const double *block) { + for (int i = 0; i < p; ++i) + if (free_blocks[i] == block) return i; + return -1; + }; + const auto sigma = [&](const double *block) { + const int i = index_of(block); + return i < 0 ? 0.0 : std::sqrt(std::max(0.0, cov[i * p + i])); + }; + const auto corr = [&](const double *a, const double *b) { + const int i = index_of(a), j = index_of(b); + if (i < 0 || j < 0) return 0.0; + const double d = std::sqrt(cov[i * p + i] * cov[j * p + j]); + return d > 0.0 ? cov[i * p + j] / d : 0.0; + }; + + unc.sigma_beam_x_pxl = sigma(centre_x); + unc.sigma_beam_y_pxl = sigma(centre_y); + unc.sigma_distance_mm = sigma(distance); + unc.sigma_rot1_rad = sigma(rot1); + unc.sigma_rot2_rad = sigma(rot2); + unc.corr_beam_x_rot1 = corr(centre_x, rot1); + unc.corr_beam_y_rot2 = corr(centre_y, rot2); + unc.valid = true; +} + +} // namespace + +DiffractionGeometry RingOptimizer::Run(const std::vector &input, + RingFitUncertainty *unc) { // Initial guess for the parameters double center_x = reference.GetBeamX_pxl(); double center_y = reference.GetBeamY_pxl(); @@ -108,7 +183,8 @@ DiffractionGeometry RingOptimizer::Run(const std::vector &in && std::all_of(input.begin(), input.end(), [&](const RingOptimizerInput &p) { return p.q_expected == input.front().q_expected; }); - if (single_ring || !refine_tilt) { + const bool tilt_free = !(single_ring || !refine_tilt); + if (!tilt_free) { problem.SetParameterBlockConstant(&rot1); problem.SetParameterBlockConstant(&rot2); } @@ -125,6 +201,16 @@ DiffractionGeometry RingOptimizer::Run(const std::vector &in // Run optimization ceres::Solve(options, &problem, &summary); + if (unc && summary.IsSolutionUsable()) { + std::vector free_blocks = {¢er_x, ¢er_y, &distance}; + if (tilt_free) { + free_blocks.push_back(&rot1); + free_blocks.push_back(&rot2); + } + ComputeUncertainty(problem, summary, free_blocks, + ¢er_x, ¢er_y, &distance, &rot1, &rot2, *unc); + } + // A failed fit must not move the detector geometry. Both callers assign the result straight // back over the geometry they passed in, so handing back the reference leaves the calibration // where it was instead of committing a diverged beam centre and distance. diff --git a/image_analysis/geom_refinement/RingOptimizer.h b/image_analysis/geom_refinement/RingOptimizer.h index 927cfd079..af134c200 100644 --- a/image_analysis/geom_refinement/RingOptimizer.h +++ b/image_analysis/geom_refinement/RingOptimizer.h @@ -12,6 +12,34 @@ struct RingOptimizerInput { double q_expected; }; +// What the ring fit knows about its own answer, from the covariance of the converged problem. +// +// The beam centre and the tilt are not independent: both displace a ring's radius as cos(phi), and +// only how that amplitude scales with the ring's radius tells them apart, which takes two well-sampled +// rings. On one ring they are exactly degenerate; on two they are merely badly correlated, and the fit +// still returns an answer - it just spends tens of pixels of beam centre to buy a tilt the data cannot +// support. The sigmas below say when that has happened and the correlations say why. +// +// Sigmas are in each parameter's own unit and are scaled by the residual variance of THIS fit +// (chi2_per_dof), so they are the usual "how far could this parameter move before the fit got visibly +// worse" and not Ceres' bare (J^T J)^-1. A tilt held fixed reports sigma 0 - it was not a parameter. +struct RingFitUncertainty { + bool valid = false; // false if the covariance could not be computed - itself a verdict: + // the problem is exactly rank-deficient at the solution + double sigma_beam_x_pxl = 0.0; + double sigma_beam_y_pxl = 0.0; + double sigma_distance_mm = 0.0; + double sigma_rot1_rad = 0.0; + double sigma_rot2_rad = 0.0; + // The two degenerate pairs. rot1 tips the detector about the vertical, so it trades against the + // beam centre in x; rot2 tips it about the horizontal and trades against y. |corr| approaching 1 + // is the signature of a tilt the rings do not separate from a beam-centre shift. + double corr_beam_x_rot1 = 0.0; + double corr_beam_y_rot2 = 0.0; + double chi2_per_dof = 0.0; // in the residual's own (q^2) units - a scale, not a goodness of fit + int free_parameters = 0; +}; + class RingOptimizer { DiffractionGeometry reference; bool refine_tilt; @@ -21,8 +49,8 @@ public: // accepts one - XDS has no place to put it - so a calibration meant for such a program is better // measured with the tilt pinned than with it refined and then dropped. RingOptimizer(const DiffractionGeometry& geom, bool refine_tilt = true); - DiffractionGeometry Run(const std::vector &input); + // unc, when given, receives the covariance of the converged fit. Computing it is a 5x5 SVD and + // costs nothing next to the solve, so there is no option to switch it off - pass nullptr instead. + DiffractionGeometry Run(const std::vector &input, + RingFitUncertainty *unc = nullptr); }; - - - diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 4fd15cffe..3e3c661be 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -2054,6 +2054,28 @@ 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; + // 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 + // trade tens of pixels of beam centre for itself. The correlations name that trade - rot1 + // against the beam in x, rot2 against y - and approach 1 as the two stop being separable. + if (const auto &u = cal.uncertainty; u.valid) { + std::cout << fmt::format("Fit sigma: PONI {:.3f}, {:.3f} px distance {:.4f} mm", + u.sigma_beam_x_pxl, u.sigma_beam_y_pxl, u.sigma_distance_mm) + << std::endl; + if (u.sigma_rot1_rad > 0.0 || u.sigma_rot2_rad > 0.0) + std::cout << fmt::format(" rot1 {:.4f} deg, rot2 {:.4f} deg " + "correlation with PONI {:+.3f} / {:+.3f}", + u.sigma_rot1_rad * RAD_TO_DEG, u.sigma_rot2_rad * RAD_TO_DEG, + u.corr_beam_x_rot1, u.corr_beam_y_rot2) + << std::endl; + else + std::cout << " tilt held fixed" << std::endl; + } else { + std::cout << "Fit sigma: not available - the fit is degenerate at its solution" + << std::endl; + } + const std::string poni_path = output_prefix + ".poni"; try { WritePoniFile(poni_path, experiment, g);