diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index c6455113..d5a1034b 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -71,7 +71,7 @@ IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool indexer_(indexer), rotation_indexer_counter(x) { if (indexer && x.IsRotationIndexing()) - rotation_indexer = std::make_unique(x, *indexer); + rotation_indexer = std::make_unique(x, *indexer, real_time); // Only retain the whole-run per-image reflections when a later scaling/merge pass will read them. if (retain_outcomes_) integration_outcome.resize(x.GetImageNum()); diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index d79c2b0a..d10e66c7 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -92,10 +92,10 @@ class IndexAndRefine { std::optional RotationAngle(int64_t image) const; // mid-exposure angle for the indexer public: - // real_time: bound the per-image geometry refinement by WALL CLOCK, as online acquisition must - - // it has a real per-image budget. Offline (rugnux, the viewer) passes false and the refinement is - // bounded by iteration count instead, so the same file reprocesses to the same answer regardless of - // what else the machine was doing. + // real_time: bound the geometry refinements - the per-image one here and the candidate-cell ones + // in the rotation indexer - by WALL CLOCK, as online acquisition must, it having a real budget. + // Offline (rugnux, the viewer) passes false and they are bounded by iteration count instead, so the + // same file reprocesses to the same answer regardless of what else the machine was doing. IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer, bool retain_outcomes = true, bool real_time = false); diff --git a/image_analysis/geom_refinement/XtalOptimizer.cpp b/image_analysis/geom_refinement/XtalOptimizer.cpp index 348cc9b0..1bd2a7a7 100644 --- a/image_analysis/geom_refinement/XtalOptimizer.cpp +++ b/image_analysis/geom_refinement/XtalOptimizer.cpp @@ -44,9 +44,10 @@ struct XtalResidualRotationOnlyPrecomp { T astar_rot[3], bstar_rot[3], cstar_rot[3]; - ceres::AngleAxisRotatePoint(rot_aa, astar_unrot, astar_rot); - ceres::AngleAxisRotatePoint(rot_aa, bstar_unrot, bstar_rot); - ceres::AngleAxisRotatePoint(rot_aa, cstar_unrot, cstar_rot); + const AngleAxisRotator rot(rot_aa); + rot.Rotate(astar_unrot, astar_rot); + rot.Rotate(bstar_unrot, bstar_rot); + rot.Rotate(cstar_unrot, cstar_rot); const Eigen::Matrix s_pred(T(h) * astar_rot[0] + T(k) * bstar_rot[0] + T(l) * cstar_rot[0], T(h) * astar_rot[1] + T(k) * bstar_rot[1] + T(l) * cstar_rot[1], diff --git a/image_analysis/geom_refinement/XtalResidual.h b/image_analysis/geom_refinement/XtalResidual.h index cf6a5a4d..93f8dcc3 100644 --- a/image_analysis/geom_refinement/XtalResidual.h +++ b/image_analysis/geom_refinement/XtalResidual.h @@ -13,6 +13,62 @@ #include "../../common/JFJochException.h" +// Rodrigues rotation with everything that depends only on the angle-axis worked out once. This is +// ceres::AngleAxisRotatePoint term for term - so the rotated point is bit-identical to it - but that +// function derives theta and its sine and cosine from the angle-axis on every call, which is pure +// repetition when one rotation is applied to several points, as it is to the three lattice columns +// below. On a Jet those trigonometric evaluations are the expensive part of the rotation. +template +struct AngleAxisRotator { + explicit AngleAxisRotator(const T *const angle_axis) + : aa{angle_axis[0], angle_axis[1], angle_axis[2]} { + using std::cos; + using std::fpclassify; + using std::hypot; + using std::sin; + + const T theta = hypot(aa[0], aa[1], aa[2]); + at_zero = (fpclassify(theta) == FP_ZERO); + if (at_zero) + return; + + costheta = cos(theta); + sintheta = sin(theta); + const T theta_inverse = T(1.0) / theta; + w[0] = aa[0] * theta_inverse; + w[1] = aa[1] * theta_inverse; + w[2] = aa[2] * theta_inverse; + } + + void Rotate(const T pt[3], T result[3]) const { + if (at_zero) { + // At zero the rotation is I + hat(angle_axis); the Taylor form keeps the derivatives + // meaningful where the angle-axis vanishes. + const T w_cross_pt[3] = {aa[1] * pt[2] - aa[2] * pt[1], + aa[2] * pt[0] - aa[0] * pt[2], + aa[0] * pt[1] - aa[1] * pt[0]}; + result[0] = pt[0] + w_cross_pt[0]; + result[1] = pt[1] + w_cross_pt[1]; + result[2] = pt[2] + w_cross_pt[2]; + return; + } + + const T w_cross_pt[3] = {w[1] * pt[2] - w[2] * pt[1], + w[2] * pt[0] - w[0] * pt[2], + w[0] * pt[1] - w[1] * pt[0]}; + const T tmp = (w[0] * pt[0] + w[1] * pt[1] + w[2] * pt[2]) * (T(1.0) - costheta); + + result[0] = pt[0] * costheta + w_cross_pt[0] * sintheta + w[0] * tmp; + result[1] = pt[1] * costheta + w_cross_pt[1] * sintheta + w[1] * tmp; + result[2] = pt[2] * costheta + w_cross_pt[2] * sintheta + w[2] * tmp; + } + + const T aa[3]; + T w[3]; + T costheta, sintheta; + bool at_zero; +}; + // Detector -> reciprocal geometry residual, shared by the per-image XtalOptimizer (one lattice, one // frame) and the offline GeometryRefiner (shared beam/distance/cell blocks, one orientation block per // frame). Parameter blocks: beam(2), distance_mm(1), detector_rot(2 = rot1,rot2), rotation_axis(3), @@ -32,7 +88,8 @@ struct XtalResidual { : obs_x(x), obs_y(y), inv_lambda(1.0/lambda), pixel_size(pixel_size), - rot3(rot3), + cos_rot3(std::cos(rot3)), + sin_rot3(std::sin(rot3)), exp_h(exp_h), exp_k(exp_k), exp_l(exp_l), @@ -67,9 +124,11 @@ struct XtalResidual { const T c2 = ceres::cos(rot2); const T s2 = ceres::sin(rot2); - // Rz(-rot3): rotation around Z (beam); constant, identity when rot3 == 0 - const T c3 = T(cos(rot3)); - const T s3 = T(sin(rot3)); + // Rz(-rot3): rotation around Z (beam); constant, identity when rot3 == 0. Its sine and cosine + // are taken in the constructor - they do not depend on any parameter, so recomputing them per + // evaluation is two libm calls per residual for nothing. + const T c3 = T(cos_rot3); + const T s3 = T(sin_rot3); // Detector coordinates in mm const T det_x = (T(obs_x) - beam[0]) * T(pixel_size); @@ -170,9 +229,10 @@ struct XtalResidual { T col2_unrot[3] = {B(0, 2) * L2, B(1, 2) * L2, B(2, 2) * L2}; T col0_rot[3], col1_rot[3], col2_rot[3]; - ceres::AngleAxisRotatePoint(p0, col0_unrot, col0_rot); - ceres::AngleAxisRotatePoint(p0, col1_unrot, col1_rot); - ceres::AngleAxisRotatePoint(p0, col2_unrot, col2_rot); + const AngleAxisRotator rot_p0(p0); + rot_p0.Rotate(col0_unrot, col0_rot); + rot_p0.Rotate(col1_unrot, col1_rot); + rot_p0.Rotate(col2_unrot, col2_rot); const Eigen::Matrix A(col0_rot[0], col0_rot[1], col0_rot[2]); const Eigen::Matrix Bv(col1_rot[0], col1_rot[1], col1_rot[2]); @@ -205,7 +265,7 @@ struct XtalResidual { const double obs_x, obs_y; const double inv_lambda; const double pixel_size; - const double rot3; + const double cos_rot3, sin_rot3; const double exp_h; const double exp_k; const double exp_l; diff --git a/image_analysis/rotation_indexer/RotationIndexer.cpp b/image_analysis/rotation_indexer/RotationIndexer.cpp index 1276ba2b..16965c50 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.cpp +++ b/image_analysis/rotation_indexer/RotationIndexer.cpp @@ -39,6 +39,12 @@ namespace { // How far from a whole number the volume ratio may sit and still count as an axis multiple. constexpr double ROT_SUPERCELL_INTEGER_TOL = 0.15; + // Iteration bound for the candidate-cell refinements when the run is not real-time. Without one, + // XtalOptimizerData falls back to its wall-clock bound and the cell a candidate refines to - and + // therefore which candidate wins - depends on how busy the machine was. Same value as Ceres' own + // default iteration limit, so it only bites where the clock was biting before. + constexpr int ROT_REFINE_ITERATIONS = 50; + // Order of the lattice point group, so "lower symmetry" is a well-defined comparison // (gemmi's enum orders Trigonal after Tetragonal, which have 6 and 8 rotations). int LatticePointGroupOrder(gemmi::CrystalSystem s) { @@ -106,9 +112,11 @@ namespace { } } -RotationIndexer::RotationIndexer(const DiffractionExperiment &x, IndexerThreadPool &indexer) +RotationIndexer::RotationIndexer(const DiffractionExperiment &x, IndexerThreadPool &indexer, + bool real_time) : experiment(x), index_ice_rings(x.GetIndexingSettings().GetIndexIceRings()), + real_time(real_time), v_(experiment.GetImageNum()), angle_deg_(experiment.GetImageNum()), axis_(x.GetGoniometer()), @@ -208,6 +216,7 @@ void RotationIndexer::RunIndexing() { .refine_detector_angles = true, .refine_rotation_axis = true, .index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings(), + .max_iterations = real_time ? 0 : ROT_REFINE_ITERATIONS, .axis = orig_axis }; if (d.crystal_system == gemmi::CrystalSystem::Trigonal) @@ -421,6 +430,7 @@ void RotationIndexer::RunIndexing() { .refine_unit_cell = false, .refine_rotation_axis = false, .index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings(), + .max_iterations = real_time ? 0 : ROT_REFINE_ITERATIONS, .axis = axis_ }; diff --git a/image_analysis/rotation_indexer/RotationIndexer.h b/image_analysis/rotation_indexer/RotationIndexer.h index 1fe5de58..d9b8dd51 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.h +++ b/image_analysis/rotation_indexer/RotationIndexer.h @@ -34,6 +34,7 @@ class RotationIndexer { constexpr static size_t max_spots_per_image = 200; const bool index_ice_rings; + const bool real_time; // see the constructor std::vector> v_; // Per-image rotation angle (mid-exposure, deg) supplied by the caller; falls back to the @@ -53,7 +54,11 @@ class RotationIndexer { std::optional indexed_lattice; public: - RotationIndexer(const DiffractionExperiment& x, IndexerThreadPool& indexer); + // real_time: bound the candidate-cell refinements by WALL CLOCK, as online acquisition must - it + // has a real budget. Offline (rugnux, the viewer) passes false and they are bounded by iteration + // count instead, so the same file reprocesses to the same cell regardless of what else the machine + // was doing at the time. + RotationIndexer(const DiffractionExperiment& x, IndexerThreadPool& indexer, bool real_time = false); // angle_deg is the image's mid-exposure rotation angle; if omitted, the goniometer angle at // `image` is used (only valid when `image` is the goniometer's own image index). void ProcessImage(int64_t image, const std::vector& spots,