From a39e18b80fd0aadff1c5a2ae428126cf26087782 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 22 Oct 2026 19:26:33 +0200 Subject: [PATCH] stills: detector tilt refined in the geometry pre-pass, kept where held-out frames agree The stills geometry bundle (RefineGlobalGeometry) can now refine the detector tilt together with the beam, the distance, the cell and the modules. With the tilt free the beam block is the direct beam, from which the PONI is worked out inside the residual: in PONI terms a tilt moves the whole pattern as a moved beam does (D tan(tilt)/pixel), more than the beam's 8 px box holds, while the direct beam does not move with it and the tilt is left with the keystone. The two tilt-like rows of ModuleGauge are gone: within a module a tilt is a gradient no rigid move makes, and the rows cost held-out accuracy. The shift, turn and outward gauges stay. The tilt is a hypothesis (DecideStillsTilt): the strong frames of the last geometry refinement are split in two halves, each fitted with the tilt held and refined (modules placed in both where the run placed them), each fit scored on the other half, and the tilt is kept where the refined fits predict the held-out spots better. No threshold. It is decided after the stills beam-centre trial and its re-refinement, on frames indexed at the settled centre: a centre error on frames indexed at it reads as a tilt. --rot1/--rot2 or --no-refine-tilt hold the tilt on stills. The stills beam-centre trial now compares direct beams. Its hypotheses (the file's centre, the background walk, the point-symmetry capture) are taken as the point the beam lands on and each trial puts the PONI where its distance and the run's tilt land the direct beam there; the walk, fitted through the geometry's tilt, returns a PONI and is converted. The pre-scan beam-centre check log compares the walk's direct beam, not its PONI, with the file's direct beam. Without a tilt nothing changes. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu --- .../geom_refinement/GeometryRefiner.cpp | 130 ++++++++++++++---- .../geom_refinement/GeometryRefiner.h | 12 +- rugnux/ModuleDiagnostics.cpp | 2 +- rugnux/Rugnux.cpp | 4 + rugnux/Rugnux.h | 18 ++- rugnux/RugnuxPreScan.cpp | 32 ++--- rugnux/RugnuxStills.cpp | 129 +++++++++++++++-- rugnux/rugnux_cli.cpp | 6 +- tests/ModuleRefinementTest.cpp | 50 +++++++ 9 files changed, 319 insertions(+), 64 deletions(-) diff --git a/image_analysis/geom_refinement/GeometryRefiner.cpp b/image_analysis/geom_refinement/GeometryRefiner.cpp index d18ddba7f..bebb5a41f 100644 --- a/image_analysis/geom_refinement/GeometryRefiner.cpp +++ b/image_analysis/geom_refinement/GeometryRefiner.cpp @@ -139,12 +139,56 @@ void PlaceOnModule(const ModuleShift &m, const double move[3], double x, double py = m.y0 + s * ux + c * uy + move[1]; } +// How a residual reads the beam block. Without the tilt refined it is the PONI, as everywhere else. With the +// tilt refined it is the DIRECT BEAM, the point the undeflected beam lands on, and the PONI is worked out from +// it: a tilt moves the whole pattern by D tan(tilt)/pixel, exactly as a moved PONI does, so in PONI terms +// the two trade tens of pixels for each other (and the beam's box cannot hold the trade), while the direct +// beam does not move with the tilt and the tilt is left with only what no beam can imitate - the keystone. +// In XtalResidual's frame the direct beam is where M (u, v) = (-D tan rot1, -D tan rot2 / cos rot1), with +// (u, v) its offset from the PONI in mm and M the image orientation (a mirror and quarter turns, so M^-1 = M^T). +struct BeamBlock { + bool direct = false; + double m00 = 1, m01 = 0, m10 = 0, m11 = 1; // DetectorOrientation::Matrix(), as XtalResidual takes it + double pixel = 1; + + template + void Poni(const T *const beam, const T &distance_mm, const T *const detector_rot, T *poni) const { + if (!direct) { + poni[0] = beam[0]; + poni[1] = beam[1]; + return; + } + using std::cos; + using std::tan; + const T dx = -distance_mm * tan(detector_rot[0]); + const T dy = -distance_mm * tan(detector_rot[1]) / cos(detector_rot[0]); + poni[0] = beam[0] - (T(m00) * dx + T(m10) * dy) / T(pixel); + poni[1] = beam[1] - (T(m01) * dx + T(m11) * dy) / T(pixel); + } +}; + +// XtalResidual with the beam block read through BeamBlock. +struct BeamXtalResidual { + XtalResidual residual; + BeamBlock beam_block; + + template + bool operator()(const T *const beam, const T *const distance_mm, const T *const detector_rot, + const T *const rotation_axis, const T *const p0, const T *const p1, const T *const p2, + T *out) const { + T poni[2]; + beam_block.Poni(beam, distance_mm[0], detector_rot, poni); + return residual(poni, distance_mm, detector_rot, rotation_axis, p0, p1, p2, out); + } +}; + // XtalResidual for a spot on a module that is being placed: block `move` = (dx, dy, rot) of the module, // (ux, uy) the spot's offset from the module centre. Moving the spot by R(rot) u - u + (dx, dy) shifts its // offset from the beam exactly as moving the beam the other way does, so the beam handed to XtalResidual // carries the move and everything else is the residual of an unmoved spot. struct ModuleXtalResidual { XtalResidual residual; + BeamBlock beam_block; double ux, uy; template @@ -153,9 +197,11 @@ struct ModuleXtalResidual { const T *const move, T *out) const { using std::cos; using std::sin; + T poni[2]; + beam_block.Poni(beam, distance_mm[0], detector_rot, poni); const T c = cos(move[2]), s = sin(move[2]); - const T moved_beam[2] = {beam[0] - ((c - T(1)) * T(ux) - s * T(uy) + move[0]), - beam[1] - (s * T(ux) + (c - T(1)) * T(uy) + move[1])}; + const T moved_beam[2] = {poni[0] - ((c - T(1)) * T(ux) - s * T(uy) + move[0]), + poni[1] - (s * T(ux) + (c - T(1)) * T(uy) + move[1])}; return residual(moved_beam, distance_mm, detector_rot, rotation_axis, p0, p1, p2, out); } }; @@ -167,36 +213,30 @@ struct ModuleXtalResidual { // weighted by the spots each module holds: a module the frames barely reach is placed poorly, and an // unweighted sum would hand its error to all the others. // -// With the tilt refined, the two tilt-like patterns are held at zero as well: a detector tilted by a small -// angle about y moves a point (x, y) (from the beam) by (x^2, x y) / D times the angle, one about x by -// (x y, y^2) / D. Each is taken as the modules' moves projected on that field, in units of the lever. +// The detector tilt is not among them. Within one module a tilt is a gradient (a scale and a shear across +// it), which no rigid move of the module makes, and with the tilt refined the placement gains no new +// direction it cannot tell apart; holding the modules' moves to no tilt-like pattern as well only cost +// accuracy on frames the fit had not seen. struct ModuleGauge { std::vector> outward; // unit vector from the beam to each module centre - std::vector> centre; // each module centre from the beam, over lever_px std::vector share; // each module's fraction of the spots double lever_px; // turns are weighed at this distance from the centre double weight; - bool tilt = false; template bool operator()(T const *const *move, T *out) const { - T sx(0), sy(0), sr(0), so(0), t1(0), t2(0); + T sx(0), sy(0), sr(0), so(0); for (size_t i = 0; i < outward.size(); ++i) { const T w(share[i]); sx += w * move[i][0]; sy += w * move[i][1]; sr += w * move[i][2] * T(lever_px); so += w * (move[i][0] * T(outward[i][0]) + move[i][1] * T(outward[i][1])); - const double x = centre[i][0], y = centre[i][1]; - t1 += w * (move[i][0] * T(x * x) + move[i][1] * T(x * y)); - t2 += w * (move[i][0] * T(x * y) + move[i][1] * T(y * y)); } out[0] = T(weight) * sx; out[1] = T(weight) * sy; out[2] = T(weight) * sr; out[3] = T(weight) * so; - out[4] = tilt ? T(weight) * t1 : T(0); - out[5] = tilt ? T(weight) * t2 : T(0); return true; } }; @@ -211,7 +251,8 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge try { const gemmi::CrystalSystem system = MapSystem(settings.crystal_system); - // Shared geometry parameter blocks (refined). Detector tilt and the rotation axis are constant. + // Shared geometry parameter blocks (refined). The rotation axis is constant, and so is the detector tilt + // unless settings.refine_tilt. double beam[2] = {nominal_geom.GetBeamX_pxl(), nominal_geom.GetBeamY_pxl()}; double distance_mm = nominal_geom.GetDetectorDistance_mm(); double detector_rot[2] = {nominal_geom.GetPoniRot1_rad(), nominal_geom.GetPoniRot2_rad()}; @@ -223,6 +264,28 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge const double lambda = nominal_geom.GetWavelength_A(); const double pixel = nominal_geom.GetPixelSize_mm(); + // With the tilt refined, the beam block is the direct beam (BeamBlock), started and bounded where the + // nominal geometry puts the direct beam. + const bool refine_tilt = settings.refine_tilt && settings.refine_beam_distance; + BeamBlock beam_block; + if (refine_tilt) { + const RotMatrix m = orientation.Matrix(); + beam_block = {true, m.Column(0).x, m.Column(1).x, m.Column(0).y, m.Column(1).y, pixel}; + double poni_at_zero[2]; + const double zero[2] = {0, 0}; + beam_block.Poni(zero, distance_mm, detector_rot, poni_at_zero); + beam[0] -= poni_at_zero[0]; + beam[1] -= poni_at_zero[1]; + } + const double beam0[2] = {beam[0], beam[1]}; + // The PONI the beam block stands for, at the current tilt and distance. + auto poni = [&](double &x, double &y) { + double p[2]; + beam_block.Poni(beam, distance_mm, detector_rot, p); + x = p[0]; + y = p[1]; + }; + // Shared cell blocks. Lengths refined (regularized); angles constant. double cell_len[3] = {input_cell.a, input_cell.b, input_cell.c}; const double cell_len0[3] = {input_cell.a, input_cell.b, input_cell.c}; @@ -304,8 +367,10 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge for (int round = 0; round < settings.rounds; ++round) { EffectiveCell(system, cell_len, cell_ang, eff_len, eff_alpha, eff_beta, eff_gamma); + double poni_x, poni_y; + poni(poni_x, poni_y); DiffractionGeometry cur = nominal_geom; - cur.BeamX_pxl(beam[0]).BeamY_pxl(beam[1]).DetectorDistance_mm(distance_mm) + cur.BeamX_pxl(poni_x).BeamY_pxl(poni_y).DetectorDistance_mm(distance_mm) .PoniRot1_rad(static_cast(detector_rot[0])).PoniRot2_rad(static_cast(detector_rot[1])); const double tol = round_tol(round); @@ -337,11 +402,18 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge new ceres::AutoDiffCostFunction( new ModuleXtalResidual{XtalResidual(s.x, s.y, lambda, pixel, cos_rot3, sin_rot3, 0.0, h, k, l, system, orientation), - s.x - m.x0, s.y - m.y0}), + beam_block, s.x - m.x0, s.y - m.y0}), new ceres::CauchyLoss(loss_scale), beam, &distance_mm, detector_rot, rot_vec, orient[i].data(), cell_len, cell_ang, move[s.module].data()); - } else + } else if (refine_tilt) + problem.AddResidualBlock( + new ceres::AutoDiffCostFunction( + new BeamXtalResidual{XtalResidual(x, y, lambda, pixel, cos_rot3, sin_rot3, 0.0, + h, k, l, system, orientation), beam_block}), + new ceres::CauchyLoss(loss_scale), + beam, &distance_mm, detector_rot, rot_vec, orient[i].data(), cell_len, cell_ang); + else problem.AddResidualBlock( new ceres::AutoDiffCostFunction( new XtalResidual(x, y, lambda, pixel, cos_rot3, sin_rot3, 0.0, @@ -396,22 +468,18 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge const double oy = modules[m].y0 - nominal_geom.GetBeamY_pxl(); const double r = std::max(1.0, std::hypot(ox, oy)); gauge->outward.push_back({ox / r, oy / r}); - gauge->centre.push_back({ox, oy}); gauge->share.push_back(module_spots[m] / total); lever += r; } gauge->lever_px = lever / static_cast(moved.size()); - for (auto &c : gauge->centre) - c = {c[0] / gauge->lever_px, c[1] / gauge->lever_px}; gauge->weight = 10.0 * recip_per_px * sqrtN; - gauge->tilt = settings.refine_tilt; auto *cost = new ceres::DynamicAutoDiffCostFunction(gauge); std::vector blocks; for (size_t m : moved) { cost->AddParameterBlock(3); blocks.push_back(move[m].data()); } - cost->SetNumResiduals(6); + cost->SetNumResiduals(4); problem.AddResidualBlock(cost, nullptr, blocks); } for (size_t m : moved) { @@ -426,10 +494,10 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge // Beam + distance bounds around the nominal value; cell lengths loosely bounded (the // regularizer does the anchoring). - problem.SetParameterLowerBound(beam, 0, nominal_geom.GetBeamX_pxl() - settings.beam_range_px); - problem.SetParameterUpperBound(beam, 0, nominal_geom.GetBeamX_pxl() + settings.beam_range_px); - problem.SetParameterLowerBound(beam, 1, nominal_geom.GetBeamY_pxl() - settings.beam_range_px); - problem.SetParameterUpperBound(beam, 1, nominal_geom.GetBeamY_pxl() + settings.beam_range_px); + problem.SetParameterLowerBound(beam, 0, beam0[0] - settings.beam_range_px); + problem.SetParameterUpperBound(beam, 0, beam0[0] + settings.beam_range_px); + problem.SetParameterLowerBound(beam, 1, beam0[1] - settings.beam_range_px); + problem.SetParameterUpperBound(beam, 1, beam0[1] + settings.beam_range_px); problem.SetParameterLowerBound(&distance_mm, 0, std::max(1.0, nominal_geom.GetDetectorDistance_mm() - settings.distance_range_mm)); problem.SetParameterUpperBound(&distance_mm, 0, @@ -441,7 +509,7 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge if (!settings.refine_distance) problem.SetParameterBlockConstant(&distance_mm); - if (settings.refine_tilt && settings.refine_beam_distance) { + if (refine_tilt) { const double tilt_range = 2.0 * PI / 180.0; for (int j = 0; j < 2; ++j) { problem.SetParameterLowerBound(detector_rot, j, (j == 0 ? nominal_geom.GetPoniRot1_rad() @@ -478,8 +546,10 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge // Final fit quality: median detector residual of the retained spots at the refined geometry. EffectiveCell(system, cell_len, cell_ang, eff_len, eff_alpha, eff_beta, eff_gamma); + double poni_x, poni_y; + poni(poni_x, poni_y); DiffractionGeometry cur = nominal_geom; - cur.BeamX_pxl(beam[0]).BeamY_pxl(beam[1]).DetectorDistance_mm(distance_mm) + cur.BeamX_pxl(poni_x).BeamY_pxl(poni_y).DetectorDistance_mm(distance_mm) .PoniRot1_rad(static_cast(detector_rot[0])).PoniRot2_rad(static_cast(detector_rot[1])); std::vector residuals; constexpr double CAP_PX = 2.0; @@ -531,8 +601,8 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge } result.ok = solve_usable && result.spots_used >= settings.min_spots_total && result.frames_used > 0; - result.beam_x_px = beam[0]; - result.beam_y_px = beam[1]; + result.beam_x_px = poni_x; + result.beam_y_px = poni_y; result.distance_mm = distance_mm; result.cell = BuildUnitCell(system, cell_len, cell_ang); result.rot1_rad = detector_rot[0]; diff --git a/image_analysis/geom_refinement/GeometryRefiner.h b/image_analysis/geom_refinement/GeometryRefiner.h index f624b4a2b..ac8ebe192 100644 --- a/image_analysis/geom_refinement/GeometryRefiner.h +++ b/image_analysis/geom_refinement/GeometryRefiner.h @@ -19,7 +19,8 @@ // re-indexing with it converts more frames. Mirrors the validated Python prototype: one Ceres problem // with shared beam/distance/cell blocks and a per-frame orientation block, robust loss, and a cell // regularizer that anchors the (known) cell to break the distance<->cell-scale degeneracy at low -// resolution. Detector tilt (rot1/rot2) is NOT refined (gauge-coupled with the beam; zero indexing gain). +// resolution. The detector tilt (rot1/rot2) is refined only on request (refine_tilt), and then with the beam +// written as the direct beam, which the tilt does not move. struct GeomRefineSpot { float x = 0, y = 0; // observed (raw) spot position, pixels @@ -59,15 +60,16 @@ struct GeometryRefinerSettings { std::vector modules; bool refine_modules = false; int min_module_spots = 500; - // Refine the detector tilt (PONI rot1, rot2) with the beam and distance. A tilt moves the modules far - // from the beam more than those near it, in a pattern module moves could otherwise take up: with it - // refined, the modules are also held to no tilt-like pattern of their own (ModuleGauge). + // Refine the detector tilt (PONI rot1, rot2), within 2 degrees of the nominal one, together with the beam, + // the distance and (with refine_modules) the modules. The beam is then fitted as the direct beam, and its + // box (beam_range_px) is about where the nominal geometry puts the direct beam: in PONI terms a tilt moves + // the whole pattern as a moved beam does, by D tan(tilt) / pixel, more than the box would hold. bool refine_tilt = false; }; struct GeometryRefinerResult { bool ok = false; - double beam_x_px = 0, beam_y_px = 0; // refined beam centre (== direct beam, tilt not refined) + double beam_x_px = 0, beam_y_px = 0; // refined beam centre, the PONI (DiffractionGeometry::BeamX_pxl) double distance_mm = 0; UnitCell cell{}; // refined cell int frames_used = 0; diff --git a/rugnux/ModuleDiagnostics.cpp b/rugnux/ModuleDiagnostics.cpp index 97fb691d6..f261e742d 100644 --- a/rugnux/ModuleDiagnostics.cpp +++ b/rugnux/ModuleDiagnostics.cpp @@ -253,7 +253,7 @@ void RunModuleDiagnostics(const ModuleDiagnosticsInput &in) { } } - // The tilt refined, and the modules held to no tilt-like pattern. + // The tilt refined as well. { GeometryRefinerSettings s2 = st; s2.refine_tilt = true; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 68b55bd77..d557b1749 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -893,6 +893,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b if (full && !experiment_.IsRotationIndexing() && config_.beam_center_check && !config_.live && DecideStillsGeometry(start_image, images_to_process) && config_.refine_geometry.has_value()) RefineStillsGeometry(start_image, end_image, images_to_process, observer); + // The detector tilt, now that the centre is settled: a centre error on frames indexed at it is exactly a + // tilt's pattern shift, so a tilt fitted before the centre is decided takes up the centre's error. + if (full) + DecideStillsTilt(); // The geometry pre-pass integrates at the radius the run started with; the width the pre-scan just // measured is applied to the canonical pass instead (RunAllPasses, between the two). See diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index bdb296d91..c54b2bbf2 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -40,6 +40,7 @@ #include "../image_analysis/scale_merge/CrystalSetting.h" // MetricViolation, MAX_METRIC_VIOLATION #include "../image_analysis/geom_refinement/PostRefine.h" // PostRefineResult #include "../image_analysis/geom_refinement/BeamCenterFromBackground.h" // BeamCenterEstimate +#include "../image_analysis/geom_refinement/GeometryRefiner.h" // GeomRefineFrame #include "../image_analysis/rotation_indexer/RotationIndexer.h" // RotationIndexerResult #include "../image_analysis/geom_refinement/PowderCalibration.h" // CalibrationMethod, CalibrationResult #include "../image_analysis/WriteReflections.h" // OutputProvenance @@ -118,7 +119,7 @@ struct ProcessConfig { std::vector calibrant_ring_q; // Fit the detector tilt (PONI rot1/rot2) along with the beam centre and the distance. False holds // the tilt where the header put it, for a calibration destined for a program that cannot express - // one - XDS. + // one - XDS. On stills it also holds the tilt (DecideStillsTilt). bool calibration_refine_tilt = true; // Measure the spot budget instead of taking experiment.GetMaxSpotCount() as given (FullAnalysis, @@ -180,6 +181,9 @@ struct ProcessConfig { // The detector distance was given on the command line: on stills it is held - not searched // (DecideStillsGeometry), not refined (RefineStillsGeometry). bool detector_distance_given = false; + // The detector tilt was given on the command line (--rot1/--rot2): on stills it is held, not refined + // (DecideStillsTilt). + bool detector_tilt_given = false; // Neither the input nor the command line states the distance (an XTC): the input's value is a // placeholder, which the stills search does not hold against the others and nothing refines from. bool detector_distance_unknown = false; @@ -1418,6 +1422,18 @@ class Rugnux { const std::vector &frames); // Returns whether the run's geometry was moved off the file's. bool DecideStillsGeometry(int start_image, int images_to_process); + // Stills: the detector tilt as a hypothesis, refined or held, kept where the frames it was not fitted to + // say so. Decided once the beam centre is (RunPipeline, after DecideStillsGeometry), on the frames the last + // RefineStillsGeometry fitted: on frames indexed at a wrong centre, the centre's error reads as a tilt. + struct StillsTiltBundle { + std::vector frames; + // Where the modules were placed on these frames, from where they were given (`given`); empty: the + // modules stay as the frames have them. + std::vector modules; + DetectorModules given; + }; + std::optional stills_tilt_bundle_; + void DecideStillsTilt(); // Stills: the cell the user gave (-C / reference), kept apart from the one the geometry refinement // writes into experiment_. See RefineStillsGeometry. std::optional stills_user_cell_; diff --git a/rugnux/RugnuxPreScan.cpp b/rugnux/RugnuxPreScan.cpp index 6786be5be..6247f61bb 100644 --- a/rugnux/RugnuxPreScan.cpp +++ b/rugnux/RugnuxPreScan.cpp @@ -1116,15 +1116,16 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru FindBeamCenter(x, mask, *projection, nthreads, &capture); const auto capture_ms = std::chrono::duration_cast( std::chrono::steady_clock::now() - t0).count(); + const auto [direct_x, direct_y] = x.GetDiffractionGeometry().GetDirectBeam_pxl(); 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; margin {:.1f}% " + "detector is ({:.1f},{:.1f}), {:.1f} px from the file's direct beam; margin {:.1f}% " "(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 - x.GetBeamX_pxl(), - capture.point[0].beam_y_pxl - x.GetBeamY_pxl()), + std::hypot(capture.point[0].beam_x_pxl - direct_x, + capture.point[0].beam_y_pxl - direct_y), 100.0f * capture.margin_point, capture_ms); const float need = BeamCenterNeed_pxl(x); if (!background_center) { @@ -1132,17 +1133,17 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru "this run - the file's ({:.2f},{:.2f}) is neither confirmed nor contradicted", x.GetBeamX_pxl(), x.GetBeamY_pxl()); } else { - // Against the DIRECT BEAM, not the PONI. Both the FFT capture and the walk read the - // centrosymmetry of the scattered background about the point the beam lands on; the - // geometry's beam_x_pxl is the PONI, which parts from it by distance*tan(rot)/pixel as - // soon as the detector is tilted. Measured on a calibrant at 110/150/200/300 mm the - // capture sits 0.5-1.5 px from the fitted direct beam and 5.0/6.8/8.4/17.9 px from the - // fitted PONI, so comparing it with the PONI reports a disagreement the data do not - // have - on the in-house files, whose header carries a 0.234 deg tilt, 3.3 px of it at - // 120 mm against a stated need of 0.7 px. - const auto [direct_x, direct_y] = x.GetDiffractionGeometry().GetDirectBeam_pxl(); - const float dx = background_center->beam_x_pxl - direct_x; - const float dy = background_center->beam_y_pxl - direct_y; + // DIRECT BEAM against direct beam. The capture reads the centrosymmetry of the image + // about the point the beam lands on; the geometry's beam_x_pxl is the PONI, which parts + // from it by distance*tan(rot)/pixel as soon as the detector is tilted. The walk fits the + // background's isotropy in 2theta through the geometry's own tilt, so what it returns is a + // PONI, and its direct beam is as far from it as the file's direct beam is from the file's + // PONI. (A run given a 1 deg tilt put the walk 12 px from the direct beam it implies, and + // that direct beam 1 px from where the run without the tilt measured it.) + const float walk_x = background_center->beam_x_pxl + (direct_x - x.GetBeamX_pxl()); + const float walk_y = background_center->beam_y_pxl + (direct_y - x.GetBeamY_pxl()); + const float dx = walk_x - direct_x; + const float dy = walk_y - direct_y; const float moved = std::hypot(dx, dy); const auto across = AcrossSpindle_pxl(x, dx, dy); // Everything a reader needs to judge the header, on one line: what the file says, what the @@ -1154,8 +1155,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru "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", - direct_x, direct_y, - background_center->beam_x_pxl, background_center->beam_y_pxl, + direct_x, direct_y, walk_x, walk_y, 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 diff --git a/rugnux/RugnuxStills.cpp b/rugnux/RugnuxStills.cpp index 9ef1e3108..f2c7b16a5 100644 --- a/rugnux/RugnuxStills.cpp +++ b/rugnux/RugnuxStills.cpp @@ -167,6 +167,7 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ return; } + stills_tilt_bundle_.reset(); const int refine_frames = std::max(1, config_.refine_geometry.value()); if (observer) observer->OnPhase("Geometry refinement (first pass)"); @@ -359,6 +360,7 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ constexpr double REFIT_ABOVE = 0.003, REFIT_CONVERGED = 0.0005; constexpr int MAX_REFITS = 3; double moved_by = std::abs(r.distance_mm / nominal.GetDetectorDistance_mm() - 1.0); + std::vector refitted; // the bundle as last indexed again, if it was for (int refit = 0; refit < MAX_REFITS && !placing_modules && moved_by > REFIT_CONVERGED && (refit > 0 || moved_by > REFIT_ABOVE) && !cancelled_; refit++) { TimingMark refit_mark(logger, "geometry refit"); @@ -407,7 +409,12 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ logger.Info("Geometry refinement: indexed again at the fitted geometry ({} frames), distance {:.4f} -> {:.4f} mm", again.size(), r.distance_mm, r2.distance_mm); r = r2; + refitted = std::move(again); } + // The frames the tilt is decided on (DecideStillsTilt): the bundle as last fitted, or, where the modules are + // placed below, the frames they are placed on, with the placement to start from. + if (config_.calibration_refine_tilt && !config_.detector_tilt_given) + stills_tilt_bundle_ = StillsTiltBundle{refitted.empty() ? strong : refitted, {}, experiment_.GetModules()}; // A detector of several modules: each placed as a rigid block, where that predicts the spots better on // frames the placement was not fitted to. The strong frames are dealt alternately into two halves; each @@ -512,6 +519,8 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ } module_shifts_ = rm.modules; if (config_.apply_modules) { + if (stills_tilt_bundle_) + stills_tilt_bundle_ = StillsTiltBundle{placing, rm.modules, given}; DetectorModules placed = given; placed.Apply(rm.modules); experiment_.Modules(placed); @@ -536,6 +545,93 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_ experiment_.SetUnitCell(r.cell); } +// The detector tilt of a stills run, as a hypothesis against the tilt the run has (the file's, normally none): +// the strong frames of the geometry refinement are dealt alternately into two halves, each half is fitted with +// the tilt held and with it refined - the beam, the distance, the cell and the modules (where they were placed) +// refined in both - and each fit is scored on the other half with its geometry frozen. The tilt is kept if the +// refined fits predict the held-out spots better. Where the tilt is real it bends the pattern in a way nothing +// else in the fit can (the keystone); where it is not, the gain on frames the fit has not seen is nothing and the +// tilt found is too small to matter either way, so no threshold has to say how large a tilt must be. +void Rugnux::DecideStillsTilt() { + if (!stills_tilt_bundle_) + return; + Logger logger("Rugnux"); + TimingMark mark(logger, "stills tilt decision"); + const StillsTiltBundle bundle = std::move(*stills_tilt_bundle_); + stills_tilt_bundle_.reset(); + + const DiffractionGeometry start = experiment_.GetDiffractionGeometry(); + const UnitCell cell = experiment_.GetUnitCell().value(); + GeometryRefinerSettings st; + st.crystal_system = experiment_.GetCrystalSystem(); + st.num_threads = config_.nthreads; + st.refine_distance = !config_.detector_distance_given; + st.cell_reg_coeff = 0.0; + st.modules = bundle.modules; + st.refine_modules = !bundle.modules.empty(); + + std::vector half[2]; + for (size_t i = 0; i < bundle.frames.size(); i++) + half[i % 2].push_back(bundle.frames[i]); + + // Held-out residual of the geometry fitted on half `fit`, with the tilt refined or held, scored on every + // spot of the other half. + auto held_out = [&](bool tilt, int fit) { + GeometryRefinerSettings sf = st; + sf.refine_tilt = tilt; + const auto f = RefineGlobalGeometry(start, cell, half[fit], sf); + if (!f.ok) + return std::numeric_limits::infinity(); + DiffractionGeometry g = start; + g.BeamX_pxl(f.beam_x_px).BeamY_pxl(f.beam_y_px).DetectorDistance_mm(f.distance_mm) + .PoniRot1_rad(f.rot1_rad).PoniRot2_rad(f.rot2_rad); + GeometryRefinerSettings ev = st; + ev.refine_beam_distance = false; + ev.modules = f.modules; + ev.refine_modules = false; + const auto e = RefineGlobalGeometry(g, f.cell, half[1 - fit], ev); + if (!e.ok) + return std::numeric_limits::infinity(); + return e.capped_rms_px; + }; + // The four held-out fits and the tilted fit on every frame side by side, as for the module placement. + std::future held_fits[2], tilt_fits[2]; + for (int fit = 0; fit < 2; fit++) { + held_fits[fit] = std::async(std::launch::async, held_out, false, fit); + tilt_fits[fit] = std::async(std::launch::async, held_out, true, fit); + } + auto tilted = std::async(std::launch::async, [&] { + GeometryRefinerSettings sf = st; + sf.refine_tilt = true; + return RefineGlobalGeometry(start, cell, bundle.frames, sf); + }); + const double held = (held_fits[0].get() + held_fits[1].get()) / 2.0; + const double refined = (tilt_fits[0].get() + tilt_fits[1].get()) / 2.0; + const GeometryRefinerResult r = tilted.get(); + const bool keep = r.ok && refined < held; + logger.Info("Detector tilt: rot1 {:+.3f} deg, rot2 {:+.3f} deg refined ({:+.3f}, {:+.3f} held); held-out spot " + "residual (capped rms) {:.3f} px held, {:.3f} px refined, {} frames - the tilt is {}", + r.rot1_rad * 180.0 / PI, r.rot2_rad * 180.0 / PI, start.GetPoniRot1_rad() * 180.0 / PI, + start.GetPoniRot2_rad() * 180.0 / PI, held, refined, bundle.frames.size(), + keep ? "refined" : "held"); + if (!keep) + return; + + experiment_.BeamX_pxl(r.beam_x_px).BeamY_pxl(r.beam_y_px).DetectorDistance_mm(r.distance_mm) + .PoniRot1_rad(r.rot1_rad).PoniRot2_rad(r.rot2_rad); + experiment_.SetUnitCell(r.cell); + if (!bundle.modules.empty()) { + DetectorModules placed = bundle.given; + placed.Apply(r.modules); + experiment_.Modules(placed); + module_shifts_ = r.modules; + } + const auto [direct_x, direct_y] = experiment_.GetDiffractionGeometry().GetDirectBeam_pxl(); + logger.Info("Detector tilt: direct beam ({:.2f}, {:.2f}) px, PONI ({:.2f}, {:.2f}) px, distance {:.4f} mm, " + "cell a/b/c {:.3f}/{:.3f}/{:.3f} A", direct_x, direct_y, r.beam_x_px, r.beam_y_px, r.distance_mm, + r.cell.a, r.cell.b, r.cell.c); +} + std::vector Rugnux::FindStillsSampleSpots(const std::vector &ordinals, int start_image) { const auto dataset = reader_.GetDataset(); AzimuthalIntegrationMapping mapping(experiment_, pixel_mask_); @@ -659,8 +755,24 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { // Where nothing indexes at any distance, the sample is this many times larger (see below). constexpr int WIDER_SAMPLE = 10; - const float file_x = experiment_.GetBeamX_pxl(), file_y = experiment_.GetBeamY_pxl(); + // Every centre below is a DIRECT BEAM, the point the undeflected beam lands on, and so is what the capture + // measures; the geometry's beam_x_pxl is the PONI, which parts from it by distance * tan(tilt) / pixel + // once the detector is tilted. A trial puts the PONI where its distance and the run's tilt land the direct + // beam on the centre tried (place). Without a tilt the two are the same point. + const auto place = [](DiffractionExperiment &e, float x, float y, float distance) { + e.DetectorDistance_mm(distance); + const auto [direct_x, direct_y] = e.GetDiffractionGeometry().GetDirectBeam_pxl(); + e.BeamX_pxl(x + (e.GetBeamX_pxl() - direct_x)).BeamY_pxl(y + (e.GetBeamY_pxl() - direct_y)); + }; + const auto file_direct = experiment_.GetDiffractionGeometry().GetDirectBeam_pxl(); + const float file_x = file_direct.first, file_y = file_direct.second; const float file_distance = experiment_.GetDetectorDistance_mm(); + // The background walk is fitted through the run's tilt, so it returns a PONI: its direct beam is as far + // from it as the file's direct beam is from the file's PONI. + std::optional> walk; + if (background_center_) + walk = std::pair{background_center_->beam_x_pxl + (file_x - experiment_.GetBeamX_pxl()), + background_center_->beam_y_pxl + (file_y - experiment_.GetBeamY_pxl())}; std::vector sample = select_equally_spaced_image_ordinals( images_to_process, std::min(images_to_process, SAMPLE_FRAMES)); int scan_count = SCAN_FRAMES; @@ -669,8 +781,7 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { // centre where there is one (it does not depend on the distance), strongest frames first. const auto find_sample = [&]() { frames = FindStillsSampleSpots(sample, start_image); - const float bx = background_center_ ? background_center_->beam_x_pxl : experiment_.GetBeamX_pxl(); - const float by = background_center_ ? background_center_->beam_y_pxl : experiment_.GetBeamY_pxl(); + const auto [bx, by] = walk ? *walk : experiment_.GetDiffractionGeometry().GetDirectBeam_pxl(); const size_t dropped = DropPowderRingSpots(frames, bx, by, experiment_.GetPixelSize_mm()); if (dropped > 0) logger.Info("Stills geometry trial: {} sampled spots on powder rings left out", dropped); @@ -691,7 +802,7 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { if (const auto it = tried.find(key); it != tried.end()) return it->second; DiffractionExperiment e = experiment_; - e.BeamX_pxl(x).BeamY_pxl(y).DetectorDistance_mm(distance); + place(e, x, y, distance); if (stills_user_cell_) e.SetUnitCell(*stills_user_cell_); // the grid asks the user's cell, not the one refined at the file's distance // The distance scan's frames are indexed with a wider acceptance: the scan only has to find the peak, and @@ -715,8 +826,8 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { // same background, which do not fail alike (a shadow across the background biases the isotropy fit, // a flat background leaves the capture without a peak). std::vector> candidates; - if (background_center_) - candidates.emplace_back(background_center_->beam_x_pxl, background_center_->beam_y_pxl); + if (walk) + candidates.push_back(*walk); if (point_symmetry_center_) candidates.push_back(*point_symmetry_center_); float x = file_x, y = file_y; @@ -851,7 +962,7 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { // asked again on them (a sample found eight times too far out put a seventh as many spots on the // lattice at the right distance as one found nearer). if (std::abs(distance / file_distance - 1.0f) > RESAMPLE_ABOVE) { - experiment_.BeamX_pxl(x).BeamY_pxl(y).DetectorDistance_mm(distance); + place(experiment_, x, y, distance); find_sample(); tried.clear(); best = -1; @@ -894,9 +1005,9 @@ bool Rugnux::DecideStillsGeometry(int start_image, int images_to_process) { } } - logger.Info("Stills beam centre: the run uses ({:.2f},{:.2f}){}", x, y, + logger.Info("Stills beam centre: the run uses the direct beam at ({:.2f},{:.2f}){}", x, y, x == file_x && y == file_y ? ", the file's" : ""); - experiment_.BeamX_pxl(x).BeamY_pxl(y).DetectorDistance_mm(distance); + place(experiment_, x, y, distance); return x != file_x || y != file_y || distance != file_distance; } diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index e628abd00..b6d2a0979 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -157,7 +157,7 @@ void print_usage() { std::cout << " Calibration (--mode calibration)" << std::endl; std::cout << " --calibrant Powder standard: lab6|agbh|ceo2|si|ice (default: lab6, case-insensitive). Overridden by -C: a unit cell given there IS the standard, with its rings enumerated from the cell and its absences from -S where that is given, or taken as primitive where it is not. -S supplies SYMMETRY absences only: a standard whose extinctions come from its atomic basis rather than from a symmetry element - silicon, whose 222 is absent in Fd-3m for that reason - is better named with --calibrant" << std::endl; std::cout << " --calibration How the rings are measured: rings|spots (default: rings). rings sums the (q x azimuth) azimuthal profile over every processed image and fits the ring arcs in it; spots pools the found spots and fits those. -s/-e/--stride select the images; rings defaults --azim-phi-bins to 32" << std::endl; - std::cout << " --no-refine-tilt Do not refine the detector tilt: hold rot1/rot2 at the header value and fit only the beam centre and the distance. For a calibration handed to a program that cannot express a tilted detector, e.g. XDS" << std::endl; + std::cout << " --no-refine-tilt Do not refine the detector tilt: hold rot1/rot2 at the header value and fit only the beam centre and the distance. For a calibration handed to a program that cannot express a tilted detector, e.g. XDS. On stills it holds the tilt in the geometry refinement too" << std::endl; std::cout << std::endl; std::cout << " Detector geometry" << std::endl; @@ -277,7 +277,7 @@ void print_usage() { std::cout << " --pixel-size Pixel size (um)" << std::endl; std::cout << " --gain Detector units per photon, for a file that states none (SACLA, CXI, XTC and --crystfel-geometry input; the frames are converted to photons with it). Default: for SACLA and CXI input, measured on the background of a few frames; XTC frames stay in detector units" << std::endl; std::cout << " --saturation Pixel value from which a pixel is overloaded, in the units the frames are read in (photons where they are converted). Replaces the file's limit and the one stills read off the frames" << std::endl; - std::cout << " --rot1 PONI rotation 1 (rad)" << std::endl; + std::cout << " --rot1 PONI rotation 1 (rad). On stills, --rot1 or --rot2 holds the tilt: it is not refined. Without them stills refine the tilt where frames the fit has not seen agree" << std::endl; std::cout << " --rot2 PONI rotation 2 (rad)" << std::endl; std::cout << " --rot3 PONI rotation 3, about the beam (rad)" << std::endl; std::cout << " --detector-mirror-y Stored image is mirrored in Y vs the detector" << std::endl; @@ -3477,6 +3477,8 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { config.beam_center_check = beam_center_check; config.detector_distance_given = detector_distance_mm.has_value(); config.detector_distance_unknown = detector_distance_unknown; + config.detector_tilt_given = rot1_rad.has_value() || rot2_rad.has_value(); + config.calibration_refine_tilt = calibration_refine_tilt; config.detect_saturation = !saturation_value && (input_is_cbf || input_is_marccd || input_is_smv || input_is_crystfel || input_is_panel_table || input_is_cxi || input_is_xtc); diff --git a/tests/ModuleRefinementTest.cpp b/tests/ModuleRefinementTest.cpp index ede657447..08494e13c 100644 --- a/tests/ModuleRefinementTest.cpp +++ b/tests/ModuleRefinementTest.cpp @@ -149,3 +149,53 @@ TEST_CASE("ScanModuleShifts_FindsShiftedModule") { CHECK(many[m].on_lattice_at_zero == one[m].on_lattice_at_zero); } } + +// A tilted detector, one panel: started from the untilted geometry with the beam where the beam really lands +// (the direct beam, which is what a beam-centre measurement gives), the refinement with the tilt must find it +// and put the PONI where the tilt has it - 10 px from the direct beam here, further than the beam's 8 px box +// would reach in PONI terms - and the fit with the tilt held must be worse. With the image mirrored and turned +// as well, the direct beam has to be taken through the image orientation. +TEST_CASE("RefineGlobalGeometry_Tilt") { + for (const auto &orientation : {DetectorOrientation(), DetectorOrientation(true, 1)}) { + DiffractionExperiment exp; + exp.IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(150) + .PoniRot1_rad(0.3 * PI / 180.0).PoniRot2_rad(-0.2 * PI / 180.0); + exp.Detector().ImageOrientation(orientation); + const DiffractionGeometry truth = exp.GetDiffractionGeometry(); + const auto [direct_x, direct_y] = truth.GetDirectBeam_pxl(); + REQUIRE(std::hypot(direct_x - 1000, direct_y - 1000) > 10); + + // One module covering the whole image: the frames' spots, each where the tilted detector shows it. + ModuleShift whole; + whole.x0 = 1000; + whole.y0 = 1000; + const CrystalLattice latt(50, 60, 70, 90, 90, 90); + std::vector frames = FramesOnMovedModules(exp, latt, {whole, whole, whole, whole}); + REQUIRE(frames.size() > 30); + for (auto &f: frames) + for (auto &s: f.spots) + s.module = -1; + + DiffractionGeometry nominal = truth; + nominal.PoniRot1_rad(0).PoniRot2_rad(0).BeamX_pxl(direct_x).BeamY_pxl(direct_y); + GeometryRefinerSettings st; + st.crystal_system = gemmi::CrystalSystem::Orthorhombic; + st.num_threads = 2; + + const auto held = RefineGlobalGeometry(nominal, latt.GetUnitCell(), frames, st); + REQUIRE(held.ok); + CHECK(held.rot1_rad == 0); + CHECK(held.rot2_rad == 0); + + st.refine_tilt = true; + const auto tilted = RefineGlobalGeometry(nominal, latt.GetUnitCell(), frames, st); + REQUIRE(tilted.ok); + CHECK(tilted.rot1_rad * 180.0 / PI == Catch::Approx(0.3).margin(0.01)); + CHECK(tilted.rot2_rad * 180.0 / PI == Catch::Approx(-0.2).margin(0.01)); + CHECK(tilted.beam_x_px == Catch::Approx(1000).margin(0.2)); + CHECK(tilted.beam_y_px == Catch::Approx(1000).margin(0.2)); + CHECK(tilted.distance_mm == Catch::Approx(150).margin(0.05)); + CHECK(tilted.median_residual_px < 0.1); + CHECK(tilted.capped_rms_px < held.capped_rms_px); + } +}