PixelRefine: env-gated orientation + cell-scale sweep (PR_SWEEP) for R-free A/B

R-free validation (full jet, 1.5 A) confirmed the r1_multiplier fix: PR x6 beats
traditional (R-free 0.2625 vs 0.2802), the multiplier optimum is ~6 (x6==x9 on
R-free; x9 only buys CC1/2 internal consistency), and per-image orientation
*refinement* is a no-op (0.2618 vs 0.2625).

Re-adds the reference-driven orientation + uniform cell-scale SWEEP behind PR_SWEEP
(off by default) to test whether the cell-scale degree of freedom - which the
gradient orientation refinement lacks, and which moves the high-res shells radially
- helps R-free at 1.5 A on serial data. CCref is a near-no-op on the jet (as for
the gradient path), but that does not certify R-free, so it is left for validation.

NOTE: PR_RMULT / PR_ORIENT / PR_SWEEP remain temporary diagnostic env knobs; once
the sweep R-free is in, bake r1_multiplier=6, drop the no-op paths, strip the knobs.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

image_analysis/pixel_refinement/PixelRefine.cpp

@@ -457,6 +457,144 @@ void PixelRefine::BuildParameterBlocks(const PixelRefineData &data,
     }
 }
 
+// Optional pre-pass (env-gated): a small GLOBAL orientation + uniform cell-scale sweep that
+// maximises CC of the box-summed intensities against the reference. Unlike the per-pixel
+// orientation refinement it also adjusts a per-image cell scale (a radial degree of freedom),
+// and makes coarse global moves the local gradient cannot. Coordinate descent over the three
+// Rodrigues axes + cell scale within geometry-derived pixel bounds (highest-resolution spot
+// moves ~1 px/step, lowest barely moves). Writes the best orientation/cell into data.latt.
+template<class T>
+void PixelRefine::SweepOrientationCell(const T *image, BraggPrediction &prediction,
+                                       PixelRefineData &data) const {
+    const int radius = data.shoebox_radius;
+    const double beam_x = data.geom.GetBeamX_pxl();
+    const double beam_y = data.geom.GetBeamY_pxl();
+    const auto qnan = std::numeric_limits<double>::quiet_NaN();
+
+    // Box-sum minus local (perimeter) background MEAN, raw counts; NaN off-detector/masked.
+    auto integrate = [&](double px, double py) -> double {
+        const int cx = static_cast<int>(std::lround(px));
+        const int cy = static_cast<int>(std::lround(py));
+        const int outer = radius + 1;
+        if (cx - outer < 0 || cy - outer < 0 ||
+            cx + outer >= static_cast<int>(xpixel) || cy + outer >= static_cast<int>(ypixel))
+            return qnan;
+        double sig = 0.0;
+        int nsig = 0;
+        std::vector<double> ring;
+        ring.reserve((2 * outer + 1) * (2 * outer + 1));
+        for (int y = cy - outer; y <= cy + outer; ++y) {
+            for (int x = cx - outer; x <= cx + outer; ++x) {
+                const T raw = image[static_cast<size_t>(xpixel) * y + x];
+                if (raw == std::numeric_limits<T>::max())
+                    return qnan;
+                if (std::is_signed_v<T> && raw == std::numeric_limits<T>::min())
+                    return qnan;
+                const double v = static_cast<double>(raw);
+                if (std::abs(x - cx) <= radius && std::abs(y - cy) <= radius) {
+                    sig += v; ++nsig;
+                } else {
+                    ring.push_back(v);
+                }
+            }
+        }
+        if (ring.size() < 5)
+            return qnan;
+        double rsum = 0.0;
+        for (const double v : ring)
+            rsum += v;
+        return sig - nsig * (rsum / static_cast<double>(ring.size()));
+    };
+
+    DiffractionExperiment exp_iter = experiment;
+    exp_iter.BeamX_pxl(beam_x).BeamY_pxl(beam_y)
+            .DetectorDistance_mm(data.geom.GetDetectorDistance_mm())
+            .PoniRot1_rad(data.geom.GetPoniRot1_rad())
+            .PoniRot2_rad(data.geom.GetPoniRot2_rad());
+    const BraggPredictionSettings settings{
+        .high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(),
+        .ewald_dist_cutoff = static_cast<float>(data.ewald_dist_cutoff),
+        .max_hkl = 100,
+        .centering = data.centering,
+        .bandwidth_sigma = static_cast<float>(data.bandwidth)
+    };
+    const int nrefl = prediction.Calc(exp_iter, data.latt, settings);
+    const auto &predicted = prediction.GetReflections();
+
+    struct Matched { int h, k, l; double refI; };
+    std::vector<Matched> matched;
+    double r_max = 0.0, r_min = std::numeric_limits<double>::max();
+    for (int i = 0; i < nrefl; ++i) {
+        const auto &r = predicted[i];
+        const auto it = reference_data.find(hkl_key_generator(r));
+        if (it == reference_data.end())
+            continue;
+        matched.push_back({r.h, r.k, r.l, it->second});
+        const double dx = r.predicted_x - beam_x, dy = r.predicted_y - beam_y;
+        const double rad = std::sqrt(dx * dx + dy * dy);
+        r_max = std::max(r_max, rad);
+        r_min = std::min(r_min, rad);
+    }
+    if (matched.size() < 20 || r_min <= 1.0 || r_max <= r_min)
+        return; // too little to anchor a meaningful sweep
+
+    auto score = [&](const CrystalLattice &L) -> double {
+        const Coord A = L.Astar(), B = L.Bstar(), C = L.Cstar();
+        double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0;
+        int n = 0;
+        for (const auto &m : matched) {
+            const Coord g = A * static_cast<float>(m.h) + B * static_cast<float>(m.k)
+                            + C * static_cast<float>(m.l);
+            const auto [x, y] = data.geom.RecipToDetector(g);
+            if (!std::isfinite(x) || !std::isfinite(y))
+                continue;
+            const double I = integrate(x, y);
+            if (!std::isfinite(I))
+                continue;
+            sx += I; sy += m.refI; sxx += I * I; syy += m.refI * m.refI; sxy += I * m.refI; ++n;
+        }
+        if (n < 10)
+            return -2.0;
+        const double nd = n;
+        const double cov = sxy - sx * sy / nd, vx = sxx - sx * sx / nd, vy = syy - sy * sy / nd;
+        return (vx > 0.0 && vy > 0.0) ? cov / std::sqrt(vx * vy) : -2.0;
+    };
+
+    const double step = 1.0 / r_max;
+    const int n_rot = std::clamp(
+        static_cast<int>(std::lround(data.sweep_max_deg * M_PI / 180.0 * r_max)), 1, 25);
+    const int n_scale = std::clamp(
+        static_cast<int>(std::lround(data.sweep_max_cell_frac * r_max)), 1, 25);
+    const Coord axes[3] = {Coord(1, 0, 0), Coord(0, 1, 0), Coord(0, 0, 1)};
+
+    CrystalLattice best = data.latt;
+    double best_cc = score(best);
+    for (int round = 0; round < 2; ++round) {
+        for (const auto &axis : axes) {
+            CrystalLattice axis_best = best;
+            double axis_cc = best_cc;
+            for (int i = -n_rot; i <= n_rot; ++i) {
+                if (i == 0) continue;
+                CrystalLattice cand = best.Multiply(RotMatrix(static_cast<float>(i * step), axis));
+                const double cc = score(cand);
+                if (cc > axis_cc) { axis_cc = cc; axis_best = cand; }
+            }
+            best = axis_best; best_cc = axis_cc;
+        }
+        CrystalLattice scale_best = best;
+        double scale_cc = best_cc;
+        for (int i = -n_scale; i <= n_scale; ++i) {
+            if (i == 0) continue;
+            const double s = 1.0 / (1.0 + i * step);
+            CrystalLattice cand = best.Multiply(gemmi::Mat33(s, 0, 0, 0, s, 0, 0, 0, s));
+            const double cc = score(cand);
+            if (cc > scale_cc) { scale_cc = cc; scale_best = cand; }
+        }
+        best = scale_best; best_cc = scale_cc;
+    }
+    data.latt = best;
+}
+
 template<class T>
 void PixelRefine::Run(const T *image,
                       BraggPrediction &prediction,
@@ -464,6 +602,10 @@ void PixelRefine::Run(const T *image,
     data.solved = false;
     data.reflections.clear();
 
+    // Optional reference-driven orientation + cell-scale sweep before prediction (env-gated).
+    if (data.sweep_orientation)
+        SweepOrientationCell(image, prediction, data);
+
     const double lambda = data.geom.GetWavelength_A();
     const double pixel_size = data.geom.GetPixelSize_mm();