beam centre: run the four consistency restarts at once

FindBeamCenterFromSpotSymmetry calls Estimate() five times - once for the answer and once from each
of four starts 25 px away - and three quarters of each of those is an 861-point brute-force grid
over the spindle. So making the uncertainty gate live was paid for by multiplying the estimator by
five, which is the whole of the pre-scan's cost.

The four restarts share nothing: each takes its own copy of the geometry and only reads the spots.
What is wanted from them is a max, which is order-independent, so running them concurrently gives
the same number. Measured on three rotation crystals, the committed centre, the reported sigma and
the fitted spindle angles are unchanged to every printed digit.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01T3yNBXk4wKdMZy1ak2NY7f
This commit is contained in:
2026-08-30 08:37:55 +02:00
co-authored by Claude Opus 5
parent 3ddfbb2da9
commit 97f57a430d
@@ -6,7 +6,9 @@
#include "../../common/JFJochMath.h" // PI
#include <algorithm>
#include <array>
#include <cmath>
#include <future>
#include <numeric>
#include <tuple>
@@ -888,19 +890,32 @@ FindBeamCenterFromSpotSymmetry(const DiffractionExperiment &experiment,
// nothing here can say which is the crystal's, so the scatter of the frame pairs - which stays
// small for either of them - is not the uncertainty and this is. The spindle is re-fitted from
// each start, so a fit that depends on where the search began is part of what is reported.
//
// The four run at once. Each takes its own copy of the geometry and only reads the spots, so they
// share nothing, and what is wanted from them is a MAX - which does not care in what order they
// finish. This is the whole of the estimator's cost: Estimate() is a brute-force grid over the
// spindle, and asking for the uncertainty runs it five times.
const std::array<std::pair<float, float>, 4> starts{{{CONSISTENCY_START_PXL, 0.0f},
{-CONSISTENCY_START_PXL, 0.0f},
{0.0f, CONSISTENCY_START_PXL},
{0.0f, -CONSISTENCY_START_PXL}}};
std::vector<std::future<float>> restarts;
restarts.reserve(starts.size());
for (const auto &[dx, dy]: starts)
restarts.push_back(std::async(std::launch::async, [&, dx = dx, dy = dy] {
DiffractionGeometry from = geom;
from.BeamX_pxl(estimate->beam_x_pxl + dx).BeamY_pxl(estimate->beam_y_pxl + dy);
SpindleEstimate again;
if (const auto other = Estimate(from, *goniometer, frame_angle_deg, spots,
spindle_estimate ? &again : nullptr))
return std::hypot(other->beam_x_pxl - estimate->beam_x_pxl,
other->beam_y_pxl - estimate->beam_y_pxl);
return 0.0f;
}));
float spread = 0.0f;
for (const auto &[dx, dy]: {std::pair<float, float>{CONSISTENCY_START_PXL, 0.0f},
{-CONSISTENCY_START_PXL, 0.0f},
{0.0f, CONSISTENCY_START_PXL},
{0.0f, -CONSISTENCY_START_PXL}}) {
DiffractionGeometry from = geom;
from.BeamX_pxl(estimate->beam_x_pxl + dx).BeamY_pxl(estimate->beam_y_pxl + dy);
SpindleEstimate again;
if (const auto other = Estimate(from, *goniometer, frame_angle_deg, spots,
spindle_estimate ? &again : nullptr))
spread = std::max(spread, std::hypot(other->beam_x_pxl - estimate->beam_x_pxl,
other->beam_y_pxl - estimate->beam_y_pxl));
}
for (auto &restart: restarts)
spread = std::max(spread, restart.get());
estimate->sigma_pxl = std::max(estimate->sigma_pxl, spread);
return estimate;
}