Per-image refinement: weight each spot by how strong it is for its resolution

`RefineGeometryIfNeeded` hands XtalOptimizer the WHOLE spot list, not the indexed subset, and the first pass admits anything within 0.3 fractional-Miller units of an integer - which is 11.3% of RANDOMLY placed spots, since the admitted volume is (4/3)*pi*t^3. Every one of them then enters an unweighted L2 fit with an arbitrary rounded index. On images with many detections the refined orientation ends up 2.3-2.8 degrees from the goniometer-consistent one and explains 14 of its own 250 spots where the undragged orientation explains 68; mosaicity and profile radius inherit the error and integration follows. Weight every spot by its intensity divided by the median intensity of its own equal-count resolution shell, applied as w^2 on the squared residual with w^2 = r/(1+r). The shell normalisation is the point: refinement needs the high-resolution spots because they carry the cell and distance, and those are LEGITIMATELY weaker, so a raw intensity weight would suppress exactly the spots the fit depends on. Measured, the weight is resolution-neutral - median exactly 0.707 in every shell, and corr(w, 1/d^2) = -0.20 / -0.11 against -0.32 / -0.34 for the same function of un-normalised intensity. This is a PRIOR: it is computed from the spot alone and never looks at the current residual, so unlike a robust loss it cannot mistake a genuine spot for an outlier while the starting geometry is still far off and leave the fit unable to move. That failure is not hypothetical - a CauchyLoss on this same residual, at the scale the multi-frame GeometryRefiner uses, collapsed one crystal's indexing rate from 99.89% to 19.83% and was rejected. It does not work by telling good spots from bad, and it does not need to. No per-spot property separates spots that index from spots that do not: measured AUC is 0.53 for peak pixel, 0.53 for total intensity, 0.51 for pixel count, 0.45 for peakedness, and a logistic regression on all twelve available features with pairwise interactions reaches only 0.64. What the weight does is halve the EFFECTIVE COUNT of every spot (mean w^2 = 0.517), and the damage scales with the absolute count of unexplained spots in the objective - 80.6 per frame here against 36.8 for the finder that was never damaged. That is also why an empirical `--max-spots 66` cap works while leaving the list no purer than before: it reaches the same operating point by discarding spots. This reaches it without discarding any, and without a tuned constant. Rotation battery, 33 crystals, both spot finders: finder A 29/33 -> 30/33 point groups (one crystal P222 -> P4212 = XDS, its high-shell CC1/2 86.0 -> 98.4) finder B 28/33 -> 29/33 point groups (one crystal I222 -> I23, its high-shell CC1/2 14.8 -> 38.0) No crystal lost its point group in either mode and no run failed. On the meta-stable multi-lattice dataset the CC1/2 spread over four frame ranges falls 19.7 -> 13.1 for finder B, and the indexing rate rises in 8 of 8 configurations. The crystal that the rejected robust loss destroyed keeps its 99.89% indexing rate exactly. The cost, stated plainly: ISa falls by 0.2-1.7 on about five crystals (and rises on two). Point-group correctness is worth more than that - merging in the wrong symmetry cannot be undone from the output, whereas ISa is a quality metric of data that remain correct - but it is a real trade and not a free win. Off by default. The indexers pass a spot list they have already selected, so their calls are unchanged; only the per-image refinement, which gets the raw list, turns it on. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
2026-07-28 09:53:13 +02:00
parent ec7a826136
commit ae126c3d5b
3 changed files with 70 additions and 3 deletions
@@ -2,6 +2,8 @@
 // SPDX-License-Identifier: GPL-3.0-only

 #include "../../common/JFJochMath.h"
+#include <algorithm>
+#include <numeric>
 #include <Eigen/Dense>

 #include "XtalOptimizer.h"
@@ -81,6 +83,45 @@ struct RotationNormRegularizer {
    const double weight;
 };

