diff --git a/common/BraggIntegrationSettings.cpp b/common/BraggIntegrationSettings.cpp index bdc1fdcb..d69a0028 100644 --- a/common/BraggIntegrationSettings.cpp +++ b/common/BraggIntegrationSettings.cpp @@ -77,15 +77,6 @@ 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 36c0473b..b8292916 100644 --- a/common/BraggIntegrationSettings.h +++ b/common/BraggIntegrationSettings.h @@ -20,7 +20,6 @@ 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); @@ -29,7 +28,6 @@ public: BraggIntegrationSettings& DMinLimit_A(float input); BraggIntegrationSettings& FixedProfileRadius_recipA(std::optional input); BraggIntegrationSettings& Integrator(IntegratorMode input); - BraggIntegrationSettings& ReciprocalProfile(bool input); [[nodiscard]] IntegratorMode GetIntegrator() const; @@ -40,5 +38,4 @@ public: [[nodiscard]] float GetDMinLimit_A() const; [[nodiscard]] float GetMinimumSigmaInRegardsToI() const; - [[nodiscard]] bool GetReciprocalProfile() const; }; diff --git a/docs/JFJOCH_PROCESS.md b/docs/JFJOCH_PROCESS.md index 0ad4bc07..136eb482 100644 --- a/docs/JFJOCH_PROCESS.md +++ b/docs/JFJOCH_PROCESS.md @@ -171,5 +171,4 @@ Integration: | --- | --- | | `--integrator ` | Spot integrator: `gaussian` (profile-fit, default) \| `empirical` \| `boxsum` (classical fallback) | | `--integration-radius ` | Signal-box radius `r1`, or `r1,r2,r3` (px). One value ⇒ `r2=r1+2`, `r3=r1+4` | -| `--reciprocal-profile` | Learn one global reciprocal-space profile width (`A+B·|q|+C·|q|²`) instead of per-shell; helps mosaic/sparse data | | `--bandwidth ` | Relative X-ray bandwidth FWHM (e.g. `0.01` for a 1% DMM); default from file or 0 (monochromatic) | diff --git a/image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md b/image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md index 0bc3acf8..9d450be2 100644 --- a/image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md +++ b/image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md @@ -322,33 +322,41 @@ 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) +### Reciprocal-space profile width — `--reciprocal-profile` (tried and REMOVED 2026-07-02) -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|)²`). +**The idea.** The per-shell width is learned in **pixels**, so it changes ~4× with resolution (the +geometric projection of a roughly-constant reciprocal-space relrod) and must be binned per shell, which +**starves** at high resolution / on sparse data. Mapping the tangential 2nd-moment into **reciprocal +space** (Jacobian `g_tan² = cos²2θ`) should remove the projection and give a transferable mosaicity law +`σ²_q,tan = A + B|q| + C|q|²` (`C|q|²` = relrod variance `~(η|q|)²`), fit once per frame — like XDS/DIALS +integrating in a reciprocal profile coordinate. It was shipped as an off-by-default dial after being +tuned on the single sharp HEWL crystal, where it was metric-neutral. -`--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. +**Why it was removed — an 11-crystal mosaicity-stratified re-test (`/data/rotation_test`, off vs on vs a +de-contaminated variant, plus a per-frame dump of the fitted widths) showed it is net-negative and cannot +work per-frame.** Three linked, proven facts: +1. **The `C|q|²` mosaicity term — the whole point — is unfittable per-frame.** The fitted curvature `a2` + came out ≈ 0 (often slightly *negative*) on **every** crystal, with **zero correlation** to the XDS + mosaicity (0.09→0.42°). Strong spots live at low `q` (q̄ ≈ 0.24–0.47) where `η²q²` is invisible; the + curvature only appears at high `q`, where there are ~0 strong spots to constrain it. The law degenerates + to a straight line. +2. With `a2 ≈ 0`, the high-res width becomes a **blind `1/cos²2θ` extrapolation → 2–4× wider than + per-shell** (e.g. px²: pding4_001 1.07 flat → 2.0–3.0; cytC 1.10 → 2.7; MyoB 0.89 → 3.4). Meanwhile the + per-shell path *starves* at high res (n_strong<30 ⇒ flat global fallback that even **shrinks** vs + mid-res) — which turns out to be **accidentally correct**: weak, crowded high-res spots want a *narrow* + aperture, not the true (wide) spot shape. +3. The over-wide profile pulls background into weak spots ⇒ **per-observation scatter (R-meas) rises**, + CC1/2 drops in reliable, high-multiplicity shells (pding4_001 mult 25, pding4_003 mult 12, EcwtCQ066, + MyoB); the mean intensities survive (CC-overall flat). A confirmatory cap (clamp the reciprocal width at + the widest well-sampled per-shell value) **recovers the regression** — direct proof the over-widening is + the harm. The only "improvements" seen (EcwtAL500, EP_cs_01-17 overall CC) were entirely in + noise-dominated shells (mult 2–3, CC<20 %, R-meas>100 %). **No crystal reliably wins.** -**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). +The per-frame paradigm structurally cannot fit mosaicity, and un-starving the high-res width is net harm, +so the dial (CLI flag, `BraggIntegrationSettings::reciprocal_profile`, and the per-frame fit block) was +deleted. A real win would need cross-frame / rolling-window pooling to constrain `a2`, or a two-pass +mosaicity seed from scaling — both against the per-frame / online + stills philosophy. Also dropped along +the way: a separate **radial** relrod fit (radial is divergence + rotation-smear, not mosaicity → no gain). ## Lineage: PixelRefine (removed 2026-06-25) diff --git a/image_analysis/bragg_integration/ProfileIntegrate2D.cpp b/image_analysis/bragg_integration/ProfileIntegrate2D.cpp index 0a61a7b1..fa0cd49b 100644 --- a/image_analysis/bragg_integration/ProfileIntegrate2D.cpp +++ b/image_analysis/bragg_integration/ProfileIntegrate2D.cpp @@ -17,10 +17,6 @@ namespace { constexpr int N_SHELL = 6; constexpr double STRONG_I_OVER_SIGMA = 5.0; constexpr int MIN_STRONG_PER_SHELL = 30; -// --reciprocal-profile width fit: gentle ridge on the curvature coefficient (sharp-crystal prior, -// so curvature must be earned) and IRLS robust-fit iterations (Huber, to reject outlier spots). -constexpr double RECIP_RIDGE = 0.05; -constexpr int RECIP_IRLS_ITERS = 3; // Radial parallax broadening as the coefficient of tan^2(2theta), i.e. Var(z)/pixel^2 [px^2]. A photon // converts at a random depth z in the sensor (exponential with attenuation length L, truncated at the @@ -288,89 +284,6 @@ 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 profile width (--reciprocal-profile): a per-frame model of the tangential - // variance in reciprocal space, y(q) = a0 + a1*t + a2*t^2 with t = (q - qbar)/qscale, replacing the - // per-shell pixel width. The Jacobian g_tan = cos(2theta) maps the pixel tangential moment into - // reciprocal space (removing the ~4x geometric growth with resolution); the t^2 term is the crystal - // MOSAICITY (relrod variance ~ (eta*|q|)^2), ~0 for a sharp crystal. Fitting in the CENTERED, - // STANDARDIZED variable t rather than raw q keeps the 3x3 normal matrix well-conditioned even when - // the strong spots span a narrow q-range (small cell / sparse still) - the raw {1,q,q^2} fit went - // near-singular there, letting tiny per-frame jitter swing the curvature into a wild over-wide - // profile. The fit is robust (IRLS / Huber: outlier spots can't drag it) with a gentle ridge on the - // curvature (sharp crystal = prior). Applied per reflection as sigma2_tan,px = y(q)/g_tan^2, with q - // clamped to the fitted strong-spot range (never extrapolated). - const bool recip_on = settings.GetReciprocalProfile(); - double rp_a0 = 0.0, rp_a1 = 0.0, rp_a2 = 0.0, rp_qbar = 0.0, rp_qscale = 1.0, rp_tmin = 0.0, rp_tmax = 0.0; - bool use_recip = false; - if (recip_on && !empirical) { - std::vector qv, yv; - qv.reserve(npredicted); yv.reserve(npredicted); - 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 y = cos2t * cos2t * (m2 / m2w); - if (!(y > 0.0) || !std::isfinite(y)) continue; - qv.push_back(1.0 / predicted[i].d); yv.push_back(y); - } - const size_t n = qv.size(); - if (n >= 30) { - double qmin = qv[0], qmax = qv[0], qsum = 0.0; - for (double q : qv) { qmin = std::min(qmin, q); qmax = std::max(qmax, q); qsum += q; } - rp_qbar = qsum / static_cast(n); - double var = 0.0; for (double q : qv) var += (q - rp_qbar) * (q - rp_qbar); - rp_qscale = std::sqrt(var / static_cast(n)); - if (!(rp_qscale > 1e-9)) rp_qscale = std::max(1e-6, 0.5 * (qmax - qmin)); - rp_tmin = (qmin - rp_qbar) / rp_qscale; rp_tmax = (qmax - rp_qbar) / rp_qscale; - 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); - }; - std::vector wob(n, 1.0); - for (int it = 0; it < RECIP_IRLS_ITERS; ++it) { - double S0=0,S1=0,S2=0,S3=0,S4=0, T0=0,T1=0,T2=0; - for (size_t i = 0; i < n; ++i) { - const double t = (qv[i] - rp_qbar) / rp_qscale, t2 = t*t, w = wob[i]; - S0 += w; S1 += w*t; S2 += w*t2; S3 += w*t2*t; S4 += w*t2*t2; - T0 += w*yv[i]; T1 += w*t*yv[i]; T2 += w*t2*yv[i]; - } - const double S4r = S4 + RECIP_RIDGE * S4; // ridge: shrink curvature toward the sharp prior - const double D = det3(S0,S1,S2, S1,S2,S3, S2,S3,S4r); - if (!(std::fabs(D) > 1e-12)) { use_recip = false; break; } - rp_a0 = det3(T0,S1,S2, T1,S2,S3, T2,S3,S4r) / D; - rp_a1 = det3(S0,T0,S2, S1,T1,S3, S2,T2,S4r) / D; - rp_a2 = det3(S0,S1,T0, S1,S2,T1, S2,S3,T2) / D; - use_recip = true; - if (it + 1 < RECIP_IRLS_ITERS) { // Huber reweight from residuals (scale = 1.4826*MAD) - std::vector res(n); - for (size_t i = 0; i < n; ++i) { - const double t = (qv[i]-rp_qbar)/rp_qscale; - res[i] = std::fabs(yv[i] - (rp_a0 + rp_a1*t + rp_a2*t*t)); - } - std::vector tmp(res); std::nth_element(tmp.begin(), tmp.begin()+n/2, tmp.end()); - const double mad = std::max(1e-6, 1.4826 * tmp[n/2]); - for (size_t i = 0; i < n; ++i) { const double r = res[i]/(1.5*mad); wob[i] = r <= 1.0 ? 1.0 : 1.0/r; } - } - } - } - } - // --- Pass B: profile-fit each reflection (Kabsch, de-biased variance v = B + I*P; iterate). --- std::vector out; out.reserve(npredicted); @@ -386,13 +299,7 @@ std::vector ProfileIntegrateInternal(const DiffractionExperiment &ex 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) { // per-frame reciprocal width instead of the per-shell pixel width - const double q = 1.0 / std::max(predicted[i].d, 1e-6f); - const double t = std::clamp((q - rp_qbar) / rp_qscale, rp_tmin, rp_tmax); - const double cos2t = 1.0 / std::sqrt(1.0 + tan2t * tan2t); - s2t = std::max(0.25, (rp_a0 + rp_a1 * t + rp_a2 * t * t) / (cos2t * cos2t)); - } + const double s2t = shell_sigma2[sh]; double s2r = s2t; double ux = 1.0, uy = 0.0; bool elong = false; diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index d5659a9b..3c493243 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -85,7 +85,6 @@ 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 { @@ -112,7 +111,6 @@ enum { OPT_CAPTURE_UNCERTAINTY, OPT_MOSAICITY, OPT_SMOOTH_G, - OPT_RECIPROCAL_PROFILE, OPT_DETECT_ICE_RINGS, OPT_NO_SCALE_FULLS, OPT_WRITE_PROCESS_H5, @@ -166,7 +164,6 @@ 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}, {"detect-ice-rings", no_argument, nullptr, OPT_DETECT_ICE_RINGS}, {"reject-outliers", required_argument, nullptr, OPT_REJECT_OUTLIERS}, {"reject-delta-cchalf", required_argument, nullptr, OPT_REJECT_DELTA_CCHALF}, @@ -345,7 +342,6 @@ 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 @@ -581,9 +577,6 @@ 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; @@ -860,13 +853,6 @@ 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;