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;