diff --git a/image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md b/image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md index 201dca40..f0e7120f 100644 --- a/image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md +++ b/image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md @@ -1,105 +1,282 @@ -# Next-generation integrator — design draft (2026-06-24) +# Next-generation integrator & the XDS-quality investigation -## The question +_Last updated 2026-06-25. Supersedes the 2026-06-24 design draft — `ProfileIntegrate2D` is now +built, wired, and the default integrator; this document records the design **and** everything the +data-quality investigation established (especially the dead ends, so they are not re-tried)._ -Can PixelRefine's **profile-fitted intensity extraction** be separated from its **reference-based -scaling**, so we get a profile-fitting integrator that needs no reference intensities — and could -PixelRefine then be rebuilt as *(profile-fit integrator) + (smarter scaling)*? +--- -**Yes.** They are already algebraically separate in the code. +## TL;DR -## Why we want it (the diagnosis) +- **The question that started this:** jfjoch's merged data was ~3–4× worse than XDS on per-observation + precision (R-meas, ISa, anomalous peak height) on the HEWL test crystals. Where is the gap, and can + we close it **without** XDS-style 3D shoeboxes (Filip's philosophy: integrate per-frame, combine at + the end, keep ~300 Hz throughput)? +- **Diagnosis (decisive):** our **mean intensities are accurate everywhere** (`CC(ours, XDS) = 0.994`, + ≥0.986 in every resolution shell) and **unbiased** (ANODE: our anomalous peaks sit on the same real + atoms as XDS). The gap is **per-observation noise**, localized to the **2D box-sum integration** of + **strong** reflections — a ~18% multiplicative floor that a fixed disk imposes by capturing a + width-dependent fraction of each spot. +- **What moved the needle:** swapping box-sum for a **profile-fitting integrator** (`ProfileIntegrate2D`, + Gaussian profile). That is the single lever that ever improved intensity quality. Everything on the + *scaling* side was neutral or negative. +- **Current best on the test crystal** (HEWL `fixed_master.h5`, anomalous S/Cl peak height as a fraction + of XDS, `1.0` = match): **0.61×**, ISa **11.0**, via Gaussian profile-fit + `-R100` + `rot3d` + + `--scale-fulls`. Box-sum baseline was **0.46×**, ISa 8.4. XDS is 1.0× / ISa ~28. +- **The remaining gap is XDS's full XSCALE/CORRECT scaling machinery**, whose individual pieces are each + neutral on this well-behaved crystal — not a quick win, and not a bias to fix. -The residual ~4× gap to XDS on R-meas (jfjoch fulls 18.9% vs XDS 4.5%) was localized to the **2D -spot integration**, not background/scaling/partiality (see `scaling-results-roadmap` memory / the -session experiments): -- the gap is in **strong** reflections (jf/XDS ratio 2.1× weak → 6.5× strong), a ~18% *multiplicative* - per-observation floor, **flat across resolution**, **per-reflection** (only 6.6% of 20.7% is - per-image/scaling), and **immune to box radius** (tight/default/wide all ≈ 18–24%). -- => it is the **box-sum method** (uniform pixel weighting over a fixed disk), `BraggIntegrate2D`. - A fixed disk captures a *width-dependent fraction* of each spot; spot widths vary (mosaicity, - divergence, resolution, detector position, rocking phase) → ~18% multiplicative scatter. -- swapping in PixelRefine's profile-fit swung ISa 6.4→16.5 (confirming the integrator is the lever, - though that path is currently broken — see below). +--- -## The decomposition (PixelRefine today) +## Philosophy & constraints (Filip) -PixelRefine does three things; only the third needs the reference: +1. **Per-frame integrate, combine at the end.** Jungfraujoch integrates each frame as a 2D image and + combines the per-frame partials of a rotation reflection into one "full" afterwards. We deliberately + do **not** batch frames into XDS-style 3D shoeboxes in the online path. Push per-frame + combine as + far as it goes. +2. **Throughput matters.** ~300 Hz end-to-end was demonstrated; that is an order of magnitude over XDS. + Judge every improvement by its **performance penalty**, not just quality (this is why global scaling + was removed — see below). +3. **Modular, not bolted-on.** The integrator is swappable (`BoxSum` ↔ `ProfileGaussian` ↔ + `ProfileEmpirical`) behind a single setting; scaling, the rot3d combine, and the merge consume a + plain `Reflection{I, σ, partiality, d}` and are untouched by the integrator choice. Both **stills and + rotation** benefit. +4. **Simple, readable code.** No speculative/defensive machinery; neutral experiments get reverted. -1. **Profile (Term 2).** Measure the tangential width `R1 = sqrt(2·)` per resolution shell - from the intensity-weighted 2nd moment of strong spots (`PixelRefine.cpp:433`). Reference-free. -2. **Extraction.** `J = Σ Pₚ(Iₚ−B)/vₚ / Σ Pₚ²/vₚ`, `var(J)=1/Σ Pₚ²/vₚ`, with `Pₚ` the - area-normalized tangential profile (Gaussian of width R1) and the **de-biased** variance - `vₚ = max(B,1)` (`PixelRefine.cpp:491-525`, METHODS §1). It also computes the **radial partiality** - `p = P_radial` and a detector-clipping **completeness**. Reference-free. -3. **Scaling (Term 1).** Fit per-image `G,B` from `J ≈ G·B_DW·p·pol·I_ref` (`PixelRefine.cpp:365`, - METHODS §4). **This is the only step that uses I_ref.** +--- -So **integrator = steps 1+2** (output `I=J/(p·pol)`, `σ`, `partiality`), **scaling = step 3** (or the -existing reference-free `ScaleOnTheFly`). A clean split. +## The pipeline today -This split also **fixes the rot3d-combine incompatibility**: PixelRefine-as-a-whole emits -`image_scale_corr=1/G` with no rotation partiality, which is why `-r pixelrefine -P rot3d` produced -garbage σ (CCref 94→84, ⟨I/σ⟩→1.3). The *decomposed* integrator emits `partiality` like -`BraggIntegrate2D`, so `ScaleOnTheFly` + `Combine3D` + merge consume it unchanged. +``` +per-frame 2D integration (BraggIntegrate2D | ProfileIntegrate2D) <- the lever + -> per-image scaling (ScaleOnTheFly: Ceres fit G, B_resolution, mosaicity, wedge) + -> rot3d combine (Combine3D: weight-sum per-frame partials -> one full per event) + -> scale fulls (optional) (ScaleOnTheFly again, Unity model, on the combined fulls) + -> merge + error model + stats (MergeOnTheFly: RefineErrorModel, CC1/2, R-meas, CCref, ISa) +``` -## The new-generation integrator +Relevant knobs (`jfjoch_process`, mirrored in the viewer's Processing-settings tabs): +- `--integrator boxsum|gaussian|empirical` — default **gaussian** (`IntegratorMode::ProfileGaussian`). +- `-P fixed|rot|rot3d` — partiality model. `rot3d` = `PartialityModel::Rotation` **+** the orthogonal + `ScalingSettings::combine_3d` bool (see "Why rot3d is a bool" below). +- `--scale-fulls` — second scaling pass on the combined fulls (XDS order). Implies `-M`, needs `rot3d`. +- `-R` — first-pass rotation-indexing images; default **100**. +- `-M` — run scaling; `-A` — anomalous (keep Friedel mates separate). -A new module `ProfileIntegrate2D`, **swappable with `BraggIntegrate2D`**, same output -(`vector` with I, σ, partiality, d, …), **no reference**, no scaling. Algorithm: +### Current integrator: `ProfileIntegrate2D` -1. **Pass A (rough + learn the profile).** Box-sum every reflection (reuse `IntegrateReflection`) - to get a rough I and the observed centroid; select strong spots (significance ≥ 5). -2. **Learn an EMPIRICAL profile**, per *(detector region × resolution shell)*, by averaging the - background-subtracted, intensity-normalized, centroid-aligned pixel grids of the strong spots. - **This is the key difference from PixelRefine**, which uses a *Gaussian* of width R1×6 — a model - that has the *same width-mismatch* as a box (and ×6 ≈ near-box, which is why PixelRefine's stills - ISa was only 1.1). A *learned* profile matches the real spot — including the 0.4-px undersampling - smear, asymmetry, and position dependence — so the weighting is genuinely optimal. (XDS learns its - profile over 9 detector regions × resolution; this is that idea.) -3. **Pass B (profile-fit extract).** For each reflection, `I = Σ Pₚ(cₚ−B)/vₚ / Σ Pₚ²/vₚ` with the - learned `P` for its region/shell, de-biased variance `vₚ = B + max(I,0)·Pₚ` (Kabsch; iterate a - couple of times), `σ = sqrt(1/Σ Pₚ²/vₚ)` corrected for the signal term. Carry the rotation - `partiality` exactly as `BraggIntegrate2D` does (so rot3d works). +Drop-in for `BraggIntegrate2D`, same `vector` output, **no reference intensities**. +Selected at `IndexAndRefine.cpp:347`. Per frame: -**Why it removes the floor:** profile-weighting integrates *out to the tail* (capturing the -width-varying part a fixed disk clips) while *downweighting the noisy edge* (avoiding the wide-box -background penalty that made radius 6,9,12 worse). A *learned* profile (not a Gaussian) removes the -model-mismatch that limited PixelRefine. The floor is multiplicative and in strong reflections → this -is the right lever. +- **Pass A — rough + pick.** Box-sum every predicted reflection (rough `I`, background, centroid); + select strong spots (significance ≥ 5). +- **Pass B — build the profile**, per resolution shell, from the strong spots: + - `ProfileGaussian` (**default, the keeper**): a Gaussian of the measured second moment. + - `ProfileEmpirical`: the averaged, centroid-aligned, background-subtracted pixel grid. +- **Pass C — profile-fit extract (Kabsch).** `I = Σ P(c−B)/v ÷ Σ P²/v`, de-biased variance + `v = B + max(I,0)·P` (iterated), `σ = √(1/Σ P²/v)`. Rotation `partiality` carried exactly as + `BraggIntegrate2D` (so `rot3d` works unchanged). -**Caution (from PixelRefine):** a *tight* profile LOST to a generous box for stills (the 0.4-px -centroid undersampling floor). An empirical profile is immune to this *because it is measured from the -undersampled spots themselves* — the template already contains the smear. The lever is *matched -weighting + de-biased variance*, never a tight aperture. +**Why it removes the floor:** profile weighting integrates out to the spot tail (the width-varying part +a fixed disk clips) while down-weighting the noisy edge (avoiding the wide-box background penalty). The +floor is multiplicative and lives in strong reflections, so matched weighting + de-biased variance is +the right lever — **not** a tighter aperture (a tight profile *lost* to a generous box for stills, the +0.4-px centroid-undersampling floor). -## Plug-in (no touch to rot3d / scaling / merge) +**Measured, do not re-try:** +- **Anisotropy / per-region profiles add nothing.** Radial/tangential anisotropic Gaussians and + per-detector-region profiles were tested → neutral-to-slightly-worse. The 2D detector-plane spots are + ~round; the real anisotropy is in the **rocking** direction, which per-frame 2D integration discards. + Keep the **isotropic per-shell** width. +- **The empirical profile under-performs the Gaussian** as built (per-frame, integer-pixel-binned, + centred on `round(predicted)`): sub-pixel smear + too few strong spots per frame/shell. It would need + sub-pixel alignment + cross-frame accumulation to compete; not currently worth it. +- _(Note: the header comment in `ProfileIntegrate2D.h` still says "DRAFT / not yet wired in" — that is + stale; it is wired and default.)_ -- New `ProfileIntegrate2D(experiment, image, predicted, …) -> vector`, output-compatible - with `BraggIntegrate2D`. -- Select it like the geometry-refinement choice (a `BraggIntegrationSettings` enum or an - `IntegrationAlgorithm`), default box-sum. Everything downstream (ScaleOnTheFly, Combine3D, Merge) - is unchanged — they only see `Reflection{I,σ,partiality,d}`. -- Reusable by **stills and rotation** (and the FPGA path later, if the profile is precomputable). +### The rot3d combine (`Combine3D`) -## Staging +Buckets per-frame partials by raw `(h,k,l)`, splits into contiguous-frame rocking events (frame gap ≤2), +inverse-variance weight-sums each event's corrected partials into **one full** +(`F = Σ wⱼ·I_corr,j / Σ wⱼ`, `σ = 1/√Σ wⱼ`; partiality = 1, `image_scale_corr` = 1), dropping events +with `Σ pⱼ < min_partiality`. -- **v1** — standalone profile-fit with the *measured-R1 Gaussian* (= PixelRefine steps 1+2 lifted out, - no reference). Cheap; reproduces PixelRefine extraction and proves the decomposition + the - rot3d-combine compatibility. Likely only a modest ISa gain (Gaussian model). -- **v2** — **empirical learned profile per resolution shell.** The expected real win. -- **v3** — empirical profile per *detector region × shell* (XDS-grade); anisotropic. +- **De-biased Poisson variance (kept):** weight each frame by `var = max(0, corr²(σ²−I)) + corr·max(0,F)` + (background noise + **model** signal, iterated 3×) instead of `σ_corr²`. Using the model signal (not + the down-fluctuating observed `I`) stops weak partials being over-weighted. ISa 5.3→6.4 (ref), + 6.7→8.4 (no-ref). +- **Validation that the grouping is structurally right:** jfjoch produces ~the **same observation count + as XDS** (936k vs 961k, multiplicity 25.1 vs 25.9). The residual gap is pure per-observation noise, + not bad event grouping. -## Validation +**Why `rot3d` is a `ScalingSettings::combine_3d` bool, not a 4th `PartialityModel` enum.** Many sites +gate on `PartialityModel::Rotation` *exactly* (e.g. `GetRotationWedgeForScaling`, rotation prediction). +A 4th enum value silently broke per-image rotation handling → NaN scales, all partials rejected. So +`-P rot3d` sets `PartialityModel::Rotation` **and** `combine_3d = true`; the pipeline only ever sees +`Rotation`, and the post-pass keys off `GetScalingSettings().GetCombine3D()`. **Do not reintroduce a +rot3d enum value.** -A/B per-observation **R-meas / ISa** vs box-sum, both with rot3d, on crystal 2 (`fixed_master.h5`) -against `XDS_ASCII.HKL` (4.5%); repeat treating it as stills (crystal 1 jet) where box-sum's hidden -floor is ISa 1.1. Use `--dump-observations`. Target: close the 18% strong-reflection floor toward -XDS's ~3%. +### Scale-on-fulls (`--scale-fulls`) — XDS order + +XDS/DIALS scale on the 3D-integrated **fulls**; jfjoch scales per-frame **partials** then combines, so +the per-frame scale `G` is fit via `G·partiality·B·lp·Itrue − Iobs` — rocking-curve/partiality **model +error leaks into G**, and the combined fulls inherit that partial scale (never scaled *as* fulls). Fix: +a second `ScaleOnTheFly` pass on the combined fulls with the **Unity** model (`G·Itrue − I_full`, no +partiality), between combine and merge. **This is the only post-integration win:** anomalous 0.53→0.57×, +ISa 9.4→10.5 — and it improves **both** precision and accuracy (unlike outlier rejection, which trades +accuracy for CC½). Kept. + +--- + +## The diagnosis — where the gap is (and is not) + +The investigation localized the gap with three independent instruments (see Tooling) and a per-reflection +comparison against `XDS_ASCII.HKL` (same 4/mmm ASU). Findings: + +**Pipeline-wide systematic-error map (ours vs XDS, per reflection):** +- `CC(ours, XDS) = 0.994` overall, ≥0.986 in **every** resolution shell → **mean intensities are accurate + everywhere.** The deficiencies are systematic *scale trends*, not random error. +- Ratio ours/XDS by resolution: 3.0 Å **0.94 (strong, −6%)** → 1.8 Å 1.02 → 1.45 Å 1.19 → 1.2 Å + **1.27 (weak, +27%)**. Two distinct problems: + 1. **Strong low-res:** mean ~right but ~6% low (mild saturation/non-linearity at high counts) **and** + per-observation scatter R-meas **13% vs XDS 2.9%** — averages out (CC fine) but kills **anomalous** + (Bijvoet differences need <1% agreement). This is the integration-quality floor. + 2. **Weak high-res:** **+20–27% high = background under-subtraction.** Code-confirmed: the integrator + subtracts a **flat ring mean** (`ProfileIntegrate2D.cpp` `out.bkg = bkg_sum/n_bkg`), so an unmodeled + background **gradient** under the peak biases weak reflections high, worst at high res. + +**Anomalous accuracy (ANODE, external benchmark):** averaged anomalous peak height at the model atoms, +ours vs XDS — SD_MET 6.05 vs 13.23, SG_CYS 5.56 vs 11.79, CL_CL 4.06 vs 9.58 (XDS ~2.2× stronger). **But +both put their strongest peaks on the same real atoms** (<0.2 Å off; non-anomalous atoms in the noise). +→ our data is **accurate / unbiased**; the gap is **precision (noise)**, ~2.2× on the anomalous ΔF, +consistent with R-meas (~3×) and ISa (~3.5×). We are not systematically wrong, just ~2.2× noisier. + +**Where the strong-reflection floor lives:** it is in **strong** reflections (jf/XDS R-meas ratio 2.1× +weak → 6.5× top-5%), a ~18% **multiplicative** floor, **flat across resolution**, **per-observation** +(per-image/scaling is only 6.6% of 20.7%), **immune to box radius** (tight/default/wide all ≈ 18–24%), +and present even in fully-captured spots (so not partiality/capture). It survives every box size → +it is the **box-sum method** (uniform pixel weighting), confirmed by the profile-fit win. + +--- + +## Experiment ledger — what was tried (the dead ends are the value) + +Scored on the HEWL test crystal, anomalous S/Cl peak height as ×XDS (`1.0` = XDS), unless noted. +Baseline progression that **worked**: + +| change | effect | verdict | +|---|---|---| +| box-sum → **2D Gaussian profile-fit** | 0.46 → **0.53×**, ISa 8.4 → 9.4, R-meas 18.9 → 14.6% | ✅ **the lever** (committed default) | +| **`rot3d` combine** (+ de-biased Poisson var) | ISa 1.7 → 5.3 → 6.4; R-meas (partials 67% → fulls 18.9%) | ✅ keep | +| **`--scale-fulls`** (XDS order) | 0.53 → **0.57×**, ISa 9.4 → 10.5 (precision **and** accuracy) | ✅ keep | +| **`-R100`** (more first-pass indexing images) | 0.57 → **0.61×**, ISa 10.5 → 11.0 (~free; sweet spot ~100) | ✅ keep (committed default) | +| masked-sentinel fix (`±1` band) | neutral on lyso; correctness (gap-clipping reflections elsewhere) | ✅ keep (committed) | + +Things that **did not work** (do not re-prototype without new evidence): + +| change | result | why it failed | +|---|---|---| +| 3D box-sum / single 3D background | R-meas identical, ISa bump is a σ-model artifact | the 3D box-sum is reconstructible from the 2D partials; recovers nothing | +| **3D profile-fit** (`Integrate3D`, voxel Kabsch) | 0.53 → **0.56× (+6%)**, diminishing (box→2D +15%, 2D→3D +6%) | integration cannot bridge the rest; offline-only, slower, dropped | +| **global (joint) scaling** | 0.54× (= 0.53 alternating), **+120 s** | block-coordinate loop already reaches the joint optimum; **removed (perf)** | +| **absorption surface** (per-image / global quadratic in detector xy) | 0.54×, ISa 9.4 → 7.3 (worse) | radially degenerate with per-HKL resolution structure → "intensity sponge"; we're unbiased, a position-scale can't help | +| **learned rocking-curve profile** (on partials) | ISa 9.4 → 6.7; with `--scale-fulls` 0.37× | partials are already divided by the model partiality; the learned `c(t)` spikes at curve edges and up-weights noise | +| **per-batch geometry refinement** (30° batches) | 0.61 → 0.56–0.57× | batch deltas confirm real smooth drift, **but** per-image refine wants the most *stable* start; a local-but-noisy per-batch start hurts | +| **sinusoidal geom-drift model** (beam/distance vs φ) | ~0 drift found; 0.57–0.58× | this crystal is well-centred (beam sub-pixel, distance 0.04%); physically sound, no drifty data to demonstrate | +| background **plane fit** (a + gₓ·dx + g_y·dy over the ring) | anomalous identical; +27% weak bias persists | the weak bias is **not** a linear gradient (it's radial / profile-vs-summation); no superfluous code | +| per-frame **B-factor** decay (`--scale-fulls -B`) | 0.57 → 0.56× | no significant radiation damage; the scalar `G` already absorbs it | +| outlier rejection (`--reject-outliers`, `--reject-delta-cchalf`) | CC½ ↑ **but anomalous ↓** (0.53 → 0.49×) | **the accuracy trap** — improves precision, degrades accuracy; not Bijvoet-safe. Do not enable for anomalous | +| raise `--min-partiality` (0.02 → 0.6) | 0.53 → 0.35× | low-partiality partials *feed* the scale fit; partiality model is good, post-refining it doesn't help | +| cut high resolution | hurts (1.2 Å > 2.0 Å) | high-res helps; not over-weighted noise | +| more scaling iterations; error-model tuning | flat | scaling converged; the error model is a **thermometer** (fits the scatter, doesn't reduce it) | +| beam-stop masking | anomalous unchanged | only ~5 low-res reflections affected on this well-behaved crystal | +| classical post-refinement | already done | cell (a=78.23, c=37.78 ≈ model), orientation (XtalOptimizer), and **mosaicity** (Ceres in `ScaleOnTheFly`) are already refined | + +**Radiation damage ruled out:** anomalous rises monotonically with dose (first 300 frames 0.24× → full +1800 0.61×) → redundancy-dominated, not damage-limited; per-frame B fits flat. + +**`-R100` is refinement, not indexing:** FFT saturates by ~33 first-pass images; the 10.5→11.0 gain is +refining the **global** geometry against more spots (30→50 populated images). **Per-image refinement is +essential** and fits real per-frame signal: `-r none` (global only) collapses to ISa 5.1 / 0.31×; +`-r orientation` 0.56×; default `BeamCenter` (per-image orient + beam) 11.0 / 0.61×. Finer + more DOF +wins — the per-image beam refit is really a per-frame position correction (crystal motion / +goniometer sphere-of-confusion). + +--- + +## Tooling (reusable instruments) + +- **`--dump-observations `** (`jfjoch_process`): `Combine3D` writes the unmerged fulls + (`h k l I σ d n_frames captured_fraction peak_frame`) for direct A/B vs `XDS_ASCII.HKL`. +- **`lyso_test/anomalous_scoreboard.sh [label]`**: one command → SHELXC + ANODE + under `qemu-x86_64-static` (the SHELX binaries hit the legacy vsyscall page; this WSL2 kernel is + `vsyscall=none`, so run them under qemu) → prints `⟨d″/σ⟩` per shell + averaged anomalous peak height + at SD_MET/SG_CYS/CL_CL vs a cached XDS reference. ~30 s/rerun. **This is the accuracy benchmark** — + anomalous ΔF is tiny, so it is exquisitely sensitive to systematic error, and it caught the outlier- + rejection accuracy trap that CC½/R-meas hid. +- **Per-reflection compare** (`ccp4-python` + gemmi, P43212 ASU): `CC(ours, XDS)` per shell and the + ours/XDS ratio-vs-resolution that produced the systematic-error map. +- **`[timing]` log lines** in the scaling/merge post-pass (per sub-step, plus the serial `MergeAll` vs + parallel `ScaleAllImages` split) — to find performance bottlenecks. + +--- + +## Performance notes + +- **Scaling is already multi-threaded.** `ScaleOnTheFly` parallelizes across images (all cores); each + per-image Ceres solve is single-threaded (correct grain). +- **Global (joint) scaling was removed** (2026-06-25): 122 s vs 3.5 s for the alternating loop (~35×) for + **no** quality gain (0.54 ≈ 0.53). Deleted everywhere (`GlobalScale.{h,cpp}`, the CLI flag, the + `ScalingSettings` field, the viewer checkbox). +- **The serial post-pass bottlenecks** (the real MT targets, if pursued): `Combine3D` (~1.8 s), + `MergeAll` (~1.4 s × scaling iterations), the final merge (~0.5 s). +- **Latent perf trap (open):** stats computation blows up to ~110 s when **no space group is set** + (huge unmerged set → `SearchSpaceGroup` auto-detection over all of it; with a SG, 0.14 s). Needs a + guard/cap — not yet touched. + +--- + +## Where the remaining gap is — honest assessment + +Box-sum → 2D profile (+15%) → 3D profile (+6%) shows integration improvements with **diminishing +returns**; integration alone cannot reach 1.0×. ANODE says we are unbiased, just ~2× noisier. Every +isolated *scaling* improvement (global fit, absorption surface, per-frame B, more iterations, error-model +tuning) was neutral or negative on this well-behaved crystal. The conclusion: **the residual ~1.8× is +XDS's full XSCALE/CORRECT machinery acting together** (multi-dimensional smooth scale + per-batch +corrections + tighter error model, fit jointly), whose individual pieces are each neutral here — a large +reimplementation with low odds of a quick win. + +**Candidate directions, with current expectation:** +1. **Weak high-res background modelling** (the +20–27% bias) — radial / profile-vs-summation, **not** a + linear gradient (the plane fit was neutral). A genuine target, but the simple fix doesn't work. +2. **Strong-reflection per-observation profile scatter** — the integration-quality floor; 2D profile-fit + already took the easy part, 3D profile gives a diminishing +6% and needs a 3D pixel shoebox (departs + from the per-frame philosophy; reader-blocked online). +3. **A learnable model.** The systematic-error map is, in effect, the "what to do" answer for an ML + integrator: the weak background bias, the strong-reflection saturation/response, and the per-observation + profile scatter are all **model-learnable** corrections — none is a simple hand-coded fix. + +**What is settled and should not be re-litigated:** we are unbiased (ANODE); 3D summation/background is +null; global scaling, absorption surface, per-batch geometry, learned rocking profile, per-frame B, +outlier rejection, partiality cuts, and beam-stop masking are all dead on this data. The lever that ever +worked is **2D profile-fit integration**; the keepers are **Gaussian profile-fit + rot3d + scale-fulls + +-R100** (≈ 0.61× anomalous, ISa 11.0). + +--- ## "Could PixelRefine be rebuilt as integrator + scaling?" -Yes: `ProfileIntegrate2D` (steps 1+2) + a scaling routine. The scaling routine can be the existing -reference-free `ScaleOnTheFly`, or a reference-driven one (Term 1) for the EM/self-seed loop — both -just consume `Reflection{I,σ,partiality}`. PixelRefine then becomes *= ProfileIntegrate2D + reference -ScaleOnTheFly*, and the experimental coupling (and its rot3d breakage) disappears. +Yes, and this is the clean decomposition that `ProfileIntegrate2D` realizes. PixelRefine does three +things; only the third needs a reference: +1. **Profile** — measure the tangential width per shell from strong spots (`PixelRefine.cpp:433`). + Reference-free. +2. **Extraction** — `J = Σ P(I−B)/v ÷ Σ P²/v`, de-biased variance, plus radial partiality + (`PixelRefine.cpp:491-525`). Reference-free. +3. **Scaling** — `J ≈ G·B_DW·p·pol·I_ref` (`PixelRefine.cpp:365`). **The only step using the reference.** + +So **integrator = steps 1+2** (= `ProfileIntegrate2D`), **scaling = step 3** (or the reference-free +`ScaleOnTheFly`). The decomposition also **fixes the rot3d incompatibility**: whole-PixelRefine emits +`image_scale_corr = 1/G` with no rotation partiality (which gave garbage σ with `-P rot3d`), whereas the +decomposed integrator emits `partiality` like `BraggIntegrate2D`, so `ScaleOnTheFly` + `Combine3D` + +merge consume it unchanged. PixelRefine then becomes `ProfileIntegrate2D + reference ScaleOnTheFly`, and +the experimental coupling disappears. diff --git a/image_analysis/bragg_integration/ProfileIntegrate2D.h b/image_analysis/bragg_integration/ProfileIntegrate2D.h index af00ffc6..14a501c3 100644 --- a/image_analysis/bragg_integration/ProfileIntegrate2D.h +++ b/image_analysis/bragg_integration/ProfileIntegrate2D.h @@ -4,7 +4,7 @@ #pragma once // ============================================================================= -// ProfileIntegrate2D — profile-fitting 2D integrator (DRAFT / not yet wired in) +// ProfileIntegrate2D — profile-fitting 2D integrator (the DEFAULT integrator) // ============================================================================= // // A drop-in alternative to BraggIntegrate2D that replaces uniform box summation with @@ -34,9 +34,9 @@ // decomposition + rot3d-combine compatibility); v2 = empirical per-shell (the expected win); // v3 = empirical per detector-region x shell (XDS-grade). // -// EXPERIMENTAL: selected by BraggIntegrationSettings::Integrator (ProfileGaussian = v1, -// ProfileEmpirical = v2); jfjoch_process exposes it as `--integrator gaussian|empirical` (default -// box-sum). For A/B vs XDS_ASCII.HKL via --dump-observations. Not yet exposed in the OpenAPI / viewer. +// Selected by BraggIntegrationSettings::Integrator: ProfileGaussian (the DEFAULT, v1) or +// ProfileEmpirical (v2); BoxSum (BraggIntegrate2D) is the fallback. jfjoch_process exposes it as +// `--integrator boxsum|gaussian|empirical`. For A/B vs XDS_ASCII.HKL via --dump-observations. // ============================================================================= #include