From 3aec235f67012deaed3e05dd174297de4913e1ba Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Tue, 30 Jun 2026 14:26:08 +0200 Subject: [PATCH] integration: add --reciprocal-profile dial (global reciprocal-space width) The per-shell profile width is learned in pixels, so it varies ~4x with resolution (mostly the geometric projection of a near-constant reciprocal-space relrod) and must be binned per shell -> it starves at high resolution / on sparse data. The new --reciprocal-profile flag instead learns ONE global width in reciprocal space, sigma2_q,tan = A + B|q| + C|q|^2: the Jacobian g_tan=cos(2theta) removes the geometric projection, and C|q|^2 is the crystal mosaicity relrod (variance ~(eta|q|)^2). Applied per reflection as sigma2_tan,px = (A + B|q| + C|q|^2)/g_tan^2 (B,C clamped >=0; quadratic->linear->constant fallback). Off by default. On the sharp HEWL test crystal (mosaicity 0.091deg, so C fits to ~0 and it reduces to the validated linear form) it is metric-neutral: ISa 16.2->16.3, anomalous 0.92x unchanged, CCref band 90.0->89.9, CC1/2 a touch lower (per-shell isn't starved at 23k spots/shell, and a global fit is less flexible). So: simpler + more transferable at a small CC1/2 cost, ISa/anomalous held. Its payoff is on MOSAIC crystals (large C|q|^2), where per-shell starves on the wide weak high-res spots and 6 shells are too coarse; both lyso test crystals are sharp, so it ships as a dial to try on mosaic data elsewhere. A separate radial relrod fit was tried and dropped (no gain). See NEXTGEN_INTEGRATOR.md. Co-Authored-By: Claude Opus 4.8 --- common/BraggIntegrationSettings.cpp | 9 +++ common/BraggIntegrationSettings.h | 3 + .../bragg_integration/NEXTGEN_INTEGRATOR.md | 28 ++++++++ .../bragg_integration/ProfileIntegrate2D.cpp | 72 +++++++++++++++++-- tools/jfjoch_process.cpp | 16 ++++- 5 files changed, 123 insertions(+), 5 deletions(-) diff --git a/common/BraggIntegrationSettings.cpp b/common/BraggIntegrationSettings.cpp index d69a0028..bdc1fdcb 100644 --- a/common/BraggIntegrationSettings.cpp +++ b/common/BraggIntegrationSettings.cpp @@ -77,6 +77,15 @@ IntegratorMode BraggIntegrationSettings::GetIntegrator() const { return integrator_mode; } +BraggIntegrationSettings &BraggIntegrationSettings::ReciprocalProfile(bool input) { + reciprocal_profile = input; + return *this; +} + +bool BraggIntegrationSettings::GetReciprocalProfile() const { + return reciprocal_profile; +} + float BraggIntegrationSettings::GetR1() const { return r_1; } diff --git a/common/BraggIntegrationSettings.h b/common/BraggIntegrationSettings.h index b8292916..36c0473b 100644 --- a/common/BraggIntegrationSettings.h +++ b/common/BraggIntegrationSettings.h @@ -20,6 +20,7 @@ class BraggIntegrationSettings { float d_min_limit_A = 1.0; std::optional fixed_profile_radius; float minimum_sigma_in_regards_to_i = 0.02; + bool reciprocal_profile = false; public: BraggIntegrationSettings& R1(float input); @@ -28,6 +29,7 @@ public: BraggIntegrationSettings& DMinLimit_A(float input); BraggIntegrationSettings& FixedProfileRadius_recipA(std::optional input); BraggIntegrationSettings& Integrator(IntegratorMode input); + BraggIntegrationSettings& ReciprocalProfile(bool input); [[nodiscard]] IntegratorMode GetIntegrator() const; @@ -38,4 +40,5 @@ public: [[nodiscard]] float GetDMinLimit_A() const; [[nodiscard]] float GetMinimumSigmaInRegardsToI() const; + [[nodiscard]] bool GetReciprocalProfile() const; }; diff --git a/image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md b/image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md index 808d6ba3..0bc3acf8 100644 --- a/image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md +++ b/image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md @@ -322,6 +322,34 @@ Open follow-up (deferred, untested): a neighbour mask on the learning accumulati alternative to the `r1` window (keeps this spot's genuine [r1,r2] tail; needs a multiplicity mask to separate self from neighbour in the overlap annulus). +### Reciprocal-space profile width — `--reciprocal-profile` (experimental dial, default OFF) + +The per-shell width is learned in **pixels**, so it changes ~4× with resolution (mostly the geometric +projection of a roughly-constant reciprocal-space relrod) and must be binned per shell to have enough +strong spots — which **starves** at high resolution / on sparse data. Mapping the spot 2nd-moment into +**reciprocal space** (multiply the pixel tangential moment by the Jacobian `g_tan² = cos²2θ`) removes the +projection and makes the width far more transferable: across resolution shells the spread collapses from +CV ≈ 44 % (pixel) to ≈ 11 % radial / 27 % tangential, and the tangential residual is a clean mosaicity +law `σ²_q,tan = A + B|q| + C|q|²` (`C|q|²` = relrod variance `~(η|q|)²`). + +`--reciprocal-profile` replaces the per-shell pixel width with **one global** fit of that law (`B`,`C` +clamped ≥ 0; falls back linear→constant). On the **sharp** HEWL test crystal (η = 0.091°, so the `C|q|²` +mosaicity term fits to noise/≈0 and it reduces to the validated linear form) it is **metric-neutral**: +ISa 16.2→16.3, anomalous 0.92×→0.91× (no drop), CCref band 90.0→89.9, CC1/2 a touch lower (1.12 Å +95→93) — because per-shell isn't starved here (~23 k spots/shell) and a 2–3-parameter global fit is just +less flexible, while the faithful pure-tangential width is slightly tighter than the (generous) isotropic +per-shell seed the de-biased fit mildly prefers. So on a sharp crystal it is **simpler + more transferable +at a small CC1/2 cost, ISa/anomalous unchanged** — kept as an off-by-default dial, not the default. + +**Where it should pay off, and why it can't be shown here:** a **mosaic** crystal makes the `C|q|²` term +large and dominant at high resolution, where (a) per-shell starves on the wide, weak high-res spots, +(b) 6 discrete shells are too coarse for the steep width growth, and (c) the wide spots overflow the +fixed grid (the ellipse path grows to hold them). This is exactly why DIALS/XDS integrate in a reciprocal +profile coordinate. Both lyso test crystals are sharp (and the jets are sharp serial stills), so the dial +can only be **demonstrated on a genuinely mosaic rotation dataset** — the reason it ships as a dial to try +elsewhere. Tested and dropped along the way: a separate **radial** relrod fit (radial is divergence + +rotation-smear, not mosaicity → no metric gain on the sharp crystal). + ## Lineage: PixelRefine (removed 2026-06-25) `ProfileIntegrate2D` is the surviving half of an earlier experimental still-integrator, **PixelRefine**, diff --git a/image_analysis/bragg_integration/ProfileIntegrate2D.cpp b/image_analysis/bragg_integration/ProfileIntegrate2D.cpp index cf2fe8c2..ae210159 100644 --- a/image_analysis/bragg_integration/ProfileIntegrate2D.cpp +++ b/image_analysis/bragg_integration/ProfileIntegrate2D.cpp @@ -284,6 +284,64 @@ std::vector ProfileIntegrateInternal(const DiffractionExperiment &ex const double F_px = geom.GetDetectorDistance_mm() / std::max(1e-6f, geom.GetPixelSize_mm()); const bool use_ellipse = !empirical && (bw_sigma > 0.0 || c_radial > 0.0); + // Reciprocal-space global profile width (--reciprocal-profile): one global model + // sigma2_q,tan = A + B*|q| + C*|q|^2 over all strong spots replaces the per-shell pixel width. The + // Jacobian g_tan = cos(2theta) maps the pixel tangential moment into reciprocal space, removing the + // geometric projection that makes the pixel width grow ~4x with resolution. The C*|q|^2 term is the + // crystal MOSAICITY (relrod variance ~ (eta*|q|)^2): ~0 for a sharp crystal (the fit collapses to + // A + B*|q|, tying the per-shell model) but dominant for a mosaic one, where per-shell starves on the + // wide, weak high-res spots. Applied per reflection as sigma2_tan,px = (A + B*|q| + C*|q|^2)/g_tan^2. + const bool recip_on = settings.GetReciprocalProfile(); + double recip_A = 0.0, recip_B = 0.0, recip_C = 0.0; + bool use_recip = false; + if (recip_on && !empirical) { + double s1 = 0, sq = 0, sq2 = 0, sq3 = 0, sq4 = 0, sy = 0, sqy = 0, sq2y = 0; + for (size_t i = 0; i < npredicted; ++i) { + const auto &rh = rough[i]; + if (!rh.ok || !rh.strong || rh.I <= 0.0 || !(predicted[i].d > 0.0f)) continue; + const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y; + const double Rpx = std::hypot(rx, ry); + if (Rpx < 1e-6) continue; + const double ux = rx / Rpx, uy = ry / Rpx; + double m2 = 0.0, m2w = 0.0; + for (int dy = -R; dy <= R; ++dy) + for (int dx = -R; dx <= R; ++dx) { + if (dx * dx + dy * dy >= r1_sq) continue; + const int64_t x = rh.cx + dx, y = rh.cy + dy; + if (x < 0 || y < 0 || x >= static_cast(xpixel) || y >= static_cast(ypixel)) continue; + const auto px = ptr[y * xpixel + x]; + if (px == special || px == special + 1 || px == saturation || px == saturation - 1) continue; + const double w = std::max(0.0, (static_cast(px) - rh.bkg) / rh.I); + const double tn = -dx * uy + dy * ux; + m2 += w * tn * tn; m2w += w; + } + if (m2w <= 0.0) continue; + const double tan2t = Rpx / F_px, cos2t = 1.0 / std::sqrt(1.0 + tan2t * tan2t); + const double q = 1.0 / predicted[i].d, yv = cos2t * cos2t * (m2 / m2w); + s1 += 1; sq += q; sq2 += q * q; sq3 += q * q * q; sq4 += q * q * q * q; + sy += yv; sqy += q * yv; sq2y += q * q * yv; + } + // Fit quadratic, falling back to linear then constant so B and C stay >= 0 (a relrod cannot shrink + // with resolution; a sharp crystal's C comes out as noise/<0 and drops to the linear branch). + auto det3 = [](double a, double b, double c, double d, double e, double f, double g, double h, double i) { + return a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g); + }; + if (s1 >= 30) { + const double D3 = det3(s1, sq, sq2, sq, sq2, sq3, sq2, sq3, sq4), D2 = s1 * sq2 - sq * sq; + const double C3 = std::fabs(D3) > 1e-9 ? det3(s1, sq, sy, sq, sq2, sqy, sq2, sq3, sq2y) / D3 : -1.0; + const double B3 = std::fabs(D3) > 1e-9 ? det3(s1, sy, sq2, sq, sqy, sq3, sq2, sq2y, sq4) / D3 : 0.0; + const double B2 = std::fabs(D2) > 1e-9 ? (s1 * sqy - sq * sy) / D2 : -1.0; + if (std::fabs(D3) > 1e-9 && C3 > 0.0 && B3 >= 0.0) { + recip_A = det3(sy, sq, sq2, sqy, sq2, sq3, sq2y, sq3, sq4) / D3; recip_B = B3; recip_C = C3; + } else if (std::fabs(D2) > 1e-9 && B2 >= 0.0) { + recip_A = (sy * sq2 - sqy * sq) / D2; recip_B = B2; + } else { + recip_A = sy / s1; // constant width fallback + } + use_recip = recip_A > 0.0; + } + } + // --- Pass B: profile-fit each reflection (Kabsch, de-biased variance v = B + I*P; iterate). --- std::vector out; out.reserve(npredicted); @@ -296,13 +354,19 @@ std::vector ProfileIntegrateInternal(const DiffractionExperiment &ex int Rf = R; const std::vector *Pvec = &shell_P[sh]; // ProfileEmpirical uses the shared learned grid - double s2t = shell_sigma2[sh], s2r = s2t; + const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y; + const double Rpx = std::hypot(rx, ry); + const double tan2t = Rpx / F_px; + double s2t = shell_sigma2[sh]; + if (use_recip) { // global reciprocal width instead of the per-shell pixel width + const double q = 1.0 / std::max(predicted[i].d, 1e-6f), cos2t = 1.0 / std::sqrt(1.0 + tan2t * tan2t); + s2t = std::max(0.25, (recip_A + recip_B * q + recip_C * q * q) / (cos2t * cos2t)); + } + double s2r = s2t; double ux = 1.0, uy = 0.0; bool elong = false; if (use_ellipse) { - const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y; - const double Rpx = std::hypot(rx, ry), sbw = bw_sigma * Rpx; - const double tan2t = Rpx / F_px; + const double sbw = bw_sigma * Rpx; const double radial_extra = sbw * sbw + c_radial * tan2t * tan2t; // bandwidth + parallax + capture // Only elongate where the radial streak adds a genuine fraction of a pixel of variance; at // low/mid resolution the smear is sub-pixel and elongating just adds noise. diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 9cee0d60..4df0671e 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -72,6 +72,7 @@ void print_usage() { std::cout << " --bandwidth Relative X-ray bandwidth FWHM (e.g. 0.01 for 1% DMM); default from file or 0" << std::endl; std::cout << " --integration-radius Signal-box radius r1, or r1,r2,r3 (px). One value => r2=r1+2, r3=r1+4" << std::endl; std::cout << " --integrator Spot integrator boxsum|gaussian|empirical (default: gaussian profile-fit; boxsum is the classical fallback)" << std::endl; + std::cout << " --reciprocal-profile Learn one global reciprocal-space profile width (A+B|q|+C|q|^2) instead of per-shell; helps mosaic/sparse data" << std::endl; } enum { @@ -97,7 +98,8 @@ enum { OPT_SCALE_FULLS, OPT_CAPTURE_UNCERTAINTY, OPT_MOSAICITY, - OPT_SMOOTH_G + OPT_SMOOTH_G, + OPT_RECIPROCAL_PROFILE }; static option long_options[] = { @@ -142,6 +144,7 @@ static option long_options[] = { {"bandwidth", required_argument, nullptr, OPT_BANDWIDTH}, {"integration-radius", required_argument, nullptr, OPT_INTEGRATION_RADIUS}, {"integrator", required_argument, nullptr, OPT_INTEGRATOR}, + {"reciprocal-profile", no_argument, nullptr, OPT_RECIPROCAL_PROFILE}, {"reject-outliers", required_argument, nullptr, OPT_REJECT_OUTLIERS}, {"reject-delta-cchalf", required_argument, nullptr, OPT_REJECT_DELTA_CCHALF}, {nullptr, 0, nullptr, 0} @@ -314,6 +317,7 @@ int main(int argc, char **argv) { std::optional d_min_scale_merge; std::optional integration_radius_arg; // "r1" or "r1,r2,r3" std::optional integrator_mode; // --integrator boxsum|gaussian|empirical + bool reciprocal_profile = false; // --reciprocal-profile std::optional outlier_reject_nsigma; // merge per-observation outlier rejection std::optional delta_cchalf_nsigma; // per-crystal CC1/2-delta rejection @@ -535,6 +539,9 @@ int main(int argc, char **argv) { else if (strcmp(optarg, "empirical") == 0) integrator_mode = IntegratorMode::ProfileEmpirical; else { logger.Error("--integrator expects boxsum|gaussian|empirical"); return 1; } break; + case OPT_RECIPROCAL_PROFILE: + reciprocal_profile = true; + break; case OPT_REJECT_OUTLIERS: outlier_reject_nsigma = std::stod(optarg); break; @@ -773,6 +780,13 @@ int main(int argc, char **argv) { : "profile (empirical)"); } + if (reciprocal_profile) { + BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings(); + bis.ReciprocalProfile(true); + experiment.ImportBraggIntegrationSettings(bis); + logger.Info("Reciprocal-space global profile width enabled (per-shell width replaced)"); + } + SpotFindingSettings spot_settings; spot_settings.enable = true; spot_settings.indexing = true;