mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 72 Wiki Activity

Ask the geometry refinement for the derivatives it uses, and solve the normal equations

Browse Source 
Four changes to the same least-squares fit, which the rotation first pass runs on every
candidate lattice and the per-image path runs on every frame.

The linear solver was DENSE_QR on a problem that is very tall and thin - thousands of spots
against at most seventeen parameters. That is the shape QR handles worst: it copies the
Jacobian out of Ceres' row-major storage into a column-major buffer on every solve, and
Eigen's blocked Householder then degenerates to the unblocked path because its block size is
the column count. Accumulating J^T J reads the Jacobian once instead. Both solve the same
damped system, so the step is the same to round-off.

Ceres sizes its dual numbers from the declared parameter blocks, not from which of them the
caller then holds constant. Nothing outside a test set refine_distance_mm - the positional
residual leaves the distance degenerate with the cell scale, which is why the rotation
post-refinement fits it in a step of its own with the cell held fixed - so the block was
declared only to be frozen, and every residual differentiated seventeen parameters to use
sixteen. It is gone, along with the test that exercised distance recovery; that test seeded
the distance off truth, which the cell would now absorb, so its seed moves to the true value.
The post-refinement's own detector step held five of its seven blocks constant and now bakes
them into the residual, leaving beam and distance.

The predicted reciprocal vector was built by rotating all three direct columns and then
crossing them. A rotation commutes with the cross product and leaves the triple product
alone, so the same vector comes out of crossing the unrotated columns and turning the result
once - three rotations become one, for every crystal system.