+// Prior confidence weight per spot: how strong the spot is FOR ITS RESOLUTION. The frame's spots are
+// ordered by resolution and cut into equal-count shells, and each intensity is divided by its shell
+// median. Refinement needs the high-resolution spots (they carry the cell and distance information) and
+// those are legitimately weaker, so a raw intensity weight would suppress exactly the wrong ones; the
+// shell normalisation makes the weight resolution-neutral by construction.
+//
+// The weight enters as w^2 on the squared residual, w^2 = r/(1+r): the shell median contributes half,
+// a 4x-median spot 0.8, a quarter-median spot 0.2. Weak spots still pull, they just do not drive. Unlike
+// a robust loss this is a PRIOR - it never looks at the current residual, so it cannot mistake a genuine
+// spot for an outlier when the starting geometry is far off and leave the fit unable to move.
+static std::vector<double> SpotConfidenceWeights(const std::vector<SpotToSave> &spots) {
+    constexpr size_t spots_per_shell = 32;
+
+    std::vector<size_t> by_res(spots.size());
+    std::iota(by_res.begin(), by_res.end(), 0);
+    std::ranges::sort(by_res, {}, [&](size_t i) { return spots[i].d_A; });
+
+    const size_t nshells = std::max<size_t>(1, spots.size() / spots_per_shell);
+    std::vector<double> weight(spots.size());
+    std::vector<float> shell_intensity;
+
+    for (size_t s = 0; s < nshells; s++) {
+        const size_t begin = s * spots.size() / nshells;
+        const size_t end = (s + 1) * spots.size() / nshells;
+
+        shell_intensity.clear();
+        for (size_t i = begin; i < end; i++)
+            shell_intensity.push_back(spots[by_res[i]].intensity);
+        std::ranges::nth_element(shell_intensity, shell_intensity.begin() + shell_intensity.size() / 2);
+        const double median = std::max(1e-3f, shell_intensity[shell_intensity.size() / 2]);
+
+        for (size_t i = begin; i < end; i++) {
+            const double r = std::max(0.0f, spots[by_res[i]].intensity) / median;
+            weight[by_res[i]] = std::sqrt(r / (1.0 + r));
+        }
+    }
+    return weight;
+}
+
 bool XtalOptimizerInternal(XtalOptimizerData &data,
                           const std::vector<std::vector<SpotToSave>> &spots,
                           const float tolerance,
@@ -144,10 +185,19 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,

        const float tolerance_sq = tolerance * tolerance;

+        // Sum of w^2 over the spots that entered - the beam prior below is scaled by it so that its
+        // strength relative to the data is the same weighted or not. Equals the residual block count
+        // when the spots are unweighted.
+        double effective_spots = 0.0;
+
        for (int i = 0; i < spots.size(); i++) {
            if (spots[i].empty())
                continue;

+            std::vector<double> weight; // empty = unweighted
+            if (data.weight_spots_by_confidence)
+                weight = SpotConfidenceWeights(spots[i]);
+
            double angle_rad = 0.0;
            std::optional<RotMatrix> rot_matr;

@@ -158,7 +208,8 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
            }

            // Add residuals for each point
-            for (const auto &pt: spots[i]) {
+            for (size_t j = 0; j < spots[i].size(); j++) {
+                const auto &pt = spots[i][j];
                if (!data.index_ice_rings && pt.ice_ring)
                    continue;

@@ -180,6 +231,9 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
                if (norm_sq > tolerance_sq)
                    continue;

+                const double weight_sq = weight.empty() ? 1.0 : weight[j] * weight[j];
+                effective_spots += weight_sq;
+
                problem.AddResidualBlock(
                    new ceres::AutoDiffCostFunction<XtalResidual, 3, 2, 1, 2, 3, 3, 3, 3>(
                        new XtalResidual(pt.x, pt.y,
@@ -189,7 +243,11 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
                                         angle_rad,
                                         h, k, l,
                                         data.crystal_system)),
-                    nullptr,
+                    // Ceres has no per-residual weight; ScaledLoss(nullptr, a) multiplies the squared
+                    // residual by the constant a, i.e. it applies a weight of sqrt(a) to the residual.
+                    weight.empty()
+                        ? nullptr
+                        : new ceres::ScaledLoss(nullptr, weight_sq, ceres::TAKE_OWNERSHIP),
                    beam,
                    &distance_mm,
                    detector_rot,
@@ -231,7 +289,7 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
            // perpendicular direction, the prior wins the gauge one.
            constexpr double sigma_px = 3.0;
            const double k = data.geom.GetPixelSize_mm() / (distance_mm * data.geom.GetWavelength_A());
-            const double w = k * std::sqrt(static_cast<double>(problem.NumResidualBlocks())) / sigma_px;
+            const double w = k * std::sqrt(effective_spots) / sigma_px;
            problem.AddResidualBlock(
                new ceres::AutoDiffCostFunction<BeamComponentPrior, 1, 2>(
                    new BeamComponentPrior(parallel, beam[parallel], w)),