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); + } +}