diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index c610aba2..cf907267 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,6 +3,7 @@ ### 1.0.0-rc.161 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. +* rugnux: New **beam-centre measurement before indexing** (`--estimate-beam-center`), which takes the direct beam from the symmetry of the spots where the sweep reaches at least half a turn, from the radial background profile where it does not, and keeps the value in the file where neither can measure it; the rotation axis' skew about the beam is fitted at the same time (`--no-fit-spindle` keeps the axis from the file). * Lattice search: a reduced cell whose β sits on the Niggli type boundary now keeps its centring instead of falling through to triclinic, so a centred lattice is no longer read as primitive according to which way its refined cell happened to reduce. * rugnux: The second rotation pass now reuses a space group that came from the intensity-based centred-lattice test, as it already did for every other group, instead of re-deciding the symmetry at the refined geometry. * Space-group search: where the search on the well-measured observations (`--search-min-zeta`) and the search on all of them disagree, the answer is now always the one the merge of all the observations supports, instead of whichever of the two found more symmetry. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index c0a049e4..bf987c50 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -151,6 +151,37 @@ Jungfraujoch uses $|\Delta_\mathrm{Ewald}|$ as an operational proxy for excitati - profile radius estimation (see §11.1), - still partiality option in scaling/merging (§10.2). +### 1.4 Measuring the direct beam before indexing + +Two exact facts about a rotation sweep let the direct beam be measured from spot positions alone, +before anything is indexed (`--estimate-beam-center`). + +**Friedel mates half a turn apart.** Rotating 180° about the spindle $\mathbf{m}$ and taking $-h$ +negates a reflection's component along $\mathbf{m}$ and leaves the rest. With the spindle +perpendicular to the beam this preserves $\mathbf{q}\cdot\mathbf{S}_0$, so $-h$ satisfies the Laue +condition at $\varphi+180°$ exactly where $h$ satisfies it at $\varphi$, and the spots recorded half a +turn apart are mirror images along the spindle direction. This gives the beam coordinate **along** +the spindle. Note these are Friedel mates, not the same reflection. + +**The second crossing.** The same reflection meets the Ewald sphere twice, at two angles that are +generally *not* 180° apart, differing only in the sign of the lab component perpendicular to both +$\mathbf{m}$ and the beam. This gives the remaining coordinate. The two crossings are separated by a +sweep angle fixed by the reflection's own position, so genuine pairs are identified without a cell or +an orientation matrix. + +Neither observable requires the reflections to be indexed: each candidate pairing votes for a beam +coordinate, and the true value accumulates while wrong pairings scatter. + +The mirror is exact in the **laboratory** frame, so it is sensitive to the spindle's orientation. A +skew of the spindle about the beam *spreads* the vote rather than shifting it, and is fitted +alongside the centre (`--no-fit-spindle` keeps the axis from the file); a tilt of the spindle towards +the beam is measured and reported but not applied. The frames read are sampled away from both ends of +the sweep, where shutter synchronisation can spoil an image. + +Where the sweep is shorter than half a turn the spot symmetry cannot be formed, and the centre is +taken instead from the centroid of the radial background profile, which needs only a few images. +Where neither method can measure the centre, the value from the file is kept. + --- ## 2. Azimuthal integration (radial profiles) diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index 804f3f84..c26598a1 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -393,6 +393,13 @@ Detector mask: | --- | --- | | `--detect-beam-stop[=N\|off]` | Find the beam stop and its holder in a projection of N images and add them to the pixel mask as bit 9, so nothing shadowed by them is integrated. **On by default** (60 images); `=off` disables. Reflections behind the stop are attenuated but not flagged, so they integrate low with a plausible sigma and no existing rejection catches them | +Geometry: + +| Option | Description | +| --- | --- | +| `--estimate-beam-center` | Measure the direct beam before indexing, from the symmetry of the spots where the sweep reaches at least half a turn and from the radial background profile where it does not; the value in the file is kept where neither can measure it. Off by default | +| `--no-fit-spindle` | With the above, keep the rotation axis given in the file instead of fitting its skew about the beam | + Spot finding: | Option | Description | diff --git a/image_analysis/beam_stop/ShadowFinder.cpp b/image_analysis/beam_stop/ShadowFinder.cpp index 93928281..59a65d9b 100644 --- a/image_analysis/beam_stop/ShadowFinder.cpp +++ b/image_analysis/beam_stop/ShadowFinder.cpp @@ -243,6 +243,16 @@ uint32_t ShadowFinder::GetFrameCount() const { return frames; } +std::vector ShadowFinder::GetMeanProjection() const { + std::unique_lock ul(m); + + std::vector mean(static_cast(width) * height, NAN); + for (size_t i = 0; i < mean.size(); i++) + if (valid_count[i] > 0 && pixel_mask[i] == 0) + mean[i] = static_cast(static_cast(sum_value[i]) / valid_count[i]); + return mean; +} + std::vector ShadowFinder::GetMask() const { std::unique_lock ul(m); diff --git a/image_analysis/beam_stop/ShadowFinder.h b/image_analysis/beam_stop/ShadowFinder.h index ef56fa7d..7426b481 100644 --- a/image_analysis/beam_stop/ShadowFinder.h +++ b/image_analysis/beam_stop/ShadowFinder.h @@ -59,5 +59,10 @@ public: // Recomputed from the accumulators on each call - meant to be called once at the end. [[nodiscard]] std::vector GetMask() const; + // Mean counts per pixel over the frames added, NAN where nothing was counted. This is the + // projection GetMask() tests, so anything else that wants the background before indexing + // gets it without reading the frames a second time. + [[nodiscard]] std::vector GetMeanProjection() const; + [[nodiscard]] uint32_t GetFrameCount() const; }; diff --git a/image_analysis/geom_refinement/BeamCenterFromBackground.cpp b/image_analysis/geom_refinement/BeamCenterFromBackground.cpp new file mode 100644 index 00000000..bba29ad9 --- /dev/null +++ b/image_analysis/geom_refinement/BeamCenterFromBackground.cpp @@ -0,0 +1,270 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "BeamCenterFromBackground.h" + +#include +#include + +#include "../../common/JFJochMath.h" + +namespace { + +// The band the background is fitted over. The low-resolution end sits outside the beam stop +// and its penumbra, the high-resolution end where the solvent ring has died away. +constexpr float BAND_LOW_RES_A = 12.0f; +constexpr float BAND_HIGH_RES_A = 2.2f; + +constexpr int SECTORS = 36; +constexpr int RADIAL_BINS = 120; + +// A cell with fewer pixels than this has no usable mean. +constexpr int MIN_PIXELS_PER_CELL = 20; + +// A radial bin missing more azimuth than this is a partial ring - it leaves the detector, or a +// module gap eats it - and a partial ring biases the profile it is compared against. +constexpr float MIN_SECTOR_COVERAGE = 0.85f; + +// Fractions of a bin's pixels are Bragg peaks. Two rounds of clipping at the upper 2 sigma take +// the mean back to the background without needing the pixel values a second time. +constexpr float CLIP_SIGMA = 2.0f; +constexpr int CLIP_ROUNDS = 2; + +// The profile is rebuilt at the trial centre every iteration, so a centre that is off smears the +// solvent ring and flattens g', which over-estimates the shift. Half steps damp that; the fixed +// point is unchanged, only the path to it. +constexpr float DAMPING = 0.5f; +constexpr int MAX_ITERATIONS = 10; +constexpr float CONVERGED_PXL = 0.02f; + +// Below these the fit has not seen enough of the detector to be believed at all. +constexpr int MIN_USABLE_SECTORS = SECTORS * 3 / 5; +constexpr int MIN_USABLE_RADIAL_BINS = 15; + +float median_of(std::vector &v) { + const size_t half = v.size() / 2; + std::nth_element(v.begin(), v.begin() + half, v.end()); + return v[half]; +} + +} // namespace + +std::optional +FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const PixelMask &mask, + const std::vector &mean) { + const auto W = static_cast(experiment.GetXPixelsNumConv()); + const auto H = static_cast(experiment.GetYPixelsNumConv()); + const size_t n_pixels = static_cast(W) * H; + if (mean.size() != n_pixels) + return {}; + + const auto &pixel_mask = mask.GetMask(experiment); + auto geom = experiment.GetDiffractionGeometry(); + + const float wavelength = geom.GetWavelength_A(); + const float sin_high = wavelength / (2.0f * BAND_HIGH_RES_A); + if (sin_high >= 1.0f) + return {}; + const float tt_lo = 2.0f * std::asin(wavelength / (2.0f * BAND_LOW_RES_A)); + const float tt_hi = 2.0f * std::asin(sin_high); + const float d_tt = (tt_hi - tt_lo) / RADIAL_BINS; + + const auto rot = geom.GetPoniRotMatrix().arr(); // row major + const float pixel_size = geom.GetPixelSize_mm(); + const float distance = geom.GetDetectorDistance_mm(); + + float beam_x = geom.GetBeamX_pxl(); + float beam_y = geom.GetBeamY_pxl(); + + constexpr int n_cells = RADIAL_BINS * SECTORS; + std::vector cell_of(n_pixels); + std::vector sum(n_cells), sum_sq(n_cells), sum_jx(n_cells), sum_jy(n_cells); + std::vector count(n_cells), count_all(n_cells); + std::vector profile(RADIAL_BINS), d_profile(RADIAL_BINS), clip_limit(n_cells); + std::vector radial_ok(RADIAL_BINS); + + float step_x = 0.0f, step_y = 0.0f, sigma_x = 0.0f, sigma_y = 0.0f; + for (int iteration = 0; iteration < MAX_ITERATIONS; iteration++) { + std::fill(sum.begin(), sum.end(), 0.0); + std::fill(sum_sq.begin(), sum_sq.end(), 0.0); + std::fill(sum_jx.begin(), sum_jx.end(), 0.0); + std::fill(sum_jy.begin(), sum_jy.end(), 0.0); + std::fill(count.begin(), count.end(), 0); + + for (int y = 0; y < H; y++) { + for (int x = 0; x < W; x++) { + const size_t i = static_cast(y) * W + x; + cell_of[i] = -1; + if (pixel_mask[i] != 0 || !std::isfinite(mean[i])) + continue; + const float u = (x - beam_x) * pixel_size; + const float v = (y - beam_y) * pixel_size; + const float lx = rot[0] * u + rot[1] * v + rot[2] * distance; + const float ly = rot[3] * u + rot[4] * v + rot[5] * distance; + const float lz = rot[6] * u + rot[7] * v + rot[8] * distance; + const float rho = std::sqrt(lx * lx + ly * ly); + const float two_theta = std::atan2(rho, lz); + if (two_theta < tt_lo || two_theta >= tt_hi || rho == 0.0f) + continue; + + const float phi = std::atan2(ly, lx); + // Both bins are clamped: a pixel one float ulp below the top of the band divides + // to exactly RADIAL_BINS, which is one cell past the end of every accumulator. + const int r_bin = std::clamp(static_cast((two_theta - tt_lo) / d_tt), 0, RADIAL_BINS - 1); + const int s_bin = std::clamp(static_cast((phi + PI) / (2 * PI) * SECTORS), 0, SECTORS - 1); + const int cell = r_bin * SECTORS + s_bin; + + // d(2theta)/d(beam), through the lab coordinate: the detector coordinate depends + // on the centre only as (x - beam_x), so moving the centre is moving the pixel. + const float denominator = rho * rho + lz * lz; + const float g_x = lz * lx / (rho * denominator); + const float g_y = lz * ly / (rho * denominator); + const float g_z = -rho / denominator; + cell_of[i] = cell; + count[cell]++; + sum[cell] += mean[i]; + sum_sq[cell] += static_cast(mean[i]) * mean[i]; + sum_jx[cell] += -pixel_size * (g_x * rot[0] + g_y * rot[3] + g_z * rot[6]); + sum_jy[cell] += -pixel_size * (g_x * rot[1] + g_y * rot[4] + g_z * rot[7]); + } + } + + count_all = count; // the Jacobian sums belong to the unclipped pixel set + + for (int round = 0; round < CLIP_ROUNDS; round++) { + for (int c = 0; c < n_cells; c++) { + if (count[c] < MIN_PIXELS_PER_CELL) { clip_limit[c] = -1.0f; continue; } + const double m = sum[c] / count[c]; + const double variance = std::max(sum_sq[c] / count[c] - m * m, 0.0); + clip_limit[c] = static_cast(m + CLIP_SIGMA * std::sqrt(variance)); + } + std::fill(sum.begin(), sum.end(), 0.0); + std::fill(sum_sq.begin(), sum_sq.end(), 0.0); + std::fill(count.begin(), count.end(), 0); + for (size_t i = 0; i < n_pixels; i++) { + const int32_t c = cell_of[i]; + if (c < 0 || clip_limit[c] < 0.0f || mean[i] > clip_limit[c]) + continue; + count[c]++; + sum[c] += mean[i]; + sum_sq[c] += static_cast(mean[i]) * mean[i]; + } + } + + // Radial profile: the median over the sectors that have a mean, on rings that are + // almost fully covered. + int usable_radial = 0; + for (int r = 0; r < RADIAL_BINS; r++) { + std::vector present; + for (int s = 0; s < SECTORS; s++) + if (count[r * SECTORS + s] >= MIN_PIXELS_PER_CELL) + present.push_back(static_cast(sum[r * SECTORS + s] / count[r * SECTORS + s])); + radial_ok[r] = static_cast(present.size()) >= MIN_SECTOR_COVERAGE * SECTORS; + profile[r] = radial_ok[r] ? median_of(present) : 0.0f; + usable_radial += radial_ok[r]; + } + if (usable_radial < MIN_USABLE_RADIAL_BINS) + return {}; + // Central difference, so a bin next to a gap in the profile drops out with it. The test + // reads the ring BEFORE it, so it has to read the covered/not-covered flags as they were, + // not as this same loop has already rewritten them. + const std::vector covered = radial_ok; + for (int r = 0; r < RADIAL_BINS; r++) { + const bool have = r > 0 && r + 1 < RADIAL_BINS && covered[r - 1] && covered[r] && covered[r + 1]; + d_profile[r] = have ? (profile[r + 1] - profile[r - 1]) / (2 * d_tt) : 0.0f; + radial_ok[r] = have; + } + + // Per sector: regress (profile of the sector - common profile) on {g, g'}. The first + // coefficient is the sector's amplitude, the second its radial shift; only the shift + // is carried on. + std::vector shift, weight, jacobian_x, jacobian_y; + for (int s = 0; s < SECTORS; s++) { + double a11 = 0, a12 = 0, a22 = 0, b1 = 0, b2 = 0; + double jx = 0, jy = 0; + int n = 0; + for (int r = 0; r < RADIAL_BINS; r++) { + const int c = r * SECTORS + s; + if (!radial_ok[r] || count[c] < MIN_PIXELS_PER_CELL) + continue; + const double g = profile[r], dg = d_profile[r]; + const double y = sum[c] / count[c] - profile[r]; + a11 += g * g; a12 += g * dg; a22 += dg * dg; + b1 += g * y; b2 += dg * y; + jx += sum_jx[c] / count_all[c]; + jy += sum_jy[c] / count_all[c]; + n++; + } + const double det = a11 * a22 - a12 * a12; + if (n < MIN_USABLE_RADIAL_BINS || det <= 0) + continue; + const double amplitude = (a22 * b1 - a12 * b2) / det; + const double this_shift = (a11 * b2 - a12 * b1) / det; + // Residual sum of squares from the normal equations, without a second pass. + double residual = 0; + for (int r = 0; r < RADIAL_BINS; r++) { + const int c = r * SECTORS + s; + if (!radial_ok[r] || count[c] < MIN_PIXELS_PER_CELL) + continue; + const double e = sum[c] / count[c] - profile[r] - amplitude * profile[r] - this_shift * d_profile[r]; + residual += e * e; + } + const double variance = residual / (n - 2) * (a11 / det); + if (!(variance > 0)) + continue; + shift.push_back(static_cast(this_shift)); + weight.push_back(static_cast(1.0 / variance)); + jacobian_x.push_back(static_cast(jx / n)); + jacobian_y.push_back(static_cast(jy / n)); + } + if (static_cast(shift.size()) < MIN_USABLE_SECTORS) + return {}; + + // shift_k = Jx_k dx + Jy_k dy, robustified so one bad sector cannot carry the answer. + std::vector w = weight; + double c11 = 0, c12 = 0, c22 = 0; + for (int round = 0; round < 3; round++) { + c11 = c12 = c22 = 0; + double r1 = 0, r2 = 0; + for (size_t k = 0; k < shift.size(); k++) { + c11 += w[k] * jacobian_x[k] * jacobian_x[k]; + c12 += w[k] * jacobian_x[k] * jacobian_y[k]; + c22 += w[k] * jacobian_y[k] * jacobian_y[k]; + r1 += w[k] * jacobian_x[k] * shift[k]; + r2 += w[k] * jacobian_y[k] * shift[k]; + } + const double det = c11 * c22 - c12 * c12; + if (det <= 0) + return {}; + step_x = static_cast((c22 * r1 - c12 * r2) / det); + step_y = static_cast((c11 * r2 - c12 * r1) / det); + std::vector residual(shift.size()); + for (size_t k = 0; k < shift.size(); k++) + residual[k] = shift[k] - jacobian_x[k] * step_x - jacobian_y[k] * step_y; + std::vector absolute(residual.size()); + for (size_t k = 0; k < residual.size(); k++) absolute[k] = std::abs(residual[k]); + const float scale = 1.4826f * median_of(absolute) + 1e-30f; + for (size_t k = 0; k < shift.size(); k++) { + const float t = residual[k] / (3 * scale); + w[k] = weight[k] / (1.0f + t * t); + } + } + + double chi2 = 0; + for (size_t k = 0; k < shift.size(); k++) { + const double e = shift[k] - jacobian_x[k] * step_x - jacobian_y[k] * step_y; + chi2 += w[k] * e * e; + } + chi2 = std::max(chi2 / (shift.size() - 2), 1.0); + const double det = c11 * c22 - c12 * c12; + sigma_x = static_cast(std::sqrt(c22 / det * chi2)); + sigma_y = static_cast(std::sqrt(c11 / det * chi2)); + + beam_x += DAMPING * step_x; + beam_y += DAMPING * step_y; + if (std::hypot(step_x, step_y) < CONVERGED_PXL) + break; + } + + return BeamCenterEstimate{beam_x, beam_y, std::max(sigma_x, sigma_y)}; +} diff --git a/image_analysis/geom_refinement/BeamCenterFromBackground.h b/image_analysis/geom_refinement/BeamCenterFromBackground.h new file mode 100644 index 00000000..07ed43cb --- /dev/null +++ b/image_analysis/geom_refinement/BeamCenterFromBackground.h @@ -0,0 +1,37 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +#include "../../common/DiffractionExperiment.h" +#include "../../common/PixelMask.h" + +struct BeamCenterEstimate { + float beam_x_pxl = 0.0f; + float beam_y_pxl = 0.0f; + float sigma_pxl = 0.0f; // 1 sigma on the fitted shift, the larger of the two axes +}; + +// Beam centre from the isotropy of the scattered background, before anything is indexed. +// +// The solvent and air scatter is isotropic in 2-theta about the beam, so a centre that is off +// shifts each azimuthal sector's radial profile by a different amount. Sector k's profile is +// m_k * g(2theta + d_k), with the shift d_k = Jx_k*dx + Jy_k*dy and m_k an amplitude that +// absorbs anything multiplicative and azimuthal - a holder arm, a cryostream shadow, a +// flat-field gradient. Fitting the amplitude alongside the shift is what makes this usable: +// a 50% shadow over one sextant otherwise reads as several tens of pixels of centre error. +// +// The leverage comes from the CURVATURE of the radial profile - the water ring - because for a +// pure exponential decay g' is proportional to g and shift and amplitude are indistinguishable. +// The sigma measures that leverage, so it grows as the curvature weakens, and the caller's gate on +// it is what keeps an ill-determined centre out. It is a precision and not an accuracy: on a +// background with NO curvature at all there is nothing to separate the two parameters, the fit +// follows the noise in g' instead, and it does so confidently. +// +// `mean` is a per-pixel projection over a few tens of frames, NAN where no frame contributed. +std::optional +FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const PixelMask &mask, + const std::vector &mean); diff --git a/image_analysis/geom_refinement/BeamCenterFromSpots.cpp b/image_analysis/geom_refinement/BeamCenterFromSpots.cpp new file mode 100644 index 00000000..6f67cc84 --- /dev/null +++ b/image_analysis/geom_refinement/BeamCenterFromSpots.cpp @@ -0,0 +1,906 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "BeamCenterFromSpots.h" + +#include "../../common/JFJochMath.h" // PI + +#include +#include +#include +#include + +namespace { + +// The vote. A mirrored pair's two coordinates sum to twice the beam's, so the true pairs pile up at +// that value on a broad pedestal of accidental ones; a running mean over VOTE_BACKGROUND_BINS +// measures the pedestal and takes it away. +-60 px about twice the current guess covers any header +// error worth correcting, and 0.2 px bins are fine enough to seed the refinement. +constexpr float VOTE_BIN_PXL = 0.2f; +constexpr float VOTE_REACH_PXL = 60.0f; +constexpr int VOTE_BACKGROUND_BINS = 61; +constexpr int MIN_VOTES = 50; + +// A reciprocal-lattice ROW makes the vote a comb: pairs one lattice step out of register vote at +// 2*beam + n*p, and such a tooth is sometimes TALLER than the true one. So every tooth above this +// fraction of the tallest, and this far from a taller one, is refined, and they are told apart +// afterwards on evidence the height does not carry. +constexpr float CANDIDATE_MIN_HEIGHT = 0.30f; +constexpr float CANDIDATE_SEPARATION_PXL = 6.0f; +constexpr int MAX_CANDIDATES = 6; + +// The Friedel mirror does not touch the other coordinate, so a genuine pair agrees in it: loosely +// while the centre is still unknown (the vote), tightly once it is (the intensity correlation). +constexpr float FRIEDEL_OTHER_TOL_PXL = 3.0f; +constexpr float FRIEDEL_MATCH_WIN_PXL = 1.2f; +constexpr float FRIEDEL_OTHER_COST_WEIGHT = 0.05f; +constexpr float FRIEDEL_TIGHT_PXL = 0.8f; +constexpr int MIN_FRIEDEL_MATCHES = 20; +constexpr int MIN_MATCHES_PER_PAIR = 5; + +// Fewer independent estimates than this - frame pairs for the Friedel fit, frames for the crossing +// - and there is nothing to take a scatter of. +constexpr int MIN_INDEPENDENT_ESTIMATES = 3; + +// A frame at phi can only pair with one at phi+180, so on a sweep of S degrees the pairs' first +// frames span S - 180 however many of them there are - and it is that span, not their number, that +// the second crossing lives on: a reflection's two crossings are an angle apart and the sweep has to +// contain it. Just past half a turn the estimator still forms its full complement of pairs and +// answers from them, and the answer is tens of pixels out on most crystals: at 185 deg of sweep six +// of six committed answers were 2.6-28 px wrong and at 190 deg five of seven, all of it on the +// crossing coordinate while the Friedel one stayed inside 0.4 px; by 200 deg not one of 65 was +// wrong. The fit cannot notice this about itself - every pair agrees with every other - so it is a +// condition on the sweep, tested before anything is fitted. Measured on complete sweeps restricted +// to a span W, the crossing fit as it stood before the timing test above commits gross answers at +// W = 10, 15 and 20 and is clean only from 25; with the timing test it is clean at every W. So this +// is a floor under a guard that is now carried elsewhere, and its value is the largest that costs +// nothing measured: the shortest sweep in the regression set, 199.8 deg, yields a span of 17.6 deg - +// a sweep of S degrees gives S - 180, less one oscillation because the paired sample stops one frame +// short of a partner it would have no room for, less the frames the shutter margin keeps back at +// each end. +constexpr float MIN_PAIR_SPAN_DEG = 17.0f; + +// A false tooth pairs unrelated reflections, and unrelated reflections do not have equal structure +// factors - |F(h)| = |F(-h)| holds for a real Friedel pair and for nothing else - so the intensity +// correlation of the matched set separates the teeth. It is a statistic and needs a population: +// below WEAK_MATCH_COUNT tight matches it is noise, and the count of them is the better evidence. +constexpr int WEAK_MATCH_COUNT = 200; +constexpr int MIN_CORRELATION_MATCHES = 30; + +// The second crossing pairs spots within the whole pool rather than between two known frames, so +// its match has to be unambiguous: exactly one spot within the radius, and the spot must sit far +// enough off the mirror line that the two crossings are genuinely different measurements. +constexpr float CROSSING_OTHER_TOL_PXL = 0.6f; +constexpr float CROSSING_MATCH_RADIUS_PXL = 1.0f; +constexpr float CROSSING_MIN_LEVER_PXL = 25.0f; + +// The crossing's own tooth test. The two crossings are the SAME reflection, so the geometry fixes +// not only where the second one is but WHEN: writing m1 = (m x s0)^ and m3 = (m1 x m)^ - the beam +// projected perpendicular to the spindle - the Laue condition holds q.m and q.m3 fixed along the +// sweep and lets only q.m1 change sign, so the two crossings are 2*atan2(q.m1, q.m3) of sweep apart. +// That angle is read off ONE spot's position, with no cell, no orientation and no indexing, and a +// false pairing has no reason to obey it: over the sampled frames it does so at the accidental rate +// of one frame in the turn. Measured over 36 datasets: true pairs keep it to a median 0.64 deg +// (worst 1.79), and the cut enriches them 29x - which is what tells the teeth apart, since the +// crossing's structure-factor correlation does not (it separates 13 of 17 against the timing's 17). +constexpr float CROSSING_PHI_TOL_DEG = 3.0f; + +// A crossing pair is ONE measurement - the partner search finds it from both ends - so these count +// pairs, not matches. The 60-frame pre-scan yields 40 to 2100 of them (median 370), because both +// crossings of a reflection have to fall on sampled frames, which is quadratic in the budget: 30 +// frames gives a median of 96 and 20 frames a median of 41. +constexpr int MIN_CROSSING_PAIRS = 10; +constexpr int MIN_CROSSING_PAIRS_PER_FRAME = 2; + +// A reflection is recorded over a few consecutive frames, so its two crossings count as one +// measurement unless the frames are this many oscillation widths apart. +constexpr float MIN_EVENT_SEPARATION_WEDGES = 4.0f; + +constexpr int REFINE_ITERATIONS = 6; +constexpr float REFINE_CONVERGED_PXL = 0.002f; + +// Where the answer is checked from. The vote reaches only +-VOTE_REACH/2 px about wherever the +// search starts, so a taller tooth just outside that is invisible from the header but plain from a +// start displaced towards it; and the crossing fit sees the Friedel fit's coordinate, so a start +// elsewhere tests the two against each other as well. Four starts half a window away either come +// back to the same answer - and then it is the crystal's, not the header's - or find a competing +// one, which the estimator has no means of telling apart from it. +constexpr float CONSISTENCY_START_PXL = 25.0f; + +// Persistent artefacts - the beam-stop halo, a hot pixel, the edge of a mask - sit at the SAME +// place on every frame, so they vote at twice their own position with one vote per frame pair, and +// that beats the real peak. A reflection never does: its Friedel mate is mirrored, not coincident. +// So a position that recurs across frames is not diffraction and is dropped. +constexpr float PERSISTENT_RADIUS_PXL = 1.5f; +constexpr float PERSISTENT_FRAME_FRACTION = 0.05f; +constexpr int PERSISTENT_MIN_FRAMES = 4; + +// How far the spindle is searched. Both reaches are a real beamline's worst case with room to +// spare: measured against XDS's refined rotation axis over 38 sweeps the largest azimuth is +// 5.4 mrad and the largest tip 9.2 mrad. The grid is coarse on purpose - it exists to pick the +// TOOTH, and a least-squares fit on that tooth's own members then places the two angles far more +// finely than any grid could. Its step has to keep the tooth visible, which means keeping the +// spread it leaves - 2*step*(detector half-height) for the azimuth, 2*step*D*(1/cos2theta - 1) for +// the tip - to about a pixel. +constexpr float SPINDLE_AZIMUTH_REACH_RAD = 0.010f; +constexpr float SPINDLE_TIP_REACH_RAD = 0.020f; +constexpr float SPINDLE_GRID_STEP_RAD = 0.001f; +constexpr float SPINDLE_PEAK_WINDOW_PXL = 1.5f; +constexpr int SPINDLE_REFINE_ITERATIONS = 4; + +float Coordinate(const BeamCenterSpot &spot, int axis) { + return axis == 0 ? spot.x : spot.y; +} + +float Coordinate(const std::pair &point, int axis) { + return axis == 0 ? point.first : point.second; +} + +float Median(std::vector v) { + if (v.empty()) + return NAN; + const size_t half = v.size() / 2; + std::nth_element(v.begin(), v.begin() + half, v.end()); + return v[half]; +} + +// The scatter of a set, as a median absolute deviation scaled to a standard deviation. +float RobustSpread(const std::vector &v) { + const float centre = Median(v); + std::vector deviation(v.size()); + for (size_t i = 0; i < v.size(); i++) + deviation[i] = std::abs(v[i] - centre); + return 1.4826f * Median(deviation); +} + +float BeamCoordinate(const DiffractionGeometry &geom, int axis) { + return axis == 0 ? geom.GetBeamX_pxl() : geom.GetBeamY_pxl(); +} + +void SetBeamCoordinate(DiffractionGeometry &geom, int axis, float value) { + if (axis == 0) + geom.BeamX_pxl(value); + else + geom.BeamY_pxl(value); +} + +// The exact lab-space mirror. Kept as an object because the trial centre changes far less often +// than the spots it is applied to. +class LabMirror { + const DiffractionGeometry &geom; + const RotMatrix inverse_rotation; + const int axis; +public: + LabMirror(const DiffractionGeometry &geometry, int mirror_axis) + : geom(geometry), inverse_rotation(geometry.GetPoniRotMatrix().transpose()), axis(mirror_axis) {} + + // The image of a detector point: negate the mirrored lab component of the ray to it and project + // the result back onto the detector. The reflecting plane contains the beam, so a point's image + // sits opposite it about the DIRECT BEAM - not about the PONI, which the detector rotations put + // up to (distance/pixel)*rot pixels away. + [[nodiscard]] std::pair operator()(float x, float y) const { + Coord lab = geom.LabCoord(x, y); + lab[axis] = -lab[axis]; + const Coord ray = inverse_rotation * lab; + if (!(ray.z > 0)) + return {NAN, NAN}; + const float scale = geom.GetDetectorDistance_mm() / (ray.z * geom.GetPixelSize_mm()); + return {geom.GetBeamX_pxl() + ray.x * scale, geom.GetBeamY_pxl() + ray.y * scale}; + } +}; + +// The distance from the PONI to the direct beam on one axis. The vote is over detector coordinates +// mirrored about the direct beam, and the fit moves the PONI, so this is what converts one to the +// other. It does not depend on where the beam is, only on the detector rotations. +float DirectBeamOffset(const DiffractionGeometry &geom, int axis) { + const auto direct = geom.GetDirectBeam_pxl(); + return Coordinate(direct, axis) - BeamCoordinate(geom, axis); +} + +// The vote histogram with its pedestal taken away, bin by bin. +std::vector VoteExcess(const std::vector &sums, float guess) { + const int n_bins = static_cast(2 * VOTE_REACH_PXL / VOTE_BIN_PXL); + const float first_bin = 2 * guess - VOTE_REACH_PXL; + std::vector histogram(n_bins, 0.0f); + for (const float sum: sums) { + // Floor, not truncation: a sum just below the first bin divides to a small negative number, + // which truncates to zero and would pile every underflow into bin 0 - a tooth of its own, + // and on real data a tall one. + const int bin = static_cast(std::floor((sum - first_bin) / VOTE_BIN_PXL)); + if (bin >= 0 && bin < n_bins) + histogram[bin] += 1.0f; + } + + std::vector excess(n_bins); + std::vector cumulative(n_bins + 1, 0.0); + for (int i = 0; i < n_bins; i++) + cumulative[i + 1] = cumulative[i] + histogram[i]; + for (int i = 0; i < n_bins; i++) { + const int lo = std::max(0, i - VOTE_BACKGROUND_BINS / 2); + const int hi = std::min(n_bins, i + VOTE_BACKGROUND_BINS / 2 + 1); + excess[i] = histogram[i] - static_cast((cumulative[hi] - cumulative[lo]) / (hi - lo)); + } + return excess; +} + +// The values of 2*beam that the vote supports, tallest first: every bin whose excess over the +// pedestal is within CANDIDATE_MIN_HEIGHT of the tallest one's, and CANDIDATE_SEPARATION_PXL away +// from an already accepted candidate. +std::vector VoteCandidates(const std::vector &sums, float guess) { + const int n_bins = static_cast(2 * VOTE_REACH_PXL / VOTE_BIN_PXL); + const float first_bin = 2 * guess - VOTE_REACH_PXL; + const std::vector excess = VoteExcess(sums, guess); + + std::vector order(n_bins); + std::iota(order.begin(), order.end(), 0); + std::ranges::sort(order, [&](int a, int b) { return excess[a] > excess[b]; }); + + std::vector candidates; + if (!(excess[order.front()] > 0.0f)) + return candidates; + const float threshold = CANDIDATE_MIN_HEIGHT * excess[order.front()]; + for (const int bin: order) { + if (excess[bin] < threshold || static_cast(candidates.size()) >= MAX_CANDIDATES) + break; + const float value = first_bin + (bin + 0.5f) * VOTE_BIN_PXL; + if (std::ranges::none_of(candidates, [&](float taken) { + return std::abs(value - taken) < CANDIDATE_SEPARATION_PXL; })) + candidates.push_back(value); + } + return candidates; +} + +// One matched Friedel pair, in the frame where the mirror has already been applied. +struct FriedelMatch { + float residual; // the partner minus the mirror image, in the coordinate the mirror flips + float other; // and in the coordinate it leaves alone + float log_intensity_a; + float log_intensity_b; +}; + +// Mutual nearest neighbours between one frame's spots, mirrored, and its partner's. A spot may be +// matched only if it is its partner's best candidate as well, so a dense frame cannot pile several +// spots onto one. +void MatchFriedelFrames(const LabMirror &mirror, int axis, + const std::vector &a, const std::vector &b, + std::vector &out) { + const int other = 1 - axis; + std::vector best_of_a(a.size(), -1), best_of_b(b.size(), -1); + std::vector cost_of_a(a.size(), INFINITY), cost_of_b(b.size(), INFINITY); + + for (size_t i = 0; i < a.size(); i++) { + const auto image = mirror(a[i].x, a[i].y); + if (!std::isfinite(image.first)) + continue; + for (size_t j = 0; j < b.size(); j++) { + const float along = Coordinate(image, axis) - Coordinate(b[j], axis); + const float across = Coordinate(image, other) - Coordinate(b[j], other); + if (std::abs(along) >= FRIEDEL_MATCH_WIN_PXL || std::abs(across) >= FRIEDEL_OTHER_TOL_PXL) + continue; + const float cost = std::abs(along) + FRIEDEL_OTHER_COST_WEIGHT * std::abs(across); + if (cost < cost_of_a[i]) { cost_of_a[i] = cost; best_of_a[i] = static_cast(j); } + if (cost < cost_of_b[j]) { cost_of_b[j] = cost; best_of_b[j] = static_cast(i); } + } + } + + for (size_t i = 0; i < a.size(); i++) { + const int j = best_of_a[i]; + if (j < 0 || best_of_b[j] != static_cast(i)) + continue; + const auto image = mirror(a[i].x, a[i].y); + out.push_back({Coordinate(b[j], axis) - Coordinate(image, axis), + Coordinate(b[j], other) - Coordinate(image, other), + std::log(std::max(a[i].intensity, 1.0f)), + std::log(std::max(b[j].intensity, 1.0f))}); + } +} + +struct FriedelCandidate { + float poni = NAN; + float sigma = NAN; + int matches = 0; + int tight_matches = 0; + float correlation = -1.0f; +}; + +// Refine one tooth of the vote: mirror every frame onto its partner, and move the centre by half +// the median residual until it stops moving. Half, because a centre that is off by e puts the +// mirror image 2e away from the spot it belongs to. +std::optional RefineFriedel(DiffractionGeometry geom, int axis, + const std::vector> &by_frame, + const std::vector> &pairs, + float start_direct_beam) { + SetBeamCoordinate(geom, axis, start_direct_beam - DirectBeamOffset(geom, axis)); + + for (int iteration = 0; iteration < REFINE_ITERATIONS; iteration++) { + const LabMirror mirror(geom, axis); + std::vector matches; + for (const auto &[first, second]: pairs) + MatchFriedelFrames(mirror, axis, by_frame[first], by_frame[second], matches); + if (static_cast(matches.size()) < MIN_FRIEDEL_MATCHES) + return {}; + + std::vector residual(matches.size()); + for (size_t i = 0; i < matches.size(); i++) + residual[i] = matches[i].residual; + const float step = Median(residual) / 2.0f; + SetBeamCoordinate(geom, axis, BeamCoordinate(geom, axis) + step); + if (std::abs(step) < REFINE_CONVERGED_PXL) + break; + } + + // The converged centre, and what each frame pair says about it on its own. The scatter of those + // is the honest uncertainty: it is what a different pair of frames would have said, not how + // finely their common median is determined, which is smaller by the square root of their number + // and is a precision rather than an accuracy. + const LabMirror mirror(geom, axis); + std::vector matches; + std::vector per_pair; + for (const auto &[first, second]: pairs) { + const size_t before = matches.size(); + MatchFriedelFrames(mirror, axis, by_frame[first], by_frame[second], matches); + if (static_cast(matches.size() - before) < MIN_MATCHES_PER_PAIR) + continue; + std::vector residual; + for (size_t i = before; i < matches.size(); i++) + residual.push_back(matches[i].residual); + per_pair.push_back(Median(residual) / 2.0f); + } + if (static_cast(matches.size()) < MIN_FRIEDEL_MATCHES + || static_cast(per_pair.size()) < MIN_INDEPENDENT_ESTIMATES) + return {}; + + // A genuine Friedel pair sits at a constant offset in the coordinate the mirror does not touch; + // a lattice-shifted one does not, so the correlation is measured over the ones that do. + std::vector across(matches.size()); + for (size_t i = 0; i < matches.size(); i++) + across[i] = matches[i].other; + const float centre = Median(across); + + std::vector log_a, log_b; + for (const auto &m: matches) + if (std::abs(m.other - centre) < FRIEDEL_TIGHT_PXL) { + log_a.push_back(m.log_intensity_a); + log_b.push_back(m.log_intensity_b); + } + + FriedelCandidate candidate; + candidate.poni = BeamCoordinate(geom, axis); + candidate.sigma = RobustSpread(per_pair); + candidate.matches = static_cast(matches.size()); + candidate.tight_matches = static_cast(log_a.size()); + if (candidate.tight_matches > MIN_CORRELATION_MATCHES) { + const float mean_a = std::accumulate(log_a.begin(), log_a.end(), 0.0f) / log_a.size(); + const float mean_b = std::accumulate(log_b.begin(), log_b.end(), 0.0f) / log_b.size(); + double covariance = 0, variance_a = 0, variance_b = 0; + for (size_t i = 0; i < log_a.size(); i++) { + covariance += (log_a[i] - mean_a) * (log_b[i] - mean_b); + variance_a += (log_a[i] - mean_a) * (log_a[i] - mean_a); + variance_b += (log_b[i] - mean_b) * (log_b[i] - mean_b); + } + if (variance_a > 0 && variance_b > 0) + candidate.correlation = static_cast(covariance / std::sqrt(variance_a * variance_b)); + } + return candidate; +} + +std::optional FitFriedel(const DiffractionGeometry &geom, int axis, + const std::vector> &by_frame, + const std::vector> &pairs) { + const int other = 1 - axis; + std::vector sums; + for (const auto &[first, second]: pairs) + for (const auto &a: by_frame[first]) + for (const auto &b: by_frame[second]) + if (std::abs(Coordinate(a, other) - Coordinate(b, other)) < FRIEDEL_OTHER_TOL_PXL) + sums.push_back(Coordinate(a, axis) + Coordinate(b, axis)); + if (static_cast(sums.size()) < MIN_VOTES) + return {}; + + std::vector refined; + for (const float candidate: VoteCandidates(sums, Coordinate(geom.GetDirectBeam_pxl(), axis))) + if (const auto fit = RefineFriedel(geom, axis, by_frame, pairs, candidate / 2.0f)) + refined.push_back(*fit); + if (refined.empty()) + return {}; + + // Which tooth is the real one. Normally the structure-factor correlation says so; where there + // are too few tight matches for it to mean anything it is noise, and the tooth that matched the + // most spots is the best evidence available. + const bool weak = std::ranges::max_element(refined, {}, &FriedelCandidate::tight_matches) + ->tight_matches < WEAK_MATCH_COUNT; + if (weak) + return *std::ranges::max_element(refined, {}, &FriedelCandidate::tight_matches); + return *std::ranges::max_element(refined, {}, &FriedelCandidate::correlation); +} + +// How far along the sweep a spot's second crossing lies, for every spot. Kabsch's frame: m2 is the +// spindle, m3 the beam with its spindle component taken out, m1 the third. Along the sweep q.m2 and +// |q| never change and the Laue condition fixes q.m3, so a reflection reaches the sphere exactly +// twice, at azimuths +-atan2(q.m1, q.m3) about m3, and the sweep between them is twice that. +std::vector CrossingSeparation_deg(const DiffractionGeometry &geom, const Coord &spindle, + const std::vector &spots) { + const Coord m2 = spindle.Normalize(); + const Coord m1 = (m2 % geom.GetScatteringVector()).Normalize(); + const Coord m3 = (m1 % m2).Normalize(); + + std::vector separation(spots.size()); + for (size_t i = 0; i < spots.size(); i++) { + const Coord q = geom.DetectorToRecip(spots[i].x, spots[i].y); + separation[i] = 2.0f * std::atan2(q * m1, q * m3) * 180.0f / static_cast(PI); + } + return separation; +} + +// One tooth of the crossing vote, refined, and how much of the pool holds it up. +struct CrossingCandidate { + float poni = NAN; + float sigma = NAN; + int pairs = 0; +}; + +// The other coordinate, from the second crossing. Both crossings of a reflection are somewhere in +// the pool rather than on two frames known in advance, so the match is over the whole pool: the +// spots are indexed on the coordinate the mirror leaves alone, and a window on it holds every +// candidate partner. +std::optional FitCrossing(DiffractionGeometry geom, int axis, + const Coord &spindle, + const std::vector &spots, + const std::vector &frame_angle_deg, + float min_separation_deg) { + const int other = 1 - axis; + std::vector order(spots.size()); + std::iota(order.begin(), order.end(), 0); + std::ranges::sort(order, [&](int a, int b) { + return Coordinate(spots[a], other) < Coordinate(spots[b], other); }); + std::vector across(spots.size()); + for (size_t i = 0; i < order.size(); i++) + across[i] = Coordinate(spots[order[i]], other); + + const std::vector separation = CrossingSeparation_deg(geom, spindle, spots); + + // Two spots can be the two crossings of one reflection when they are far enough apart in the + // sweep to be two events at all, and when the sweep between them is the one the first spot's + // own position asks for. + const auto crossing_pair = [&](int a, int b) { + const float delta = frame_angle_deg[spots[b].frame] - frame_angle_deg[spots[a].frame]; + return std::abs(delta) > min_separation_deg + && std::abs(std::remainder(delta - separation[a], 360.0f)) < CROSSING_PHI_TOL_DEG; + }; + + std::vector sums; + for (size_t i = 0; i < order.size(); i++) { + const auto last = std::upper_bound(across.begin() + i + 1, across.end(), + across[i] + CROSSING_OTHER_TOL_PXL); + for (auto it = across.begin() + i + 1; it != last; ++it) { + const int j = order[it - across.begin()]; + if (crossing_pair(order[i], j)) + sums.push_back(Coordinate(spots[order[i]], axis) + Coordinate(spots[j], axis)); + } + } + if (static_cast(sums.size()) < MIN_VOTES) + return {}; + + // Refine one tooth: mirror every spot, take the partner the mirror lands on, and move the centre + // by half the median residual until it stops moving. + const auto refine = [&](float start_2beam) -> std::optional { + DiffractionGeometry trial = geom; + SetBeamCoordinate(trial, axis, start_2beam / 2.0f - DirectBeamOffset(trial, axis)); + + std::vector residual; + std::vector frame; + for (int iteration = 0; iteration <= REFINE_ITERATIONS; iteration++) { + const LabMirror mirror(trial, axis); + residual.clear(); + frame.clear(); + for (size_t q = 0; q < spots.size(); q++) { + const auto image = mirror(spots[q].x, spots[q].y); + if (!std::isfinite(image.first)) + continue; + // Both crossings measure the same |q|, so a spot close to the mirror line is its own + // image and says nothing about where that line is. + if (std::abs(Coordinate(spots[q], axis) - Coordinate(image, axis)) / 2.0f < CROSSING_MIN_LEVER_PXL) + continue; + const float key = Coordinate(image, other); + const auto first = std::lower_bound(across.begin(), across.end(), key - CROSSING_MATCH_RADIUS_PXL); + const auto last = std::upper_bound(across.begin(), across.end(), key + CROSSING_MATCH_RADIUS_PXL); + int found = -1; + int count = 0; + for (auto it = first; it != last; ++it) { + const int j = order[it - across.begin()]; + if (!crossing_pair(static_cast(q), j)) + continue; + if (std::hypot(Coordinate(spots[j], axis) - Coordinate(image, axis), + Coordinate(spots[j], other) - key) > CROSSING_MATCH_RADIUS_PXL) + continue; + found = j; + count++; + } + // The mirror is its own inverse, so the pair is found from both ends and gives the + // same residual twice. It is one measurement and is counted once. + if (count != 1 || found < static_cast(q)) + continue; + residual.push_back(Coordinate(spots[found], axis) - Coordinate(image, axis)); + frame.push_back(spots[q].frame); + } + if (static_cast(residual.size()) < MIN_CROSSING_PAIRS) + return {}; + if (iteration == REFINE_ITERATIONS) + break; + + const float step = Median(residual) / 2.0f; + SetBeamCoordinate(trial, axis, BeamCoordinate(trial, axis) + step); + if (std::abs(step) < REFINE_CONVERGED_PXL) + break; + } + + // As for the Friedel fit: the scatter of what the individual frames say, not of their mean. + std::vector> of_frame(frame_angle_deg.size()); + for (size_t i = 0; i < residual.size(); i++) + of_frame[frame[i]].push_back(residual[i] / 2.0f); + std::vector per_frame; + for (const auto &one: of_frame) + if (static_cast(one.size()) >= MIN_CROSSING_PAIRS_PER_FRAME) + per_frame.push_back(Median(one)); + if (static_cast(per_frame.size()) < MIN_INDEPENDENT_ESTIMATES) + return {}; + return CrossingCandidate{BeamCoordinate(trial, axis), RobustSpread(per_frame), + static_cast(residual.size())}; + }; + + // Which tooth is the real one. A false tooth pairs unrelated reflections, and unrelated + // reflections do not keep the sweep angle their positions ask for, so the pairs that survive + // the timing test are the evidence the height does not carry - the crossing's counterpart of + // the Friedel side's |F(h)| = |F(-h)|. + std::optional best; + for (const float candidate: VoteCandidates(sums, Coordinate(geom.GetDirectBeam_pxl(), axis))) + if (const auto fit = refine(candidate)) + if (!best || fit->pairs > best->pairs) + best = fit; + return best; +} + +// Drop the spots that are not diffraction: a position that appears on frame after frame. +std::vector DropPersistentSpots(const std::vector &spots, int n_frames) { + std::vector order(spots.size()); + std::iota(order.begin(), order.end(), 0); + std::ranges::sort(order, [&](int a, int b) { return spots[a].x < spots[b].x; }); + std::vector sorted_x(spots.size()); + for (size_t i = 0; i < order.size(); i++) + sorted_x[i] = spots[order[i]].x; + + const int limit = std::max(PERSISTENT_MIN_FRAMES, + static_cast(PERSISTENT_FRAME_FRACTION * n_frames)); + std::vector kept; + std::vector frames; + for (size_t i = 0; i < order.size(); i++) { + const BeamCenterSpot &spot = spots[order[i]]; + const auto first = std::lower_bound(sorted_x.begin(), sorted_x.end(), spot.x - PERSISTENT_RADIUS_PXL); + const auto last = std::upper_bound(sorted_x.begin(), sorted_x.end(), spot.x + PERSISTENT_RADIUS_PXL); + frames.clear(); + for (auto it = first; it != last; ++it) { + const BeamCenterSpot &other = spots[order[it - sorted_x.begin()]]; + if (std::hypot(other.x - spot.x, other.y - spot.y) <= PERSISTENT_RADIUS_PXL) + frames.push_back(other.frame); + } + std::ranges::sort(frames); + if (std::unique(frames.begin(), frames.end()) - frames.begin() <= limit) + kept.push_back(spot); + } + return kept; +} + +// THE SPINDLE IS NEVER PERPENDICULAR TO THE BEAM, AND THE FILE NEVER SAYS SO +// +// Everything above assumes the spindle m satisfies m.s0 = 0 and lies on a lab axis. Neither holds +// at a beamline, and the header cannot be asked: every master of the regression set writes the axis +// as exactly (-1,0,0) or (0,1,0), while XDS's refined axis for the same 38 sweeps departs from it +// by up to 5.4 mrad in the plane and 9.2 mrad out of it. So the deviation has to be measured here +// or not at all. It splits into two components that behave completely differently. +// +// AZIMUTH - the spindle turned about the beam by `a`. Both mirror planes turn with it (the +// Friedel one is normal to the spindle, the crossing one contains the spindle and the beam), so a +// mirror taken about the nominal lab axis leaves, along the mirrored coordinate, a residual of +// exactly 2*a*(the other coordinate). It does not bias the fit - that residual is odd in the other +// coordinate, so its median over a symmetric set is zero - it SMEARS the vote, by up to +-3 px at +// 2 mrad across a 1500 px detector, until a neighbouring comb tooth outvotes the true one. On a +// synthetic sweep that is 9.1 px of centre error at 2 mrad and 17.7 px at 5 mrad, reported at a +// sigma of 0.09-0.14 - a silent, confident, wrong answer. This is what is fitted below. +// +// TILT OUT OF THE PLANE - the spindle tipped towards the beam by `e`. Then the Friedel mate is no +// longer in diffracting condition at phi+180: with q' = q - 2(q.m)m the Ewald residual is +// -2(q.m)(m.s0), and the mate diffracts at phi+180+dphi with dphi = 2(q.m)e / (q.(m x s0_hat)). +// Carrying that rotation through to the detector leaves, along the mirrored coordinate, exactly +// +// 2 * e * D_pxl * (1/cos(2theta) - 1) +// +// (verified against the forward geometry to 1e-4 px), a purely RADIAL term - and the crossing axis +// does not see it at all, because its mirror plane's normal is m x s0, which is independent of the +// component of m along the beam. Uncorrected it costs the Friedel coordinate one tip times the +// median of that lever over the matched pairs: on the regression set a median 0.024 px per mrad, +// worst 0.060, i.e. a median 0.03 px and a worst case of 0.24 px. It is fitted below but NOT +// applied; the note at the call site says why. +// +// One thing this must NOT do is mirror about the spindle itself. The plane normal to m does not +// contain the beam once e is non-zero, so mirroring a ray about it moves the direct beam bodily by +// 2*e*D_pxl - 4.8 px, i.e. a 2.4 px centre error, at e = 2 mrad and 2400 px of distance. Only the +// azimuth of the mirror plane follows the spindle; the plane itself always contains the beam. That +// is why the correction below is a rotation of the whole problem about the beam and nothing else. + +// A point turned about a centre, used to bring the sweep into the frame where the spindle does lie +// on a lab axis - which is the frame the mirror and both votes are written for. +std::pair Turn(float x, float y, float cx, float cy, float sin_a, float cos_a) { + const float dx = x - cx, dy = y - cy; + return {cx + dx * cos_a - dy * sin_a, cy + dx * sin_a + dy * cos_a}; +} + +// How far the tip moves the mirror image of a spot, per radian of tip - half the pair residual +// above, which is what a shift of the spot itself has to be for the two to cancel. +float TipLever_pxl(float radius_pxl, float distance_pxl) { + return std::hypot(radius_pxl, distance_pxl) - distance_pxl; +} + +// One candidate Friedel pair, reduced to what the vote needs at any spindle: the sum of the two +// positions, their difference, and the two spots' tip levers together. +struct SpindlePair { + float sum_x, sum_y; + float diff_x, diff_y; + float tip_lever; +}; + +// The spindle, from the Friedel vote. +// +// Neither component moves the true tooth, both spread it - the azimuth over 2*a*(detector height), +// the tip over 2*e*D*(1/cos2theta - 1) - so the spindle the spots were taken at is the one that +// makes the vote tallest. That is how the tooth is found; where it is found is then a least-squares +// question, because on the tooth's own members the vote value is linear in both angles: +// +// value = 2*beam_along + (azimuth error)*(the pair's sum ACROSS the mirror line) +// - (tip error)*(the pair's tip lever) +// +// and those two columns are well separated - the first is odd across the mirror line, the second is +// even and radial. What is NOT in the fit is the third column a wrong detector rot1 would need, the +// mirrored coordinate squared; it is separable from the tip in principle (the two correlate at 0.66 +// over a real spot distribution, condition number 6.1) but it is not separated here, which is the +// whole reason the tip is reported rather than used. +// +// A grid rather than a descent, because the two are not separable when both are far out: a 9 mrad +// tip - which the regression set has - spreads the tooth enough on its own that a spurious azimuth +// sharpens it, and a descent that meets the azimuth first never leaves that minimum. +// +// The pairs are collected once with a window wide enough for the whole search - the vote's own +// window, opened by the largest displacement an azimuth in range can produce - and every trial is +// then one pass over that list. +std::pair FitSpindle(const std::vector> &by_frame, + const std::vector> &pairs, + int axis, float guess_x, float guess_y, float distance_pxl, + float &excess_at_fit, float &excess_at_nominal) { + const int other = 1 - axis; + std::vector candidates; + for (const auto &[first, second]: pairs) + for (const auto &a: by_frame[first]) + for (const auto &b: by_frame[second]) { + const float delta_other = Coordinate(a, other) - Coordinate(b, other); + const float delta_axis = Coordinate(a, axis) - Coordinate(b, axis); + if (std::abs(delta_other) + >= FRIEDEL_OTHER_TOL_PXL + SPINDLE_AZIMUTH_REACH_RAD * std::abs(delta_axis)) + continue; + candidates.push_back( + {a.x + b.x, a.y + b.y, a.x - b.x, a.y - b.y, + TipLever_pxl(std::hypot(a.x - guess_x, a.y - guess_y), distance_pxl) + + TipLever_pxl(std::hypot(b.x - guess_x, b.y - guess_y), distance_pxl)}); + } + + const float guess = axis == 0 ? guess_x : guess_y; + const float centre_x = 2 * guess_x, centre_y = 2 * guess_y; + + // One trial: the pairs that still look like Friedel pairs at this spindle, and what each of + // them votes for. `along` is the vote value, `across` the sum along the mirror line - the + // azimuth's column in the fit below. + std::vector along, across, lever; + const auto evaluate = [&](float trial_azimuth, float trial_tip) { + const float sin_a = std::sin(trial_azimuth), cos_a = std::cos(trial_azimuth); + along.clear(); across.clear(); lever.clear(); + for (const auto &c: candidates) { + const float difference = other == 0 ? c.diff_x * cos_a + c.diff_y * sin_a + : -c.diff_x * sin_a + c.diff_y * cos_a; + if (std::abs(difference) >= FRIEDEL_OTHER_TOL_PXL) + continue; + const auto turned = Turn(c.sum_x, c.sum_y, centre_x, centre_y, -sin_a, cos_a); + along.push_back(Coordinate(turned, axis) - trial_tip * c.tip_lever); + across.push_back(Coordinate(turned, other)); + lever.push_back(c.tip_lever); + } + }; + const auto excess_of = [&](float trial_azimuth, float trial_tip) { + evaluate(trial_azimuth, trial_tip); + if (static_cast(along.size()) < MIN_VOTES) + return -INFINITY; + const std::vector excess = VoteExcess(along, guess); + return *std::ranges::max_element(excess); + }; + + excess_at_nominal = excess_of(0.0f, 0.0f); + float azimuth = 0.0f, tip = 0.0f, best = excess_at_nominal; + const int azimuth_steps = static_cast(SPINDLE_AZIMUTH_REACH_RAD / SPINDLE_GRID_STEP_RAD); + const int tip_steps = static_cast(SPINDLE_TIP_REACH_RAD / SPINDLE_GRID_STEP_RAD); + for (int i = -azimuth_steps; i <= azimuth_steps; i++) + for (int j = -tip_steps; j <= tip_steps; j++) { + const float trial_azimuth = i * SPINDLE_GRID_STEP_RAD; + const float trial_tip = j * SPINDLE_GRID_STEP_RAD; + const float excess = excess_of(trial_azimuth, trial_tip); + if (excess > best) { best = excess; azimuth = trial_azimuth; tip = trial_tip; } + } + excess_at_fit = best; + + // The tooth is now known to a grid step; where it sits is a straight three-parameter fit over + // its own members, iterated because moving the angles changes which pairs land on it. + for (int iteration = 0; iteration < SPINDLE_REFINE_ITERATIONS; iteration++) { + evaluate(azimuth, tip); + if (static_cast(along.size()) < MIN_VOTES) + return {azimuth, tip}; + const float peak = Median(along); + double n = 0, s_a = 0, s_l = 0, s_aa = 0, s_al = 0, s_ll = 0, s_y = 0, s_ay = 0, s_ly = 0; + for (size_t i = 0; i < along.size(); i++) { + const double y = along[i] - peak; + if (std::abs(y) >= SPINDLE_PEAK_WINDOW_PXL) + continue; + const double a = across[i], l = lever[i]; + n += 1; s_a += a; s_l += l; s_aa += a * a; s_al += a * l; s_ll += l * l; + s_y += y; s_ay += a * y; s_ly += l * y; + } + if (n < MIN_VOTES) + return {azimuth, tip}; + // The 2x2 system in (azimuth error, tip error) after the constant is projected out. + const double caa = s_aa - s_a * s_a / n, cal = s_al - s_a * s_l / n, cll = s_ll - s_l * s_l / n; + const double cay = s_ay - s_a * s_y / n, cly = s_ly - s_l * s_y / n; + const double determinant = caa * cll - cal * cal; + if (!(std::abs(determinant) > 0)) + return {azimuth, tip}; + azimuth -= static_cast((cay * cll - cly * cal) / determinant); + tip += static_cast((cly * caa - cay * cal) / determinant); + if (std::abs(azimuth) > SPINDLE_AZIMUTH_REACH_RAD || std::abs(tip) > SPINDLE_TIP_REACH_RAD) + return {0.0f, 0.0f}; + } + return {azimuth, tip}; +} + +// One run of the two fits, from wherever `geom` says the beam is. With `spindle_estimate` non-null +// the spindle is fitted here as well, and the sweep turned into its frame before anything else - +// so a run started elsewhere re-fits it, and a spindle that depends on where the search began shows +// up in the answer's spread like everything else. +std::optional +Estimate(DiffractionGeometry geom, const GoniometerAxis &goniometer, + const std::vector &frame_angle_deg, const std::vector &spots, + SpindleEstimate *spindle_estimate) { + // Which detector coordinate the Friedel mirror flips: the one the spindle lies along. The lab + // frame's x and y are the detector's own, so the spindle's larger lab component names the axis + // - a vertical spindle simply swaps the two estimators over. + const Coord spindle = goniometer.GetAxis().Normalize(); + const int friedel_axis = std::abs(spindle.x) >= std::abs(spindle.y) ? 0 : 1; + + // Frames half a turn apart. The pairing has to be as exact as the frames allow: a mate that + // sits half an oscillation away is a different reflection. + const float wedge = std::max(std::abs(goniometer.GetWedge_deg()), 1e-3f); + std::vector> pairs; + for (size_t i = 0; i < frame_angle_deg.size(); i++) + for (size_t j = i + 1; j < frame_angle_deg.size(); j++) + if (std::abs(std::abs(frame_angle_deg[j] - frame_angle_deg[i]) - 180.0f) < wedge) + pairs.emplace_back(i, j); + if (static_cast(pairs.size()) < MIN_INDEPENDENT_ESTIMATES) + return {}; + + // The earlier frame of each pair, taken by ANGLE and not by position: the frames arrive sorted + // on the first pass, but a second pass that reads more of the sweep appends them, and an index + // order that is not an angle order would make this span the whole turn and the test vacuous. + float earliest = INFINITY, latest = -INFINITY; + for (const auto &pair: pairs) { + const float first = std::min(frame_angle_deg[pair.first], frame_angle_deg[pair.second]); + earliest = std::min(earliest, first); + latest = std::max(latest, first); + } + if (latest - earliest < MIN_PAIR_SPAN_DEG) + return {}; + + auto kept = DropPersistentSpots(spots, static_cast(frame_angle_deg.size())); + std::vector> by_frame(frame_angle_deg.size()); + for (const auto &spot: kept) + by_frame[spot.frame].push_back(spot); + + // The spindle, and the sweep turned by its azimuth about the beam. In that frame the spindle + // does lie on a lab axis, so everything below is the estimator as written, and the answer is + // turned back out of it at the end. Turning the SPOTS rather than the mirror is what carries the + // correction into the two VOTES as well - they work in raw detector coordinates, where a mirror + // line that is not parallel to a pixel axis breaks them just as thoroughly as it breaks the + // mirror. + // + // The TIP is fitted alongside and then NOT applied. It has to be in the fit - left out, the + // 9 mrad the regression set has spreads the vote enough on its own that a spurious azimuth + // sharpens it - but it is a nuisance parameter, not a correction, and for two reasons. Its + // column is radial, and so is what a wrong detector rot1 leaves behind, so the fitted tip is + // the spindle's plus about 0.9 of the header's tilt error and is not the spindle alone + // (measured over 38 sweeps: slope 1.03 on XDS's refined axis, offset +4.0 mrad). And its + // scatter, about 2.6 mrad, times the lever it acts on is larger than the 0.07 px it would + // correct - applied, it moves the centre a median 0.08 px the wrong way. + const auto pivot = geom.GetDirectBeam_pxl(); + const float distance_pxl = geom.GetDetectorDistance_mm() / geom.GetPixelSize_mm(); + float azimuth = 0.0f; + if (spindle_estimate) { + float excess = NAN, excess_nominal = NAN, tip = NAN; + std::tie(azimuth, tip) = FitSpindle(by_frame, pairs, friedel_axis, pivot.first, pivot.second, + distance_pxl, excess, excess_nominal); + *spindle_estimate = SpindleEstimate{azimuth, tip, excess, excess_nominal}; + } + if (azimuth != 0.0f) { + const float sin_a = std::sin(azimuth), cos_a = std::cos(azimuth); + for (auto &spot: kept) + std::tie(spot.x, spot.y) = Turn(spot.x, spot.y, pivot.first, pivot.second, -sin_a, cos_a); + for (auto &frame: by_frame) + frame.clear(); + for (const auto &spot: kept) + by_frame[spot.frame].push_back(spot); + } + + const auto friedel = FitFriedel(geom, friedel_axis, by_frame, pairs); + if (!friedel) + return {}; + SetBeamCoordinate(geom, friedel_axis, friedel->poni); + + const auto crossing = FitCrossing(geom, 1 - friedel_axis, spindle, kept, frame_angle_deg, + MIN_EVENT_SEPARATION_WEDGES * wedge); + if (!crossing) + return {}; + SetBeamCoordinate(geom, 1 - friedel_axis, crossing->poni); + + if (azimuth != 0.0f) { + // Back out of the spindle frame. The turn is about the beam, so it acts on the DIRECT beam + // and the PONI follows from it: the offset between the two is a property of the detector + // rotations alone and does not move. + const float sin_a = std::sin(azimuth), cos_a = std::cos(azimuth); + const auto direct = geom.GetDirectBeam_pxl(); + const auto turned = Turn(direct.first, direct.second, pivot.first, pivot.second, sin_a, cos_a); + geom.BeamX_pxl(turned.first - DirectBeamOffset(geom, 0)); + geom.BeamY_pxl(turned.second - DirectBeamOffset(geom, 1)); + } + + return BeamCenterEstimate{geom.GetBeamX_pxl(), geom.GetBeamY_pxl(), + std::max(friedel->sigma, crossing->sigma)}; +} + +} // namespace + +std::optional +FindBeamCenterFromSpotSymmetry(const DiffractionExperiment &experiment, + const std::vector &frame_angle_deg, + const std::vector &spots, + SpindleEstimate *spindle_estimate) { + const auto goniometer = experiment.GetGoniometer(); + if (!goniometer || spots.empty()) + return {}; + + const auto geom = experiment.GetDiffractionGeometry(); + auto estimate = Estimate(geom, *goniometer, frame_angle_deg, spots, spindle_estimate); + if (!estimate) + return {}; + + // How far the answer moves when the search is started somewhere else. On a sound measurement it + // does not move at all; where it does, the two answers are both consistent with the spots and + // 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. + float spread = 0.0f; + for (const auto &[dx, dy]: {std::pair{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)); + } + estimate->sigma_pxl = std::max(estimate->sigma_pxl, spread); + return estimate; +} diff --git a/image_analysis/geom_refinement/BeamCenterFromSpots.h b/image_analysis/geom_refinement/BeamCenterFromSpots.h new file mode 100644 index 00000000..25d8482c --- /dev/null +++ b/image_analysis/geom_refinement/BeamCenterFromSpots.h @@ -0,0 +1,77 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include + +#include "BeamCenterFromBackground.h" +#include "../../common/DiffractionExperiment.h" + +// One spot of the pre-scan sample: its position in image pixels, its photon count, and which of the +// sampled frames it came from (an index into the angles the estimator is given). +struct BeamCenterSpot { + float x = 0.0f; + float y = 0.0f; + float intensity = 0.0f; + int frame = 0; +}; + +// Where the spindle really points. The goniometer axis in the file is a nominal - every master of +// the regression set writes it as exactly a lab axis, and no beamline can hold one there - so the +// two angles the beam-centre estimator cannot do without are measured from the spots instead. +// +// AZIMUTH: the spindle turned about the beam. Both mirror lines turn with it, and a mirror taken +// about the nominal one instead leaves a residual of twice the azimuth times the other detector +// coordinate - which does not bias the vote but spreads it until the wrong tooth wins. +// +// TIP: the spindle tipped out of the plane perpendicular to the beam, i.e. the part that makes it +// not perpendicular at all. The Friedel mate is then no longer in diffracting condition at exactly +// phi+180, and the compensating rotation leaves a purely radial residual along the mirrored +// coordinate. The crossing axis does not see it; only the Friedel one does. +// +// Both are angles of the spindle, not of the detector, and neither is ever written to the geometry. +struct SpindleEstimate { + float azimuth_rad = NAN; // fitted rotation of the spindle about the beam + float tip_rad = NAN; // fitted rotation of the spindle towards the beam + float vote_excess = NAN; // the Friedel vote's tallest tooth there + float vote_excess_nominal = NAN; // and at the nominal spindle, for comparison +}; + +// Beam centre from the symmetry of a rotation sweep's spot positions, before anything is indexed. +// Nothing here uses the background, a ring, a unit cell or an orientation matrix - only where the +// spots are, and two facts that hold exactly for any crystal on a single-axis goniometer. +// +// Friedel. Rotating the crystal 180 deg about the spindle m and taking -h negates the scattering +// vector's component ALONG m and leaves the rest untouched, so the spots on the frame at phi+180 are +// the exact mirror image, in that lab component, of the spots on the frame at phi. Mirroring one +// frame's spots onto its partner's therefore measures the beam coordinate along m. +// +// Second crossing. Along the sweep |q| and q.m are both constant and the Ewald condition fixes the +// component along the beam, so the remaining one can only change sign: every reflection is recorded +// twice, at positions mirrored in the coordinate PERPENDICULAR to m. That gives the other coordinate. +// +// The estimate is of the physical direct beam; it is returned in the PONI parametrisation the rest of +// the geometry uses, so the caller can put it straight into DiffractionExperiment. +// +// `sigma_pxl` is what the caller gates on, and it is the larger of two things: the scatter of what +// the individual frame pairs said, and how far the answer moves when the search is started +// elsewhere. The second is what matters. The vote sees only +-30 px about wherever it starts, so a +// header far enough out leaves the true peak outside the window and a false one wins it with a +// perfectly small scatter - a precision, not an accuracy. Only a start somewhere else finds that. +// +// The sweep has to reach 180 deg for the first and about a full turn for the second, which is why +// this cannot serve stills, screening wedges or short sweeps - the background estimator does, on +// three frames. `frame_angle_deg` is the rotation angle of each sampled frame, and the estimator +// pairs them itself; `spots` is the pooled spot list of those frames. +// +// std::nullopt when the sweep, or the diffraction on it, does not support the measurement. +// A non-null `spindle` is what asks for the spindle to be fitted and applied; with a null pointer +// this is exactly the estimator that was here before. +std::optional +FindBeamCenterFromSpotSymmetry(const DiffractionExperiment &experiment, + const std::vector &frame_angle_deg, + const std::vector &spots, + SpindleEstimate *spindle = nullptr); diff --git a/image_analysis/geom_refinement/CMakeLists.txt b/image_analysis/geom_refinement/CMakeLists.txt index c364dad5..08029708 100644 --- a/image_analysis/geom_refinement/CMakeLists.txt +++ b/image_analysis/geom_refinement/CMakeLists.txt @@ -1,5 +1,9 @@ ADD_LIBRARY(JFJochGeomRefinement STATIC + BeamCenterFromBackground.cpp + BeamCenterFromBackground.h + BeamCenterFromSpots.cpp + BeamCenterFromSpots.h RingOptimizer.cpp RingOptimizer.h AssignSpotsToRings.cpp diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 611626bb..c06c303d 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -26,8 +26,11 @@ #include "../image_analysis/MXAnalysisWithoutFPGA.h" #include "../image_analysis/beam_stop/ShadowFinder.h" #include "../image_analysis/IndexAndRefine.h" +#include "../image_analysis/geom_refinement/BeamCenterFromBackground.h" +#include "../image_analysis/geom_refinement/BeamCenterFromSpots.h" #include "../image_analysis/geom_refinement/GeometryRefiner.h" #include "../image_analysis/indexing/IndexerThreadPool.h" +#include "../image_analysis/spot_finding/ImageSpotFinderCPU.h" #include "../image_analysis/azint/AzIntEngineCPU.h" #include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h" #include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h" @@ -52,6 +55,46 @@ #include namespace { + // How precisely the beam centre has to be placed before the estimate is used at all. An estimate + // that fails it is treated as no estimate: the spot symmetry falls through to the background and + // the background to the header. A fit's own scatter will not serve as the gate - the ways either + // estimator goes badly wrong all leave a small one - so the spot symmetry reports instead how + // far its answer moves when the search is started elsewhere, which does separate them. On the 38 + // rotation regression crystals that stays below 0.35 px against errors of at most 0.74; over a + // ring of injected geometries it rises past 1.0 px on 96 % of the answers that are 2 px or more + // out, and on none of the 38 as they are. + constexpr float MAX_BEAM_CENTER_SIGMA_PXL = 1.0f; + + // Images the projection is built from when the beam centre is wanted but the beam-stop pre-pass + // is off; with it on, that pre-pass's image count is used and one projection serves both. + constexpr int BEAM_CENTER_PROJECTION_IMAGES = 60; + + // Spots per pre-scan frame kept for the beam centre, the strongest first. + constexpr size_t BEAM_CENTER_SPOTS_PER_IMAGE = 500; + + // Frames read for the spot symmetry alone when the first sample could not measure it. Both + // estimators are counting statistics over PAIRS of frames - the Friedel one over frames half a + // turn apart, the second crossing over any two frames a reflection is recorded on - so a sweep + // whose frames carry few spots can only be measured by reading more of them. On the two + // rotation regression crystals that decline, this is what turns the decline into an answer + // 0.04 px from the post-refined centre; the sample the beam stop is built from is left alone. + constexpr int BEAM_CENTER_SPARSE_IMAGES = 400; + + // Frames the spot symmetry reads, as a multiple of the beam stop's. Its two estimators are + // counts of matched PAIRS and the pairs grow faster than the frames do, so this is the one knob + // that buys independence from WHICH frames were drawn - and that dependence is the estimator's + // largest remaining failure. Asking the same code for one frame more or less re-draws the whole + // sample: over five draws of 38 sweeps, 60 frames leaves three answers more than 1 px out, none + // of which reports an unusual sigma, and 120 leaves none in 186. It is also what makes the end + // margin below free: at 60 frames excluding the ends moved three centres, at 120 it moves none. + constexpr int BEAM_CENTER_IMAGE_FACTOR = 2; + + // Images left out at each end of the sweep. Shutter synchronisation disturbs the first and last + // frame or two of a run, and a partly exposed frame is dark over the whole detector - which is + // exactly the kind of error a projection cannot average away and a background comparison reads + // as signal. The margin is capped at a tenth of the sweep so a short run still has a sample. + constexpr int PRESCAN_END_MARGIN_IMAGES = 5; + // Pick up to requested_images ordinals spread evenly across [0, images_to_process) for the // first pass of two-pass rotation indexing. std::vector select_equally_spaced_image_ordinals(int images_to_process, int requested_images) { @@ -75,6 +118,29 @@ namespace { return ret; } + // A sample arranged as pairs of images half a turn apart: half the ordinals spread over the + // part of the sweep that has a partner 180 deg later, and those partners. That is what the + // Friedel mirror needs. It is NOT the plain selection rearranged - even on a full turn the two + // sets are all but disjoint, sharing 2 of 60 ordinals on a 3600 image sweep - so the caller + // keeps them apart and gives each consumer the one it needs. Empty when the sweep is shorter + // than half a turn, which is the caller's signal that the spot symmetry cannot be measured here. + std::vector select_half_turn_paired_ordinals(int images_to_process, int requested_images, + float degrees_per_ordinal) { + if (!(degrees_per_ordinal > 0.0f)) + return {}; + const int half_turn = static_cast(std::lround(180.0 / degrees_per_ordinal)); + if (half_turn < 1 || images_to_process - half_turn < 1) + return {}; + + std::set ordinals; + for (const int ordinal : select_equally_spaced_image_ordinals(images_to_process - half_turn, + requested_images / 2)) { + ordinals.insert(ordinal); + ordinals.insert(ordinal + half_turn); + } + return {ordinals.begin(), ordinals.end()}; + } + // Apply a goniometer rotation SCALE: the stage turned k times the angle the file records. What is // wrong is the SWEEP, not where it began, so the per-frame increment and the per-frame oscillation // width both scale by k while the starting angle is left alone - the stage reached that position @@ -109,29 +175,111 @@ Rugnux::Rugnux(JFJochHDF5Reader &reader, DiffractionExperiment experiment, } } -void Rugnux::FindBeamStop(int start_image, int images_to_process, int frame_count) { +void Rugnux::PreScan(int start_image, int images_to_process, int frame_count) { Logger logger("Rugnux"); // Two-pass rotation runs this twice. The shadow does not move, and re-detecting would find - // nothing (its pixels are masked by now) and so clear the mask the first pass established. + // nothing (its pixels are masked by now) and so clear the mask the first pass established; + // the beam centre must not be re-read either, because by then the geometry is post-refined. const auto ¤t = pixel_mask_.GetMask(); - if (std::any_of(current.begin(), current.end(), - [](uint32_t v) { return (v & (1u << PixelMask::BeamStopPixelBit)) != 0; })) + const bool have_shadow = std::any_of(current.begin(), current.end(), + [](uint32_t v) { return (v & (1u << PixelMask::BeamStopPixelBit)) != 0; }); + const bool want_shadow = config_.detect_beam_stop.has_value() && !have_shadow; + const bool want_beam_center = config_.estimate_beam_center && !beam_center_placed_; + if (!want_shadow && !want_beam_center) return; - const auto sample = select_equally_spaced_image_ordinals(images_to_process, frame_count); + // The two consumers want different frames and each gets its own set. The shadow is built from + // the even spread over the sweep it has always been built from, so that asking for a beam + // centre cannot move the mask; the spot symmetry needs frames paired half a turn apart, where + // the sweep reaches that far. The two sets barely overlap, so the loop below walks their union + // and every frame is read once. + // + // Both samples are taken from the middle of the sweep, both ends left out: an even spread over + // [0, images_to_process) always picks the very first and the very last image, and those are the + // ones a shutter out of step spoils - the last carries under 90 % of the sample's median counts + // on 12 of 38 regression sweeps and is essentially unexposed on 5. + const auto goniometer = experiment_.GetGoniometer(); + const int margin = std::min(PRESCAN_END_MARGIN_IMAGES, images_to_process / 10); + auto shadow_sample = select_equally_spaced_image_ordinals(images_to_process - 2 * margin, + frame_count); + auto spot_sample = want_beam_center && goniometer + ? select_half_turn_paired_ordinals( + images_to_process - 2 * margin, + BEAM_CENTER_IMAGE_FACTOR * frame_count, + std::abs(goniometer->GetIncrement_deg()) * config_.stride) + : std::vector(); + for (int &ordinal : shadow_sample) + ordinal += margin; + for (int &ordinal : spot_sample) + ordinal += margin; + const bool want_spot_symmetry = !spot_sample.empty(); + const std::set shadow_set(shadow_sample.begin(), shadow_sample.end()); + const std::set spot_set(spot_sample.begin(), spot_sample.end()); + std::set sample = shadow_set; + sample.insert(spot_set.begin(), spot_set.end()); if (sample.empty()) return; ShadowFinder finder(experiment_, pixel_mask_); std::vector buffer; + + // The spot symmetry is read in the same pass as the projection rather than in a pre-pass of its + // own. The image is preprocessed a second time for it: the shadow is accumulated from raw counts + // and the spot finder works on the converted image. + std::unique_ptr prescan_mapping; + std::unique_ptr preprocessor; + std::unique_ptr preprocessed; + std::unique_ptr spot_finder; + std::vector decompression_buffer; + std::vector frame_angle_deg; + std::vector beam_center_spots; + if (want_spot_symmetry) { + prescan_mapping = std::make_unique(experiment_, pixel_mask_); + preprocessor = std::make_unique(experiment_, pixel_mask_); + preprocessed = std::make_unique(experiment_.GetPixelsNum()); + if (config_.spot_finding.adaptive_threshold) + spot_finder = std::make_unique(*prescan_mapping); + else + spot_finder = std::make_unique(experiment_.GetXPixelsNumConv(), + experiment_.GetYPixelsNumConv()); + } + + // The spots of one image, added to the pool as a frame of its own. Taken as a step of its own + // because the sparse-sweep pass below reads images the shadow does not. + const auto add_spots = [&](CompressedImage &image, int image_idx) { + const int frame = static_cast(frame_angle_deg.size()); + frame_angle_deg.push_back(goniometer->GetAngle_deg(static_cast(image_idx))); + try { + preprocessor->Analyze(*preprocessed, image.GetUncompressedPtr(decompression_buffer), + image.GetMode()); + } catch (const std::exception &e) { + logger.Warning("Pre-scan: failed to preprocess image {}: {}", image_idx, e.what()); + return; + } + auto spots = spot_finder->Run(*preprocessed, config_.spot_finding); + // The strongest of a crowded frame: the symmetry is over-determined either way, and the + // matching is quadratic in the spots of one frame. + if (spots.size() > BEAM_CENTER_SPOTS_PER_IMAGE) { + std::partial_sort(spots.begin(), spots.begin() + BEAM_CENTER_SPOTS_PER_IMAGE, spots.end(), + [](const DiffractionSpot &a, const DiffractionSpot &b) { + return a.Count() > b.Count(); }); + spots.resize(BEAM_CENTER_SPOTS_PER_IMAGE); + } + for (const auto &spot : spots) { + const Coord centroid = spot.RawCoord(); + beam_center_spots.push_back({centroid.x, centroid.y, + static_cast(spot.Count()), frame}); + } + }; + for (const int ordinal : sample) { const int image_idx = start_image + ordinal * config_.stride; std::shared_ptr img; try { img = reader_.GetRawImage(image_idx); } catch (const std::exception &e) { - logger.Warning("Beam stop detection: failed to load image {}: {}", image_idx, e.what()); + logger.Warning("Pre-scan: failed to load image {}: {}", image_idx, e.what()); continue; } if (!img) continue; @@ -140,20 +288,111 @@ void Rugnux::FindBeamStop(int start_image, int images_to_process, int frame_coun msg.image = img->image; msg.number = ordinal; msg.original_number = image_idx; - finder.AddImage(msg, buffer); + if (shadow_set.contains(ordinal)) + finder.AddImage(msg, buffer); + + if (spot_set.contains(ordinal)) + add_spots(msg.image, image_idx); } if (finder.GetFrameCount() == 0) { - logger.Warning("Beam stop detection: no image could be read. Skipping."); + logger.Warning("Pre-scan: no image could be read. Skipping."); return; } - const auto shadow = finder.GetMask(); - const auto shadowed = std::count(shadow.begin(), shadow.end(), 1u); - pixel_mask_.LoadBeamStopMask(experiment_, shadow); - logger.Info("Beam stop shadow: {} pixels ({:.2f}% of the detector) found in {} images", - shadowed, 100.0 * static_cast(shadowed) / static_cast(shadow.size()), - finder.GetFrameCount()); + if (want_shadow) { + const auto shadow = finder.GetMask(); + const auto shadowed = std::count(shadow.begin(), shadow.end(), 1u); + pixel_mask_.LoadBeamStopMask(experiment_, shadow); + logger.Info("Beam stop shadow: {} pixels ({:.2f}% of the detector) found in {} images", + shadowed, 100.0 * static_cast(shadowed) / static_cast(shadow.size()), + finder.GetFrameCount()); + } + + // The beam centre comes after, so that the shadow is in the mask by the time it is placed and + // the holder arm - the largest azimuthal asymmetry on the detector - is out of the way of the + // background the centre is read from. + if (!want_beam_center) + return; + beam_center_placed_ = true; + + // The spot symmetry first where the sweep supports it: it is exact geometry rather than a + // property of the scattering, and it measures the sets the background cannot - a flat radial + // profile at long wavelength leaves that fit nothing to bite on. The background answers + // wherever the sweep does not reach half a turn, which is every still and every screening + // wedge, so the two cover between them what neither covers alone. + std::string source = "spot symmetry"; + std::optional estimate; + SpindleEstimate spindle; + if (want_spot_symmetry) + estimate = FindBeamCenterFromSpotSymmetry(experiment_, frame_angle_deg, beam_center_spots, + config_.fit_spindle ? &spindle : nullptr); + + // Where it does not come out - no answer at all, or one that moves when the search is started + // elsewhere - read more of the sweep and ask again. Both estimators are counts of matched spot + // PAIRS, the Friedel one over frames half a turn apart and the second crossing over any two + // frames one reflection is recorded on, so a sweep whose frames carry few spots, or whose pairs + // match rarely, is measurable only from more of them. More data first and a weaker estimator + // second: this runs before the fall-through, and only what survives it falls through. The extra + // frames go to the spot finder alone, so the beam-stop projection keeps the images it was + // validated on. + if (want_spot_symmetry && (!estimate || estimate->sigma_pxl > MAX_BEAM_CENTER_SIGMA_PXL)) { + // What has already been read FOR THE SPOTS, which is not the whole sample: the frames the + // shadow was built from carry none of them, and testing against the union would silently + // skip the ones this pass exists to read. + std::vector extra; + for (int ordinal : select_half_turn_paired_ordinals( + images_to_process - 2 * margin, BEAM_CENTER_SPARSE_IMAGES, + std::abs(goniometer->GetIncrement_deg()) * config_.stride)) { + ordinal += margin; + if (!spot_set.contains(ordinal)) + extra.push_back(ordinal); + } + if (!extra.empty()) { + logger.Info("Beam centre: the spot symmetry does not come out on {} frames " + "({} spots); reading {} more frames", + frame_angle_deg.size(), beam_center_spots.size(), extra.size()); + for (const int ordinal : extra) { + const int image_idx = start_image + ordinal * config_.stride; + std::shared_ptr img; + try { + img = reader_.GetRawImage(image_idx); + } catch (const std::exception &e) { + logger.Warning("Pre-scan: failed to load image {}: {}", image_idx, e.what()); + continue; + } + if (img) + add_spots(img->image, image_idx); + } + estimate = FindBeamCenterFromSpotSymmetry(experiment_, frame_angle_deg, + beam_center_spots, + config_.fit_spindle ? &spindle : nullptr); + } + } + if (config_.fit_spindle && std::isfinite(spindle.azimuth_rad)) + logger.Info("Spindle: azimuth about the beam {:.3f} mrad, tip towards it {:.3f} mrad " + "(Friedel vote {:.0f} -> {:.0f})", + 1e3 * spindle.azimuth_rad, 1e3 * spindle.tip_rad, + spindle.vote_excess_nominal, spindle.vote_excess); + if (!estimate || estimate->sigma_pxl > MAX_BEAM_CENTER_SIGMA_PXL) { + source = "background"; + estimate = FindBeamCenterFromBackground(experiment_, pixel_mask_, finder.GetMeanProjection()); + } + if (!estimate) { + logger.Info("Beam centre: not measurable, keeping ({:.2f},{:.2f})", + experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl()); + return; + } + const float moved = std::hypot(estimate->beam_x_pxl - experiment_.GetBeamX_pxl(), + estimate->beam_y_pxl - experiment_.GetBeamY_pxl()); + const bool commit = estimate->sigma_pxl <= MAX_BEAM_CENTER_SIGMA_PXL; + logger.Info("Beam centre from {}: ({:.2f},{:.2f}) -> ({:.2f},{:.2f}), moved {:.2f} px, " + "sigma {:.2f} px => {}", + source, experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(), + estimate->beam_x_pxl, estimate->beam_y_pxl, moved, estimate->sigma_pxl, + commit ? "COMMIT" : "reject (kept header)"); + if (commit) + experiment_.BeamX_pxl(estimate->beam_x_pxl).BeamY_pxl(estimate->beam_y_pxl); } void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_process, @@ -328,6 +567,8 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ experiment_.BeamX_pxl(r.beam_x_px).BeamY_pxl(r.beam_y_px).DetectorDistance_mm(r.distance_mm); experiment_.SetUnitCell(r.cell); + // The beam centre now comes from indexed spots, which the background cannot better. + beam_center_placed_ = true; } ProcessResult Rugnux::Run(RugnuxObserver *observer) { @@ -639,8 +880,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // After any geometry refinement, so the shadow is found about the beam centre actually used, // and before the azimuthal mapping and the output mask below, which both read pixel_mask_. - if (config_.detect_beam_stop.has_value()) - FindBeamStop(start_image, images_to_process, config_.detect_beam_stop.value()); + if (config_.detect_beam_stop.has_value() || config_.estimate_beam_center) + PreScan(start_image, images_to_process, + config_.detect_beam_stop.value_or(BEAM_CENTER_PROJECTION_IMAGES)); AzimuthalIntegrationMapping mapping(experiment_, pixel_mask_); diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index a2d2e543..4c76ed35 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -81,6 +81,14 @@ struct ProcessConfig { // pixel mask, so nothing behind them is integrated. The value is the number of frames projected. std::optional detect_beam_stop; + // Beam centre before anything is indexed (--estimate-beam-center). Runs in the pre-scan, on the + // frames the beam-stop pre-pass reads when that is on and on its own sample when it is off, and + // replaces the header beam centre when it is measured precisely enough. From the symmetry of the + // spot positions where the sweep reaches half a turn, and from the isotropy of the background + // where it does not. + bool estimate_beam_center = false; + bool fit_spindle = false; + // Rotation two-pass geometry post-refinement (FullAnalysis, rotation only; on by default in the rugnux // CLI, --rotation-no-postrefine disables it). When set, a first pass integrates and post-refines the // detector distance + beam (from the observed spot positions) and the cell scale + rotation axis (from the @@ -249,14 +257,22 @@ class Rugnux { bool prepass_promoted_point_group_ = false; + // Whether the beam centre has already been placed by a measurement the pre-scan cannot better: + // a stills bundle adjustment, or an earlier pre-scan of this same run. The pre-scan estimate is + // a starting point for indexing, so it must not overwrite either of those. + bool beam_center_placed_ = false; + + // Pre-scan: read a spread sample of frames and take two things off them - the shadow of the beam + // stop and its holder, added to the pixel mask (config_.detect_beam_stop), and the beam centre + // (config_.estimate_beam_center), from the symmetry of the spots where the sweep reaches half a + // turn and from the isotropy of the scattered background where it does not. Either may be asked + // for without the other. + void PreScan(int start_image, int images_to_process, int frame_count); + // Stills global geometry-refinement first pass (config_.refine_geometry): index a spread sample of // frames, bundle-adjust the shared beam/distance/cell from the strongest ones, and apply the result // to experiment_ so the main pass re-indexes with it. No-op (leaves experiment_ unchanged) if it // cannot run. Stills only; rotation has its own two-pass. - // Beam-stop shadow pre-pass (config_.detect_beam_stop): project a spread sample of frames and - // add the shadow of the stop and its holder to the pixel mask. - void FindBeamStop(int start_image, int images_to_process, int frame_count); - void RefineStillsGeometry(int start_image, int end_image, int images_to_process, RugnuxObserver *observer); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 21b26669..7513d25f 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -97,6 +97,8 @@ void print_usage() { std::cout << " Detector mask" << std::endl; std::cout << " --detect-beam-stop[=N|off] Find the beam stop and its holder in a projection of N images and add them to the pixel mask (bit 9), so nothing shadowed by them is integrated. ON by default (60 images); =off disables. Reflections behind the stop are attenuated but not flagged, so they are integrated low with a plausible sigma and no existing rejection catches them" << std::endl; + std::cout << " --estimate-beam-center Place the beam centre before anything is indexed, and use it in place of the header value when it is measured precisely enough. On a sweep that reaches half a turn it comes from the symmetry of the spot positions - the frames 180 deg apart are each other's mirror image, and every reflection is recorded twice - and where the sweep does not reach that far, from the isotropy of the scattered background, which needs only a few frames. It reads frames of its own, chosen as pairs half a turn apart, so it does not change the mask --detect-beam-stop finds. Ignored when --beam-x/--beam-y are given, and when a stills geometry refinement has already placed the centre from indexed spots" << std::endl; + std::cout << " --no-fit-spindle Take the goniometer axis from the file rather than measuring the spindle's rotation about the beam from the spots. Measuring it is the DEFAULT and the file is never right: every master writes an exact lab axis and no goniometer is one. Both mirror lines of the beam-centre estimator turn with the spindle, so an axis a few tenths of a milliradian out smears the vote until a neighbouring tooth wins. Only has an effect with --estimate-beam-center" << std::endl; std::cout << std::endl; std::cout << " Spot finding" << std::endl; @@ -217,6 +219,9 @@ enum { OPT_BACKGROUND_RADIAL, OPT_REFINE_GEOMETRY, OPT_DETECT_BEAM_STOP, + OPT_ESTIMATE_BEAM_CENTER, + OPT_FIT_SPINDLE, + OPT_NO_FIT_SPINDLE, OPT_BANDWIDTH, OPT_INTEGRATION_RADIUS, OPT_INTEGRATION_STENCIL, @@ -314,6 +319,9 @@ static option long_options[] = { {"rotation-scale", required_argument, nullptr, OPT_ROTATION_SCALE}, {"refine-geometry", optional_argument, nullptr, OPT_REFINE_GEOMETRY}, {"detect-beam-stop", optional_argument, nullptr, OPT_DETECT_BEAM_STOP}, + {"estimate-beam-center", no_argument, nullptr, OPT_ESTIMATE_BEAM_CENTER}, + {"fit-spindle", no_argument, nullptr, OPT_FIT_SPINDLE}, + {"no-fit-spindle", no_argument, nullptr, OPT_NO_FIT_SPINDLE}, {"spot-sigma", required_argument, nullptr, OPT_SPOT_SIGMA}, @@ -626,6 +634,8 @@ static int RunRugnux(int argc, char **argv) { bool background_radial_given = false; // --background-radial seen at all (unset => auto) std::optional background_radial_arg; // when given: set = force on/off, unset = auto std::optional detect_beam_stop = 60; // --detect-beam-stop[=N|off]; on by default + bool estimate_beam_center = false; // --estimate-beam-center + bool fit_spindle = true; // --fit-spindle / --no-fit-spindle std::optional refine_geometry; // --refine-geometry[=N]: stills global geometry-refinement pass bool refine_geometry_disabled = false; // --refine-geometry=off: opt out of the stills default-on @@ -722,6 +732,16 @@ static int RunRugnux(int argc, char **argv) { ? parse_number_arg(optarg, "--detect-beam-stop", logger, 1, 1000000) : 60; break; + case OPT_ESTIMATE_BEAM_CENTER: + estimate_beam_center = true; + break; + case OPT_FIT_SPINDLE: + fit_spindle = true; + estimate_beam_center = true; + break; + case OPT_NO_FIT_SPINDLE: + fit_spindle = false; + break; case OPT_REFINE_GEOMETRY: { if (optarg && std::string(optarg) == "off") { refine_geometry = std::nullopt; @@ -1570,6 +1590,8 @@ static int RunRugnux(int argc, char **argv) { // azimuthal-integration settings are derived, which depend on the geometry. if (beam_x) experiment.BeamX_pxl(beam_x.value()); if (beam_y) experiment.BeamY_pxl(beam_y.value()); + // A centre the user typed in is the answer, not a starting point. + if (beam_x || beam_y) estimate_beam_center = false; if (detector_distance_mm) experiment.DetectorDistance_mm(detector_distance_mm.value()); if (wavelength_A) experiment.IncidentEnergy_keV(WVL_1A_IN_KEV / wavelength_A.value()); if (rot1_rad) experiment.PoniRot1_rad(rot1_rad.value()); @@ -1622,6 +1644,8 @@ static int RunRugnux(int argc, char **argv) { config.nthreads = nthreads; config.output_prefix = output_prefix; config.detect_beam_stop = detect_beam_stop; + config.estimate_beam_center = estimate_beam_center; + config.fit_spindle = fit_spindle; Rugnux process(reader, experiment, *dataset->pixel_mask, config); g_active_process = &process; @@ -1664,6 +1688,8 @@ static int RunRugnux(int argc, char **argv) { config.nthreads = nthreads; config.output_prefix = output_prefix; config.detect_beam_stop = detect_beam_stop; + config.estimate_beam_center = estimate_beam_center; + config.fit_spindle = fit_spindle; config.write_process_h5 = false; // the .poni below is the output of this mode // Spot finding for --calibration spots. Indexing is off: a calibration wants the spot positions @@ -2051,6 +2077,8 @@ static int RunRugnux(int argc, char **argv) { config.rotation_indexing = rotation_indexing; config.two_pass_rotation = two_pass_rotation; config.detect_beam_stop = detect_beam_stop; + config.estimate_beam_center = estimate_beam_center; + config.fit_spindle = fit_spindle; config.rotation_postrefine_geometry = rotation_postrefine_geometry; config.rotation_scale = rotation_scale; config.rotation_indexing_image_count = rotation_indexing_image_count; diff --git a/tests/BeamCenterFromBackgroundTest.cpp b/tests/BeamCenterFromBackgroundTest.cpp new file mode 100644 index 00000000..dc0007b9 --- /dev/null +++ b/tests/BeamCenterFromBackgroundTest.cpp @@ -0,0 +1,126 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include +#include + +#include "../image_analysis/geom_refinement/BeamCenterFromBackground.h" +#include "../common/DetectorSetup.h" +#include "../common/JFJochMath.h" + +namespace { + +// Solvent and air scatter: a decaying continuum with the water ring on it. The ring is where the +// leverage comes from - the continuum here is a pure exponential, on which g' is proportional to g +// and a shift and an amplitude are the same thing - so `ring` is how much there is to fit. +float background(float two_theta_rad, float ring) { + const float ring_two_theta = 0.3239f; // ~3.1 A at 1 A + const float t = (two_theta_rad - ring_two_theta) / 0.035f; + return 140.0f * std::exp(-two_theta_rad / 0.25f) + ring * std::exp(-0.5f * t * t); +} + +// The projection the pre-scan hands over: the mean of a few tens of frames, laid out about +// geom_true, NAN where the detector has nothing. `shadow_sector` multiplies one sextant, the way a +// holder arm or a cryostream does. +std::vector SynthesiseProjection(const DiffractionExperiment &experiment, + const PixelMask &mask, + const DiffractionGeometry &geom_true, + float ring, float shadow_sector) { + const auto W = static_cast(experiment.GetXPixelsNumConv()); + const auto H = static_cast(experiment.GetYPixelsNumConv()); + const auto &pixel_mask = mask.GetMask(experiment); + + std::vector mean(static_cast(W) * H, NAN); + std::mt19937 rng(20260812); + std::normal_distribution gauss(0.0f, 1.0f); + constexpr float FRAMES = 60.0f; // the mean of this many frames, so the noise is that far down + + for (int y = 0; y < H; y++) { + for (int x = 0; x < W; x++) { + const size_t i = static_cast(y) * W + x; + if (pixel_mask[i] != 0) + continue; + float value = background(geom_true.TwoTheta_rad(static_cast(x), static_cast(y)), ring); + const float phi = geom_true.Phi_rad(static_cast(x), static_cast(y)); + if (phi > 0.0f && phi < static_cast(PI) / 3.0f) + value *= shadow_sector; + mean[i] = value + gauss(rng) * std::sqrt(value / FRAMES); + } + } + return mean; +} + +DiffractionExperiment TestExperiment() { + DiffractionExperiment x(DetJF4M()); + x.IncidentEnergy_keV(WVL_1A_IN_KEV).DetectorDistance_mm(100.0f); + // The band the estimator fits, 12-2.2 A, has to be on the detector, so start from its centre. + x.BeamX_pxl(static_cast(x.GetXPixelsNumConv()) / 2.0f) + .BeamY_pxl(static_cast(x.GetYPixelsNumConv()) / 2.0f); + return x; +} + +DiffractionGeometry OffsetBy(const DiffractionGeometry &geom, float dx, float dy) { + DiffractionGeometry out = geom; + out.BeamX_pxl(geom.GetBeamX_pxl() + dx).BeamY_pxl(geom.GetBeamY_pxl() + dy); + return out; +} + +} // namespace + +// The measurement: the background is isotropic in 2-theta about the beam, so a centre that is off +// shifts each azimuthal sector's radial profile by a different amount, and the shifts give the +// centre back. Nothing here is indexed, so this is what a de-novo run has to start from - and the +// estimate has to arrive, because a routine that quietly returns "not measurable" is +// indistinguishable from a careful refusal in every log line and every merging statistic. +TEST_CASE("BeamCenterFromBackground_RecoversAnInjectedOffset", "[BeamCenter]") { + DiffractionExperiment x = TestExperiment(); + PixelMask pixel_mask(x); + + const DiffractionGeometry geom_true = OffsetBy(x.GetDiffractionGeometry(), 3.0f, -2.5f); + const auto projection = SynthesiseProjection(x, pixel_mask, geom_true, 60.0f, 1.0f); + const auto estimate = FindBeamCenterFromBackground(x, pixel_mask, projection); + + REQUIRE(estimate.has_value()); + CHECK(estimate->beam_x_pxl == Catch::Approx(geom_true.GetBeamX_pxl()).margin(0.5)); + CHECK(estimate->beam_y_pxl == Catch::Approx(geom_true.GetBeamY_pxl()).margin(0.5)); + // And it has to say so precisely enough to be used: the caller commits at 1 px. + CHECK(estimate->sigma_pxl < 1.0f); +} + +// The sigma is the only thing standing between a bad background and a wrong geometry, so it has to +// grow when the ring it is fitting does not. With the ring at 1.4% of the continuum the centre is +// still found, but the fit says it is an order of magnitude less sure of it. +TEST_CASE("BeamCenterFromBackground_SigmaTracksTheLeverage", "[BeamCenter]") { + DiffractionExperiment x = TestExperiment(); + PixelMask pixel_mask(x); + + const DiffractionGeometry geom_true = OffsetBy(x.GetDiffractionGeometry(), 3.0f, -2.5f); + const auto strong = FindBeamCenterFromBackground( + x, pixel_mask, SynthesiseProjection(x, pixel_mask, geom_true, 60.0f, 1.0f)); + const auto weak = FindBeamCenterFromBackground( + x, pixel_mask, SynthesiseProjection(x, pixel_mask, geom_true, 2.0f, 1.0f)); + + REQUIRE(strong.has_value()); + REQUIRE(weak.has_value()); + CHECK(weak->beam_x_pxl == Catch::Approx(geom_true.GetBeamX_pxl()).margin(1.0)); + CHECK(weak->beam_y_pxl == Catch::Approx(geom_true.GetBeamY_pxl()).margin(1.0)); + CHECK(weak->sigma_pxl > 5.0f * strong->sigma_pxl); +} + +// A holder arm or a cryostream is multiplicative and azimuthal, and a sector that is simply darker +// looks exactly like a sector whose profile has moved. The per-sector amplitude is what tells them +// apart: without it half a sextant of shadow reads as tens of pixels of centre error. +TEST_CASE("BeamCenterFromBackground_AnAzimuthalShadowIsNotACentreError", "[BeamCenter]") { + DiffractionExperiment x = TestExperiment(); + PixelMask pixel_mask(x); + + const DiffractionGeometry geom_true = x.GetDiffractionGeometry(); // the centre is already right + const auto projection = SynthesiseProjection(x, pixel_mask, geom_true, 60.0f, 0.5f); + const auto estimate = FindBeamCenterFromBackground(x, pixel_mask, projection); + + REQUIRE(estimate.has_value()); + CHECK(estimate->beam_x_pxl == Catch::Approx(geom_true.GetBeamX_pxl()).margin(1.0)); + CHECK(estimate->beam_y_pxl == Catch::Approx(geom_true.GetBeamY_pxl()).margin(1.0)); +} diff --git a/tests/BeamCenterFromSpotsTest.cpp b/tests/BeamCenterFromSpotsTest.cpp new file mode 100644 index 00000000..e5ee73d9 --- /dev/null +++ b/tests/BeamCenterFromSpotsTest.cpp @@ -0,0 +1,284 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include +#include + +#include "../image_analysis/geom_refinement/BeamCenterFromSpots.h" +#include "../common/DetectorSetup.h" +#include "../common/JFJochMath.h" + +namespace { + +constexpr float WEDGE_DEG = 1.0f; +constexpr int PAIRS = 30; +constexpr float PAIR_SPACING_DEG = 180.0f / PAIRS; +constexpr float POSITION_NOISE_PXL = 0.15f; + +struct Sweep { + std::vector angle_deg; + std::vector spots; +}; + +DiffractionExperiment TestExperiment(const Coord &spindle) { + DiffractionExperiment x(DetDECTRIS(1500, 1600, "Test detector", "")); + x.IncidentEnergy_keV(WVL_1A_IN_KEV).DetectorDistance_mm(180.0f); + x.BeamX_pxl(static_cast(x.GetXPixelsNumConv()) / 2.0f) + .BeamY_pxl(static_cast(x.GetYPixelsNumConv()) / 2.0f); + x.Goniometer(GoniometerAxis("omega", 0.0f, WEDGE_DEG, spindle, {})); + return x; +} + +DiffractionGeometry OffsetBy(const DiffractionGeometry &geom, float dx, float dy) { + DiffractionGeometry out = geom; + out.BeamX_pxl(geom.GetBeamX_pxl() + dx).BeamY_pxl(geom.GetBeamY_pxl() + dy); + return out; +} + +// |F(h)| = |F(-h)|, which is what tells a true Friedel match from a lattice-shifted one, so the +// intensity has to depend on the reflection through |h| |k| |l| alone. +float StructureFactor(int h, int k, int l) { + const uint32_t key = 2654435761u * (std::abs(h) + 41u * std::abs(k) + 1723u * std::abs(l) + 1u); + return 50.0f + static_cast(key % 4000u); +} + +// A synthetic rotation sweep. An orthorhombic lattice in a fixed orientation is turned about the +// spindle; every reflection meets the Ewald sphere exactly twice, and each meeting is put on the +// detector at the geometry the estimator is asked to find. Only the frames of the pre-scan sample +// are kept, so what comes out is what the pre-pass sees: pairs half a turn apart, and about a +// thirtieth of the sweep. +Sweep MakeSweep(const DiffractionExperiment &experiment, const DiffractionGeometry &truth, int pairs) { + const Coord spindle = experiment.GetGoniometer()->GetAxis().Normalize(); + const Coord S0 = truth.GetScatteringVector(); + + Sweep sweep; + for (int i = 0; i < pairs; i++) { + sweep.angle_deg.push_back(i * PAIR_SPACING_DEG); + sweep.angle_deg.push_back(i * PAIR_SPACING_DEG + 180.0f); + } + + const RotMatrix orientation(0.7f, Coord(0.3f, 0.5f, 0.8f)); + const Coord astar = orientation * Coord(1.0f / 55.0f, 0, 0); + const Coord bstar = orientation * Coord(0, 1.0f / 65.0f, 0); + const Coord cstar = orientation * Coord(0, 0, 1.0f / 75.0f); + const float q_max = 1.0f / 2.4f; + + std::mt19937 rng(20260812); + std::normal_distribution noise(0.0f, POSITION_NOISE_PXL); + + for (int h = -30; h <= 30; h++) + for (int k = -35; k <= 35; k++) + for (int l = -40; l <= 40; l++) { + if (h == 0 && k == 0 && l == 0) + continue; + const Coord q0 = astar * static_cast(h) + bstar * static_cast(k) + + cstar * static_cast(l); + const float q_sq = q0 * q0; + if (q_sq > q_max * q_max) + continue; + + // Turning about the spindle leaves q.m alone and rotates the rest, so the Ewald + // condition q.S0 = -|q|^2/2 is one sinusoid in the rotation angle: two solutions, + // which are the reflection's two crossings. + // + // Which WAY the angle turns the crystal is rugnux's convention, not a free choice: + // BraggPredictionRot takes the offset to the diffracting condition as + // -atan2(sin, cos) in this same frame, i.e. dq/dphi = -m x q. The sweep this + // generator makes is otherwise time-reversed with respect to every real dataset, + // which no estimator that reads only positions can notice. + const Coord along = spindle * (q0 * spindle); + const Coord across = q0 - along; + const Coord turned = across % spindle; + const float a = across * S0, b = turned * S0; + const float radius = std::hypot(a, b); + const float target = -0.5f * q_sq - along * S0; + if (radius <= 0.0f || std::abs(target) > radius) + continue; + const float phase = std::atan2(b, a); + const float offset = std::acos(target / radius); + + for (const float sign: {-1.0f, 1.0f}) { + float phi_deg = (phase + sign * offset) * 180.0f / static_cast(PI); + while (phi_deg < 0.0f) phi_deg += 360.0f; + while (phi_deg >= 360.0f) phi_deg -= 360.0f; + for (size_t frame = 0; frame < sweep.angle_deg.size(); frame++) { + if (phi_deg < sweep.angle_deg[frame] || phi_deg >= sweep.angle_deg[frame] + WEDGE_DEG) + continue; + const float phi = phi_deg * static_cast(PI) / 180.0f; + const Coord q = along + across * std::cos(phi) + turned * std::sin(phi); + const auto [x, y] = truth.RecipToDetector(q); + if (!std::isfinite(x) || x < 0 || y < 0 + || x >= experiment.GetXPixelsNumConv() || y >= experiment.GetYPixelsNumConv()) + continue; + sweep.spots.push_back({x + noise(rng), y + noise(rng), StructureFactor(h, k, l), + static_cast(frame)}); + } + } + } + return sweep; +} + +} // namespace + +// The measurement: the sweep at phi+180 is the mirror image of the sweep at phi in the coordinate +// along the spindle, and every reflection crosses the Ewald sphere twice, mirrored in the other +// one. Between them the two put the beam centre where nothing has been indexed yet. Both a +// horizontal spindle and a vertical one, because which estimator supplies which coordinate is +// decided by the goniometer axis and one of the two is a real detector's arrangement. +TEST_CASE("BeamCenterFromSpots_RecoversAnInjectedOffset", "[BeamCenter]") { + const Coord spindle = GENERATE(Coord(1, 0, 0), Coord(0, 1, 0)); + const DiffractionExperiment x = TestExperiment(spindle); + const DiffractionGeometry truth = OffsetBy(x.GetDiffractionGeometry(), 3.0f, -2.5f); + const Sweep sweep = MakeSweep(x, truth, PAIRS); + + const auto estimate = FindBeamCenterFromSpotSymmetry(x, sweep.angle_deg, sweep.spots); + + REQUIRE(estimate.has_value()); + CHECK(estimate->beam_x_pxl == Catch::Approx(truth.GetBeamX_pxl()).margin(0.5)); + CHECK(estimate->beam_y_pxl == Catch::Approx(truth.GetBeamY_pxl()).margin(0.5)); + // And it has to say so precisely enough to be used: the caller commits at 1 px. + CHECK(estimate->sigma_pxl < 1.0f); + + // The header is where the search starts and nothing more. A file whose centre is 15 px out has + // to give the same answer as one whose centre is right, or the estimate is partly the header's. + DiffractionExperiment moved = x; + moved.BeamX_pxl(x.GetBeamX_pxl() + 15.0f).BeamY_pxl(x.GetBeamY_pxl() - 15.0f); + const auto from_elsewhere = FindBeamCenterFromSpotSymmetry(moved, sweep.angle_deg, sweep.spots); + REQUIRE(from_elsewhere.has_value()); + CHECK(from_elsewhere->beam_x_pxl == Catch::Approx(estimate->beam_x_pxl).margin(0.1)); + CHECK(from_elsewhere->beam_y_pxl == Catch::Approx(estimate->beam_y_pxl).margin(0.1)); +} + +// A hot pixel, the beam-stop halo or the edge of a mask sits at the SAME place on every frame, so +// it pairs with itself and votes at twice its own position - one vote per frame pair, which is +// enough to out-vote the real peak. On real data that was two crystals wrong by 12 and 1.9 px. +TEST_CASE("BeamCenterFromSpots_APersistentArtefactIsNotABeamCentre", "[BeamCenter]") { + const DiffractionExperiment x = TestExperiment(Coord(1, 0, 0)); + const DiffractionGeometry truth = OffsetBy(x.GetDiffractionGeometry(), 3.0f, -2.5f); + Sweep sweep = MakeSweep(x, truth, PAIRS); + + for (size_t frame = 0; frame < sweep.angle_deg.size(); frame++) + for (int i = 0; i < 12; i++) + sweep.spots.push_back({truth.GetBeamX_pxl() + 40.0f + 3.0f * i, + truth.GetBeamY_pxl() - 25.0f - 2.0f * i, + 9000.0f, static_cast(frame)}); + + const auto estimate = FindBeamCenterFromSpotSymmetry(x, sweep.angle_deg, sweep.spots); + + REQUIRE(estimate.has_value()); + CHECK(estimate->beam_x_pxl == Catch::Approx(truth.GetBeamX_pxl()).margin(0.5)); + CHECK(estimate->beam_y_pxl == Catch::Approx(truth.GetBeamY_pxl()).margin(0.5)); +} + +// The vote reaches 30 px about wherever the search starts, so a header further out than that leaves +// the true peak outside the window and a lattice-shifted one wins it - and wins it cleanly, with a +// frame-pair scatter as small as a correct answer's. What separates them is that the wrong one +// depends on where the search began: started elsewhere, the estimator finds something else. The +// reported sigma has to carry that, or the caller commits a centre tens of pixels out. +TEST_CASE("BeamCenterFromSpots_AnAnswerThatDependsOnTheHeaderIsNotACentre", "[BeamCenter]") { + const DiffractionExperiment x = TestExperiment(Coord(1, 0, 0)); + const DiffractionGeometry truth = OffsetBy(x.GetDiffractionGeometry(), 35.0f, 35.0f); + const Sweep sweep = MakeSweep(x, truth, PAIRS); + + const auto estimate = FindBeamCenterFromSpotSymmetry(x, sweep.angle_deg, sweep.spots); + + REQUIRE(estimate.has_value()); + // It answers, and the answer is nowhere near the truth - that is the failure mode, not a bug. + CHECK(std::hypot(estimate->beam_x_pxl - truth.GetBeamX_pxl(), + estimate->beam_y_pxl - truth.GetBeamY_pxl()) > 10.0f); + // What must not happen is that it says so quietly: the caller commits at 1 px. + CHECK(estimate->sigma_pxl > 1.0f); +} + +// Not every spot is one of a reflection's two crossings. A second read-out, a satellite, a family +// of detector artefacts - anything that shadows the real spots at a fixed displacement pairs with +// them at that displacement and votes as a tooth of its own, and here that tooth is TWICE the true +// one's height. What the decoys cannot copy is when their partner appears: the sweep angle between +// two crossings is fixed by where the first one is, and a decoy's is not. +TEST_CASE("BeamCenterFromSpots_ATallerDecoyToothDoesNotWin", "[BeamCenter]") { + const DiffractionExperiment x = TestExperiment(Coord(1, 0, 0)); + const DiffractionGeometry truth = OffsetBy(x.GetDiffractionGeometry(), 3.0f, -2.5f); + Sweep sweep = MakeSweep(x, truth, PAIRS); + + std::mt19937 rng(20260813); + const size_t recorded = sweep.spots.size(); + for (size_t i = 0; i < recorded; i++) { + const BeamCenterSpot &spot = sweep.spots[i]; + sweep.spots.push_back({spot.x, spot.y + 9.0f, spot.intensity, + static_cast(rng() % sweep.angle_deg.size())}); + } + + const auto estimate = FindBeamCenterFromSpotSymmetry(x, sweep.angle_deg, sweep.spots); + + REQUIRE(estimate.has_value()); + // The decoy tooth is 4.5 px away in the coordinate the second crossing supplies. + CHECK(estimate->beam_y_pxl == Catch::Approx(truth.GetBeamY_pxl()).margin(0.5)); +} + +// The Friedel mirror needs half a turn. On a screening wedge, on stills, on anything short there is +// no partner frame to mirror onto, and the estimator has to say so rather than answer from the +// pairs it has not got - the background estimator is what covers those. +TEST_CASE("BeamCenterFromSpots_DeclinesOnAShortSweep", "[BeamCenter]") { + const DiffractionExperiment x = TestExperiment(Coord(1, 0, 0)); + const DiffractionGeometry truth = OffsetBy(x.GetDiffractionGeometry(), 3.0f, -2.5f); + const Sweep sweep = MakeSweep(x, truth, PAIRS); + + // The first half of the sample only: 30 frames spread over 174 deg, none of them a pair. + std::vector angle_deg; + std::vector spots; + for (size_t frame = 0; frame < sweep.angle_deg.size(); frame += 2) + angle_deg.push_back(sweep.angle_deg[frame]); + for (const auto &spot: sweep.spots) + if (spot.frame % 2 == 0) + spots.push_back({spot.x, spot.y, spot.intensity, spot.frame / 2}); + + CHECK_FALSE(FindBeamCenterFromSpotSymmetry(x, angle_deg, spots).has_value()); +} + +// No beamline can hold the spindle exactly perpendicular to the beam, or exactly on a lab axis - +// and the file always claims it does: every master of the regression set writes the goniometer axis +// as an exact lab vector, so the deviation has to come from the spots or from nowhere. What the +// estimator cannot survive is the part turned ABOUT THE BEAM: both mirror lines turn with it, and a +// mirror taken about the lab axis instead smears the vote by twice that angle times the other +// coordinate until a neighbouring comb tooth outvotes the true one. It does not decline when that +// happens - it answers, at the same sigma as a good fit. +// +// The other component, the spindle tipped towards the beam, is here too and must not be corrected +// in the same way: mirroring about the spindle itself would move the direct beam by twice that +// angle times the distance. Both mirror planes contain the beam whatever the spindle does. +TEST_CASE("BeamCenterFromSpots_SurvivesASpindleOffPerpendicular", "[BeamCenter]") { + const float azimuth = 0.002f; // turned about the beam + const float tipped = 0.002f; // and out of the plane perpendicular to it + const bool vertical = GENERATE(false, true); + const Coord in_plane = vertical ? Coord(-std::sin(azimuth), std::cos(azimuth), 0) + : Coord(std::cos(azimuth), std::sin(azimuth), 0); + const Coord spindle = in_plane * std::cos(tipped) + Coord(0, 0, 1) * std::sin(tipped); + + const DiffractionExperiment x = TestExperiment(spindle.Normalize()); + const DiffractionGeometry truth = OffsetBy(x.GetDiffractionGeometry(), 3.0f, -2.5f); + const Sweep sweep = MakeSweep(x, truth, PAIRS); + + SpindleEstimate fitted; + const auto estimate = FindBeamCenterFromSpotSymmetry(x, sweep.angle_deg, sweep.spots, &fitted); + + REQUIRE(estimate.has_value()); + // The azimuth is measured, not assumed, and the vote is taller for it. + CHECK(fitted.azimuth_rad == Catch::Approx(azimuth).margin(2e-4)); + CHECK(fitted.vote_excess > fitted.vote_excess_nominal); + CHECK(estimate->beam_x_pxl == Catch::Approx(truth.GetBeamX_pxl()).margin(0.5)); + CHECK(estimate->beam_y_pxl == Catch::Approx(truth.GetBeamY_pxl()).margin(0.5)); + + // And this is the failure it exists to remove: without it the same sweep is answered, with a + // sigma that says nothing is wrong, further away. How much further depends on what else is in + // the estimator - once the second crossing is guarded by its timing test the comb mis-pick that + // made this several pixels is already gone, and what is left is the azimuth's own smearing of + // the Friedel vote. So the claim is a ratio, which is what the correction actually owns. + const auto uncorrected = FindBeamCenterFromSpotSymmetry(x, sweep.angle_deg, sweep.spots); + REQUIRE(uncorrected.has_value()); + const float corrected_error = std::hypot(estimate->beam_x_pxl - truth.GetBeamX_pxl(), + estimate->beam_y_pxl - truth.GetBeamY_pxl()); + CHECK(std::hypot(uncorrected->beam_x_pxl - truth.GetBeamX_pxl(), + uncorrected->beam_y_pxl - truth.GetBeamY_pxl()) > 3.0f * corrected_error); +} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 65eb3eb4..4439900f 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -95,6 +95,8 @@ ADD_EXECUTABLE(jfjoch_test CCTest.cpp MultiLatticeSearchTest.cpp LatticeReductionTest.cpp + BeamCenterFromBackgroundTest.cpp + BeamCenterFromSpotsTest.cpp ) target_link_libraries(jfjoch_test Catch2WithMain JFJochBroker JFJochReceiver JFJochReader JFJochStreamWriter