The documentation described the arrangement before all this, and had drifted in a second way:
the first-pass rotation indexing has been refining the detector tilt and the rotation axis by
default, which the text said were held fixed.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
jungfrau
2026-08-24 05:31:59 -04:00
co-authored by Claude Opus 5
parent 4b48d19064
commit 0d07141b5c
8 changed files with 168 additions and 81 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/geom_refinement/PostRefine.cpp
+70 -34
View File
@@ -475,32 +475,47 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
const int cost_nt = static_cast<int>(std::max<size_t>(1, std::min(nthreads,
static_cast<size_t>(n_terms))));
const int cost_chunk = (n_terms + cost_nt - 1) / cost_nt;
const auto cost_by_fifth = [&](double k) {
std::vector<std::array<double, 5>> per_chunk(cost_nt);
// Several k are always wanted at once (the grid below asks for 101), and they all sweep the
// same event list, so sweep it ONCE and evaluate every k on each event while it is still in
// registers. The per-thread accumulator is one slot per (k, fifth) - 4 kB for the grid, small
// enough to stay in L1 - against re-reading the whole term array once per k. Each (k, fifth)
// still receives its events in the same order and the chunks are still summed in chunk order,
// so the sums are the ones a k-at-a-time loop produced, bit for bit.
const auto cost_grid = [&](const std::vector<double> &ks) {
const int nk = static_cast<int>(ks.size());
std::vector<std::vector<std::array<double, 5>>> per_chunk(
cost_nt, std::vector<std::array<double, 5>>(nk));
ParallelChunks(n_terms, nthreads, [&](int lo, int hi) {
std::array<double, 5> acc{};
std::vector<std::array<double, 5>> acc(nk);
for (int e = lo; e < hi; ++e) {
const double r = residual_at(terms[e], k);
const double s2 = r * r;
acc[fifth_of[e]] += (s2 <= huber_d2) ? s2
: (2.0 * huber_delta * std::sqrt(s2) - huber_d2);
const ScaleTerm &term = terms[e];
const int fifth = fifth_of[e];
for (int g = 0; g < nk; ++g) {
const double r = residual_at(term, ks[g]);
const double s2 = r * r;
acc[g][fifth] += (s2 <= huber_d2) ? s2
: (2.0 * huber_delta * std::sqrt(s2) - huber_d2);
}
}
per_chunk[lo / cost_chunk] = acc;
per_chunk[lo / cost_chunk] = std::move(acc);
});
std::array<double, 5> total{};
std::vector<std::array<double, 5>> total(nk);
for (const auto &acc : per_chunk)
for (int j = 0; j < 5; ++j) total[j] += acc[j];
for (int g = 0; g < nk; ++g)
for (int j = 0; j < 5; ++j) total[g][j] += acc[g][j];
return total;
};
const auto cost_by_fifth = [&](double k) { return cost_grid({k})[0]; };
// Scan the interval Ceres was bounded to, then close in. No event's phase can move by more than
// a fraction of a period over an interval this narrow, so the objective has no structure the
// grid could step over; the refinement is only there to place the minimum precisely.
constexpr int SCALE_GRID = 101;
constexpr double SCALE_K_LO = 0.95, SCALE_K_HI = 1.05;
std::vector<std::array<double, 5>> grid(SCALE_GRID);
std::vector<double> grid_k(SCALE_GRID);
for (int g = 0; g < SCALE_GRID; ++g)
grid[g] = cost_by_fifth(SCALE_K_LO + (SCALE_K_HI - SCALE_K_LO) * g / (SCALE_GRID - 1));
grid_k[g] = SCALE_K_LO + (SCALE_K_HI - SCALE_K_LO) * g / (SCALE_GRID - 1);
const std::vector<std::array<double, 5>> grid = cost_grid(grid_k);
auto solve_scale = [&](int drop_fifth) {
const auto total = [&](const std::array<double, 5> &f) {
double t = 0.0;
@@ -601,23 +616,45 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
double rot_vec[3] = {axv[0], axv[1], axv[2]}; // committed (or nominal) axis
// ===== Step B: detector distance + beam from the observed positions, cell fixed =====
std::vector<const Partial *> obs;
for (size_t i = 0; i < n_pts; ++i)
if (std::isfinite(pts[i].obs_x) && std::isfinite(pts[i].obs_y)) obs.push_back(&pts[i]);
// Count first, then fill, exactly as the partial gather above does and for the same reason:
// this walks the same tens of millions of partials, and a pointer vector grown by push_back
// copies itself every time it doubles. new[] rather than a sized vector so the array is not
// zeroed on one thread before the parallel fill overwrites it. The fill lands in the order
// the serial loop produced, so the selection below sees the same sequence it always did.
const int n_obs_chunks = static_cast<int>(std::clamp<size_t>(nthreads, 1,
std::max<size_t>(1, n_pts)));
const size_t obs_chunk = (n_pts + n_obs_chunks - 1) / n_obs_chunks;
std::vector<size_t> obs_offset(n_obs_chunks + 1, 0);
const auto keep_obs = [&](size_t i) {
return std::isfinite(pts[i].obs_x) && std::isfinite(pts[i].obs_y);
};
ParallelChunks(static_cast<int>(n_pts), nthreads, [&](int lo, int hi) {
size_t keep = 0;
for (int i = lo; i < hi; ++i) if (keep_obs(i)) keep++;
obs_offset[static_cast<size_t>(lo) / obs_chunk + 1] = keep;
});
for (int c = 0; c < n_obs_chunks; ++c) obs_offset[c + 1] += obs_offset[c];
size_t n_obs = obs_offset[n_obs_chunks];
std::unique_ptr<const Partial *[]> obs(new const Partial *[n_obs]);
ParallelChunks(static_cast<int>(n_pts), nthreads, [&](int lo, int hi) {
size_t at = obs_offset[static_cast<size_t>(lo) / obs_chunk];
for (int i = lo; i < hi; ++i) if (keep_obs(i)) obs[at++] = &pts[i];
});
constexpr size_t MAX_OBS = 20000;
if (obs.size() > MAX_OBS) {
std::nth_element(obs.begin(), obs.begin() + MAX_OBS, obs.end(),
if (n_obs > MAX_OBS) {
std::nth_element(obs.get(), obs.get() + MAX_OBS, obs.get() + n_obs,
[](const Partial *a, const Partial *b) {
return a->I / std::max(1e-9, static_cast<double>(a->sigma))
> b->I / std::max(1e-9, static_cast<double>(b->sigma)); });
obs.resize(MAX_OBS);
n_obs = MAX_OBS;
}
result.obs_used = static_cast<int>(obs.size());
result.obs_used = static_cast<int>(n_obs);
const double beam_x0 = nominal_geom.GetBeamX_pxl(), beam_y0 = nominal_geom.GetBeamY_pxl();
const double dist0 = nominal_geom.GetDetectorDistance_mm();
auto pos_cost = [&](Subset s, const double beam[2], const double dist[1]) {
double c = 0.0; int n = 0;
for (const Partial *pp : obs) {
for (size_t oi = 0; oi < n_obs; ++oi) {
const Partial *pp = obs[oi];
if (!in(pp->h, pp->k, pp->l, s)) continue;
XtalResidual r(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, rot3, angle_rad(pp->img),
pp->h, pp->k, pp->l, sys);
@@ -630,21 +667,20 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
auto solve_detector = [&](Subset s, double beam_out[2], double &dist_out) {
double beam[2] = {beam_x0, beam_y0}, dist[1] = {dist0};
ceres::Problem p;
for (const Partial *pp : obs) {
for (size_t oi = 0; oi < n_obs; ++oi) {
const Partial *pp = obs[oi];
if (!in(pp->h, pp->k, pp->l, s)) continue;
p.AddResidualBlock(new ceres::AutoDiffCostFunction<XtalResidual, 3, 2, 1, 2, 3, 3, 3, 3>(
new XtalResidual(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, rot3, angle_rad(pp->img),
pp->h, pp->k, pp->l, sys)),
new ceres::CauchyLoss(0.02), beam, dist,
const_cast<double *>(det_rot), rot_vec, p0, p1, p2);
p.AddResidualBlock(new ceres::AutoDiffCostFunction<XtalResidualBeamDistance, 3, 2, 1>(
new XtalResidualBeamDistance(
XtalResidual(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, rot3, angle_rad(pp->img),
pp->h, pp->k, pp->l, sys),
det_rot, rot_vec, p0, p1, p2)),
new ceres::CauchyLoss(0.02), beam, dist);
}
if (p.NumResidualBlocks() == 0) { beam_out[0] = beam_x0; beam_out[1] = beam_y0; dist_out = dist0; return false; }
// Everything but the beam and the distance is held at its step-A value. Once per problem:
// the block only has to exist, and repeating it per observation was up to MAX_OBS times
// five calls for the same five blocks.
p.SetParameterBlockConstant(const_cast<double *>(det_rot));
p.SetParameterBlockConstant(rot_vec);
p.SetParameterBlockConstant(p0); p.SetParameterBlockConstant(p1); p.SetParameterBlockConstant(p2);
// Everything but the beam and the distance is held at its step-A value, and the residual
// above bakes those five blocks in rather than declaring them and freezing them, so there
// is nothing left to hold constant here.
p.SetParameterLowerBound(dist, 0, dist0 * 0.95); p.SetParameterUpperBound(dist, 0, dist0 * 1.05);
for (int j = 0; j < 2; ++j) { p.SetParameterLowerBound(beam, j, beam[j] - 15.0);
p.SetParameterUpperBound(beam, j, beam[j] + 15.0); }
@@ -655,7 +691,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
return sum.IsSolutionUsable();
};
double beam[2] = {beam_x0, beam_y0}, dist = dist0;
if (obs.size() >= static_cast<size_t>(settings.min_events)) {
if (n_obs >= static_cast<size_t>(settings.min_events)) {
double beam_fit[2], dist_fit;
const bool convB = solve_detector(FIT, beam_fit, dist_fit);
const double b_nom[2] = {beam_x0, beam_y0}, d_nom[1] = {dist0};
@@ -679,7 +715,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector<IntegrationOutcome
result.detector_refined ? beam[0] : beam_x0, result.detector_refined ? beam[1] : beam_y0,
cvB_nom, cvB_ref, result.detector_refined ? "COMMIT" : "reject (kept nominal detector)");
} else {
logger.Info("Post-refine GEOM step B: only {} positional observations - skipped", obs.size());
logger.Info("Post-refine GEOM step B: only {} positional observations - skipped", n_obs);
}
// Assemble the committed geometry.