From d75102b9b43ef985ecf6dbfa840ae9faa8d6c019 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 12 Sep 2026 20:36:15 +0200 Subject: [PATCH] beam centre: the FFT captures the whole detector, the walk refines what it captured --- .../BeamCenterFromBackground.cpp | 42 +++++++++++- .../BeamCenterFromBackground.h | 41 +++++++++++- rugnux/Rugnux.cpp | 64 ++++++++----------- 3 files changed, 106 insertions(+), 41 deletions(-) diff --git a/image_analysis/geom_refinement/BeamCenterFromBackground.cpp b/image_analysis/geom_refinement/BeamCenterFromBackground.cpp index 7c6227960..9129589bc 100644 --- a/image_analysis/geom_refinement/BeamCenterFromBackground.cpp +++ b/image_analysis/geom_refinement/BeamCenterFromBackground.cpp @@ -82,7 +82,8 @@ float median_of(std::vector &v) { std::optional FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const PixelMask &mask, - const std::vector &mean, size_t nthreads) { + const std::vector &mean, size_t nthreads, + std::optional> start) { if (nthreads == 0) nthreads = std::max(1u, std::thread::hardware_concurrency()); const auto W = static_cast(experiment.GetXPixelsNumConv()); @@ -106,8 +107,8 @@ FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const Pixe 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(); + float beam_x = start ? start->first : geom.GetBeamX_pxl(); + float beam_y = start ? start->second : geom.GetBeamY_pxl(); constexpr int n_cells = RADIAL_BINS * SECTORS; std::vector cell_of(n_pixels); @@ -359,3 +360,38 @@ FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const Pixe return BeamCenterEstimate{beam_x, beam_y, std::max({sigma_x, sigma_y, std::hypot(step_x, step_y)})}; } + +std::optional +FindBeamCenter(const DiffractionExperiment &experiment, const PixelMask &mask, + const std::vector &mean, size_t nthreads, BeamCenterFFTResult *capture) { + // The image the capture scores: the projection with everything the mask excludes taken out. + // The beam stop is in that mask by the time this runs, and it is the part that matters - an + // opaque region on one side of the beam is centrosymmetric about ITS own centre, and where it + // is large enough that preference beats the background's. + std::vector scored = mean; + const auto &pixel_mask = mask.GetMask(experiment); + if (pixel_mask.size() == scored.size()) + for (size_t i = 0; i < scored.size(); i++) + if (pixel_mask[i] != 0) + scored[i] = NAN; + + const BeamCenterFFTResult fft = FindBeamCenterFFT(experiment, scored); + if (capture) + *capture = fft; + if (fft.point.empty()) + return FindBeamCenterFromBackground(experiment, mask, mean, nthreads); + + const std::pair seed{fft.point[0].beam_x_pxl, fft.point[0].beam_y_pxl}; + if (auto refined = FindBeamCenterFromBackground(experiment, mask, mean, nthreads, seed)) + return refined; + + // The walk has nothing to fit AT THE CAPTURE. That is a statement about the capture as much as + // about the background: a centre far enough out leaves too few fully covered rings inside the + // fitted band for the profile to exist at all. Measured on one small-molecule set whose capture + // lands 536 px out - the walk from the file's centre answers there to 2.9 px, and taking the + // capture instead would have been the composition's only regression on the corpus. So ask the + // walk again where it used to be asked, and keep the capture only where neither start works. + if (auto from_file = FindBeamCenterFromBackground(experiment, mask, mean, nthreads)) + return from_file; + return BeamCenterEstimate{seed.first, seed.second, BEAM_CENTER_CAPTURE_SIGMA_PXL}; +} diff --git a/image_analysis/geom_refinement/BeamCenterFromBackground.h b/image_analysis/geom_refinement/BeamCenterFromBackground.h index 7d9451818..350cb6ddd 100644 --- a/image_analysis/geom_refinement/BeamCenterFromBackground.h +++ b/image_analysis/geom_refinement/BeamCenterFromBackground.h @@ -4,8 +4,10 @@ #pragma once #include +#include #include +#include "BeamCenterFFT.h" #include "../../common/DiffractionExperiment.h" #include "../../common/PixelMask.h" @@ -34,6 +36,43 @@ struct BeamCenterEstimate { // `mean` is a per-pixel projection over a few tens of frames, NAN where no frame contributed. // nthreads = 0 asks for all hardware threads. The pixels are split into a fixed number of row // blocks whatever that count is, so the answer does not depend on it. +// +// `start` is where the walk begins; the centre in the file when it is not given. The walk advances +// by a bounded distance per iteration, so where it starts decides how much of its budget is spent +// travelling and - on a surface with more than one basin - which fixed point it can reach at all. std::optional FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const PixelMask &mask, - const std::vector &mean, size_t nthreads = 0); + const std::vector &mean, size_t nthreads = 0, + std::optional> start = {}); + +// The precision of a centre that is the FFT capture alone, with no walk behind it. The capture is +// a half-pixel grid position read off a surface, measured over 75 rotation datasets at a median +// 3.2 px and a 90th percentile of 11 px from the truth, so this is what it knows the centre to - +// a capture precision and not a fit precision. It is deliberately far above the ceiling the +// callers adopt a centre on: a capture is evidence about where the beam is, not a measurement of +// where it is. +constexpr float BEAM_CENTER_CAPTURE_SIGMA_PXL = 5.0f; + +// The beam centre from the background, captured globally and then refined locally. +// +// The walk in FindBeamCenterFromBackground is a good local refiner and a poor searcher: it moves +// one or two pixels per iteration, it is seeded at the centre in the file, and on a background +// whose isotropy is broken it has a second basin to fall into. The FFT score is the opposite - +// it evaluates EVERY candidate centre on the detector in one transform set, at a cost that does +// not depend on how wrong the file is, but it returns a half-pixel grid position and no sigma. +// Composing them takes the reach from one and the precision from the other: the capture chooses +// the basin, the walk finishes inside it and reports what it knows the answer to. +// +// The shadow of the beam stop is blanked out of the image the capture scores. A one-sided opaque +// region imposes a centrosymmetry of its own that can beat the background's - measured, an umbra +// over 9.6 % of the detector put the capture 48 px out, and masking it put it back to 1.1 px, +// while a random mask of the same area changed nothing. +// +// Where the walk declines at the capture the walk is asked again from the centre in the file, and +// only where neither start gives it something to fit does the capture stand alone, at +// BEAM_CENTER_CAPTURE_SIGMA_PXL - which is what lets a caller that only needs a hypothesis to test +// still get one. Returns nothing only where the capture has no candidate either. +std::optional +FindBeamCenter(const DiffractionExperiment &experiment, const PixelMask &mask, + const std::vector &mean, size_t nthreads = 0, + BeamCenterFFTResult *capture = nullptr); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 138c6ab16..75801bd35 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1215,13 +1215,33 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // second pass of a two-pass run the shadow is already in the mask, so this does not run again - // by then the centre has been post-refined and the file's value is no longer what is in it. if (want_shadow && config_.beam_center_check) { - background_center_ = FindBeamCenterFromBackground(experiment_, pixel_mask_, - finder.GetMeanProjection(), config_.nthreads); + const auto t0 = std::chrono::steady_clock::now(); + BeamCenterFFTResult capture; + background_center_ = FindBeamCenter(experiment_, pixel_mask_, finder.GetMeanProjection(), + config_.nthreads, &capture); measured_beam_center_ = background_center_; + const auto capture_ms = std::chrono::duration_cast( + std::chrono::steady_clock::now() - t0).count(); + if (!capture.point.empty()) + logger.Info("Beam centre capture: the strongest point-symmetry centre on the WHOLE " + "detector is ({:.1f},{:.1f}), {:.1f} px from the file's; line mirrors put x " + "at {:.1f} and y at {:.1f}; margins {:.1f}%/{:.1f}%/{:.1f}% (point/x/y, " + "under ~3% the surface is flat). One transform set scores every candidate " + "centre, so this costs the same however wrong the file is - {} ms including " + "the walk that refines it", + capture.point[0].beam_x_pxl, capture.point[0].beam_y_pxl, + std::hypot(capture.point[0].beam_x_pxl - experiment_.GetBeamX_pxl(), + capture.point[0].beam_y_pxl - experiment_.GetBeamY_pxl()), + capture.line_x.empty() ? std::numeric_limits::quiet_NaN() + : capture.line_x[0].beam_x_pxl, + capture.line_y.empty() ? std::numeric_limits::quiet_NaN() + : capture.line_y[0].beam_y_pxl, + 100.0f * capture.margin_point, 100.0f * capture.margin_line_x, + 100.0f * capture.margin_line_y, capture_ms); const float need = BeamCenterNeed_pxl(experiment_); if (!background_center_) { - logger.Info("Beam centre check: the background is too flat to place the centre on this " - "run - the file's ({:.2f},{:.2f}) is neither confirmed nor contradicted", + logger.Info("Beam centre check: neither the capture nor the walk places the centre on " + "this run - the file's ({:.2f},{:.2f}) is neither confirmed nor contradicted", experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl()); } else { const float dx = background_center_->beam_x_pxl - experiment_.GetBeamX_pxl(); @@ -1233,7 +1253,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // and what this particular geometry asks of the centre. The last is what turns a bare pixel // count into a verdict - 3 px is nothing at 100 mm and a lost lattice at 500 mm. logger.Info("Beam centre check: the file says ({:.2f},{:.2f}), the scattered background " - "says ({:.2f},{:.2f}) +- {:.2f} px - a difference of {:.2f} px, against the " + "puts it at ({:.2f},{:.2f}) +- {:.2f} px - a difference of {:.2f} px, against the " "{:.2f} px this geometry asks the centre to be right to. Measured only; the " "run keeps the centre in the file", experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(), @@ -1263,36 +1283,6 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru across ? "it is mostly along the spindle, where nothing announces itself" : "this is not a rotation sweep, so there is no spindle to resolve it on"); } - - // The same projection scored globally: one FFT set evaluates the centrosymmetry of the - // background about EVERY candidate centre on the detector, so its cost is the same whether - // the header is right, 2 px out or 351 px out - where the walk above the first pass would - // otherwise pay per pixel of error. Measured on every run and committed on none: a - // shortlist and a margin, never a centre (see BeamCenterFFT.h - the surface can be locally - // flat, and a margin under a few percent means the argmax is a coin toss among the listed - // candidates, not a measurement). - const auto t0 = std::chrono::steady_clock::now(); - const auto fft = FindBeamCenterFFT(experiment_, finder.GetMeanProjection()); - const auto ms = std::chrono::duration_cast( - std::chrono::steady_clock::now() - t0) - .count(); - if (!fft.point.empty()) { - const float fx = fft.point[0].beam_x_pxl, fy = fft.point[0].beam_y_pxl; - const float from_header = std::hypot(fx - experiment_.GetBeamX_pxl(), - fy - experiment_.GetBeamY_pxl()); - logger.Info("Beam centre FFT: strongest point-symmetry centre ({:.1f},{:.1f}), {:.1f} px " - "from the file's; line mirrors put x at {:.1f} and y at {:.1f}; margins " - "{:.1f}%/{:.1f}%/{:.1f}% (point/x/y, under ~3% the surface is flat and the " - "shortlist is a set of hypotheses, not an answer). Whole detector in {} ms; " - "measured only, nothing consumes it", - fx, fy, from_header, - fft.line_x.empty() ? std::numeric_limits::quiet_NaN() - : fft.line_x[0].beam_x_pxl, - fft.line_y.empty() ? std::numeric_limits::quiet_NaN() - : fft.line_y[0].beam_y_pxl, - 100.0f * fft.margin_point, 100.0f * fft.margin_line_x, - 100.0f * fft.margin_line_y, ms); - } } if (!want_beam_center) @@ -1415,8 +1405,8 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // same geometry, same projection - so this asks for it again only where it was not run. estimate = background_center_ ? background_center_ - : FindBeamCenterFromBackground(experiment_, pixel_mask_, - finder.GetMeanProjection(), config_.nthreads); + : FindBeamCenter(experiment_, pixel_mask_, finder.GetMeanProjection(), + config_.nthreads); } if (!estimate) { logger.Info("Beam centre: not measurable, keeping ({:.2f},{:.2f})",