From 45ee8c2b40204ceaa4126f004e5024830b17e01b Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sun, 14 Jun 2026 19:59:50 +0200 Subject: [PATCH] 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 --- image_analysis/IndexAndRefine.cpp | 1 + .../pixel_refinement/PixelRefine.cpp | 142 ++++++++++++++++++ image_analysis/pixel_refinement/PixelRefine.h | 12 ++ 3 files changed, 155 insertions(+) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 34244a4a..a37a6962 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -450,6 +450,7 @@ bool IndexAndRefine::PixelRefineIntegrate(DataMessage &msg, prd.shoebox_radius = static_cast(std::lround(experiment.GetBraggIntegrationSettings().GetR1())); if (const char *m = std::getenv("PR_RMULT")) prd.r1_multiplier = std::stod(m); // TEMP: Term-2 R1 multiplier sweep if (std::getenv("PR_ORIENT")) prd.refine_orientation = true; // A/B: per-image orientation refinement + if (std::getenv("PR_SWEEP")) prd.sweep_orientation = true; // A/B: orientation + cell-scale sweep std::vector buffer; const uint8_t *ptr = image.GetUncompressedPtr(buffer); diff --git a/image_analysis/pixel_refinement/PixelRefine.cpp b/image_analysis/pixel_refinement/PixelRefine.cpp index 843d6360..b2f2161f 100644 --- a/image_analysis/pixel_refinement/PixelRefine.cpp +++ b/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 +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::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(std::lround(px)); + const int cy = static_cast(std::lround(py)); + const int outer = radius + 1; + if (cx - outer < 0 || cy - outer < 0 || + cx + outer >= static_cast(xpixel) || cy + outer >= static_cast(ypixel)) + return qnan; + double sig = 0.0; + int nsig = 0; + std::vector 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(xpixel) * y + x]; + if (raw == std::numeric_limits::max()) + return qnan; + if (std::is_signed_v && raw == std::numeric_limits::min()) + return qnan; + const double v = static_cast(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(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(data.ewald_dist_cutoff), + .max_hkl = 100, + .centering = data.centering, + .bandwidth_sigma = static_cast(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; + double r_max = 0.0, r_min = std::numeric_limits::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(m.h) + B * static_cast(m.k) + + C * static_cast(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(std::lround(data.sweep_max_deg * M_PI / 180.0 * r_max)), 1, 25); + const int n_scale = std::clamp( + static_cast(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(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 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(); diff --git a/image_analysis/pixel_refinement/PixelRefine.h b/image_analysis/pixel_refinement/PixelRefine.h index 1e7a2dac..40dba239 100644 --- a/image_analysis/pixel_refinement/PixelRefine.h +++ b/image_analysis/pixel_refinement/PixelRefine.h @@ -83,6 +83,13 @@ struct PixelRefineData { double orient_reg_sigma_deg = 1.0; // orientation-prior strength (deg) double fit_signal_sigma_pix = 1.5; // signal-weighting sigma for the orientation fit (px) + // Optional reference-driven orientation + per-image cell-scale sweep before prediction + // (env-gated, off by default). Coarse global moves the gradient refinement cannot make, + // plus a radial (cell-scale) DOF the gradient path lacks. See SweepOrientationCell. + bool sweep_orientation = false; + double sweep_max_deg = 0.15; + double sweep_max_cell_frac = 0.003; + // Relative X-ray bandwidth (sigma of dlambda/lambda), e.g. ~0.004 for a 1% FWHM // DMM, ~1e-4 for Si(111). Adds a resolution-dependent radial broadening to R[0]. // 0 = monochromatic (the term switches off entirely). @@ -129,6 +136,11 @@ class PixelRefine { double beam[2], double &dist_mm, double detector_rot[2], double latt_vec0[3], double latt_vec1[3], double latt_vec2[3]) const; + + // Reference-driven global orientation + cell-scale sweep (see PixelRefineData::sweep_*). + template + void SweepOrientationCell(const T *image, BraggPrediction &prediction, + PixelRefineData &data) const; public: PixelRefine(const DiffractionExperiment &experiment, const std::vector &reference);