diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 2b331a909..03ff33cde 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -8,6 +8,7 @@ * `rugnux` and `jfjoch_viewer` read PILATUS miniCBF sweeps natively, and open masters written at other facilities, including Eiger 1.x and third-party NXmx. * `rugnux` determines the lattice and the space group more reliably - the true cell where the first pass offers a whole-number multiple of it, so a pseudo-translated crystal keeps its full-length axis and a small molecule is indexed on its own cell rather than a protein-sized one, and the point group, the setting and the systematic absences - and `-S` refuses or re-seats a fixed space group whose symmetry axes the indexed cell does not carry. * `rugnux` asks the space-group search again on the cell the lattice metric supports, whenever that metric carries more rotational symmetry than the Bravais class the indexer named, so a lattice whose reduction landed in a lower-symmetry sub-cell can still reach its true point group; the higher symmetry is adopted only where the intensities confirm it. +* `rugnux` post-refines the rotation geometry with one joint fit of the crystal and the detector against both the observed spot positions and the observed rocking angles, in place of scaling the cell against the angles and then reading the detector distance off the scaled cell; the refined distance no longer depends on how wrong the file's distance was. * `rugnux` measures the beam centre on every run and indexes with it when the file's value indexes nothing, refines only the detector-tilt component the data determine - a beam-centre error is no longer reported as a tilt - and places a detector swung out on a 2theta arm where the file says it stands. * `rugnux` writes the unmerged MTZ by default, with a P1 merge beside it, a batch header for every image the observations span, and events kept to the same `--min-captured-fraction` as the merge, so a wrong space group can be re-merged in a scaling program without reprocessing. * `rugnux` writes reflection files in the conventions downstream programs read: `FreeR_flag` is 0 for the test set and 1 for the working set - it was the other way round - the merged and P1 MTZ carry the reserved `HKL_base` dataset so a CCP4 program reads the wavelength instead of falling back to 1.54187 A, and the merged mmCIF marks the free set as `_refln.status` = `f`. diff --git a/docs/CPU_DATA_ANALYSIS_INDEXING.md b/docs/CPU_DATA_ANALYSIS_INDEXING.md index b6da4f498..9037bc34e 100644 --- a/docs/CPU_DATA_ANALYSIS_INDEXING.md +++ b/docs/CPU_DATA_ANALYSIS_INDEXING.md @@ -165,7 +165,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. -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. +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 elsewhere - by the rotation post-refinement, which frees it alongside the whole crystal and adds the distance-independent rocking-angle residual that breaks the degeneracy, 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$, @@ -206,18 +206,20 @@ The loose first stage necessarily admits some spots that are not reflections of The refinement above (§7.2) runs per image against that image's spots. For rotation data an additional **post-refinement** (on by default; `--rotation-no-postrefine` disables it) improves the detector distance, beam centre and crystal cell/axis using **all** frames at once, then re-integrates: 1. **Pass 1** integrates, scales and merges at the header geometry. -2. From pass-1's integrated reflections, the geometry is refined over all frames (Ceres, robust loss) in **two separate steps** rather than one joint fit: - - **Step A**: crystal cell scale + goniometer-axis direction, from the observed rotation angles (a distance-independent excitation residual). - - **Step B**: shared detector distance + beam centre, from the observed spot positions, with the cell held at step A — so the positional residual is no longer degenerate with the cell scale. +2. From pass-1's integrated reflections, the crystal and the detector are refined together over all frames (Ceres, robust loss) in **one joint fit**, against both residuals at once: + - the **positional** detector↔reciprocal residual at each partial's observed spot, and + - a distance-independent **Ewald excitation** residual at each reflection's observed rocking centroid $\phi_\mathrm{obs}$. - Each step is **cross-validated** on a deterministic split of the *reflections* (an avalanche-mixed $hkl$ hash, not a frame split and not an $h+k+l$ parity, which would collide with a centering condition and leave the held-out half empty): fitted on one half, committed only if it lowers the held-out residual, otherwise left at nominal. The solver bounds the move — distance within ±5 %, beam centre within ±15 px — and detector tilt is held fixed, being gauge-coupled to the crystal orientation on a single crystal. -3. **Pass 2** re-indexes de novo and re-integrates at the committed geometry. Only the **detector distance and beam centre** carry over: the refined cell and axis are used to make step B well-posed, but pass 2 re-indexes from scratch, so they are not propagated. + Free: the crystal orientation, the unit cell (every parameter the crystal system leaves free, not one overall scale), the goniometer-axis direction, the detector distance and the beam centre. The positional residual on its own *is* degenerate with the cell scale — that is why this used to be split into a cell-scale step and a distance step — but the excitation residual does not involve the detector at all, so it fixes the absolute size of the reciprocal lattice and breaks the degeneracy inside the same problem. Splitting it instead cost accuracy twice over: pass 1 frees the whole lattice against a frozen distance, so the distortion it absorbs is *anisotropic* and no single scale can undo it; and whatever bias is left in that scale goes straight into the distance, which is only ever determined relative to the cell. + + The fit is **cross-validated** on a deterministic split of the *reflections* (an avalanche-mixed $hkl$ hash, not a frame split and not an $h+k+l$ parity, which would collide with a centering condition and leave the held-out half empty): fitted on one half, committed only if it lowers the held-out residual and the move stays small — distance within 1 %, every cell length within 1 %, beam centre within 15 px of the nearest centre anything already believes — otherwise left at nominal. Detector tilt is held fixed, being gauge-coupled to the crystal orientation on a single crystal. +3. **Pass 2** re-indexes de novo and re-integrates at the committed geometry. Only the **detector distance and beam centre** carry over: the refined cell, orientation and axis are what make the distance identifiable, but pass 2 re-indexes from scratch, so they are not propagated. The space group is determined **after** pass 2, on the geometry the run refined, and pass 1 does not search at all: a decision taken on the worse of the two passes and then carried forward is a constraint on the better one, and had to be reconciled with what pass 2 later found. The guard that chooses which pass is written compares each pass's **first** merge — $P1$ on both sides, full resolution range, before the correction surfaces — which both passes produce anyway, so it never compares statistics computed in two different space groups. One index-time veto remains and is keyed to pass 1's **lattice** rather than its group: a centred pass-1 lattice against a primitive pass-2 one. Only pass 2 is written, as the canonical `_*` output. Pass 1's merge exists to give the guard something to judge pass 2 against, so it stops short of the parts of the merge that only fill in a file — the correction surfaces, the twinning and radiation-damage analyses, the R-free flags and the amplitudes — and writes no merged files of its own. -**Goniometer rotation scale (report only).** A stage that turns further than it was commanded to leaves no trace in the file, because the stored $\omega$ values *are* the commanded ones; the excess then presents as the crystal drifting, in this program and in others. Step A already measures it without a new degree of freedom: its residual rotates by $-\phi\,\mathbf{u}$ with $\mathbf{u}$ an **unnormalised** 3-vector, so $|\mathbf{u}|$ is the factor by which the stage actually turned, and normalising the axis throws it away. It is reported, and warned about beyond 0.5 %, under the same cross-validation that gates the cell move — a fold that merely soaked up noise cannot raise the flag. It is a detector, not a calibration: nothing corrects the data, and it **under-reads** the true magnitude, because the fit only sees reflections that indexed at the nominal angle and per-frame orientation refinement has already absorbed part of the error. +**Goniometer rotation scale (report only).** A stage that turns further than it was commanded to leaves no trace in the file, because the stored $\omega$ values *are* the commanded ones; the excess then presents as the crystal drifting, in this program and in others. The excitation residual already measures it without a new degree of freedom: it rotates by $-\phi\,\mathbf{u}$ with $\mathbf{u}$ an **unnormalised** 3-vector, so $|\mathbf{u}|$ is the factor by which the stage actually turned, and normalising the axis throws it away. It is reported, and warned about beyond 0.5 %, under its own leave-a-fifth-of-the-sweep-out check — a fold that merely soaked up noise cannot raise the flag. It is a detector, not a calibration: nothing corrects the data, and it **under-reads** the true magnitude, because the fit only sees reflections that indexed at the nominal angle and per-frame orientation refinement has already absorbed part of the error. ### 7.6 Detector geometry from powder rings diff --git a/image_analysis/geom_refinement/PostRefine.cpp b/image_analysis/geom_refinement/PostRefine.cpp index a900a9499..f9b0a9a05 100644 --- a/image_analysis/geom_refinement/PostRefine.cpp +++ b/image_analysis/geom_refinement/PostRefine.cpp @@ -13,7 +13,7 @@ #include #include "../../common/JFJochMath.h" // PI -#include "XtalResidual.h" // XtalResidual (the positional detector<->reciprocal residual, step B) +#include "XtalResidual.h" // the positional detector<->reciprocal residual, and the cell it is parameterised by #include "LatticeReduction.h" #include "ceres/ceres.h" #include "ceres/rotation.h" @@ -45,26 +45,36 @@ struct Partial { float obs_x, obs_y; // observed spot centroid (pixels); NAN if the box sum found no centroid }; -// A rocking event and its precomputed reference reciprocal vector (phi=0 frame, from the indexed lattice). +// A rocking event: one reflection's intensity-weighted centroid over the frames it spans. struct Event { double phi_obs; // rad, intensity-weighted rocking centroid double weight; // sqrt(sum I / sum sigma) - double e_ref[3]; // h*a* + k*b* + l*c* at the reference (unrefined) cell/orientation int h, k, l; }; -// Distance-INDEPENDENT Ewald excitation residual for a uniform cell-scale parameter s and a refined -// goniometer axis (3-vector). The header-distance miscalibration leaves a uniform cell scale; the axis is -// the other phi_obs lever. Both are phi_obs-constrained (distance-independent). e_ref is the reference -// reciprocal (h*a* + k*b* + l*c* at the indexed cell). On the Ewald sphere <=> |p|^2 + 2 p_z/lambda == 0. -struct ScaleAxisExcitationResidual { - ScaleAxisExcitationResidual(double lambda, double angle_rad, double weight, const double e_ref[3]) +// Distance-INDEPENDENT Ewald excitation residual with the whole crystal free. The observed rocking +// centroid phi_obs says where the reflection actually crossed the Ewald sphere, which fixes the +// ABSOLUTE size of the reciprocal lattice with the detector never entering - that is what breaks the +// distance <-> cell-scale degeneracy the spot positions alone leave open, and it is why the two can be +// fitted together. The parameter blocks are XtalResidual's crystal half (rotation axis, orientation, +// cell lengths, cell angles), so one problem can share them between the two residuals. +// On the Ewald sphere <=> |p|^2 + 2 p_z/lambda == 0. +struct JointExcitationResidual { + JointExcitationResidual(double lambda, double angle_rad, double weight, + int h, int k, int l, gemmi::CrystalSystem symmetry) : inv_lambda(1.0 / lambda), angle_rad(angle_rad), weight(weight), - ex(e_ref[0]), ey(e_ref[1]), ez(e_ref[2]) {} + h(h), k(k), l(l), symmetry(symmetry) {} template - bool operator()(const T *const s, const T *const axis, T *residual) const { - const T inv_s = T(1) / s[0]; - const T p_ref[3] = {T(ex) * inv_s, T(ey) * inv_s, T(ez) * inv_s}; + bool operator()(const T *const axis, const T *const p0, const T *const p1, const T *const p2, + T *residual) const { + Eigen::Matrix bxc, cxa, axb; + T invV; + XtalResidual::ReciprocalBasis(p1, p2, symmetry, bxc, cxa, axb, invV); + const Eigen::Matrix unrot = (bxc * T(h) + cxa * T(k) + axb * T(l)) * invV; + const T recip_unrot[3] = {unrot[0], unrot[1], unrot[2]}; + T p_ref[3]; + const AngleAxisRotator rot_p0(p0); + rot_p0.Rotate(recip_unrot, p_ref); const T aa[3] = {T(-angle_rad) * axis[0], T(-angle_rad) * axis[1], T(-angle_rad) * axis[2]}; T p_lab[3]; ceres::AngleAxisRotatePoint(aa, p_ref, p_lab); @@ -73,21 +83,23 @@ struct ScaleAxisExcitationResidual { residual[0] = T(weight) * zeta * T(0.5) / T(inv_lambda); return true; } - const double inv_lambda, angle_rad, weight, ex, ey, ez; + const double inv_lambda, angle_rad, weight; + const double h, k, l; + const gemmi::CrystalSystem symmetry; }; -// GONIOMETER ROTATION SCALE k: the same Ewald excitation residual, but with the cell scale and the axis -// DIRECTION already committed by step A, so the single free quantity is how far the stage actually turned -// per unit of commanded angle. Two differences from step A matter: +// GONIOMETER ROTATION SCALE k: the same Ewald excitation residual, but with the crystal and the axis +// DIRECTION already committed by the joint fit, so the single free quantity is how far the stage actually +// turned per unit of commanded angle. Two differences from that fit matter: // * the angle is measured from the CENTRE of the sweep, not from the goniometer's zero. The reference // orientation is the one rotation indexing fitted against the commanded angles, so it has already // absorbed the MEAN angle error; only the part that varies across the sweep is left to fit. Scaling the -// absolute angle instead - which is what reading k off the length of step A's axis vector does - asks +// absolute angle instead - which is what reading k off the length of the fitted axis vector does - asks // the fit to also produce a constant offset it has no parameter for, and the least-squares compromise // shrinks k towards 1 by var(phi) / (var(phi) + phi_centre^2): exactly a factor of four for the common // case of a sweep starting at zero. -// * e_mid is the reference reciprocal vector already turned to the sweep centre and divided by the -// committed cell scale, so nothing but k is free. +// * e_mid is the committed reciprocal vector already turned to the sweep centre, so nothing but k is +// free. struct RotationScaleResidual { RotationScaleResidual(double lambda, double dangle_rad, const double u[3], const double e_mid[3]) : inv_lambda(1.0 / lambda), dangle_rad(dangle_rad), @@ -124,7 +136,6 @@ PostRefineResult PostRefineRotationGeometry(const std::vector=2-frame events carry an - // unbiased phi_obs (a single-frame centroid is just the frame centre); precompute e_ref per event. + // unbiased phi_obs (a single-frame centroid is just the frame centre). constexpr float MAX_FRAME_GAP = 2.0f; const auto run_end = [&](size_t i, size_t end) { size_t j = i + 1; @@ -258,10 +269,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector 0.0 && sumSig > 0.0)) return false; - const Coord e = Astar * static_cast(pts[i].h) + Bstar * static_cast(pts[i].k) - + Cstar * static_cast(pts[i].l); - out = Event{sumIphi / sumI, std::sqrt(sumI / sumSig), {e.x, e.y, e.z}, - pts[i].h, pts[i].k, pts[i].l}; + out = Event{sumIphi / sumI, std::sqrt(sumI / sumSig), pts[i].h, pts[i].k, pts[i].l}; return true; }; // An event never crosses an h boundary - h is the leading sort key - so the buckets can be @@ -325,15 +333,24 @@ PostRefineResult PostRefineRotationGeometry(const std::vectorreciprocal residual at each partial's observed spot, and + // * the distance-INDEPENDENT Ewald excitation residual at each rocking centroid phi_obs. + // It replaces a two-step fit that scaled the whole cell by ONE scalar against phi_obs and then + // read the distance off that scaled cell. The two were separated because the positional residual + // is degenerate with the cell scale - which is true of the positions ALONE, and is exactly what + // the excitation residual breaks, so the degeneracy that motivated the split is already resolved + // inside the same problem. Splitting it cost accuracy twice over: the first pass frees the whole + // lattice against a frozen distance, so the distortion it absorbs is ANISOTROPIC and no single + // scale can undo it; and whatever bias is left in that scale goes straight into the distance, + // which is only ever determined relative to the cell. Measured over three wavelengths of one + // crystal, that scale's bias changed SIGN with the wavelength and the distance error followed it + // with a slope of 1.3 and a correlation of 0.99, while the rocking centroids on their own fixed + // the cell volume to 0.06 %. + // Detector tilt is held fixed (gauge-coupled to the orientation on a single crystal). + // Committed only if it lowers a HELD-OUT (deterministic split-half) residual and the move stays + // within the bounds below - otherwise the geometry is left at nominal ("quit when things go wrong"). if (settings.refine_geometry) { const gemmi::CrystalSystem sys = (settings.crystal_system == gemmi::CrystalSystem::Trigonal) ? gemmi::CrystalSystem::Hexagonal @@ -345,13 +362,16 @@ PostRefineResult PostRefineRotationGeometry(const std::vector(h) * 2654435761u + static_cast(k) * 2246822519u + static_cast(l) * 3266489917u; @@ -362,53 +382,246 @@ PostRefineResult PostRefineRotationGeometry(const std::vector(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); }; - auto solve_scale_axis = [&](Subset s, double &s_out, double ax_out[3]) { - double sc = 1.0, axv[3] = {ax0[0], ax0[1], ax0[2]}; - ceres::Problem p; + 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 (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)); }); + n_obs = MAX_OBS; + } + result.obs_used = static_cast(n_obs); + + // ===== The joint fit ===== + // Cost of a whole geometry over one subset, the two residual families kept apart so the log + // can say which of them moved. Each is a mean per residual VALUE, and the number they are + // summed over below is the same weighting the solver itself applies. + auto joint_cost = [&](Subset s, const double bm[2], const double ds[1], const double rv[3], + const double q0[3], const double q1[3], const double q2[3], + double &pos_out, double &exc_out) { + double cp = 0.0, ce = 0.0; + size_t np = 0, ne = 0; + 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, cos_rot3, sin_rot3, + angle_rad(pp->img), pp->h, pp->k, pp->l, sys, orientation); + double res[3] = {0, 0, 0}; + r(bm, ds, det_rot, rv, q0, q1, q2, res); + cp += res[0] * res[0] + res[1] * res[1] + res[2] * res[2]; + np += 3; + } for (const int32_t i : selected) { const Event &ev = events[i]; if (!in(ev.h, ev.k, ev.l, s)) continue; - p.AddResidualBlock(new ceres::AutoDiffCostFunction( - new ScaleAxisExcitationResidual(lambda_l, ev.phi_obs, settings.excitation_weight, ev.e_ref)), - new ceres::CauchyLoss(0.02), &sc, axv); + JointExcitationResidual r(lambda_l, ev.phi_obs, settings.excitation_weight, + ev.h, ev.k, ev.l, sys); + double res = 0.0; + r(rv, q0, q1, q2, &res); + ce += res * res; + ++ne; } - p.SetParameterLowerBound(&sc, 0, 0.9); p.SetParameterUpperBound(&sc, 0, 1.1); - for (int j = 0; j < 3; ++j) { p.SetParameterLowerBound(axv, j, ax0[j] - 0.05); - p.SetParameterUpperBound(axv, j, ax0[j] + 0.05); } - ceres::Solver::Options o; o.linear_solver_type = ceres::DENSE_QR; o.max_num_iterations = 50; + pos_out = np ? cp / np : 0.0; + exc_out = ne ? ce / ne : 0.0; + return (np + ne) ? (cp + ce) / static_cast(np + ne) : 0.0; + }; + + // One problem, both residual families, every block seeded at nominal. The detector tilt is + // declared and held constant so that a non-zero rot1/rot2 still acts on the observed side. + auto solve_joint = [&](Subset s, double bm[2], double ds[1], double rv[3], + double q0[3], double q1[3], double q2[3]) { + bm[0] = beam_x0; bm[1] = beam_y0; ds[0] = dist0; + for (int j = 0; j < 3; ++j) { + rv[j] = ax0[j]; q0[j] = p0_0[j]; q1[j] = p1_0[j]; q2[j] = p2_0[j]; + } + ceres::Problem p; + for (size_t oi = 0; oi < n_obs; ++oi) { + const Partial *pp = obs[oi]; + if (!in(pp->h, pp->k, pp->l, s)) continue; + p.AddResidualBlock(new ceres::AutoDiffCostFunction( + new XtalResidual(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, cos_rot3, sin_rot3, + angle_rad(pp->img), pp->h, pp->k, pp->l, sys, orientation)), + new ceres::CauchyLoss(0.02), bm, ds, det_rot, rv, q0, q1, q2); + } + for (const int32_t i : selected) { + const Event &ev = events[i]; + if (!in(ev.h, ev.k, ev.l, s)) continue; + p.AddResidualBlock(new ceres::AutoDiffCostFunction( + new JointExcitationResidual(lambda_l, ev.phi_obs, settings.excitation_weight, + ev.h, ev.k, ev.l, sys)), + new ceres::CauchyLoss(0.02), rv, q0, q1, q2); + } + if (p.NumResidualBlocks() == 0) return false; + p.SetParameterBlockConstant(det_rot); + if (!p2_free) p.SetParameterBlockConstant(q2); + p.SetParameterLowerBound(ds, 0, dist0 * 0.95); p.SetParameterUpperBound(ds, 0, dist0 * 1.05); + // The box has to reach everywhere the gate below would accept, or the gate is never the + // thing that decides: a fit pinned at a box face lands exactly ON the bound and is then + // refused for being there. + for (int j = 0; j < 2; ++j) { + const double m = settings.measured_beam_px ? (*settings.measured_beam_px)[j] : bm[j]; + p.SetParameterLowerBound(bm, j, std::min(bm[j], m) - BEAM_BOUND_PXL); + p.SetParameterUpperBound(bm, j, std::max(bm[j], m) + BEAM_BOUND_PXL); + } + for (int j = 0; j < 3; ++j) { + p.SetParameterLowerBound(rv, j, ax0[j] - 0.05); + p.SetParameterUpperBound(rv, j, ax0[j] + 0.05); + p.SetParameterLowerBound(q1, j, 0.97 * p1_0[j]); + p.SetParameterUpperBound(q1, j, 1.03 * p1_0[j]); + } + if (p2_free) + for (int j = 0; j < 3; ++j) { + p.SetParameterLowerBound(q2, j, p2_0[j] - 0.05); + p.SetParameterUpperBound(q2, j, p2_0[j] + 0.05); + } + ceres::Solver::Options o; o.linear_solver_type = ceres::DENSE_QR; o.max_num_iterations = 60; o.num_threads = std::max(1, settings.num_threads); o.logging_type = ceres::LoggingType::SILENT; ceres::Solver::Summary sum; ceres::Solve(o, &p, &sum); - s_out = sc; ax_out[0] = axv[0]; ax_out[1] = axv[1]; ax_out[2] = axv[2]; return sum.IsSolutionUsable(); }; - double s_fit = 1.0, ax_fit[3]; - const bool convA = solve_scale_axis(FIT, s_fit, ax_fit); - const double cvA_nom = excit_cost(VAL, 1.0, ax0); - const double cvA_ref = excit_cost(VAL, s_fit, ax_fit); - // Commit the cell scale only for a small, credible move: a well-calibrated header needs < ~0.6 %, - // so a > 1 % scale is a red flag (on multi-lattice / noisy data the excitation fit is biased the - // same way in every cross-validation fold, so the relative-improvement gate cannot catch it). - result.cell_refined = convA && cvA_ref < 0.98 * cvA_nom && std::fabs(s_fit - 1.0) < 0.01; - double s = 1.0, axv[3] = {ax0[0], ax0[1], ax0[2]}; - if (result.cell_refined) solve_scale_axis(ALL, s, axv); // commit: re-fit on all data + + double beam[2] = {beam_x0, beam_y0}, dist[1] = {dist0}; + double axv[3] = {ax0[0], ax0[1], ax0[2]}; + bool commit = false; + if (n_obs >= static_cast(settings.min_events)) { + double pos_nom = 0.0, exc_nom = 0.0, pos_ref = 0.0, exc_ref = 0.0; + const double cv_nom = joint_cost(VAL, beam, dist, axv, p0, p1, p2, pos_nom, exc_nom); + double bm_f[2], ds_f[1], rv_f[3], q0_f[3], q1_f[3], q2_f[3]; + const bool convJ = solve_joint(FIT, bm_f, ds_f, rv_f, q0_f, q1_f, q2_f); + const double cv_ref = joint_cost(VAL, bm_f, ds_f, rv_f, q0_f, q1_f, q2_f, pos_ref, exc_ref); + // Commit only for a small, credible move: distance < 1 % (a calibrated header needs + // < ~0.6 %), each cell length < 1 %, and the beam inside the bound measured from + // whichever centre anything already believes is nearer. A larger move is the red flag for + // an unreliable fit - typically a second lattice whose spots bias every cross-validation + // fold identically, so the relative "it improved" gate is blind to it. The absolute size + // of the move discriminates a genuine header correction from that failure far better than + // the absolute residual, which real marginal (noisy / iced) data shares with it. + const double from_nominal = std::hypot(bm_f[0] - beam_x0, bm_f[1] - beam_y0); + const double from_measured = settings.measured_beam_px + ? std::hypot(bm_f[0] - (*settings.measured_beam_px)[0], + bm_f[1] - (*settings.measured_beam_px)[1]) + : std::numeric_limits::infinity(); + double len_shift = 0.0; + for (int j = 0; j < 3; ++j) + len_shift = std::max(len_shift, std::fabs(q1_f[j] - p1_0[j]) / std::max(1e-9, p1_0[j])); + const bool in_bounds = std::fabs(ds_f[0] - dist0) < 0.01 * dist0 + && std::min(from_nominal, from_measured) < BEAM_BOUND_PXL + && len_shift < 0.01; + commit = convJ && cv_ref < 0.98 * cv_nom && in_bounds; + if (commit) solve_joint(ALL, beam, dist, axv, p0, p1, p2); // commit: re-fit on all data + // Name which test refused it. Four different things reject here and the geometry that + // comes out is the same in all four, so a run that silently keeps its header geometry + // says nothing about whether the fit was bad, the improvement too small, or the move too + // large for the bound - which is the one case where the number worth reading is the one + // that was thrown away. + const char *verdict = + commit ? "COMMIT" + : !convJ ? "reject (the fit did not converge)" + : !(cv_ref < 0.98 * cv_nom) ? "reject (the held-out residual did not improve enough)" + : std::fabs(ds_f[0] - dist0) >= 0.01 * dist0 + ? "reject (the distance moved more than 1 %)" + : len_shift >= 0.01 ? "reject (a cell length moved more than 1 %)" + : "reject (the beam moved further than the bound from " + "every centre anything believes)"; + // The cell the log names is the EFFECTIVE one - what the residual's B matrix builds + // from the blocks - not the blocks themselves, whose unused components a high-symmetry + // system leaves at whatever the seed happened to put there. + double len_nom[3], ang_nom[3], len_fit[3], ang_fit[3]; + EffectiveCellFromParams(sys, p1_0, p2_0, len_nom, ang_nom); + EffectiveCellFromParams(sys, q1_f, q2_f, len_fit, ang_fit); + logger.Info("Post-refine GEOM (joint crystal + detector): dist {:.3f} -> {:.3f} mm, beam " + "({:.2f},{:.2f}) -> ({:.2f},{:.2f}), cell {:.3f} {:.3f} {:.3f} -> " + "{:.3f} {:.3f} {:.3f}, held-out positional {:.3e} -> {:.3e}, excitation " + "{:.3e} -> {:.3e} => {}", + dist0, ds_f[0], beam_x0, beam_y0, bm_f[0], bm_f[1], + len_nom[0], len_nom[1], len_nom[2], len_fit[0], len_fit[1], len_fit[2], + pos_nom, pos_ref, exc_nom, exc_ref, verdict); + } else { + logger.Info("Post-refine GEOM: only {} positional observations - the joint fit needs the " + "spot positions as well as the rocking angles, so nothing is refined", n_obs); + } + result.cell_refined = commit; + result.detector_refined = commit; + const double axlen = std::sqrt(axv[0]*axv[0] + axv[1]*axv[1] + axv[2]*axv[2]); const double axdev = std::acos(std::clamp((axv[0]*ax0[0]+axv[1]*ax0[1]+axv[2]*ax0[2]) / std::max(1e-9, axlen), -1.0, 1.0)) * 180.0 / PI; - logger.Info("Post-refine GEOM step A (cell/axis): s = {:.5f}, rot-axis {:.3f} deg, held-out excit " - "{:.3e} -> {:.3e} => {}", s, axdev, cvA_nom, cvA_ref, - result.cell_refined ? "COMMIT" : "reject (kept nominal cell)"); + // The committed crystal, as a lattice again. Where nothing was committed this is the lattice + // rotation indexing handed in, unchanged - not its round trip through the symmetry-constrained + // parameterisation, which would move the cell for a fit that was refused. + double eff_len[3], eff_ang[3]; + EffectiveCellFromParams(sys, p1, p2, eff_len, eff_ang); + const CrystalLattice committed_latt = commit + ? AngleAxisAndCellToLattice(p0, eff_len, eff_ang[0], eff_ang[1], eff_ang[2]) + : reference_latt; + const Coord As = committed_latt.Astar(), Bs = committed_latt.Bstar(), Cs = committed_latt.Cstar(); + if (commit) { + result.cell = UnitCell{static_cast(eff_len[0]), static_cast(eff_len[1]), + static_cast(eff_len[2]), + static_cast(eff_ang[0] * 180.0 / PI), + static_cast(eff_ang[1] * 180.0 / PI), + static_cast(eff_ang[2] * 180.0 / PI)}; + logger.Info("Post-refine GEOM: committed cell {:.3f} {:.3f} {:.3f} {:.2f} {:.2f} {:.2f}, " + "rotation axis moved {:.3f} deg (the cell is the fit's own; the second pass " + "re-indexes at the refined detector geometry)", result.cell.a, result.cell.b, + result.cell.c, result.cell.alpha, result.cell.beta, result.cell.gamma, axdev); + } // ===== Goniometer rotation SCALE k, its own one-parameter fit on the same rocking events ===== // The angles stored in the file are the COMMANDED ones, so a stage that turned k times as far @@ -425,7 +638,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector(n_events); const double sweep_deg = (phi_hi - phi_lo) * 180.0 / PI; - // The reference reciprocal vector turned to the sweep centre, at the committed cell scale. The + // The committed reciprocal vector turned to the sweep centre. The // angle then enters the fit measured FROM that centre. A constant crystal missetting about the // spindle is k with a slope in phi, so measuring the angle from the goniometer's zero instead // lets a missetting leak into k with gain / - which depends only on where the sweep @@ -450,8 +663,10 @@ PostRefineResult PostRefineRotationGeometry(const std::vector(n_events), nthreads, [&](int lo, int hi) { for (int e = lo; e < hi; ++e) { - const double p[3] = {events[e].e_ref[0] / s, events[e].e_ref[1] / s, - events[e].e_ref[2] / s}; + const Coord ec = As * static_cast(events[e].h) + + Bs * static_cast(events[e].k) + + Cs * static_cast(events[e].l); + const double p[3] = {ec.x, ec.y, ec.z}; double em[3]; ceres::AngleAxisRotatePoint(aa_c, p, em); const double ue = u[0] * em[0] + u[1] * em[1] + u[2] * em[2]; @@ -609,177 +824,15 @@ PostRefineResult PostRefineRotationGeometry(const std::vector 1.0 ? "further" : "less far"); - // Cell (scale s, shape fixed) as the XtalResidual parameter blocks p0/p1/p2, held CONSTANT in step B. - double p0[3] = {0, 0, 0}, p1[3] = {0, 0, 0}, p2[3] = {0, 0, 0}; - double beta = r0.beta; - switch (sys) { - case gemmi::CrystalSystem::Tetragonal: - LatticeToRodriguesAndLengths_GS(reference_latt, p0, p1); - p1[0] = (p1[0] + p1[1]) / 2.0; break; - case gemmi::CrystalSystem::Cubic: - LatticeToRodriguesAndLengths_GS(reference_latt, p0, p1); - p1[0] = (p1[0] + p1[1] + p1[2]) / 3.0; break; - case gemmi::CrystalSystem::Hexagonal: - LatticeToRodriguesAndLengths_Hex(reference_latt, p0, p1); break; - case gemmi::CrystalSystem::Monoclinic: - LatticeToRodriguesLengthsBeta_Mono(reference_latt, p0, p1, beta); - p2[0] = beta; break; - case gemmi::CrystalSystem::Orthorhombic: - LatticeToRodriguesAndLengths_GS(reference_latt, p0, p1); break; - default: - LatticeToRodriguesAndLengths_GS(reference_latt, p0, p1); - p2[0] = r0.alpha * PI / 180.0; p2[1] = r0.beta * PI / 180.0; p2[2] = r0.gamma * PI / 180.0; break; - } - for (int j = 0; j < 3; ++j) p1[j] *= s; // apply the committed cell scale - double rot_vec[3] = {axv[0], axv[1], axv[2]}; // committed (or nominal) axis - - // ===== Step B: detector distance + beam from the observed positions, cell fixed ===== - // 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 (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)); }); - n_obs = MAX_OBS; - } - 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 (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, cos_rot3, sin_rot3, - angle_rad(pp->img), - pp->h, pp->k, pp->l, sys, orientation); - double resid[3] = {0, 0, 0}; - r(beam, dist, det_rot, rot_vec, p0, p1, p2, resid); - c += resid[0]*resid[0] + resid[1]*resid[1] + resid[2]*resid[2]; ++n; - } - return n ? c / n : 0.0; - }; - auto solve_detector = [&](Subset s, double beam_out[2], double &dist_out) { - double beam[2] = {beam_x0, beam_y0}, dist[1] = {dist0}; - // What the reduced residual holds fixed varies with the frame alone, so build one set of - // constants per image rather than one per observation. A map node keeps its address as - // the map grows, and the map outlives the problem that points into it. - std::map frame_const; - ceres::Problem p; - for (size_t oi = 0; oi < n_obs; ++oi) { - const Partial *pp = obs[oi]; - if (!in(pp->h, pp->k, pp->l, s)) continue; - const XtalFrameConstants &fc = frame_const.try_emplace( - pp->img, det_rot, rot_vec, angle_rad(pp->img), p1, p2, sys).first->second; - p.AddResidualBlock(new ceres::AutoDiffCostFunction( - new XtalResidualBeamDistance( - XtalResidual(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, cos_rot3, sin_rot3, - angle_rad(pp->img), pp->h, pp->k, pp->l, sys, orientation), - fc, p0)), - 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, 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); - // The box has to reach everywhere the gate below would accept, or the gate is never the - // thing that decides: a fit pinned at a box face lands exactly ON the bound and is then - // refused for being there. - for (int j = 0; j < 2; ++j) { - const double m = settings.measured_beam_px ? (*settings.measured_beam_px)[j] : beam[j]; - p.SetParameterLowerBound(beam, j, std::min(beam[j], m) - BEAM_BOUND_PXL); - p.SetParameterUpperBound(beam, j, std::max(beam[j], m) + BEAM_BOUND_PXL); - } - ceres::Solver::Options o; o.linear_solver_type = ceres::DENSE_QR; o.max_num_iterations = 60; - o.num_threads = std::max(1, settings.num_threads); o.logging_type = ceres::LoggingType::SILENT; - ceres::Solver::Summary sum; ceres::Solve(o, &p, &sum); - beam_out[0] = beam[0]; beam_out[1] = beam[1]; dist_out = dist[0]; - return sum.IsSolutionUsable(); - }; - double beam[2] = {beam_x0, beam_y0}, dist = dist0; - 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}; - const double b_ref[2] = {beam_fit[0], beam_fit[1]}, d_ref[1] = {dist_fit}; - const double cvB_nom = pos_cost(VAL, b_nom, d_nom); - const double cvB_ref = pos_cost(VAL, b_ref, d_ref); - // Commit the detector geometry only for a small, credible move: distance < 1 % (a calibrated - // header needs < ~0.6 %). A larger move is the red flag for an unreliable fit - typically a - // second lattice whose spots bias every cross-validation fold identically, so the relative - // "it improved" gate is blind to it and pulls a spurious distance<->cell pair (the radial - // degeneracy) far off. The absolute size of the move discriminates a genuine header correction - // from that failure far better than the absolute residual, which real marginal (noisy / iced) - // data shares with the multi-lattice case. - const double from_nominal = std::hypot(beam_fit[0] - beam_x0, beam_fit[1] - beam_y0); - const double from_measured = settings.measured_beam_px - ? std::hypot(beam_fit[0] - (*settings.measured_beam_px)[0], - beam_fit[1] - (*settings.measured_beam_px)[1]) - : std::numeric_limits::infinity(); - const bool in_bounds = std::fabs(dist_fit - dist0) < 0.01 * dist0 - && std::min(from_nominal, from_measured) < BEAM_BOUND_PXL; - result.detector_refined = convB && cvB_ref < 0.98 * cvB_nom && in_bounds; - if (result.detector_refined) { double bo[2]; solve_detector(ALL, bo, dist); beam[0] = bo[0]; beam[1] = bo[1]; } - // Name which test refused it. Three different things reject here and the geometry that - // comes out is the same in all three, so a run that silently keeps its header geometry - // says nothing about whether the fit was bad, the improvement too small, or the move - // too large for the bound - which is the one case where the number worth reading is the - // one that was thrown away. - const char *verdict = - result.detector_refined ? "COMMIT" - : !convB ? "reject (the fit did not converge; kept nominal detector)" - : !(cvB_ref < 0.98 * cvB_nom) ? "reject (the held-out residual did not improve enough; " - "kept nominal detector)" - : std::fabs(dist_fit - dist0) >= 0.01 * dist0 - ? "reject (the distance moved more than 1 %; kept nominal detector)" - : "reject (the beam moved further than the bound from every " - "centre anything believes; kept nominal detector)"; - logger.Info("Post-refine GEOM step B (distance/beam): dist {:.3f} -> {:.3f} mm, beam " - "({:.2f},{:.2f}) -> ({:.2f},{:.2f}), held-out pos {:.3e} -> {:.3e} => {}", - dist0, result.detector_refined ? dist : dist0, beam_x0, beam_y0, - result.detector_refined ? beam[0] : beam_x0, result.detector_refined ? beam[1] : beam_y0, - cvB_nom, cvB_ref, verdict); - } else { - logger.Info("Post-refine GEOM step B: only {} positional observations - skipped", n_obs); - } - // Assemble the committed geometry. - UnitCell cellA = r0; - if (result.cell_refined) { cellA.a = static_cast(r0.a * s); cellA.b = static_cast(r0.b * s); - cellA.c = static_cast(r0.c * s); } - result.cell = cellA; - result.distance_after_mm = dist; + result.distance_after_mm = dist[0]; result.beam_x_before_px = beam_x0; result.beam_x_after_px = beam[0]; result.beam_y_before_px = beam_y0; result.beam_y_after_px = beam[1]; result.events_used = static_cast(selected.size()); - result.ok = result.cell_refined || result.detector_refined; + result.ok = commit; if (!result.ok) - logger.Info("Post-refine GEOM: neither step passed cross-validation - geometry left at nominal"); + logger.Info("Post-refine GEOM: the joint fit did not pass cross-validation - geometry " + "left at nominal"); return result; } diff --git a/image_analysis/geom_refinement/PostRefine.h b/image_analysis/geom_refinement/PostRefine.h index ab29c6c6c..b3cf0a707 100644 --- a/image_analysis/geom_refinement/PostRefine.h +++ b/image_analysis/geom_refinement/PostRefine.h @@ -17,13 +17,15 @@ // Post-integration geometry refinement for rotation data. Unlike the at-indexing XtalOptimizer, this runs // AFTER integration/merge, where each reflection has an OBSERVED rocking centroid phi_obs (the intensity- -// weighted mean goniometer angle over the frames it spans) and an observed spot position. It refines one -// shared crystal orientation + cell (+ optionally the detector distance) against two residuals: +// weighted mean goniometer angle over the frames it spans) and an observed spot position. One JOINT fit +// refines the crystal (orientation, cell, rotation axis) and the detector (distance, beam centre) together +// against two residuals: // * an Ewald excitation residual evaluated at phi_obs (distance-independent) -> pins the absolute cell // scale that the positional residual leaves degenerate with the distance. Because phi_obs is the real // rocking angle (not a frame centre) it is unbiased. // * the positional detector<->reciprocal residual at each partial's observed spot -> pins the distance. -// Reflections are weighted by their merged I/sigma (strong, well-measured reflections dominate). +// Where there are more reflections than the fit's caps, the strongest by I/sigma are the ones kept; +// inside the fit every reflection carries the same weight. struct PostRefineResult { bool ok = false; @@ -34,11 +36,13 @@ struct PostRefineResult { double distance_before_mm = 0.0, distance_after_mm = 0.0; double beam_x_before_px = 0.0, beam_x_after_px = 0.0; // refined beam centre (GEOM mode) double beam_y_before_px = 0.0, beam_y_after_px = 0.0; - bool cell_refined = false; // GEOM step A (cell scale + axis) passed cross-validation - bool detector_refined = false; // GEOM step B (distance + beam) passed cross-validation + // The joint fit is one decision - the crystal and the detector are refined together and committed + // together - so these two are always equal. Both are kept because the report names them separately. + bool cell_refined = false; // the joint fit passed cross-validation (crystal half) + bool detector_refined = false; // the joint fit passed cross-validation (detector half) // GONIOMETER ROTATION SCALE: the factor by which the stage actually turned relative to the angle // stored in the file (which is the COMMANDED value, hence a stage calibration error is invisible in - // the header). Fitted after step A as a single free parameter, with the cell scale and the axis + // the header). Fitted after the joint fit as a single free parameter, with the crystal and the axis // direction held at their committed values. Always the fitted value; 1.0 = header and stage agree. double rotation_scale = 1.0; // Whether the fit passed every test needed to ACT on it: enough sweep and events, a significant and @@ -48,14 +52,14 @@ struct PostRefineResult { struct PostRefineSettings { gemmi::CrystalSystem crystal_system = gemmi::CrystalSystem::Triclinic; - bool refine_geometry = false; // XtalOptimizer-equivalent: cell scale + axis (from phi_obs) and detector - // distance + beam centre (from the observed spot positions X,Y), as two - // separate cross-validated steps. The only supported refinement mode. + bool refine_geometry = false; // XtalOptimizer-equivalent: the crystal (from phi_obs and the observed + // spot positions) and the detector distance + beam centre, in one + // cross-validated fit. The only supported refinement mode. double excitation_weight = 1.0; // weight of the phi/excitation residual vs the positional one int min_events = 50; int num_threads = 1; // An independent measurement of the beam centre in pixels, where the run has one (rugnux's pre-scan - // fit to the isotropy of the scattered background, which runs on every rotation run). Step B bounds + // fit to the isotropy of the scattered background, which runs on every rotation run). The fit bounds // how far it may move the beam from whichever of this and the nominal centre is NEARER; see // BEAM_BOUND_PXL in PostRefine.cpp for why the nominal centre alone is not enough to bound it. std::optional> measured_beam_px; diff --git a/image_analysis/geom_refinement/XtalResidual.h b/image_analysis/geom_refinement/XtalResidual.h index 3d206210f..9abcda661 100644 --- a/image_analysis/geom_refinement/XtalResidual.h +++ b/image_analysis/geom_refinement/XtalResidual.h @@ -76,6 +76,39 @@ struct AngleAxisRotator { bool at_zero; }; +// The cell that XtalResidual::ReciprocalBasis builds out of the parameter blocks p1 (lengths) and p2 +// (angles, radians): the components a system does not read are filled in from its symmetry, so a caller +// holding a refined parameter set can turn it back into a cell or a lattice. The inverse of what that +// function reads, and kept beside it so the two conventions cannot drift apart. +inline void EffectiveCellFromParams(gemmi::CrystalSystem symmetry, const double p1[3], const double p2[3], + double lengths[3], double angles_rad[3]) { + const double right = M_PI / 2.0; + lengths[0] = p1[0]; + lengths[1] = p1[1]; + lengths[2] = p1[2]; + angles_rad[0] = right; angles_rad[1] = right; angles_rad[2] = right; + switch (symmetry) { + case gemmi::CrystalSystem::Hexagonal: + lengths[1] = p1[0]; + angles_rad[2] = 2.0 * M_PI / 3.0; + break; + case gemmi::CrystalSystem::Tetragonal: + lengths[1] = p1[0]; + break; + case gemmi::CrystalSystem::Cubic: + lengths[1] = p1[0]; lengths[2] = p1[0]; + break; + case gemmi::CrystalSystem::Monoclinic: + angles_rad[1] = p2[0]; + break; + case gemmi::CrystalSystem::Triclinic: + angles_rad[0] = p2[0]; angles_rad[1] = p2[1]; angles_rad[2] = p2[2]; + break; + default: // orthorhombic + break; + } +} + // 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), diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 624783fac..c77dffd06 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -231,9 +231,9 @@ std::string RenderResultReport(const std::string &output_prefix, if (result.post_refine.has_value()) { const auto &pr = *result.post_refine; Section(os, "3. GEOMETRY POST-REFINEMENT"); - os << " The rotation two-pass fits the detector distance and beam centre from the observed spot\n" - << " positions, and the cell scale and rotation axis from the observed rocking angles. Each\n" - << " step is committed only if it improves a held-out residual.\n\n"; + os << " The rotation two-pass fits the crystal (orientation, cell, rotation axis) and the detector\n" + << " (distance, beam centre) together, against the observed spot positions and the observed\n" + << " rocking angles at once. It is committed only if it improves a held-out residual.\n\n"; Key(os, "POSTREFINE_EVENTS_USED", pr.events_used); Key(os, "POSTREFINE_OBS_USED", pr.obs_used); Key(os, "POSTREFINE_CELL_COMMITTED", pr.cell_refined ? "TRUE" : "FALSE");