From 0d07141b5c8234c78d1496687278e3350fe906ce Mon Sep 17 00:00:00 2001 From: jungfrau Date: Mon, 24 Aug 2026 05:31:59 -0400 Subject: [PATCH] Ask the geometry refinement for the derivatives it uses, and solve the normal equations Four changes to the same least-squares fit, which the rotation first pass runs on every candidate lattice and the per-image path runs on every frame. The linear solver was DENSE_QR on a problem that is very tall and thin - thousands of spots against at most seventeen parameters. That is the shape QR handles worst: it copies the Jacobian out of Ceres' row-major storage into a column-major buffer on every solve, and Eigen's blocked Householder then degenerates to the unblocked path because its block size is the column count. Accumulating J^T J reads the Jacobian once instead. Both solve the same damped system, so the step is the same to round-off. Ceres sizes its dual numbers from the declared parameter blocks, not from which of them the caller then holds constant. Nothing outside a test set refine_distance_mm - the positional residual leaves the distance degenerate with the cell scale, which is why the rotation post-refinement fits it in a step of its own with the cell held fixed - so the block was declared only to be frozen, and every residual differentiated seventeen parameters to use sixteen. It is gone, along with the test that exercised distance recovery; that test seeded the distance off truth, which the cell would now absorb, so its seed moves to the true value. The post-refinement's own detector step held five of its seven blocks constant and now bakes them into the residual, leaving beam and distance. The predicted reciprocal vector was built by rotating all three direct columns and then crossing them. A rotation commutes with the cross product and leaves the triple product alone, so the same vector comes out of crossing the unrotated columns and turning the result once - three rotations become one, for every crystal system. The documentation described the arrangement before all this, and had drifted in a second way: the first-pass rotation indexing has been refining the detector tilt and the rotation axis by default, which the text said were held fixed. Co-Authored-By: Claude Opus 5 (1M context) --- docs/CPU_DATA_ANALYSIS.md | 2 +- image_analysis/IndexAndRefine.cpp | 2 - image_analysis/geom_refinement/PostRefine.cpp | 104 ++++++++++++------ .../geom_refinement/XtalOptimizer.cpp | 27 ++--- .../geom_refinement/XtalOptimizer.h | 1 - image_analysis/geom_refinement/XtalResidual.h | 101 +++++++++++++---- .../rotation_indexer/RotationIndexer.cpp | 2 - tests/XtalOptimizerTest.cpp | 10 +- 8 files changed, 168 insertions(+), 81 deletions(-) diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index de41f64e..9391fe29 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -468,7 +468,7 @@ The refinement jointly optimizes, depending on mode and constraints: - crystal orientation (a global rotation), - unit-cell parameters, with constraints determined by inferred crystal system. -By default only the beam center, unit cell and crystal orientation are refined; the detector distance, tilt angles and rotation-axis direction are held fixed unless explicitly enabled. A lighter **orientation-only** mode refines just the crystal orientation, for stills whose geometry is already trusted. It carries a weak small-rotation prior penalising the whole angle-axis vector (all three components, at a low weight); what it is there for is the poorly-determined out-of-plane component, which is the one the data barely constrain. +The detector distance is not refined against one crystal's spots at all: the positional residual leaves it degenerate with the cell scale, so it is fitted separately - by the rotation post-refinement, which holds the cell at the value its own cell/axis step settled on, and by the stills `--refine-geometry` bundle. Per image, the beam centre and the crystal orientation are refined, and the unit cell as well for stills. The first-pass rotation indexing refines the detector tilt and the rotation-axis direction too, against the spots accumulated across the sweep; everywhere else both are held fixed, because on a single crystal a tilt is absorbed almost exactly by the crystal orientation. A lighter **orientation-only** mode refines just the crystal orientation, for stills whose geometry is already trusted. It carries a weak small-rotation prior penalising the whole angle-axis vector (all three components, at a low weight); what it is there for is the poorly-determined out-of-plane component, which is the one the data barely constrain. For higher symmetries, constraints are enforced, e.g. - cubic: $a=b=c,\ \alpha=\beta=\gamma=90^\circ$, diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index d5a1034b..7e43d049 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -308,7 +308,6 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In .min_angle_deg = 30.0f, .max_angle_deg = 150.0f, .refine_beam_center = true, - .refine_distance_mm = false, .refine_detector_angles = false, .refine_unit_cell = !experiment.IsRotationIndexing(), // The whole spot list is passed below, not the indexed subset, so on weak images most of what @@ -402,7 +401,6 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In .crystal_system = data.crystal_system, .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(), .refine_beam_center = false, - .refine_distance_mm = false, .refine_detector_angles = false, .refine_unit_cell = false, .refine_rotation_axis = false, diff --git a/image_analysis/geom_refinement/PostRefine.cpp b/image_analysis/geom_refinement/PostRefine.cpp index cedc68a8..11b438a3 100644 --- a/image_analysis/geom_refinement/PostRefine.cpp +++ b/image_analysis/geom_refinement/PostRefine.cpp @@ -475,32 +475,47 @@ PostRefineResult PostRefineRotationGeometry(const std::vector(std::max(1, std::min(nthreads, static_cast(n_terms)))); const int cost_chunk = (n_terms + cost_nt - 1) / cost_nt; - const auto cost_by_fifth = [&](double k) { - std::vector> per_chunk(cost_nt); + // Several k are always wanted at once (the grid below asks for 101), and they all sweep the + // same event list, so sweep it ONCE and evaluate every k on each event while it is still in + // registers. The per-thread accumulator is one slot per (k, fifth) - 4 kB for the grid, small + // enough to stay in L1 - against re-reading the whole term array once per k. Each (k, fifth) + // still receives its events in the same order and the chunks are still summed in chunk order, + // so the sums are the ones a k-at-a-time loop produced, bit for bit. + const auto cost_grid = [&](const std::vector &ks) { + const int nk = static_cast(ks.size()); + std::vector>> per_chunk( + cost_nt, std::vector>(nk)); ParallelChunks(n_terms, nthreads, [&](int lo, int hi) { - std::array acc{}; + std::vector> acc(nk); for (int e = lo; e < hi; ++e) { - const double r = residual_at(terms[e], k); - const double s2 = r * r; - acc[fifth_of[e]] += (s2 <= huber_d2) ? s2 - : (2.0 * huber_delta * std::sqrt(s2) - huber_d2); + const ScaleTerm &term = terms[e]; + const int fifth = fifth_of[e]; + for (int g = 0; g < nk; ++g) { + const double r = residual_at(term, ks[g]); + const double s2 = r * r; + acc[g][fifth] += (s2 <= huber_d2) ? s2 + : (2.0 * huber_delta * std::sqrt(s2) - huber_d2); + } } - per_chunk[lo / cost_chunk] = acc; + per_chunk[lo / cost_chunk] = std::move(acc); }); - std::array total{}; + std::vector> total(nk); for (const auto &acc : per_chunk) - for (int j = 0; j < 5; ++j) total[j] += acc[j]; + for (int g = 0; g < nk; ++g) + for (int j = 0; j < 5; ++j) total[g][j] += acc[g][j]; return total; }; + const auto cost_by_fifth = [&](double k) { return cost_grid({k})[0]; }; // Scan the interval Ceres was bounded to, then close in. No event's phase can move by more than // a fraction of a period over an interval this narrow, so the objective has no structure the // grid could step over; the refinement is only there to place the minimum precisely. constexpr int SCALE_GRID = 101; constexpr double SCALE_K_LO = 0.95, SCALE_K_HI = 1.05; - std::vector> grid(SCALE_GRID); + std::vector grid_k(SCALE_GRID); for (int g = 0; g < SCALE_GRID; ++g) - grid[g] = cost_by_fifth(SCALE_K_LO + (SCALE_K_HI - SCALE_K_LO) * g / (SCALE_GRID - 1)); + grid_k[g] = SCALE_K_LO + (SCALE_K_HI - SCALE_K_LO) * g / (SCALE_GRID - 1); + const std::vector> grid = cost_grid(grid_k); auto solve_scale = [&](int drop_fifth) { const auto total = [&](const std::array &f) { double t = 0.0; @@ -601,23 +616,45 @@ PostRefineResult PostRefineRotationGeometry(const std::vector obs; - for (size_t i = 0; i < n_pts; ++i) - if (std::isfinite(pts[i].obs_x) && std::isfinite(pts[i].obs_y)) obs.push_back(&pts[i]); + // Count first, then fill, exactly as the partial gather above does and for the same reason: + // this walks the same tens of millions of partials, and a pointer vector grown by push_back + // copies itself every time it doubles. new[] rather than a sized vector so the array is not + // zeroed on one thread before the parallel fill overwrites it. The fill lands in the order + // the serial loop produced, so the selection below sees the same sequence it always did. + const int n_obs_chunks = static_cast(std::clamp(nthreads, 1, + std::max(1, n_pts))); + const size_t obs_chunk = (n_pts + n_obs_chunks - 1) / n_obs_chunks; + std::vector obs_offset(n_obs_chunks + 1, 0); + const auto keep_obs = [&](size_t i) { + return std::isfinite(pts[i].obs_x) && std::isfinite(pts[i].obs_y); + }; + ParallelChunks(static_cast(n_pts), nthreads, [&](int lo, int hi) { + size_t keep = 0; + for (int i = lo; i < hi; ++i) if (keep_obs(i)) keep++; + obs_offset[static_cast(lo) / obs_chunk + 1] = keep; + }); + for (int c = 0; c < n_obs_chunks; ++c) obs_offset[c + 1] += obs_offset[c]; + size_t n_obs = obs_offset[n_obs_chunks]; + std::unique_ptr obs(new const Partial *[n_obs]); + ParallelChunks(static_cast(n_pts), nthreads, [&](int lo, int hi) { + size_t at = obs_offset[static_cast(lo) / obs_chunk]; + for (int i = lo; i < hi; ++i) if (keep_obs(i)) obs[at++] = &pts[i]; + }); constexpr size_t MAX_OBS = 20000; - if (obs.size() > MAX_OBS) { - std::nth_element(obs.begin(), obs.begin() + MAX_OBS, obs.end(), + if (n_obs > MAX_OBS) { + std::nth_element(obs.get(), obs.get() + MAX_OBS, obs.get() + n_obs, [](const Partial *a, const Partial *b) { return a->I / std::max(1e-9, static_cast(a->sigma)) > b->I / std::max(1e-9, static_cast(b->sigma)); }); - obs.resize(MAX_OBS); + n_obs = MAX_OBS; } - result.obs_used = static_cast(obs.size()); + result.obs_used = static_cast(n_obs); const double beam_x0 = nominal_geom.GetBeamX_pxl(), beam_y0 = nominal_geom.GetBeamY_pxl(); const double dist0 = nominal_geom.GetDetectorDistance_mm(); auto pos_cost = [&](Subset s, const double beam[2], const double dist[1]) { double c = 0.0; int n = 0; - for (const Partial *pp : obs) { + for (size_t oi = 0; oi < n_obs; ++oi) { + const Partial *pp = obs[oi]; if (!in(pp->h, pp->k, pp->l, s)) continue; XtalResidual r(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, rot3, angle_rad(pp->img), pp->h, pp->k, pp->l, sys); @@ -630,21 +667,20 @@ PostRefineResult PostRefineRotationGeometry(const std::vectorh, pp->k, pp->l, s)) continue; - p.AddResidualBlock(new ceres::AutoDiffCostFunction( - new XtalResidual(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, rot3, angle_rad(pp->img), - pp->h, pp->k, pp->l, sys)), - new ceres::CauchyLoss(0.02), beam, dist, - const_cast(det_rot), rot_vec, p0, p1, p2); + p.AddResidualBlock(new ceres::AutoDiffCostFunction( + new XtalResidualBeamDistance( + XtalResidual(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, rot3, angle_rad(pp->img), + pp->h, pp->k, pp->l, sys), + det_rot, rot_vec, p0, p1, p2)), + new ceres::CauchyLoss(0.02), beam, dist); } if (p.NumResidualBlocks() == 0) { beam_out[0] = beam_x0; beam_out[1] = beam_y0; dist_out = dist0; return false; } - // Everything but the beam and the distance is held at its step-A value. Once per problem: - // the block only has to exist, and repeating it per observation was up to MAX_OBS times - // five calls for the same five blocks. - p.SetParameterBlockConstant(const_cast(det_rot)); - p.SetParameterBlockConstant(rot_vec); - p.SetParameterBlockConstant(p0); p.SetParameterBlockConstant(p1); p.SetParameterBlockConstant(p2); + // Everything but the beam and the distance is held at its step-A value, and the residual + // above bakes those five blocks in rather than declaring them and freezing them, so there + // is nothing left to hold constant here. p.SetParameterLowerBound(dist, 0, dist0 * 0.95); p.SetParameterUpperBound(dist, 0, dist0 * 1.05); for (int j = 0; j < 2; ++j) { p.SetParameterLowerBound(beam, j, beam[j] - 15.0); p.SetParameterUpperBound(beam, j, beam[j] + 15.0); } @@ -655,7 +691,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector= static_cast(settings.min_events)) { + if (n_obs >= static_cast(settings.min_events)) { double beam_fit[2], dist_fit; const bool convB = solve_detector(FIT, beam_fit, dist_fit); const double b_nom[2] = {beam_x0, beam_y0}, d_nom[1] = {dist0}; @@ -679,7 +715,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector( - new XtalResidual(residual)), + new ceres::AutoDiffCostFunction( + new XtalResidualFixedDistance(residual, distance_mm)), loss, beam, - &distance_mm, detector_rot, rot_vec, latt_vec0, @@ -318,14 +316,6 @@ bool XtalOptimizerInternal(XtalOptimizerData &data, // Distance, detector angles, rotation axis and cell are parameter blocks only in the general // seven-block residual; the reduced one bakes them in, so there is nothing left to configure. if (!beam_and_orientation_only) { - if (!data.refine_distance_mm) - problem.SetParameterBlockConstant(&distance_mm); - else { - const double dist_range = 0.1; - problem.SetParameterLowerBound(&distance_mm, 0, distance_mm * (1.0 - dist_range)); - problem.SetParameterUpperBound(&distance_mm, 0, distance_mm * (1.0 + dist_range)); - } - if (!data.refine_detector_angles) { problem.SetParameterBlockConstant(detector_rot); } else { @@ -374,7 +364,14 @@ bool XtalOptimizerInternal(XtalOptimizerData &data, // Configure solver ceres::Solver::Options options; - options.linear_solver_type = ceres::DENSE_QR; + // Normal equations, not QR. The problem is very tall and thin - thousands of spots against at + // most 17 parameters - and that is the shape DENSE_QR handles worst: it copies the Jacobian out + // of Ceres' row-major storage into a column-major buffer on every solve, and Eigen's blocked + // Householder then degenerates to the unblocked path because its block size is min(48, columns). + // Accumulating J^T J reads the Jacobian once instead. Both solve the same damped system, so the + // step is the same to round-off; the column scaling Ceres applies by default and the LM diagonal + // keep the squared condition number in hand. + options.linear_solver_type = ceres::DENSE_NORMAL_CHOLESKY; options.minimizer_progress_to_stdout = false; if (data.max_iterations > 0) options.max_num_iterations = data.max_iterations; @@ -400,8 +397,6 @@ bool XtalOptimizerInternal(XtalOptimizerData &data, data.geom.BeamX_pxl(beam[0]).BeamY_pxl(beam[1]); } - if (data.refine_distance_mm) - data.geom.DetectorDistance_mm(distance_mm); if (data.refine_detector_angles) data.geom.PoniRot1_rad(detector_rot[0]).PoniRot2_rad(detector_rot[1]); @@ -516,7 +511,7 @@ bool XtalOptimizerRotationOnly(XtalOptimizerData &data, } ceres::Solver::Options options; - options.linear_solver_type = ceres::DENSE_QR; + options.linear_solver_type = ceres::DENSE_NORMAL_CHOLESKY; // tall and thin, as above options.minimizer_progress_to_stdout = false; if (data.max_iterations > 0) options.max_num_iterations = data.max_iterations; diff --git a/image_analysis/geom_refinement/XtalOptimizer.h b/image_analysis/geom_refinement/XtalOptimizer.h index 35375be6..6d107a63 100644 --- a/image_analysis/geom_refinement/XtalOptimizer.h +++ b/image_analysis/geom_refinement/XtalOptimizer.h @@ -23,7 +23,6 @@ struct XtalOptimizerData { float max_angle_deg = 120.0f; bool refine_beam_center = true; - bool refine_distance_mm = false; bool refine_detector_angles = false; bool refine_unit_cell = true; // This refines unit cell size + angles - orientation is always refined bool refine_rotation_axis = false; diff --git a/image_analysis/geom_refinement/XtalResidual.h b/image_analysis/geom_refinement/XtalResidual.h index 93f8dcc3..72c28145 100644 --- a/image_analysis/geom_refinement/XtalResidual.h +++ b/image_analysis/geom_refinement/XtalResidual.h @@ -228,32 +228,32 @@ struct XtalResidual { T col1_unrot[3] = {B(0, 1) * L1, B(1, 1) * L1, B(2, 1) * L1}; 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]; - 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); + // Build the reciprocal vector in the UNROTATED frame and turn it once at the end, rather than + // rotating all three direct columns first. A rotation commutes with the cross product and leaves + // the triple product alone - (Rb)x(Rc) = R(bxc), and (Ra).(R(bxc)) = a.(bxc) - so the reciprocal + // basis of the rotated cell is the rotation of the unrotated one, and h a* + k b* + l c* is the + // rotation of the same combination built from the unrotated columns. Three rotations become one, + // and the cross products now run on columns that are mostly zero for every system but triclinic. + const Eigen::Matrix a_unrot(col0_unrot[0], col0_unrot[1], col0_unrot[2]); + const Eigen::Matrix b_unrot(col1_unrot[0], col1_unrot[1], col1_unrot[2]); + const Eigen::Matrix c_unrot(col2_unrot[0], col2_unrot[1], col2_unrot[2]); - 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]); - const Eigen::Matrix C(col2_rot[0], col2_rot[1], col2_rot[2]); + const Eigen::Matrix bxc = b_unrot.cross(c_unrot); + const Eigen::Matrix cxa = c_unrot.cross(a_unrot); + const Eigen::Matrix axb = a_unrot.cross(b_unrot); - const Eigen::Matrix BxC = Bv.cross(C); - const Eigen::Matrix CxA = C.cross(A); - const Eigen::Matrix AxB = A.cross(Bv); - - const T V = A.dot(BxC); - const T invV = T(1) / V; - - const Eigen::Matrix Astar = BxC * invV; - const Eigen::Matrix Bstar = CxA * invV; - const Eigen::Matrix Cstar = AxB * invV; + const T invV = T(1) / a_unrot.dot(bxc); const T h = T(exp_h); const T k = T(exp_k); const T l = T(exp_l); - const Eigen::Matrix e_pred_recip = Astar * h + Bstar * k + Cstar * l; + const Eigen::Matrix recip_unrot = (bxc * h + cxa * k + axb * l) * invV; + + T recip_in[3] = {recip_unrot[0], recip_unrot[1], recip_unrot[2]}; + T e_pred_recip[3]; + const AngleAxisRotator rot_p0(p0); + rot_p0.Rotate(recip_in, e_pred_recip); residual[0] = e_obs_recip[0] - e_pred_recip[0]; residual[1] = e_obs_recip[1] - e_pred_recip[1]; @@ -273,6 +273,69 @@ struct XtalResidual { gemmi::CrystalSystem symmetry; }; +// Same residual with the detector angles, the rotation axis and the whole crystal - orientation, cell +// lengths, cell angles - baked in as constants, so that only beam(2) and distance(1) remain parameter +// blocks. That is the rotation post-refinement's detector step, which deliberately holds the crystal at +// the value its cell/axis step settled on: the seven-block form differentiates 17 parameters in order to +// use 3, this one runs on Jet. It forwards to XtalResidual with the very same values, so the +// residual and the beam/distance columns of its Jacobian are unchanged. +struct XtalResidualBeamDistance { + XtalResidualBeamDistance(const XtalResidual &residual, + const double *detector_rot, + const double *rotation_axis, + const double *orientation, + const double *uc_len, + const double *uc_angle) + : residual(residual), + detector_rot{detector_rot[0], detector_rot[1]}, + rotation_axis{rotation_axis[0], rotation_axis[1], rotation_axis[2]}, + orientation{orientation[0], orientation[1], orientation[2]}, + uc_len{uc_len[0], uc_len[1], uc_len[2]}, + uc_angle{uc_angle[0], uc_angle[1], uc_angle[2]} { + } + + template + bool operator()(const T *const beam, const T *const distance, T *residual_out) const { + const T rot[2] = {T(detector_rot[0]), T(detector_rot[1])}; + const T axis[3] = {T(rotation_axis[0]), T(rotation_axis[1]), T(rotation_axis[2])}; + const T p0[3] = {T(orientation[0]), T(orientation[1]), T(orientation[2])}; + const T p1[3] = {T(uc_len[0]), T(uc_len[1]), T(uc_len[2])}; + const T p2[3] = {T(uc_angle[0]), T(uc_angle[1]), T(uc_angle[2])}; + return residual(beam, distance, rot, axis, p0, p1, p2, residual_out); + } + + const XtalResidual residual; + const double detector_rot[2]; + const double rotation_axis[3]; + const double orientation[3]; + const double uc_len[3]; + const double uc_angle[3]; +}; + +// Same residual with only the detector distance baked in, leaving the other six blocks free: 16 +// declared parameters instead of 17. Nothing in rugnux refines the distance from a single crystal - the +// rotation post-refinement fits it globally, with the cell held fixed so the two are no longer +// degenerate - so the block was only ever declared in order to be frozen. Dropping it is worth more +// than one parameter in seventeen suggests: Jet keeps its derivatives in an Eigen vector, and +// 17 doubles is one lane past the four-wide boundary, so it costs five packet operations and cannot be +// over-aligned where 16 costs four and can. The distance is forwarded unchanged, so the residual and +// every remaining Jacobian column are the same. +struct XtalResidualFixedDistance { + XtalResidualFixedDistance(const XtalResidual &residual, double distance_mm) + : residual(residual), distance_mm(distance_mm) { + } + + template + bool operator()(const T *const beam, const T *const detector_rot, const T *const rotation_axis, + const T *const p0, const T *const p1, const T *const p2, T *residual_out) const { + const T distance[1] = {T(distance_mm)}; + return residual(beam, distance, detector_rot, rotation_axis, p0, p1, p2, residual_out); + } + + const XtalResidual residual; + const double distance_mm; +}; + // Same residual with distance, detector angles, rotation axis and unit cell baked in as constants, so // that only beam(2) and orientation(3) remain parameter blocks. Ceres sizes its autodiff dual numbers // from the DECLARED blocks, not from which of them the caller then holds constant, so the seven-block diff --git a/image_analysis/rotation_indexer/RotationIndexer.cpp b/image_analysis/rotation_indexer/RotationIndexer.cpp index 16965c50..8b2c62e3 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.cpp +++ b/image_analysis/rotation_indexer/RotationIndexer.cpp @@ -212,7 +212,6 @@ void RotationIndexer::RunIndexing() { .min_angle_deg = 30.0f, .max_angle_deg = 150.0f, .refine_beam_center = true, - .refine_distance_mm = false, .refine_detector_angles = true, .refine_rotation_axis = true, .index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings(), @@ -425,7 +424,6 @@ void RotationIndexer::RunIndexing() { .crystal_system = search_result_.system, .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(), .refine_beam_center = false, - .refine_distance_mm = false, .refine_detector_angles = false, .refine_unit_cell = false, .refine_rotation_axis = false, diff --git a/tests/XtalOptimizerTest.cpp b/tests/XtalOptimizerTest.cpp index f200693a..3e4a69d7 100644 --- a/tests/XtalOptimizerTest.cpp +++ b/tests/XtalOptimizerTest.cpp @@ -543,15 +543,17 @@ TEST_CASE("XtalOptimizer_rotation") { // direction that XtalOptimizer now deliberately restrains toward the header rather than refining // freely (see the BeamComponentPrior in XtalOptimizer.cpp). So only the perpendicular component // (Y) is seeded off-truth to exercise beam-centre recovery; X is seeded at its true value. + // The distance is seeded at its true value too: XtalOptimizer does not refine it - the rotation + // post-refinement fits it globally instead, where the cell is held fixed and the two are no longer + // degenerate - so a distance seeded off-truth here would simply be absorbed by the cell. XtalOptimizerData xtal_opt{}; xtal_opt.latt = CrystalLattice(39.7f, 50.6f, 79.6f, 90.0f, 94.5f, 90.5f); - xtal_opt.geom.BeamX_pxl(1000).BeamY_pxl(997).DetectorDistance_mm(203.0) + xtal_opt.geom.BeamX_pxl(1000).BeamY_pxl(997).DetectorDistance_mm(200.0) .PoniRot1_rad(0.01).PoniRot2_rad(0.02); xtal_opt.crystal_system = gemmi::CrystalSystem::Monoclinic; xtal_opt.axis = axis; xtal_opt.min_spots = 200; xtal_opt.refine_beam_center = true; - xtal_opt.refine_distance_mm = true; xtal_opt.refine_detector_angles = false; xtal_opt.max_time = 30.0; @@ -568,9 +570,6 @@ TEST_CASE("XtalOptimizer_rotation") { CHECK(fabsf(xtal_opt.geom.GetBeamX_pxl() - exp_i.GetBeamX_pxl()) < 0.2f); CHECK(fabsf(xtal_opt.geom.GetBeamY_pxl() - exp_i.GetBeamY_pxl()) < 0.2f); - // 700 micron is a lot, but for now it is important to see that refinement improves in the right direction - CHECK(fabsf(xtal_opt.geom.GetDetectorDistance_mm() - exp_i.GetDetectorDistance_mm()) < 0.7f); - // Lattice checks CHECK(fabsf(uc_ref.a - uc_out.a) < 0.2f); CHECK(fabsf(uc_ref.b - uc_out.b) < 0.2f); @@ -641,7 +640,6 @@ TEST_CASE("XtalOptimizer_refine_rotation_axis") { std::nullopt); xtal_opt.min_spots = 200; xtal_opt.refine_beam_center = true; - xtal_opt.refine_distance_mm = false; xtal_opt.refine_detector_angles = false; xtal_opt.refine_rotation_axis = true; xtal_opt.max_time = 30.0;