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

beam centre: measure it on every run, and try it when the header fails

Browse Source 
The header beam centre is wrong by more than the geometry absorbs on
two thirds of foreign depositions, and nothing measured it. The solvent
ring already gives it away: the background projection the beam-stop
pre-scan builds is enough to fit the centre, so the measurement is a fit
over an array the run has already paid for.

It is reported on every run and committed on none. The first pass runs
again at the measured centre and the two lattices are compared; the
measured centre is taken only where the header indexes nothing, and a
disagreement is reported rather than resolved, because the only arbiter
available at that stage is the frame count and it is inverted.

Recovers two depositions whose header is the geometric detector centre.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01T3yNBXk4wKdMZy1ak2NY7f
This commit is contained in:
leonarski_f
2026-08-30 14:07:13 +02:00
co-authored by Claude Opus 5
parent 65dfdb1324 bf41b1fab1
commit f75cf72d1a
6 changed files with 633 additions and 26 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/geom_refinement/BeamCenterFromSpots.cpp
+27 -12
View File
@@ -6,7 +6,9 @@
#include "../../common/JFJochMath.h" // PI
#include <algorithm>
#include <array>
#include <cmath>
#include <future>
#include <numeric>
#include <tuple>
@@ -888,19 +890,32 @@ FindBeamCenterFromSpotSymmetry(const DiffractionExperiment &experiment,
// nothing here can say which is the crystal's, so the scatter of the frame pairs - which stays
// small for either of them - is not the uncertainty and this is. The spindle is re-fitted from
// each start, so a fit that depends on where the search began is part of what is reported.
//
// The four run at once. Each takes its own copy of the geometry and only reads the spots, so they
// share nothing, and what is wanted from them is a MAX - which does not care in what order they
// finish. This is the whole of the estimator's cost: Estimate() is a brute-force grid over the
// spindle, and asking for the uncertainty runs it five times.
const std::array<std::pair<float, float>, 4> starts{{{CONSISTENCY_START_PXL, 0.0f},
{-CONSISTENCY_START_PXL, 0.0f},
{0.0f, CONSISTENCY_START_PXL},
{0.0f, -CONSISTENCY_START_PXL}}};
std::vector<std::future<float>> restarts;
restarts.reserve(starts.size());
for (const auto &[dx, dy]: starts)
restarts.push_back(std::async(std::launch::async, [&, dx = dx, dy = dy] {
DiffractionGeometry from = geom;
from.BeamX_pxl(estimate->beam_x_pxl + dx).BeamY_pxl(estimate->beam_y_pxl + dy);
SpindleEstimate again;
if (const auto other = Estimate(from, *goniometer, frame_angle_deg, spots,
spindle_estimate ? &again : nullptr))
return std::hypot(other->beam_x_pxl - estimate->beam_x_pxl,
other->beam_y_pxl - estimate->beam_y_pxl);
return 0.0f;
}));
float spread = 0.0f;
for (const auto &[dx, dy]: {std::pair<float, float>{CONSISTENCY_START_PXL, 0.0f},
{-CONSISTENCY_START_PXL, 0.0f},
{0.0f, CONSISTENCY_START_PXL},
{0.0f, -CONSISTENCY_START_PXL}}) {
DiffractionGeometry from = geom;
from.BeamX_pxl(estimate->beam_x_pxl + dx).BeamY_pxl(estimate->beam_y_pxl + dy);
SpindleEstimate again;
if (const auto other = Estimate(from, *goniometer, frame_angle_deg, spots,
spindle_estimate ? &again : nullptr))
spread = std::max(spread, std::hypot(other->beam_x_pxl - estimate->beam_x_pxl,
other->beam_y_pxl - estimate->beam_y_pxl));
}
for (auto &restart: restarts)
spread = std::max(spread, restart.get());
estimate->sigma_pxl = std::max(estimate->sigma_pxl, spread);
return estimate;
}
rugnux/Rugnux.cpp
+500 -14
View File
@@ -85,6 +85,54 @@ namespace {
// out, and on none of the 38 as they are.
constexpr float MAX_BEAM_CENTER_SIGMA_PXL = 1.0f;
// The longest cell axis assumed where no cell is known yet. It has to be LARGE: a long axis makes
// the tolerance below small, and the ceiling is the larger of it and the flat bound above, so a
// large default is the one that changes nothing.
constexpr float UNKNOWN_CELL_A_MAX_A = 200.0f;
// How far the beam centre may be wrong before the FIRST PASS stops indexing. A centre error is
// fixed in the LAB frame, so accumulating a sweep smears every reciprocal-lattice point around a
// circle and the FFT amplitude at an axis of length a is multiplied by J0(2 pi delta p a/(D
// lambda)); this is the displacement at which that has fallen to 0.70, and the first zero is
// 2.09x further out. Over the in-house and non-SLS corpora it runs from 0.6 px to 9 px - a
// fifteen-fold spread, which is why the flat bound above is only ever used as a FLOOR under it.
//
// It is printed, and it is deliberately not a relevance test. It says how wrong the header may
// be; the estimator's sigma says how well the estimator knows its own answer, and gating one on
// the other rejects a centre correct to 0.03 px because its sigma was 0.98. Nor is a move that is
// small against it a move that can be ignored: measured on real data, a 0.12 px change of centre
// - 0.03x of this - is the difference between the deposited cell and a halved axis.
float BeamCenterNeed_pxl(const DiffractionExperiment &experiment) {
const auto cell = experiment.GetUnitCell();
const float a_max = cell ? std::max({cell->a, cell->b, cell->c}) : UNKNOWN_CELL_A_MAX_A;
return 0.183f * experiment.GetDetectorDistance_mm() * experiment.GetWavelength_A()
/ (experiment.GetPixelSize_mm() * a_max);
}
// A detector-plane displacement resolved ACROSS the spindle - the only component the tolerance
// above is about. Along the spindle the error translates the derotated cloud rigidly and the FFT
// amplitude cannot see it at all; across it the peaks smear and the amplitude falls. (The lattice
// fit that follows the FFT is not translation-invariant, so "along the spindle" is free for the
// transform and not for the run - which is why this is reported and not used to decide anything.)
// The detector direction wanted is the one whose lab image is perpendicular to the spindle n:
// (ux e1 + uy e2).n = 0, so u is proportional to (-(e2.n), e1.n).
std::optional<float> AcrossSpindle_pxl(const DiffractionExperiment &experiment, float dx, float dy) {
const auto goniometer = experiment.GetGoniometer();
if (!goniometer || !goniometer->IsScanning())
return {};
const auto geom = experiment.GetDiffractionGeometry();
const float beam_x = geom.GetBeamX_pxl(), beam_y = geom.GetBeamY_pxl();
const Coord n = goniometer->GetAxis().Normalize();
const Coord origin = geom.LabCoord(beam_x, beam_y);
const Coord e1 = geom.LabCoord(beam_x + 1.0f, beam_y) - origin;
const Coord e2 = geom.LabCoord(beam_x, beam_y + 1.0f) - origin;
const float ux = -(e2 * n), uy = e1 * n;
const float length = std::hypot(ux, uy);
if (length < 1e-3f)
return {};
return std::abs(dx * ux + dy * uy) / length;
}
// 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;
@@ -303,6 +351,49 @@ namespace {
}
// What the file says about the beam centre, before a frame is read. Nothing here measures anything -
// it reads the value that is there and says where it looks like it came from - and it is free, which
// is why it runs on every run whether or not a centre was asked for.
//
// A centre ON A MASKED PIXEL is the one of these worth acting on: a real beam does not sit on a dead
// pixel or in a module gap, so the value is wrong about something. Over the 29 non-SLS datasets whose
// outcome is on record (base rate of a bad outcome 48%) it fires on 4 and 3 of those 4 went wrong,
// against 2 false alarms in 42 in-house masters. What it does NOT promise is that correcting the
// centre repairs the run - on one of the three it does not.
//
// The other two are PROVENANCE, not prediction, and are said rather than warned about. A value equal
// to the geometric detector centre, or a whole number of pixels in both coordinates, was typed rather
// than measured - true of 29 of 42 in-house masters, and those are out by a median 3.5 px against
// 0.7 px for the rest - but a typed value is often close enough that as predictors of harm both sit
// at or below the base rate. They are still worth saying, because the spot-symmetry search reaches
// only about +-55 px from wherever it starts: a placeholder tens of pixels out is unreachable by
// construction, all five of its searches agree, and the small sigma they report says nothing is
// wrong. Only the background estimator recovers that case, and only a human reading this line knows
// to look.
static void LogHeaderBeamCenter(const DiffractionExperiment &experiment, const PixelMask &pixel_mask) {
Logger logger("Rugnux");
const float beam_x = experiment.GetBeamX_pxl();
const float beam_y = experiment.GetBeamY_pxl();
const auto width = static_cast<int>(experiment.GetXPixelsNumConv());
const auto height = static_cast<int>(experiment.GetYPixelsNumConv());
if (std::hypot(beam_x - static_cast<float>(width) / 2.0f,
beam_y - static_cast<float>(height) / 2.0f) < 1.0f)
logger.Info("Beam centre ({:.2f},{:.2f}) is the geometric centre of the detector - a value "
"written rather than measured", beam_x, beam_y);
else if (beam_x == std::round(beam_x) && beam_y == std::round(beam_y))
logger.Info("Beam centre ({:.0f},{:.0f}) is a whole number of pixels in both coordinates - a "
"value written rather than measured", beam_x, beam_y);
const int x = static_cast<int>(std::lround(beam_x));
const int y = static_cast<int>(std::lround(beam_y));
if (x >= 0 && y >= 0 && x < width && y < height
&& pixel_mask.GetMask()[static_cast<size_t>(y) * width + x] != 0)
logger.Warning("Beam centre ({:.2f},{:.2f}) lands on a masked pixel. A real beam does not, so "
"the value in the file is suspect; --estimate-beam-center measures it",
beam_x, beam_y);
}
Rugnux::Rugnux(JFJochReader &reader, DiffractionExperiment experiment,
PixelMask pixel_mask, ProcessConfig config)
: reader_(reader), experiment_(std::move(experiment)),
@@ -312,6 +403,8 @@ Rugnux::Rugnux(JFJochReader &reader, DiffractionExperiment experiment,
// already carries one starts clear; the user mask (bit 8) is left as it was loaded.
pixel_mask_.ClearBeamStopMask(experiment_);
LogHeaderBeamCenter(experiment_, pixel_mask_);
// A manually asserted stage calibration applies to everything, before anything reads an angle.
if (config_.rotation_scale.has_value())
if (const auto g = experiment_.GetGoniometer()) {
@@ -740,14 +833,88 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
const auto shadow = finder.GetMask(config_.nthreads);
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<double>(shadowed) / static_cast<double>(shadow.size()),
finder.GetFrameCount());
// Not "this beamline has no beam stop". The search floods outwards from seeds within 4 px of
// the ASSUMED beam centre, so a centre far enough out puts every seed in a module gap and the
// flood finds nothing - which makes an empty shadow a statement about the centre. It is the
// sharpest free predicate there is: over 29 non-SLS datasets it fires on 3, all 3 of which had
// a bad outcome, with no false alarm in 39 in-house sweeps (which report 0.44-1.51 % masked).
if (shadowed == 0)
logger.Warning("Beam stop: no shadow found in {} images. The search starts from the beam "
"centre in the file, so this points at that centre rather than at the "
"experiment; --estimate-beam-center measures it", finder.GetFrameCount());
else
logger.Info("Beam stop shadow: {} pixels ({:.2f}% of the detector) found in {} images",
shadowed, 100.0 * static_cast<double>(shadowed) / static_cast<double>(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.
// What the scattered background makes the centre, measured on EVERY run and committed on none.
// The fit reads the projection the beam-stop pre-pass has just built and no frame of its own, so
// it costs a fraction of a second on top of a pass that has already run - which is what makes it
// affordable to do always, and doing it always is the point. A header beam centre is the metadata
// field the field jokes about, and nothing in any package says how wrong a given one is: on a run
// that indexes and merges perfectly well this line still reports that the file claims a centre the
// data put 73 px away. The run is not touched by it. What consumes the number is the second first
// pass in RunPipeline, and --estimate-beam-center's fall-through below, which is the same fit.
//
// Only where the projection exists: with --detect-beam-stop=off nothing accumulates it, and
// projecting every frame for this alone is the largest per-frame cost of the pre-scan. On the
// 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());
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",
experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl());
} else {
const float dx = background_center_->beam_x_pxl - experiment_.GetBeamX_pxl();
const float dy = background_center_->beam_y_pxl - experiment_.GetBeamY_pxl();
const float moved = std::hypot(dx, dy);
const auto across = AcrossSpindle_pxl(experiment_, dx, dy);
// Everything a reader needs to judge the header, on one line: what the file says, what the
// background says, how far apart they are, how well the background knows its own answer,
// 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 "
"{:.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(),
background_center_->beam_x_pxl, background_center_->beam_y_pxl,
background_center_->sigma_pxl, moved, need);
// Three sigma, because the fit is a PRECISION and not an accuracy: on a background with no
// curvature there is nothing to separate shift from amplitude and it follows the noise in
// g', confidently. A move it cannot resolve is not evidence about the header either way.
if (moved <= 3.0f * background_center_->sigma_pxl)
logger.Info("Beam centre check: that is under three times the fit's own sigma, so the "
"file's centre is as good as this measurement can tell - the second first "
"pass below indexes it anyway, because a fraction of a pixel is enough to "
"decide a cell and no sigma small enough to gate on would notice");
else if (across && *across > need)
// Info and not a warning: `need` is where the FFT amplitude at a 200 A axis has fallen
// to 0.70, which is a floor under what matters and not a threshold for failure - it is
// exceeded on two of three ordinary in-house sweeps whose lattice comes out right. What
// is worth a warning is what the second pass below then finds, not this.
logger.Info("Beam centre check: {:.2f} px of the difference is ACROSS the spindle, more "
"than the {:.2f} px this geometry absorbs, and a centre error is fixed in "
"the lab frame - it smears the accumulated reciprocal-space cloud, so the "
"first pass can lose peaks or take an axis harmonic. The second first pass "
"below indexes both centres and reports whether it did", *across, need);
else
logger.Info("Beam centre check: that is a real difference, and {} - the second first "
"pass below indexes both centres and compares what they give",
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");
}
}
if (!want_beam_center)
return;
beam_center_placed_ = true;
@@ -760,6 +927,12 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
std::string source = "spot symmetry";
std::optional<BeamCenterEstimate> estimate;
SpindleEstimate spindle;
// The bound an estimate has to meet. One-sided by construction: the flat constant was calibrated
// on a 244 mm / 0.95 A / 0.15 mm geometry and is far too tight for a loose one, so where the
// geometry asks less of the centre the bound is relaxed to what it asks and never below the
// constant. Nothing it accepts today is lost.
const float need = BeamCenterNeed_pxl(experiment_);
const float ceiling = std::max(MAX_BEAM_CENTER_SIGMA_PXL, need);
if (want_spot_symmetry)
estimate = FindBeamCenterFromSpotSymmetry(experiment_, frame_angle_deg, beam_center_spots,
config_.fit_spindle ? &spindle : nullptr);
@@ -772,7 +945,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
// 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)) {
if (want_spot_symmetry && (!estimate || estimate->sigma_pxl > ceiling)) {
// 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.
@@ -844,23 +1017,56 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
"(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) {
if (!estimate || estimate->sigma_pxl > ceiling) {
// Say so. Below about 220 deg of sweep the spot symmetry never clears this bound and the
// background answers instead - measured within 0.23 px of the full-sweep centre at every span
// from 120 to 220 deg - so the fall-through is what covers short sweeps, and until this line
// existed nothing in the log said the spot arm had even been tried.
if (want_spot_symmetry && estimate)
logger.Info("Beam centre: the spot symmetry answers ({:.2f},{:.2f}) at sigma {:.2f} px, "
"over the {:.2f} px ceiling - falling through to the background",
estimate->beam_x_pxl, estimate->beam_y_pxl, estimate->sigma_pxl, ceiling);
else if (want_spot_symmetry)
logger.Info("Beam centre: the spot symmetry does not come out on this sweep - falling "
"through to the background");
source = "background";
estimate = FindBeamCenterFromBackground(experiment_, pixel_mask_, finder.GetMeanProjection());
// The same fit the check above has already done wherever the projection exists - same mask,
// 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());
}
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;
const float dx = estimate->beam_x_pxl - experiment_.GetBeamX_pxl();
const float dy = estimate->beam_y_pxl - experiment_.GetBeamY_pxl();
const float moved = std::hypot(dx, dy);
const auto across = AcrossSpindle_pxl(experiment_, dx, dy);
// Two conditions, and the second is the new one. A move under three times the estimator's own
// sigma is not a measurement of anything, and a centre that is not moved cannot move the
// two-pass loop off its fixed point either. On the 39 rotation regression crystals this adopts
// 38 and keeps the header on one, a 0.74 px move at sigma 0.32.
const bool commit = estimate->sigma_pxl <= ceiling && moved > 3.0f * estimate->sigma_pxl;
const char *verdict = "COMMIT";
if (!commit)
verdict = estimate->sigma_pxl > ceiling ? "reject: sigma over the ceiling (kept header)"
: "reject: the move is under 3 sigma (kept header)";
logger.Info("Beam centre from {}: ({:.2f},{:.2f}) -> ({:.2f},{:.2f}), moved {:.2f} px, "
"sigma {:.2f} px => {}",
"sigma {:.2f} px against a {:.2f} px ceiling => {}",
source, experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(),
estimate->beam_x_pxl, estimate->beam_y_pxl, moved, estimate->sigma_pxl,
commit ? "COMMIT" : "reject (kept header)");
ceiling, verdict);
if (across)
logger.Info("Beam centre: {:.2f} px of that move is across the spindle, where this geometry "
"asks for {:.2f} px. That is a floor and not a tolerance - a move 0.03x of it has "
"decided a halved axis", *across, need);
if (across && *across > need)
logger.Warning("Beam centre: the file's centre is {:.2f} px out across the spindle, more than "
"the {:.2f} px this geometry absorbs - expect the first pass to lose peaks or "
"to take an axis harmonic", *across, need);
if (commit)
experiment_.BeamX_pxl(estimate->beam_x_pxl).BeamY_pxl(estimate->beam_y_pxl);
}
@@ -1879,14 +2085,39 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
// ratio separates a real axis multiplication from a centering coincidence. Volumes are
// primitive (see above), so a pure setting difference is a ratio of 1 and never fires.
bool integer_subcell = false;
if (bp.result.has_value() && !clearly_more && vol > 1.0 && bp.vol > 1.0) {
double harmonic_ratio = 0.0;
if (bp.result.has_value() && vol > 1.0 && bp.vol > 1.0) {
const bool tied = static_cast<float>(score) >= bp.score * 0.9f - 0.5f;
const double ratio = (vol < bp.vol) ? bp.vol / vol : vol / bp.vol;
const double nearest = std::round(ratio);
const bool integer_multiple = nearest >= 2.0 && std::abs(ratio - nearest) < 0.15;
integer_subcell = tied && integer_multiple && vol < bp.vol;
integer_subcell = tied && !clearly_more && integer_multiple && vol < bp.vol;
// sqrt(3) as well as an integer: a hexagonal supercell that lengthens one axis by
// sqrt(3) is the other harmonic seen on real data and its volume ratio is not a
// whole number. This is a REPORT and not a decision - it is deliberately taken
// over the whole pair, `clearly_more` included, because the harmonic is often the
// one that indexes more frames (see below).
if (integer_multiple || std::abs(ratio - std::sqrt(3.0)) < 0.15)
harmonic_ratio = ratio;
}
// Two candidates whose primitive volumes differ by a small integer or by sqrt(3) are
// one lattice and its axis harmonic, and what decides between them is the BEAM CENTRE.
// A centre error along the spindle translates the derotated cloud rigidly, which the
// FFT amplitude cannot see at all, so the peaks stay sharp and the indexer commits -
// with confidence - to a sub-multiple. Measured on real data: a 0.12 px change of
// centre, 0.03x of what the J0 law asks across the spindle, is the whole difference
// between the deposited cell and a halved axis, and the halved one is the answer that
// indexed MORE frames (100.00 % against 99.23 %). So the frame count cannot arbitrate
// this pair, and until this line nothing in the run said what the pair meant.
if (harmonic_ratio > 0.0)
logger.Warning("Scheme '{}' (primitive cell volume {:.0f} A^3) and '{}' ({:.0f} A^3) "
"differ by {:.2f}x - one is an axis harmonic of the other. Which of "
"the two is right is settled by the beam centre to a fraction of a "
"pixel, and the harmonic is often the one that indexes MORE frames, "
"so the frame counts above cannot decide it. Check the beam centre "
"(--estimate-beam-center)", name, vol, bp.name, bp.vol, harmonic_ratio);
if (!bp.result.has_value() || clearly_more || integer_subcell) {
if (integer_subcell)
logger.Info("Scheme '{}' cell (vol {:.0f}) is a {:.0f}x sub-cell of '{}' (vol {:.0f}) at "
@@ -1937,7 +2168,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
GoniometerAxis flipped = *gon;
flipped.Axis(-gon->GetAxis());
experiment_.Goniometer(flipped);
const FirstPass alt = pick_best(*indexer_pool, *indexer);
FirstPass alt = pick_best(*indexer_pool, *indexer);
if (alt.result.has_value() && alt.score > best.score) {
logger.Info("Rotation axis sign: the file's axis indexes {}/{} validation frames, "
"the opposite sign {}/{} - adopting the opposite sign for this run",
@@ -1950,6 +2181,261 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
}
}
// The background centre, tried as a SECOND FIRST PASS.
//
// The estimate itself is free - it is read off the projection --detect-beam-stop has already
// built - and a first pass re-uses the spots it has already found, so indexing the run a
// second time at the measured centre costs about what one rung of the ladder below costs.
// That buys the comparison the on-failure ladder structurally cannot make: a centre wrong
// ALONG the spindle does not fail. It keeps 96-100 % of frames indexed and quietly returns a
// 2x, 3x or sqrt(3) axis harmonic, so no failure trigger ever fires on it. Indexing both
// centres and comparing the two ANSWERS is what can see it.
//
// What must NOT arbitrate the two is the frame count. Acceptance is a fractional-Miller test,
// so a cell twice as long has to place every spot twice as accurately to score the same:
// measured on real data, two centres 0.12 px apart gave the deposited cell at 99.23 % and a
// halved axis at 100.00 %, and the WRONG answer indexed better. A rule of the form "take the
// centre that indexes more" would therefore pick wrong in exactly the case this comparison
// exists for. So only what needs no arbiter is decided here:
//
// * the file's centre does not index a majority -> nothing defends it, so try the measured
// one (and the axis sign with it, see below). A run that indexes nothing has nothing to
// lose; where the measured one is no better, the file's centre is put back;
// * the file's centre indexes and the two agree on the lattice -> keep the file's centre and
// say they agree. This is the overwhelming majority, and saying it is worth the pass: it
// is a free statement that the header is good enough for this crystal;
// * they index and disagree -> say so, with both cells and their primitive-volume ratio,
// and change nothing. Naming a disagreement is useful; resolving it needs an arbiter
// nobody has.
//
// Where the file's centre indexes nothing and the measured one does not either, the axis-sign
// question is asked AGAIN at the measured centre before giving up. The two unknowns are
// discrete and coupled, and the rescue above asked its question at a centre that may be 73 px
// wrong: on the two depositions here whose header is the geometric centre of the detector,
// each error hides the other, both runs sit at 0/60 until both are moved, and moving both
// takes them to 60/60 on the deposited lattice. Asking the pair costs one more first pass on
// a run that has already failed twice.
if (!cancelled_ && config_.beam_center_check && background_center_) {
// Consumed, not re-read: a two-pass rotation run comes back through here with the geometry
// the first pass post-refined, and a pre-scan estimate of the FILE's centre says nothing
// about that one. The pre-scan does not re-measure it either - by then the beam-stop shadow
// is already in the mask, so the projection is not rebuilt.
const BeamCenterEstimate measured = *background_center_;
background_center_.reset();
const float header_x = experiment_.GetBeamX_pxl(), header_y = experiment_.GetBeamY_pxl();
const float measured_x = measured.beam_x_pxl;
const float measured_y = measured.beam_y_pxl;
const float moved = std::hypot(measured_x - header_x, measured_y - header_y);
// No gate on how significant the move is. It is tempting to skip the second pass when the
// two centres are within the fit's own noise - it would save it on a quarter of the runs
// here - but that is deciding at a threshold before the evidence is in, and the evidence
// is cheap. It would also gate away exactly the case worth asking about: measured on real
// data, a 0.12 px change of centre is the whole difference between the deposited cell and
// a halved axis, and 0.12 px is well inside the sigma of any fit that could gate it.
{
const int majority = static_cast<int>(validation.size()) / 2;
experiment_.BeamX_pxl(measured_x).BeamY_pxl(measured_y);
FirstPass alt = pick_best(*indexer_pool, *indexer);
const bool header_indexes = best.result.has_value() && best.score > majority;
bool measured_indexes = alt.result.has_value() && alt.score > majority;
if (!header_indexes) {
// Nothing defends the file's centre: the data do not index at it, and the
// background says it is `moved` px away. So try the measured one, and if that
// does not index either, ask the AXIS-SIGN question again there. The two unknowns
// are discrete and coupled, and the sign rescue above asked its question at a
// centre that may itself be wrong - each error hides the other, and only moving
// both reveals anything.
bool joint = false;
const auto gonio = experiment_.GetGoniometer();
if (!measured_indexes && gonio.has_value() && gonio->IsScanning()) {
GoniometerAxis flipped = *gonio;
flipped.Axis(-gonio->GetAxis());
experiment_.Goniometer(flipped);
const FirstPass both = pick_best(*indexer_pool, *indexer);
if (both.result.has_value() && both.score > majority) {
logger.Warning("Beam centre from the background: the file's centre "
"({:.2f},{:.2f}) indexes {}/{} validation frames and the "
"measured ({:.2f},{:.2f}) {}/{}, but the measured centre "
"WITH the opposite rotation-axis sign indexes {}/{}. The "
"file has both wrong, and neither shows while the other "
"stands; the run takes the measured centre and the opposite "
"sign", header_x, header_y, best.score,
static_cast<int>(validation.size()), measured_x, measured_y,
alt.score, static_cast<int>(validation.size()), both.score,
static_cast<int>(validation.size()));
alt = both;
measured_indexes = true;
joint = true;
} else {
experiment_.Goniometer(*gonio);
}
}
if (measured_indexes) {
if (!joint)
logger.Warning("Beam centre from the background: the file's ({:.2f},{:.2f}) "
"indexes {}/{} validation frames and the measured "
"({:.2f},{:.2f}) indexes {}/{} - the beam centre in this "
"file is wrong, and the run adopts the measured one",
header_x, header_y, best.score,
static_cast<int>(validation.size()), measured_x, measured_y,
alt.score, static_cast<int>(validation.size()));
best = alt;
} else {
// Nothing worked. Put the file's centre back, so the ladder below searches
// around the value the run started with rather than around a fit that has just
// been tested and found no better - and so a run that fails for some other
// reason fails at the geometry it was given.
experiment_.BeamX_pxl(header_x).BeamY_pxl(header_y);
logger.Info("Beam centre check: neither the file's centre ({}/{} frames) nor "
"the measured one {:.2f} px away ({}/{}) indexes this run - the "
"file's centre is kept and this is not a beam-centre failure",
best.score, static_cast<int>(validation.size()), moved, alt.score,
static_cast<int>(validation.size()));
}
} else {
experiment_.BeamX_pxl(header_x).BeamY_pxl(header_y);
if (measured_indexes) {
// Primitive volumes, as everywhere else here: two passes can find the same
// lattice in different settings, and a centred setting's cell is an exact
// integer multiple of its primitive one.
const double ratio = alt.vol < best.vol ? best.vol / alt.vol : alt.vol / best.vol;
const auto &hc = best.result->search_result.conventional.GetUnitCell();
const auto &mc = alt.result->search_result.conventional.GetUnitCell();
const bool same_class =
best.result->search_result.system == alt.result->search_result.system
&& best.result->search_result.centering == alt.result->search_result.centering;
// The primitive volume is what the beam centre decides. A centre error along
// the spindle makes the FFT take an axis harmonic, which changes the volume by
// an integer factor or by sqrt(3); the Bravais class can differ for a reason
// that has nothing to do with the centre, so the two are reported separately.
if (ratio < 1.02 && same_class) {
logger.Info("Beam centre check: indexing at the measured centre "
"({:.2f},{:.2f}) returns the same lattice as the file's - "
"{}-centred {}, cells agreeing to {:.2f} % in volume. The {:.2f} px "
"difference does not decide this crystal's cell",
measured_x, measured_y, alt.result->search_result.centering,
gemmi::crystal_system_str(alt.result->search_result.system),
100.0 * (ratio - 1.0), moved);
} else if (ratio < 1.02) {
// Same cell, different metric symmetry read off it. That is a promotion
// decided on an angular tolerance against the unrefined FFT candidate, not
// a lattice the beam centre chose, so it is worth naming and is not the
// alarming case: the cell is not in question. Both cases in this corpus
// are of this kind, and in both the file's centre is the one that finds
// the symmetry.
logger.Info("Beam centre check: indexing at the measured centre "
"({:.2f},{:.2f}) returns the same CELL as the file's - the "
"primitive volumes agree to {:.2f} % - but a different metric "
"symmetry: {}-centred {} from the file's centre against "
"{}-centred {} from the measured one. The {:.2f} px difference "
"does not decide this crystal's cell; the run keeps the file's "
"centre and the symmetry it found", measured_x, measured_y,
100.0 * (ratio - 1.0), best.result->search_result.centering,
gemmi::crystal_system_str(best.result->search_result.system),
alt.result->search_result.centering,
gemmi::crystal_system_str(alt.result->search_result.system),
moved);
} else {
const double nearest = std::round(ratio);
const bool harmonic = (nearest >= 2.0 && std::abs(ratio - nearest) < 0.15)
|| std::abs(ratio - std::sqrt(3.0)) < 0.15;
logger.Warning("Beam centre check: the two centres DISAGREE about the "
"lattice. The file's ({:.2f},{:.2f}) gives {}-centred {}, "
"{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} (primitive volume "
"{:.0f} A^3, {}/{} frames); the measured ({:.2f},{:.2f}) gives "
"{}-centred {}, {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} "
"(primitive volume {:.0f} A^3, {}/{} frames){}",
header_x, header_y, best.result->search_result.centering,
gemmi::crystal_system_str(best.result->search_result.system),
hc.a, hc.b, hc.c, hc.alpha, hc.beta, hc.gamma, best.vol,
best.score, static_cast<int>(validation.size()),
measured_x, measured_y, alt.result->search_result.centering,
gemmi::crystal_system_str(alt.result->search_result.system),
mc.a, mc.b, mc.c, mc.alpha, mc.beta, mc.gamma, alt.vol,
alt.score, static_cast<int>(validation.size()),
harmonic ? fmt::format(" - the volumes differ by {:.2f}x, so one "
"is an axis harmonic of the other and the "
"beam centre is what settles which",
ratio)
: "");
// Deliberately no decision. The frame counts above are the one statistic
// that must not choose: on the measured case the harmonic indexed MORE.
logger.Warning("Beam centre check: the run keeps the file's centre. The "
"frame counts cannot arbitrate this - a cell twice as long "
"must place every spot twice as accurately to score the "
"same, so the harmonic is often the one that indexes more. "
"Re-run with --estimate-beam-center to measure the centre "
"from the spot symmetry, which is exact geometry");
}
} else {
logger.Info("Beam centre check: the measured centre ({:.2f},{:.2f}) indexes "
"{}/{} validation frames against the file's {}/{} - the file's "
"centre stands", measured_x, measured_y, alt.score,
static_cast<int>(validation.size()), best.score,
static_cast<int>(validation.size()));
}
}
}
}
// Beam centre as an INDEXING HYPOTHESIS. A beam-centre error is not repairable downstream: it
// is fixed in the LAB frame, so accumulating a sweep smears every reciprocal-lattice point
// around a circle and the FFT amplitude at an axis of length a is multiplied by
// J0(2 pi delta p a/(D lambda)). Past the first zero the true axis is gone and its harmonic
// wins, which is what a first-pass axis doubling usually is. It is decidable from the data on
// exactly the count the scheme choice and the axis-sign rescue already use: the right centre
// indexes and the wrong one does not.
//
// BOTH detector directions are searched, and that is the one thing to keep. The J0 law is
// about the FFT AMPLITUDE, which is translation-invariant, so a centre error along the spindle
// is free for the transform - but the step after it, fitting a lattice whose origin is the
// beam, is not, and a rigidly shifted lattice is fitted by a finer one. Measured by injection
// on two crystals: across the spindle the cell stays right and the indexed fraction collapses
// (99 % to 25 %), while along it the run keeps 96-100 % indexed and quietly adopts a 2x, 3x or
// sqrt(3) supercell. The direction that announces itself is not the dangerous one.
//
// Only after a poor pass, so a correct header costs nothing; and spot finding is not repeated,
// the cache being keyed by image and the spot POSITIONS not depending on the centre.
if (!cancelled_ && !beam_center_searched_ && config_.beam_center_search_pxl > 0
&& best.score < 0.5 * static_cast<double>(validation.size())) {
beam_center_searched_ = true;
const float beam_x = experiment_.GetBeamX_pxl(), beam_y = experiment_.GetBeamY_pxl();
const int reach = config_.beam_center_search_pxl;
logger.Info("Beam-centre hypothesis: {}/{} validation frames at ({:.2f},{:.2f}); trying "
"the centre +-{} px along each detector axis, a pixel at a time",
best.score, static_cast<int>(validation.size()), beam_x, beam_y, reach);
// Adopt only a hypothesis that INDEXES - the same majority the guards above test - rather
// than whichever rung scores best. Taking any improvement was tried and is wrong: on a
// dataset where no centre works the ladder wandered to the far end of its range on 10/60
// against 5/60 and reported a 30 px move, which is the count's noise floor.
const int adopt_bar = static_cast<int>(validation.size()) / 2;
bool adopted = false;
float won_x = beam_x, won_y = beam_y;
for (int k = 1; k <= reach && !cancelled_ && !adopted; k++) {
for (const auto &[dx, dy]: {std::pair<int, int>{-k, 0}, {k, 0}, {0, -k}, {0, k}}) {
const float trial_x = beam_x + static_cast<float>(dx);
const float trial_y = beam_y + static_cast<float>(dy);
experiment_.BeamX_pxl(trial_x).BeamY_pxl(trial_y);
FirstPass alt = pick_best(*indexer_pool, *indexer);
logger.Info("Beam-centre hypothesis ({:.2f},{:.2f}): {}/{} frames", trial_x,
trial_y, alt.score, static_cast<int>(validation.size()));
if (alt.result.has_value() && alt.score > adopt_bar && alt.score > best.score) {
best = alt;
won_x = trial_x;
won_y = trial_y;
adopted = true;
break;
}
}
}
experiment_.BeamX_pxl(won_x).BeamY_pxl(won_y);
if (adopted)
logger.Info("Beam centre from indexing: ({:.2f},{:.2f}) -> ({:.2f},{:.2f}), moved "
"{:.2f} px, {}/{} validation frames", beam_x, beam_y, won_x, won_y,
std::hypot(won_x - beam_x, won_y - beam_y), best.score,
static_cast<int>(validation.size()));
}
// Long-axis rescue. When the de-novo cell indexes few validation frames, a long, finely-spaced
// axis was likely lost: the unconstrained FFT either collapsed it to a short sub-multiple or let
// a denser supercell over-fit the accumulated cloud (a small global-orientation error throws the
rugnux/Rugnux.h
+26
View File
@@ -21,6 +21,7 @@
#include "ModelValidation.h" // ModelValidationResult
#include "../image_analysis/scale_merge/SearchSpaceGroup.h" // SearchSpaceGroupResult
#include "../image_analysis/geom_refinement/PostRefine.h" // PostRefineResult
#include "../image_analysis/geom_refinement/BeamCenterFromBackground.h" // BeamCenterEstimate
#include "../image_analysis/rotation_indexer/RotationIndexer.h" // RotationIndexerResult
#include "RugnuxCalibration.h" // CalibrationMethod, CalibrationResult
@@ -101,6 +102,21 @@ struct ProcessConfig {
bool estimate_beam_center = false;
bool fit_spindle = false;
// Beam centre as an INDEXING HYPOTHESIS (--beam-center-search). Runs only after a first pass that
// indexes under half the validation frames, so it costs nothing on a run that works, and it is a
// repair path rather than a measurement: it steps the centre a pixel at a time and keeps the first
// rung that indexes a majority. The value is how many pixels it reaches in each direction; 0 is off.
int beam_center_search_pxl = 12;
// Beam centre from the scattered background, measured on EVERY run and reported (--beam-center-check,
// on by default). The fit reads the projection --detect-beam-stop has already built and no frame of
// its own, so it is close to free, and its result is worth having on a run that comes out right:
// it is what lets a job reprocessing a database say a deposited header is 73 px out. It also
// indexes a SECOND first pass at that centre, where the move is significant, and compares the
// two answers - the one comparison an on-failure trigger cannot make, because a centre wrong
// along the spindle does not fail.
bool beam_center_check = true;
// Adaptive integration radius (--adaptive-integration-radius). When set, the pre-scan measures how
// wide the recorded spots are - r80, the radius holding 80 % of a spot's flux, over a fixed
// aperture that owes nothing to the integrator (SpotWidth.h) - and sets r1 from it. r1 is the
@@ -319,6 +335,16 @@ class Rugnux {
// 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;
// Whether the beam-centre hypothesis ladder has already run in this RUN, not in this pass: the
// second pass starts from the geometry the first one post-refined, so searching again there
// re-asks a question already answered, and where the answer was "no centre works" it doubles the
// cost of the failure.
bool beam_center_searched_ = false;
// What the scattered background makes the beam centre, measured in the pre-scan (--beam-center-check)
// and never committed there: the run keeps the centre it was given. RunPipeline consumes it for
// the second first pass, and clears it, so the two-pass loop does not re-ask on a post-refined
// geometry the pre-scan estimate says nothing about.
std::optional<BeamCenterEstimate> background_center_;
// The recorded spot width the pre-scan measured (config_.adaptive_integration_radius), kept so the
// two-pass rotation run measures it once and both passes integrate at the same radius.
rugnux/rugnux_cli.cpp
+30
View File
@@ -105,6 +105,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 << " --beam-center-search[=N|off] After a first pass that indexes fewer than half the validation frames, try the beam centre a pixel at a time out to N px along each detector axis and keep the first one that indexes a majority. A beam-centre error is fixed in the lab frame, so it smears the accumulated reciprocal-space cloud and the FFT takes an axis harmonic instead of the true axis; nothing downstream repairs that. ON by default (12 px); =off disables. It runs only after a pass that has already failed, so a run that indexes never pays for it" << std::endl;
std::cout << " --beam-center-check[=off] Measure the beam centre from the isotropy of the scattered background on EVERY run, report how far the file's value is from it, and index a second first pass with it to see whether the two centres give the same lattice. The fit reads the projection --detect-beam-stop already builds, so it costs no extra frames. Nothing is committed on a run that indexes: the measured centre is adopted only where the file's indexes nothing and the measured one indexes a majority. ON by default; =off disables" << 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;
@@ -235,6 +237,8 @@ enum {
OPT_BACKGROUND_RADIAL,
OPT_REFINE_GEOMETRY,
OPT_DETECT_BEAM_STOP,
OPT_BEAM_CENTER_SEARCH,
OPT_BEAM_CENTER_CHECK,
OPT_ESTIMATE_BEAM_CENTER,
OPT_FIT_SPINDLE,
OPT_NO_FIT_SPINDLE,
@@ -351,6 +355,8 @@ 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},
{"beam-center-search", optional_argument, nullptr, OPT_BEAM_CENTER_SEARCH},
{"beam-center-check", optional_argument, nullptr, OPT_BEAM_CENTER_CHECK},
{"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},
@@ -712,6 +718,8 @@ static int RunRugnux(int argc, char **argv) {
std::optional<bool> background_radial_arg; // when given: set = force on/off, unset = auto
std::optional<int> detect_beam_stop = 60; // --detect-beam-stop[=N|off]; on by default
bool estimate_beam_center = false; // --estimate-beam-center
int beam_center_search = 12; // --beam-center-search[=N|off]; on by default
bool beam_center_check = true; // --beam-center-check[=off]; on by default
bool fit_spindle = true; // --fit-spindle / --no-fit-spindle
std::optional<int> 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
@@ -814,6 +822,24 @@ static int RunRugnux(int argc, char **argv) {
? parse_number_arg<int>(optarg, "--detect-beam-stop", logger, 1, 1000000)
: 60;
break;
case OPT_BEAM_CENTER_SEARCH:
// How far the ladder reaches, in pixels, in each direction. The step is one pixel, flat:
// a header beam centre is written in pixels and is wrong by pixels. Deriving the step from
// the J0 law was tried and is wrong, because the only cell available at that point is the
// one the FAILED pass returned - on a dataset whose failed cell was a small spurious
// sub-cell the formula asked for a 6 px step, which steps clean over the lobe it is
// looking for.
if (optarg && std::string(optarg) == "off") {
beam_center_search = 0;
break;
}
beam_center_search = optarg
? parse_number_arg<int>(optarg, "--beam-center-search", logger, 1, 1000)
: 12;
break;
case OPT_BEAM_CENTER_CHECK:
beam_center_check = !optarg || std::string(optarg) != "off";
break;
case OPT_ESTIMATE_BEAM_CENTER:
estimate_beam_center = true;
break;
@@ -1902,6 +1928,7 @@ static int RunRugnux(int argc, char **argv) {
config.output_prefix = output_prefix;
config.detect_beam_stop = detect_beam_stop;
config.estimate_beam_center = estimate_beam_center;
config.beam_center_check = beam_center_check;
config.fit_spindle = fit_spindle;
Rugnux process(reader, experiment, *dataset->pixel_mask, config);
@@ -1947,6 +1974,7 @@ static int RunRugnux(int argc, char **argv) {
config.output_prefix = output_prefix;
config.detect_beam_stop = detect_beam_stop;
config.estimate_beam_center = estimate_beam_center;
config.beam_center_check = beam_center_check;
config.fit_spindle = fit_spindle;
config.write_process_h5 = false; // the .poni below is the output of this mode
@@ -2371,6 +2399,8 @@ static int RunRugnux(int argc, char **argv) {
config.two_pass_rotation = two_pass_rotation;
config.detect_beam_stop = detect_beam_stop;
config.estimate_beam_center = estimate_beam_center;
config.beam_center_search_pxl = beam_center_search;
config.beam_center_check = beam_center_check;
config.fit_spindle = fit_spindle;
config.adaptive_integration_radius = adaptive_integration_radius;
config.rotation_postrefine_geometry = rotation_postrefine_geometry;
tests/BeamCenterFromBackgroundTest.cpp
+21
View File
@@ -124,3 +124,24 @@ TEST_CASE("BeamCenterFromBackground_AnAzimuthalShadowIsNotACentreError", "[BeamC
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));
}
// The same, with the detector tilted. Every pixel's 2-theta and azimuth, and the derivative of
// 2-theta with respect to the centre that the fit is built on, go through the detector rotation
// matrix, so the tilt is not a detail of the geometry here - it is in the Jacobian.
TEST_CASE("BeamCenterFromBackground_SurvivesADetectorTilt", "[BeamCenter]") {
DiffractionExperiment x = TestExperiment();
x.PoniRot1_rad(0.005f).PoniRot2_rad(-0.003f);
PixelMask pixel_mask(x);
const DiffractionGeometry geom_true = OffsetBy(x.GetDiffractionGeometry(), 3.0f, -2.5f);
const auto [direct_x, direct_y] = geom_true.GetDirectBeam_pxl();
REQUIRE(std::hypot(direct_x - geom_true.GetBeamX_pxl(), direct_y - geom_true.GetBeamY_pxl()) > 5.0f);
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));
CHECK(estimate->sigma_pxl < 1.0f);
}
tests/BeamCenterFromSpotsTest.cpp
+29
View File
@@ -282,3 +282,32 @@ TEST_CASE("BeamCenterFromSpots_SurvivesASpindleOffPerpendicular", "[BeamCenter]"
CHECK(std::hypot(uncorrected->beam_x_pxl - truth.GetBeamX_pxl(),
uncorrected->beam_y_pxl - truth.GetBeamY_pxl()) > 3.0f * corrected_error);
}
// No detector is mounted exactly square to the beam, and a tilt separates the two centres this
// estimator works with: the PONI it reports and the DIRECT BEAM both mirror lines are taken about.
// The conversion between them is exact - it is where the vote is centred, where each tooth's
// refinement starts, and how the answer is turned back out of the spindle frame - so a sign error
// in it is a sign error of twice the offset in the answer. Here that offset is about 12 px, which
// is 24 times the tolerance below.
TEST_CASE("BeamCenterFromSpots_SurvivesADetectorTilt", "[BeamCenter]") {
DiffractionExperiment x = TestExperiment(Coord(1, 0, 0));
x.PoniRot1_rad(0.005f).PoniRot2_rad(-0.003f);
const DiffractionGeometry truth = OffsetBy(x.GetDiffractionGeometry(), 3.0f, -2.5f);
const Sweep sweep = MakeSweep(x, truth, PAIRS);
// The tilt has to be worth testing: with the two centres on top of each other the conversion
// is the identity and its sign is unobservable.
const auto [direct_x, direct_y] = truth.GetDirectBeam_pxl();
REQUIRE(std::hypot(direct_x - truth.GetBeamX_pxl(), direct_y - truth.GetBeamY_pxl()) > 5.0f);
SpindleEstimate fitted;
const auto estimate = FindBeamCenterFromSpotSymmetry(x, sweep.angle_deg, sweep.spots, &fitted);
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));
CHECK(estimate->sigma_pxl < 1.0f);
// The spindle is perpendicular to the beam and on a lab axis here, so the tilt must not be
// read as an azimuth: what the fit sees of the detector belongs to the detector.
CHECK(fitted.azimuth_rad == Catch::Approx(0.0f).margin(5e-4));
}