leonarski_fandClaude Opus 5
1d1d7fd686
The centrosymmetry score of the pre-scan projection is a self-convolution, so one FFT set scores every candidate centre on the detector at half-pixel spacing - the cost is O(N log N) and independent of how far the header centre is from the truth, where every existing search pays per pixel of error. BeamCenterFFT computes the 2D point-inversion score (a background measurement, the physics FindBeamCenterFromBackground fits locally) and the two 1D line-mirror scores (for the along-spindle coordinate this is exact Friedel physics, sharp exactly where the indexing count is blind), and returns a non-maximum-suppressed shortlist plus a peak-to-runner-up margin per surface - a shortlist and a margin, never a centre: measured on 17 datasets the surface can be locally flat (~25 px), and the margin is what says so. Measured on the gross-header cases that motivate it (offline, float64 reference): a header 351 px wrong is captured at rank 1 within 1.0 px from 30 frames, and still from a 30 deg wedge; two 73 px placeholder headers at rank 1 within 0.5 px; a dataset whose header is ~170 px wrong but which no centre can index is flagged by the lowest margin of the set (0.8 %) instead of being answered confidently. Capture survives 30-120 deg wedges on all four sets tried (union of 8 candidates within 6 px everywhere). Numerics: fftwf (the tree's FFTW is single precision) with the valid-pixel mean subtracted before the transform - the masked Pearson is exactly invariant under a global shift, and the subtraction removes the large-term cancellation - and the per-element combination done in double. Against the float64 numpy reference the shortlist positions are identical and the margins agree to 3e-6 absolute on the two controls whose margins are 0.3 %. An overlap with no image variance (a mirrored empty region) carries no evidence and minted r values of 4-274 in float64 as much as float32; such centres are now not scored (VARIANCE_FLOOR). Wired report-only into the pre-scan next to the background estimate: one log line with the strongest candidates, the margins and the wall time. Nothing consumes it. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_011GxZqDiFP3KqriBhNdcR56
72 lines
3.0 KiB
C++
72 lines
3.0 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include <catch2/catch_all.hpp>
|
|
|
|
#include <cmath>
|
|
#include <random>
|
|
|
|
#include "../image_analysis/geom_refinement/BeamCenterFFT.h"
|
|
|
|
namespace {
|
|
|
|
// An isotropic scattered background about (cx, cy): a decaying continuum with a solvent-ring bump,
|
|
// the mean of ~60 frames' worth of counting noise on top, NAN outside the detector's usable area.
|
|
// A horizontal dead stripe stands in for a module gap.
|
|
std::vector<float> Synthesise(int64_t w, int64_t h, float cx, float cy, float shadow_sector) {
|
|
std::vector<float> mean(static_cast<size_t>(w * h), NAN);
|
|
std::mt19937 rng(20260912);
|
|
std::normal_distribution<float> gauss(0.0f, 1.0f);
|
|
for (int64_t y = 0; y < h; y++) {
|
|
if (y >= 300 && y < 320)
|
|
continue; // module gap
|
|
for (int64_t x = 0; x < w; x++) {
|
|
const float r = std::hypot(static_cast<float>(x) - cx, static_cast<float>(y) - cy);
|
|
const float t = (r - 260.0f) / 30.0f;
|
|
float value = 120.0f * std::exp(-r / 220.0f) + 45.0f * std::exp(-0.5f * t * t) + 5.0f;
|
|
const float phi = std::atan2(static_cast<float>(y) - cy, static_cast<float>(x) - cx);
|
|
if (phi > 0.0f && phi < 1.0f)
|
|
value *= shadow_sector;
|
|
mean[static_cast<size_t>(y * w + x)] = value + gauss(rng) * std::sqrt(value / 60.0f);
|
|
}
|
|
}
|
|
return mean;
|
|
}
|
|
|
|
} // namespace
|
|
|
|
TEST_CASE("BeamCenterFFT_SyntheticBackground", "[BeamCenter]") {
|
|
// The centre is nowhere near the image middle - the score must not care.
|
|
const int64_t w = 900, h = 700;
|
|
const float cx = 297.25f, cy = 411.5f;
|
|
|
|
SECTION("clean isotropic background") {
|
|
const auto mean = Synthesise(w, h, cx, cy, 1.0f);
|
|
const auto r = BeamCenterFFTScore(w, h, mean);
|
|
|
|
REQUIRE(!r.point.empty());
|
|
REQUIRE(!r.line_x.empty());
|
|
REQUIRE(!r.line_y.empty());
|
|
// Whole-detector capture: the strongest point-symmetry centre is the true one.
|
|
CHECK(std::hypot(r.point[0].beam_x_pxl - cx, r.point[0].beam_y_pxl - cy) < 2.0f);
|
|
CHECK(std::abs(r.line_x[0].beam_x_pxl - cx) < 2.0f);
|
|
CHECK(std::abs(r.line_y[0].beam_y_pxl - cy) < 2.0f);
|
|
// A clean unimodal surface separates its argmax from the runner-up by a wide margin.
|
|
CHECK(r.margin_point > 0.05f);
|
|
}
|
|
|
|
SECTION("a shadowed sextant keeps the true centre in the shortlist") {
|
|
const auto mean = Synthesise(w, h, cx, cy, 0.5f);
|
|
const auto r = BeamCenterFFTScore(w, h, mean);
|
|
|
|
REQUIRE(!r.point.empty());
|
|
float best = 1e9f;
|
|
for (const auto &c : r.point)
|
|
best = std::min(best, std::hypot(c.beam_x_pxl - cx, c.beam_y_pxl - cy));
|
|
// The argmax is pulled by the asymmetry (measured ~9 px here - the same shadow bias the
|
|
// background estimator carries); the shortlist must still hold the answer within the
|
|
// 12 px a re-indexing verifier accepts.
|
|
CHECK(best < 12.0f);
|
|
}
|
|
}
|