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) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01NfuDvf5ipV3Hi8TiCUKD27

image_analysis/geom_refinement/AssignSpotsToRings.cpp

@@ -296,7 +296,7 @@ std::vector<RingOptimizerInput> AssignSpotsToRings(const DiffractionGeometry &ge
 }
 
 void OptimizeGeometry(DiffractionGeometry &geom, const std::vector<SpotToSave> &v, const std::vector<float> &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);
 }

image_analysis/geom_refinement/AssignSpotsToRings.h

@@ -54,7 +54,7 @@ std::vector<RingClusters> GuessInitialGeometry(DiffractionGeometry &geom,  const
 void GuessGeometry(DiffractionGeometry &geom,  const std::vector<SpotToSave> &v, const std::vector<float> &ring_q,
                    bool refine_tilt = true);
 void OptimizeGeometry(DiffractionGeometry &geom, const std::vector<SpotToSave> &v,  const std::vector<float> &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.

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<RingOptimizerInput> &points) {
+                            const std::vector<RingOptimizerInput> &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<float> &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<SpotToSave> &spots,
@@ -73,9 +77,10 @@ CalibrationResult CalibrateFromSpots(const std::vector<SpotToSave> &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,

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

image_analysis/geom_refinement/RingOptimizer.cpp

@@ -1,6 +1,7 @@
 // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
 // SPDX-License-Identifier: GPL-3.0-only
 
+#include <cmath>
 #include <algorithm>
 
 #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<RingOptimizerInput> &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<const double *> &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<int>(free_blocks.size());
+    const int n = static_cast<int>(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<std::pair<const double *, const double *>> 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<double> cov(static_cast<size_t>(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<double>(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<RingOptimizerInput> &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<RingOptimizerInput> &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<RingOptimizerInput> &in
     // Run optimization
     ceres::Solve(options, &problem, &summary);
 
+    if (unc && summary.IsSolutionUsable()) {
+        std::vector<const double *> free_blocks = {&center_x, &center_y, &distance};
+        if (tilt_free) {
+            free_blocks.push_back(&rot1);
+            free_blocks.push_back(&rot2);
+        }
+        ComputeUncertainty(problem, summary, free_blocks,
+                           &center_x, &center_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.

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<RingOptimizerInput> &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<RingOptimizerInput> &input,
+                            RingFitUncertainty *unc = nullptr);
 };
-
-
-

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);