rugnux: add -r multi geometry refinement (per-image best of three)

Per-image geometry refinement is a tradeoff that flips per dataset: beam+cell refinement
extends OCP merge quality but DIVERGES on sparse-spot stills (e.g. KR2's ~20 spots around a
232 A axis), pushing good lattices out of tolerance so they fail the acceptance floor.
`-r multi` runs all three (none / orientation / beam_and_lattice) on a copy per image and
keeps whichever indexes the most spots (scorer = fractional-Miller-within-tolerance count,
mirroring AnalyzeIndexing); ties prefer less refinement, to avoid overfitting the sparse list.

Validated: KR2 index 7.08% (-r beam_and_lattice default) -> 10.05% (matches the -r none best),
while OCP R-free stays ~equal to beam_and_lattice. New GeomRefinementAlgorithmEnum::Multi handled
in the CLI/HDF5/command-line echoes; the API convert maps it to beam_and_lattice (offline only).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This commit is contained in:
2026-07-14 12:04:32 +02:00
co-authored by Claude Opus 4.8
parent fb399480e5
commit 4a5028d98d
6 changed files with 62 additions and 2 deletions
+50
View File
@@ -176,6 +176,23 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data
return outcome;
}
namespace {
// Count spots whose fractional Miller index falls within the indexing tolerance of an integer for a
// given lattice + geometry - the "how well does this model explain the spots" score used by -r multi.
int CountIndexedSpots(const DiffractionGeometry &geom, const CrystalLattice &latt,
const std::vector<SpotToSave> &spots, float tol_sq) {
const Coord a = latt.Vec0(), b = latt.Vec1(), c = latt.Vec2();
int n = 0;
for (const auto &s : spots) {
const Coord recip = s.ReciprocalCoord(geom);
const float hf = recip * a, kf = recip * b, lf = recip * c;
const float dh = hf - std::round(hf), dk = kf - std::round(kf), dl = lf - std::round(lf);
if (dh * dh + dk * dk + dl * dl < tol_sq) ++n;
}
return n;
}
} // namespace
void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::IndexingOutcome &outcome) {
if (!outcome.lattice_candidate)
return;
@@ -213,6 +230,39 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In
outcome.beam_center_updated = true;
}
break;
case GeomRefinementAlgorithmEnum::Multi: {
// Try all three refinements per image and keep whichever indexes the most spots. Beam+cell
// refinement helps some stills but diverges on sparse spot lists (few spots, long axes),
// where it pushes a good lattice out of tolerance; scoring by indexed-spot count lets each
// image fall back to orientation-only or no refinement when refinement would hurt. Ties
// prefer less refinement (strict >, not >=) to avoid overfitting.
const float tol = experiment.GetIndexingSettings().GetTolerance();
const float tol_sq = tol * tol;
XtalOptimizerData d_none = data;
XtalOptimizerData d_orient = data;
XtalOptimizerRotationOnly(d_orient, msg.spots, 0.2);
XtalOptimizerRotationOnly(d_orient, msg.spots, 0.1);
XtalOptimizerRotationOnly(d_orient, msg.spots, 0.05);
XtalOptimizerData d_beam = data;
const bool beam_ok = XtalOptimizer(d_beam, {msg.spots});
const int s_none = CountIndexedSpots(d_none.geom, d_none.latt, msg.spots, tol_sq);
const int s_orient = CountIndexedSpots(d_orient.geom, d_orient.latt, msg.spots, tol_sq);
const int s_beam = beam_ok ? CountIndexedSpots(d_beam.geom, d_beam.latt, msg.spots, tol_sq) : -1;
if (s_beam > s_none && s_beam > s_orient) {
data = d_beam;
outcome.experiment.BeamX_pxl(data.geom.GetBeamX_pxl())
.BeamY_pxl(data.geom.GetBeamY_pxl());
outcome.beam_center_updated = true;
} else if (s_orient > s_none) {
data = d_orient;
} else {
data = d_none;
}
break;
}
}
outcome.lattice_candidate = data.latt;