From ae126c3d5b02ec3cbe3fa76daa700a51f0b871a1 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 28 Jul 2026 09:53:13 +0200 Subject: [PATCH] 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) --- image_analysis/IndexAndRefine.cpp | 4 ++ .../geom_refinement/XtalOptimizer.cpp | 64 ++++++++++++++++++- .../geom_refinement/XtalOptimizer.h | 5 ++ 3 files changed, 70 insertions(+), 3 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 3e5e6885..a90a7819 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -254,6 +254,10 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In .refine_distance_mm = false, .refine_detector_angles = false, .refine_unit_cell = !experiment.IsRotationIndexing(), + // The whole spot list is passed below, not the indexed subset, so on weak images most of what + // enters the fit at the loose first tolerance is arbitrarily indexed noise. Weight every spot by + // how strong it is for its resolution so those contribute without dragging the orientation. + .weight_spots_by_confidence = true, .max_time = 0.04 // 40 ms is max allowed time for the operation }; diff --git a/image_analysis/geom_refinement/XtalOptimizer.cpp b/image_analysis/geom_refinement/XtalOptimizer.cpp index 45774851..e35d4e8e 100644 --- a/image_analysis/geom_refinement/XtalOptimizer.cpp +++ b/image_analysis/geom_refinement/XtalOptimizer.cpp @@ -2,6 +2,8 @@ // SPDX-License-Identifier: GPL-3.0-only #include "../../common/JFJochMath.h" +#include +#include #include #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 SpotConfidenceWeights(const std::vector &spots) { + constexpr size_t spots_per_shell = 32; + + std::vector 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(1, spots.size() / spots_per_shell); + std::vector weight(spots.size()); + std::vector 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> &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 weight; // empty = unweighted + if (data.weight_spots_by_confidence) + weight = SpotConfidenceWeights(spots[i]); + double angle_rad = 0.0; std::optional 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( 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(problem.NumResidualBlocks())) / sigma_px; + const double w = k * std::sqrt(effective_spots) / sigma_px; problem.AddResidualBlock( new ceres::AutoDiffCostFunction( new BeamComponentPrior(parallel, beam[parallel], w)), diff --git a/image_analysis/geom_refinement/XtalOptimizer.h b/image_analysis/geom_refinement/XtalOptimizer.h index d28a884a..105a2faf 100644 --- a/image_analysis/geom_refinement/XtalOptimizer.h +++ b/image_analysis/geom_refinement/XtalOptimizer.h @@ -30,6 +30,11 @@ struct XtalOptimizerData { bool index_ice_rings = true; + // Weight each spot by how strong it is for its resolution, so that low-confidence spots contribute + // without driving the fit (see SpotConfidenceWeights). Off by default: the indexers call this with a + // spot list they have already selected, it is the per-image refinement that gets the raw list. + bool weight_spots_by_confidence = false; + float max_time = 1.0; std::optional axis;