PixelRefine: make factored Terms 1+2 the model, remove old wiring
Build Packages / build:rpm (ubuntu2404_nocuda) (push) Successful in 25m0s
Build Packages / build:rpm (rocky8_nocuda) (push) Successful in 26m42s
Build Packages / build:rpm (rocky8_sls9) (push) Successful in 27m7s
Build Packages / build:rpm (ubuntu2204_nocuda) (push) Successful in 28m25s
Build Packages / build:rpm (rocky9_nocuda) (push) Successful in 29m44s
Build Packages / build:rpm (rocky9_sls9) (push) Successful in 32m14s
Build Packages / build:rpm (rocky8) (push) Successful in 24m39s
Build Packages / build:rpm (ubuntu2404) (push) Successful in 23m52s
Build Packages / build:rpm (ubuntu2204) (push) Successful in 25m19s
Build Packages / Generate python client (push) Successful in 23s
Build Packages / XDS test (durin plugin) (push) Successful in 20m42s
Build Packages / Create release (push) Skipped
Build Packages / build:rpm (rocky9) (push) Successful in 27m2s
Build Packages / Build documentation (push) Successful in 1m23s
Build Packages / DIALS test (push) Successful in 31m5s
Build Packages / XDS test (JFJoch plugin) (push) Successful in 14m55s
Build Packages / XDS test (neggia plugin) (push) Successful in 13m7s
Build Packages / Unit tests (push) Successful in 2h14m40s
Build Packages / build:rpm (ubuntu2404_nocuda) (push) Successful in 25m0s
Build Packages / build:rpm (rocky8_nocuda) (push) Successful in 26m42s
Build Packages / build:rpm (rocky8_sls9) (push) Successful in 27m7s
Build Packages / build:rpm (ubuntu2204_nocuda) (push) Successful in 28m25s
Build Packages / build:rpm (rocky9_nocuda) (push) Successful in 29m44s
Build Packages / build:rpm (rocky9_sls9) (push) Successful in 32m14s
Build Packages / build:rpm (rocky8) (push) Successful in 24m39s
Build Packages / build:rpm (ubuntu2404) (push) Successful in 23m52s
Build Packages / build:rpm (ubuntu2204) (push) Successful in 25m19s
Build Packages / Generate python client (push) Successful in 23s
Build Packages / XDS test (durin plugin) (push) Successful in 20m42s
Build Packages / Create release (push) Skipped
Build Packages / build:rpm (rocky9) (push) Successful in 27m2s
Build Packages / Build documentation (push) Successful in 1m23s
Build Packages / DIALS test (push) Successful in 31m5s
Build Packages / XDS test (JFJoch plugin) (push) Successful in 14m55s
Build Packages / XDS test (neggia plugin) (push) Successful in 13m7s
Build Packages / Unit tests (push) Successful in 2h14m40s
PixelRefine is now an intensity-only operation: geometry is fixed (refined upstream by XtalOptimizer) and the only objective is the factored per-reflection likelihood (FACTORED_MODEL.md Terms 1+2) - measured per-resolution profile width R1 plus one Fisher-weighted intensity/scaling residual per reflection, fitting the per-image scale G and B. Validated on crystal 2 (fixed_master.h5 as stills, 1.7 A): CC1/2 84-92%, CCref 77-92%, flat - reproduces the env-flag prototype and matches the rotation path from the stills path. Removed: - the per-pixel ShoeboxResidual loss and PixelResidual cost functor; - all in-PixelRefine geometry refinement (orientation/cell/beam/distance/R), the regularised-orientation LSQ, signal-weighting, and the global sweep; - Term 3 (per-spot recentring) - a confirmed no-op on both crystals; - the diagnostic scaffolding (covariance, centroid, adaptive_R1) and the PR_* env knobs + stderr dumps in IndexAndRefine; - the PredictImage/ChiSquaredImage renderers and the entire viewer PixelRefine window/table/params + worker bindings + shoebox overlay. The sweep box-integrator background median became mean (consistency) by virtue of removing the sweep. METHODS.md rewritten for the current model; findings recorded in FINDINGS-2026-06.md. Net -2200 lines. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
Claude Opus 4.8
@@ -446,28 +446,6 @@ bool IndexAndRefine::PixelRefineIntegrate(DataMessage &msg,
|
||||
if (const auto bw = experiment.GetBandwidthFWHM())
|
||||
prd.bandwidth = bw.value() / 2.3548; // FWHM -> sigma
|
||||
|
||||
// TEMPORARY diagnostic knobs to probe the effect of the (currently fixed) spot
|
||||
// widths R[0] (radial/partiality) and R[1] (tangential/profile). Remove after.
|
||||
if (const char *r0 = std::getenv("PR_R0")) prd.R[0] = std::stod(r0);
|
||||
if (const char *r1 = std::getenv("PR_R1")) prd.R[1] = std::stod(r1);
|
||||
// TEMPORARY: PR_COV refines G,B,R and dumps their per-image correlation matrix.
|
||||
if (std::getenv("PR_COV")) {
|
||||
prd.refine_scale = true; prd.refine_B = true; prd.refine_R = true;
|
||||
prd.compute_covariance = true;
|
||||
}
|
||||
if (std::getenv("PR_FIX_R0")) prd.fix_R0 = true; // hold R0, refine R1 only
|
||||
if (std::getenv("PR_FIX_R")) prd.refine_R = false; // hold R0 and R1: G-B correlation only
|
||||
if (std::getenv("PR_ADAPT_R1")) prd.adaptive_R1 = true; // measure R1 from spot moments
|
||||
if (std::getenv("PR_CENTROID")) prd.measure_centroid = true; // observed-vs-predicted offset
|
||||
if (std::getenv("PR_RECENTER")) prd.recenter_profile = true; // recentre profile on centroid
|
||||
if (const char *s = std::getenv("PR_RECENTER_SIGNIF")) prd.recenter_min_signif = std::stod(s);
|
||||
if (std::getenv("PR_INTENSITY")) { // factored-likelihood Term 1: per-reflection intensity residual
|
||||
prd.intensity_residual = true;
|
||||
prd.refine_orientation = false; prd.refine_R = false;
|
||||
prd.refine_scale = true; prd.refine_B = true;
|
||||
}
|
||||
if (std::getenv("PR_SHAPE")) prd.shape_R1 = true; // Term 2: per-resolution R1 from spot moments
|
||||
|
||||
std::vector<uint8_t> buffer;
|
||||
const uint8_t *ptr = image.GetUncompressedPtr(buffer);
|
||||
switch (image.GetMode()) {
|
||||
@@ -487,18 +465,6 @@ bool IndexAndRefine::PixelRefineIntegrate(DataMessage &msg,
|
||||
return false;
|
||||
}
|
||||
|
||||
if (prd.covariance_valid)
|
||||
fprintf(stderr, "[cov] GB=%.3f GR0=%.3f GR1=%.3f BR0=%.3f BR1=%.3f R0R1=%.3f\n",
|
||||
prd.corr_GB, prd.corr_GR0, prd.corr_GR1, prd.corr_BR0, prd.corr_BR1, prd.corr_R0R1);
|
||||
if (prd.adaptive_R1)
|
||||
fprintf(stderr, "[R1] %.5f\n", prd.R[1]);
|
||||
if (prd.shape_R1 && std::isfinite(prd.shape_R1_lores))
|
||||
fprintf(stderr, "[shapeR1] lores=%.5f hires=%.5f\n", prd.shape_R1_lores, prd.shape_R1_hires);
|
||||
if (prd.measure_centroid && std::isfinite(prd.centroid_lo_tang_c))
|
||||
fprintf(stderr, "[res] lo_s=%.1f lo_tc=%.3f lo_tp=%.3f lo_rc=%.3f hi_s=%.1f hi_tc=%.3f hi_tp=%.3f hi_rc=%.3f\n",
|
||||
prd.centroid_lo_signif, prd.centroid_lo_tang_c, prd.centroid_lo_tang_p, prd.centroid_lo_rad_c,
|
||||
prd.centroid_hi_signif, prd.centroid_hi_tang_c, prd.centroid_hi_tang_p, prd.centroid_hi_rad_c);
|
||||
|
||||
// PixelRefine output flows into the normal save/merge path: the refined
|
||||
// geometry/lattice and the already-scaled reflections become the outcome.
|
||||
i_outcome.reflections = std::move(prd.reflections);
|
||||
|
||||
@@ -1,6 +1,8 @@
|
||||
# A factored likelihood for joint integration + scaling + geometry
|
||||
|
||||
**Status: design spec (not implemented).** Goal: replace the per-pixel least-squares
|
||||
**Status: Terms 1+2 implemented and shipping as the PixelRefine default** (see
|
||||
`PixelRefine.cpp`, `METHODS.md`, `FINDINGS-2026-06.md`); Term 3 (geometry) and the
|
||||
priors/NN extensions of §4 remain future work. Goal: replace the per-pixel least-squares
|
||||
of PixelRefine with a per-*reflection* likelihood that fuses profile-fit integration,
|
||||
scaling against the reference, and geometry refinement into one differentiable
|
||||
objective — the foundation for priors (Bayesian) and learned components (NN), and the
|
||||
|
||||
@@ -0,0 +1,122 @@
|
||||
# PixelRefine — findings, June 2026
|
||||
|
||||
A record of the main results from the still-integration investigation, validated on two
|
||||
lysozyme P4₃2₁2 datasets indexed against the same reference (`6G8A_refine_001.mtz`):
|
||||
|
||||
* **Crystal 1 — serial jet stills** (`LysozymeJet5-…`, MicroMAX DMM ~1 % bandwidth).
|
||||
* **Crystal 2 — rotation crystal** (`fixed_master.h5`, 1800 × 0.2°, Si mono), run *as
|
||||
stills* as a stress test, and as rotation for reference. XDS output of this crystal
|
||||
(`CORRECT.LP`, `XDS_ASCII.HKL`) is the external ground truth.
|
||||
|
||||
---
|
||||
|
||||
## 1. The factored model is a qualitative win (crystal 2, stills, 1.7 Å)
|
||||
|
||||
Replacing the per-pixel least squares with the factored Terms 1+2 (intensity residual +
|
||||
measured per-resolution $R_1$) turns an erratic, high-res-collapsing result into a flat one:
|
||||
|
||||
| | N_obs | ⟨I/σ⟩ | CC₁/₂ (per shell) | CCref (per shell) |
|
||||
|---|---:|---:|---|---|
|
||||
| baseline per-pixel loss | 799 k | 7.2 | 75→81→56→32→…→0 % | erratic, →0 |
|
||||
| **factored Terms 1+2** | 1.22 M | 10.7 | **84–92 %, flat to 1.7 Å** | **77–92 %, flat** |
|
||||
|
||||
This *reaches the proper rotation (`-R -P rot`) path's quality from the stills path.* The
|
||||
mechanism is clean and on-thesis: **Term 1 lifts the per-image scale CC vs the reference
|
||||
from 0.09 → 0.40** (median 0.01 → 0.39) — the per-pixel loss produced near-random per-image
|
||||
scales; the factored objective makes each image's scale agree with the reference, so the
|
||||
merge coheres. The reference enters at maximum leverage (it sets both the target and the
|
||||
Fisher weight).
|
||||
|
||||
**Term 3 (per-spot recentring) is a no-op** on both crystals (93.5 vs 93.6 % CC₁/₂): the
|
||||
generous integration aperture already contains the spot, so shifting the spotlight by the
|
||||
sub-pixel prediction offset does not change the integrated total. It was removed; the
|
||||
position residual's value can only come as a *geometry* term (refine orientation/distance
|
||||
from the centroid), which is XtalOptimizer's job, not a per-spot mask shift.
|
||||
|
||||
---
|
||||
|
||||
## 2. The residual centroid offset is a sampling floor, not a recoverable error
|
||||
|
||||
On the jet, predictions land on the spot to within a ~0.4 px tangential scatter. Three
|
||||
independent cuts show this is the **centroid undersampling floor of the ~2×2 spots**, not a
|
||||
geometry error or a shape effect:
|
||||
|
||||
* **Flat with intensity.** Binned by significance (≈6σ vs ≈39σ, a 6× span) the tangential
|
||||
centroid offset is 0.41 → 0.36 px (median). A background-limited centroid would shrink
|
||||
∝ 1/significance (≈6×); it does not, so it is a real ~0.35 px floor, not counting noise.
|
||||
* **Peak is a *worse* predictor than the centroid** (sub-pixel parabolic mode: 0.52 vs 0.36
|
||||
px at high signif). If an asymmetric tail were dragging the centroid off a correct peak,
|
||||
the peak would be *better* — it is not. So there is **no coherent shape asymmetry** to
|
||||
model; the offset is a zero-mean per-spot scatter from sub-pixel phase aliasing of a spot
|
||||
only ~2 px wide.
|
||||
* **Radial centroid offset ≈ 0.01–0.02 px, flat** → no distance/parallax error (parallax is
|
||||
radial and grows outward).
|
||||
|
||||
Consequences: recentring cannot rescue weak reflections (their prediction is already
|
||||
centred; the 0.4 px is irreducible sampling scatter), which is *why* a generous box beats a
|
||||
tight profile mask. A *symmetric* non-Gaussian shape would not move the centroid at all — it
|
||||
would show up only in Term 2 ($R_1$), not in a position term.
|
||||
|
||||
---
|
||||
|
||||
## 3. The σ gap to XDS is fulls-vs-partials, not intensities
|
||||
|
||||
XDS reports **ISa = 28.3** (asymptotic relative error 3.5 %); our stills path reports
|
||||
**ISa ≈ 1.1** and the rotation path **≈ 1.6**. The decisive clue is in `XDS_ASCII.HKL`: the
|
||||
`PEAK` column (fraction of each reflection captured) is **≈ 100 % for 97 % of reflections**,
|
||||
and the header gives mosaicity 0.091° < 0.2° oscillation. So **these are full reflections**,
|
||||
recorded over the 1–3 frames each rocking curve spans — not partials.
|
||||
|
||||
* **XDS** profile-fits in 3D (the third axis is the rotation/rocking direction) and *sums*
|
||||
the rocking curve with profile weights → one full per reflection, counting-limited σ.
|
||||
* **jfjoch `-P rot`** integrates a 2D shoebox *per frame* and recovers each full by
|
||||
*dividing* by the rocking-curve fraction $R_j$ — a cheap approximation of 3D integration.
|
||||
That division injects per-observation noise (it amplifies each frame's background noise,
|
||||
pays N independent backgrounds, and carries a random per-observation partiality error).
|
||||
|
||||
Crucially, **ISa and merged-intensity accuracy are different axes, decoupled by
|
||||
multiplicity.** Our merged intensities are correct (§4), because ~60–240× multiplicity
|
||||
averages the per-observation noise down; ISa measures the per-*observation* precision, which
|
||||
multiplicity cannot improve. So "right intensities, wrong σ" is not a contradiction.
|
||||
|
||||
**A cheap probe confirms it.** Raising `--min-partiality` from 0.02 → 0.5 on crystal 2
|
||||
(rotation) lifts ISa **1.6 → 3.8** at *zero* completeness cost (high multiplicity) and with
|
||||
CCref flat — and the high-res shells *improve*. The default 0.02 keeps deep-tail partials
|
||||
(2 % of a reflection, scaled back ×50) that were over-weighted and polluting the merge. So
|
||||
of the ~17× rotation gap, **~2.8× is tunable tail-weighting** and **~6× is structural** (the
|
||||
2D-divide vs 3D-sum difference) — the structural part needs a real rocking-curve sum, not a
|
||||
knob.
|
||||
|
||||
---
|
||||
|
||||
## 4. Our intensities are XDS-grade — the limit is σ, not I
|
||||
|
||||
Direct CC of merged intensities (both reduced to the 4/mmm asymmetric unit; 98.8 % of
|
||||
reflections matched, no manual reindex):
|
||||
|
||||
| vs XDS (CC on merged I, to 1.2 Å) | overall | at 1.2 Å |
|
||||
|---|---|---|
|
||||
| **PixelRefine stills** (factored) | **95.9 %** | 96.9 % |
|
||||
| rotation `-P rot` (min_partiality 0.3) | 98.8 % | 98.5 % |
|
||||
|
||||
PixelRefine-stills intensities track XDS at **95–98 %, flat to the 1.2 Å diffraction
|
||||
limit**, only ~3 % behind full rotation integration. This is exactly what the CC₁/₂→CC_true
|
||||
relation predicts (stills CC₁/₂≈0.88 ⇒ 0.967), so it is real, carried by the huge stills
|
||||
multiplicity averaging per-observation noise down. **Conclusion: the intensity estimation is
|
||||
right; the remaining gap to XDS is entirely the per-observation σ / partiality axis** — the
|
||||
3D rocking-curve integration frontier (`FACTORED_MODEL.md` §4), deliberately parked.
|
||||
|
||||
---
|
||||
|
||||
## What is fixed vs. parked
|
||||
|
||||
**Landed (and now the default model):** mean (not median) background in both integrators;
|
||||
de-biased (background-limited) variance; widened Ewald prediction band; the factored
|
||||
Terms 1+2 as the only PixelRefine objective; the global XDS-form merge error model with ISa.
|
||||
|
||||
**Diagnosed dead-ends (do not re-litigate):** per-image $R$ refinement (degenerate with
|
||||
scale); per-spot recentring (no-op); chasing the 0.4 px centroid floor (sampling, not
|
||||
recoverable). **Parked (rotation-side):** the 3D rocking-curve sum / one-shot
|
||||
post-refinement, a `min_partiality` default for high-multiplicity data, and
|
||||
partiality-aware variance weighting — the path to XDS-grade ISa, but a separate axis from
|
||||
the intensity work.
|
||||
@@ -1,31 +1,30 @@
|
||||
# PixelRefine — methods and improvements
|
||||
# PixelRefine — methods
|
||||
|
||||
This note documents the changes made to the still-image *pixel-refinement* integrator
|
||||
(`PixelRefine`) and, more importantly, **why** each one was needed. It is written from a
|
||||
methods point of view; the equations are the load-bearing part.
|
||||
`PixelRefine` is the still-image integrator. It integrates the Bragg reflections of one
|
||||
image by **profile fitting against a reference intensity set** $I^\mathrm{ref}$ (e.g.
|
||||
`F_calc` from a deposited model, or the current merged estimate in an EM-style outer loop)
|
||||
and returns already-scaled intensities. It is an **intensity-wise** operation: the
|
||||
detector geometry (orientation, cell, beam, distance) is taken as fixed — it was refined
|
||||
upstream by `XtalOptimizer` (`IndexAndRefine::RefineGeometryIfNeeded`) — and PixelRefine
|
||||
only measures the spot shape and fits the per-image scale.
|
||||
|
||||
`PixelRefine` integrates Bragg reflections on a **still** image by fitting a per-pixel
|
||||
forward model against a known *reference* intensity set (e.g. `F_calc` from a deposited
|
||||
model). Unlike a rotation experiment, a still samples only a thin slice of each
|
||||
reflection, so the integrator must (a) model the partiality of that slice, (b) refine the
|
||||
per-image geometry well enough that the high-resolution shoeboxes land on signal, and
|
||||
(c) scale each image onto the reference. Each of those three is where the original code
|
||||
went wrong.
|
||||
The objective is the factored per-reflection likelihood of `FACTORED_MODEL.md`, **Terms 1
|
||||
and 2**. This note records the equations and the reasons behind each design choice.
|
||||
|
||||
Throughout, a reflection's shoebox is a small box of raw detector pixels $I_p$ with a
|
||||
local flat background $B$; the (area-normalised) model spot profile at pixel $p$ is
|
||||
$P_p$, and $v_p$ is the variance used to weight pixel $p$.
|
||||
local flat background $B$; the area-normalised tangential profile at pixel $p$ is $P_p$,
|
||||
and $v_p$ is the variance used to weight pixel $p$.
|
||||
|
||||
---
|
||||
|
||||
## 0. The forward model
|
||||
|
||||
The recorded amplitude of a still reflection is estimated by profile fitting:
|
||||
The recorded amplitude of a still reflection is the profile-fit amplitude
|
||||
|
||||
$$
|
||||
J \;=\; \frac{\sum_p w_p\,P_p\,(I_p - B)}{\sum_p w_p\,P_p^{2}},
|
||||
\qquad w_p = \frac{1}{v_p},
|
||||
\qquad \operatorname{var}(J) = \frac{1}{\sum_p P_p^{2}/v_p}.
|
||||
J \;=\; \frac{\sum_p P_p\,(I_p - B)/v_p}{\sum_p P_p^{2}/v_p},
|
||||
\qquad
|
||||
\operatorname{var}(J) = \frac{1}{\sum_p P_p^{2}/v_p}.
|
||||
$$
|
||||
|
||||
The full (rotation-equivalent) intensity is recovered by dividing out the factors a still
|
||||
@@ -33,296 +32,183 @@ does not record,
|
||||
|
||||
$$
|
||||
I \;=\; \frac{J}{p\,B_\mathrm{DW}\,\mathrm{pol}},\qquad
|
||||
p = \exp\!\left(-\frac{\epsilon_r^{2}}{R_0^{2}}\right),\quad
|
||||
p = \exp\!\left(-\frac{\epsilon_r^{2}}{R_{0,\mathrm{eff}}^{2}}\right),\quad
|
||||
B_\mathrm{DW}=\exp\!\left(-\frac{B_\mathrm{fac}}{4 d^{2}}\right),
|
||||
$$
|
||||
|
||||
where the **partiality** $p$ is the fraction of the mosaic block crossing the Ewald
|
||||
sphere, $\epsilon_r = 1/\lambda - |S_{hkl}|$ is the excitation error, $R_0$ the radial
|
||||
(rocking) width, and $\mathrm{pol}$ the polarisation correction. The per-image model
|
||||
intensity for a pixel is
|
||||
sphere, $\epsilon_r$ is the radial excitation error, $R_0$ the radial (rocking) width, and
|
||||
$\mathrm{pol}$ the polarisation correction. The tangential profile is a separable,
|
||||
area-normalised Gaussian of width $R_1$:
|
||||
|
||||
$$
|
||||
I_p^\mathrm{model} = G\,I^\mathrm{ref}\,B_\mathrm{DW}\,p\,P^\mathrm{tang}_p\,\mathrm{pol} + B,
|
||||
P_p = \frac{1}{\pi R_1^{2}}\exp\!\left(-\frac{\epsilon_{t,p}^{2}}{R_1^{2}}\right).
|
||||
$$
|
||||
|
||||
with $G$ the per-image scale. The per-image least squares minimises
|
||||
$\chi^2 = \sum_p w_p\,(I_p^\mathrm{model}-I_p)^2$ over geometry, orientation, $G$, $R$.
|
||||
A finite **X-ray bandwidth** thickens the Ewald shell radially, adding a fixed,
|
||||
resolution-dependent term to the radial width, $R_{0,\mathrm{eff}}^2 = R_0^2 +
|
||||
(b\lambda)^2/(2d^4)$ ($b$ = relative bandwidth; the pink-beam / DMM signature). $b=0$ is a
|
||||
monochromatic no-op.
|
||||
|
||||
---
|
||||
|
||||
## 1. De-biased variance (the load-bearing fix)
|
||||
|
||||
**Symptom.** Mean intensities went **negative** in the high-resolution shells
|
||||
($\langle I/\sigma\rangle$ down to $-12$), which a box-sum integrator never does. The
|
||||
per-image scale $G$ also collapsed to $0$ on most images, dropping ~80 % of observations.
|
||||
($\langle I/\sigma\rangle$ down to $-12$), and the per-image scale $G$ collapsed to $0$ on
|
||||
most images, dropping ~80 % of observations.
|
||||
|
||||
**Cause.** Both the extraction and the fit weighted each pixel by its **observed** count,
|
||||
$v_p = I_p$. For a background pixel that fluctuated *down* ($I_p < B$), $v_p$ is small, so
|
||||
$w_p = 1/v_p$ is *large*, and its contribution $P_p (I_p-B)/v_p < 0$ is large in
|
||||
magnitude. Summed over the many (mostly empty) shoebox pixels, this drags $J$ below zero —
|
||||
the classic *inverse-observed-count* (Poisson-on-data) bias. It bites hardest where the
|
||||
true signal is weakest, i.e. at high resolution. In the fit it manifests differently but
|
||||
identically in origin: the weighted empty pixels make "no signal" ($G=0$) the cheapest
|
||||
solution, so $G\to 0$.
|
||||
**Cause.** The extraction weighted each pixel by its **observed** count, $v_p = I_p$. A
|
||||
down-fluctuated background pixel ($I_p < B$) then gets a small $v_p$, hence a large
|
||||
$w_p=1/v_p$, and its contribution $P_p(I_p-B)/v_p < 0$ is large in magnitude. Summed over
|
||||
the many empty shoebox pixels this drags $J$ below zero — the *inverse-observed-count*
|
||||
(Poisson-on-data) bias, worst where the true signal is weakest (high resolution).
|
||||
|
||||
**Fix.** For background-limited (weak) reflections the correct variance is the local
|
||||
background, **constant over the shoebox**:
|
||||
|
||||
$$
|
||||
v_p = \max(B,\,1)\quad\Longrightarrow\quad
|
||||
J = \frac{\sum_p P_p\,(I_p-B)}{\sum_p P_p^{2}},
|
||||
$$
|
||||
|
||||
the unbiased uniform-variance estimator. This single change turned $\langle I/\sigma\rangle$
|
||||
positive at all resolutions and stopped the scale collapse. It is applied to both the
|
||||
extraction weight and the fit weight.
|
||||
**Fix.** For background-limited reflections the correct variance is the local background,
|
||||
**constant over the shoebox**, $v_p = \max(B,1)$, giving the unbiased uniform-variance
|
||||
estimator $J = \sum_p P_p (I_p-B)/\sum_p P_p^2$. This turned $\langle I/\sigma\rangle$
|
||||
positive at all resolutions and stopped the scale collapse.
|
||||
|
||||
---
|
||||
|
||||
## 2. Prediction band and multiplicity
|
||||
|
||||
**Symptom.** PixelRefine recorded ~4× fewer observations per unique reflection than the
|
||||
classical integrator — completeness was fine, *redundancy* was not.
|
||||
|
||||
**Cause.** A reflection is given a shoebox only when it lies within a radial band of the
|
||||
Ewald sphere,
|
||||
|
||||
$$
|
||||
\bigl|\,|S_{hkl}| - 1/\lambda\,\bigr| \le \delta .
|
||||
$$
|
||||
|
||||
For randomly oriented stills the number of images on which a given $hkl$ satisfies this is
|
||||
$\propto \delta$. The default $\delta = 5\times10^{-4}\,\text{Å}^{-1}$ was 4–6× tighter
|
||||
than the classical integrator's $\delta = 2\text{–}3\times r_\mathrm{profile}$, so each
|
||||
reflection was recorded on 4–6× fewer images.
|
||||
|
||||
**Fix.** Widen to $\delta = 2\times10^{-3}\,\text{Å}^{-1}$. Multiplicity rose from
|
||||
~240 k to ~950 k observations and CC$_\mathrm{ref}$ from 49.7 % to 55.9 %. Widening is
|
||||
only safe once the fit is well-behaved (Sections 1, 3, 4); with the original
|
||||
unconstrained fit it caused divergence.
|
||||
A reflection is given a shoebox only when it lies within a radial band of the Ewald
|
||||
sphere, $\bigl|\,|S_{hkl}| - 1/\lambda\,\bigr| \le \delta$. For randomly oriented stills
|
||||
the number of images on which a given $hkl$ qualifies is $\propto \delta$. The original
|
||||
$\delta = 5\times10^{-4}\,\text{Å}^{-1}$ was 4–6× tighter than a box integrator, giving 4×
|
||||
fewer observations per reflection. Widening to $\delta = 2\times10^{-3}\,\text{Å}^{-1}$
|
||||
(`ewald_dist_cutoff`) restores the multiplicity; the partiality $p$ downweights the
|
||||
slightly-off-Ewald tails it admits. (Widening is only safe with the de-biased variance of
|
||||
§1 and the factored objective of §§3–4 — with a per-pixel geometry fit it diverged.)
|
||||
|
||||
---
|
||||
|
||||
## 3. Regularising the per-image fit
|
||||
## 3. Term 2 — measured per-resolution profile width $R_1$
|
||||
|
||||
**Symptom.** Freeing *any* per-image parameter (orientation, $R$, even the scalar scale
|
||||
$G$) collapsed the merged data (CC$_{1/2}$ from 90 % to a few %). The predict↔refine loop
|
||||
with all parameters frozen was, by contrast, byte-identical to extraction-only — proving
|
||||
the *fit*, not the loop, was at fault.
|
||||
|
||||
**Cause.** The per-image problem regularised **nothing**: orientation had no prior, $R$
|
||||
and $G$ only a lower bound. Three orientation DOF (plus $R$, $G$) against a handful of
|
||||
signal pixels per still overfit the noise, and an unconstrained $1/G$ then scrambled the
|
||||
cross-image merge.
|
||||
|
||||
**Fix.** Anchor each refined parameter to its prior with a *data-scaled* weight, as
|
||||
`ScaleOnTheFly` already does for rotation data. The data term has one residual **per
|
||||
pixel**, so the prior weight must scale with the pixel count:
|
||||
$R_1$ is **measured, not fitted**. Fitting $R_1$ inside a per-image least squares is
|
||||
degenerate with the scale $G$ (a narrower profile and a larger scale trade off), and that
|
||||
degeneracy slides the per-image scale and wrecks the merge. But a **second moment** is
|
||||
normalised by the total intensity, so it carries shape information *decoupled from scale*:
|
||||
|
||||
$$
|
||||
w_\theta = \sqrt{\frac{N_\mathrm{pix}}{\sigma_\theta^{2}}}\,,
|
||||
\qquad \chi^2_\mathrm{reg} = w_\theta^{2}\,(\theta-\theta_0)^2 .
|
||||
$$
|
||||
|
||||
Using $\sqrt{N_\mathrm{refl}}$ instead (as in the rotation scaler) is a factor
|
||||
$\sqrt{N_\mathrm{pix}/N_\mathrm{refl}}\approx\sqrt{49}\approx 7$ too weak and is simply not
|
||||
felt. Applied to:
|
||||
|
||||
* **Scale**: $\theta=G,\ \theta_0=1$. Prevents $1/G$ from wandering; restored CC$_{1/2}$
|
||||
from 4 % back to 87 %.
|
||||
* **Orientation**: $\theta$ = Rodrigues vector, $\theta_0$ = spot-centroid orientation,
|
||||
$\sigma_\theta$ in radians. At $\sigma_\theta\!\sim\!1^\circ$ this gives the best
|
||||
CC$_\mathrm{ref}$; beyond $\sim 2^\circ$ the fit overfits and the merge collapses, so
|
||||
$\sigma_\theta$ is the safety knob.
|
||||
|
||||
---
|
||||
|
||||
## 4. Signal-weighting the fit
|
||||
|
||||
**Cause.** Even de-biased, a shoebox is ~80 % empty pixels (≈ 40 of 49 for a radius-3
|
||||
box). They carry no information on $G$, $R$ or orientation but add noise and, near the
|
||||
overfitting edge, destabilise the fit.
|
||||
|
||||
**Fix.** Multiply the fit weight by a detector-space Gaussian centred on the predicted
|
||||
spot,
|
||||
|
||||
$$
|
||||
w_p \;\to\; w_p \cdot \exp\!\left(-\frac{r_p^{2}}{2\sigma_s^{2}}\right),
|
||||
\qquad r_p = \lVert (x_p,y_p) - (x_\mathrm{pred},y_\mathrm{pred})\rVert,
|
||||
$$
|
||||
|
||||
so the signal-bearing core drives the refined parameters ($\sigma_s\approx 1.5$ px). With
|
||||
the regularised orientation this lifted CC$_\mathrm{ref}$ from 60.5 % to **62.6 %**.
|
||||
|
||||
---
|
||||
|
||||
## 5. Global orientation + cell-scale sweep
|
||||
|
||||
**Motivation.** The classical refinement (`XtalOptimizer`) works on spot *centroids* from
|
||||
spot-finding, which are poor value at high resolution; the local LSQ above is a heavy
|
||||
gradient step that cannot make the ~degree-scale global moves needed to pull a
|
||||
mis-indexed crystal's high-resolution reflections onto their shoeboxes. A small,
|
||||
**global** sweep that simply asks *"at which orientation does the most high-resolution
|
||||
signal appear where the reference says it should?"* is structurally better suited.
|
||||
|
||||
**Score.** Pearson CC of the box-summed intensities against the reference, over **all**
|
||||
matched reflections (not just the strong ones):
|
||||
|
||||
$$
|
||||
\mathrm{CC}_\mathrm{ref}(L) = \operatorname{corr}_{hkl}\bigl(\,I^\mathrm{box}_{hkl}(L),\ I^\mathrm{ref}_{hkl}\,\bigr).
|
||||
$$
|
||||
|
||||
The strong low-resolution reflections **anchor** the CC (moving them off their boxes
|
||||
collapses it), so the *change* in CC across the sweep is driven almost entirely by weak
|
||||
high-resolution reflections falling onto — or off — real signal. This is the "appearing
|
||||
out of the void" behaviour.
|
||||
|
||||
**Geometry-derived bounds (parameter-free).** A spot at resolution $d$ sits at detector
|
||||
radius
|
||||
|
||||
$$
|
||||
r(d) = \frac{L}{p}\,\frac{\lambda}{d}\quad[\text{px}],
|
||||
$$
|
||||
|
||||
($L$ = detector distance, $p$ = pixel size). Both a crystal rotation $\delta\theta$ and a
|
||||
fractional cell-scale $\epsilon$ displace that spot by an amount proportional to its
|
||||
radius:
|
||||
|
||||
$$
|
||||
\Delta_\mathrm{px} = r(d)\,\delta\theta \quad(\text{rotation}),\qquad
|
||||
\Delta_\mathrm{px} = r(d)\,\epsilon \quad(\text{cell scale}).
|
||||
$$
|
||||
|
||||
So high-resolution spots (large $r$) move most. The two natural constraints fix the grid:
|
||||
|
||||
* **Step** = 1 px at the highest resolution: $\;\delta\theta_\mathrm{step} = 1/r_\mathrm{max}$
|
||||
(finer is below the detector's resolving power).
|
||||
* **Range** = the orientation uncertainty $\Delta\theta_u$ (a few px at high res). The
|
||||
number of steps per axis is then
|
||||
|
||||
$$
|
||||
n = \frac{\Delta\theta_u}{\delta\theta_\mathrm{step}} = \Delta\theta_u\, r_\mathrm{max},
|
||||
$$
|
||||
|
||||
a handful of steps, and the lowest-resolution spots move only
|
||||
$n\,\delta\theta_\mathrm{step}\,r_\mathrm{min} = \Delta\theta_u\, r_\mathrm{min} \ll 2$ px
|
||||
— i.e. the strong anchors stay put by construction.
|
||||
|
||||
> **Lesson learned.** Setting the range from "2 px at low resolution" instead gives
|
||||
> $n = 2\,r_\mathrm{max}/r_\mathrm{min} = 2\,d_\mathrm{low}/d_\mathrm{high}$ steps — tens
|
||||
> of pixels of high-resolution freedom — which lets the per-image CC overfit and *degrades*
|
||||
> the merge. The range must be tied to the (small) orientation uncertainty, not the
|
||||
> low-resolution cap.
|
||||
|
||||
The sweep is a coordinate descent over the three Rodrigues axes and the cell scale, run
|
||||
**before** the LSQ. It improves the high-resolution shells (CC$_\mathrm{ref}$ +1 to +5
|
||||
per shell) while preserving CC$_{1/2}\approx 90$ %. It is an *alternative* to the
|
||||
loose-$\sigma$ orientation LSQ, not a complement — stacking the two double-moves the
|
||||
orientation and overfits. It is therefore available but **off by default**.
|
||||
|
||||
---
|
||||
|
||||
## 6. Background-estimator bias (the largest σ error; both integrators)
|
||||
|
||||
**Symptom.** Pushed past the true resolution limit (e.g. a 1.8 Å crystal merged to 1.3 Å),
|
||||
the no-signal shells still reported $\langle I/\sigma\rangle \approx 4\text{–}6$ with
|
||||
$\mathrm{CC}_{1/2}\approx 0$ — i.e. confident "data" where there is none. Present in *both*
|
||||
`PixelRefine` and the classical `BraggIntegrate2D`.
|
||||
|
||||
**Cause.** Both estimated the local background as the **median** of the surrounding ring.
|
||||
For a Poisson / right-skewed background the median sits *below* the mean,
|
||||
|
||||
$$
|
||||
\operatorname{median}(B) < \mathbb{E}[B],
|
||||
$$
|
||||
|
||||
so subtracting it under-subtracts on every pixel. The leftover positive offset is tiny per
|
||||
pixel but coherent, so over an $n_\mathrm{pix}$-pixel peak and a multiplicity-$m$ merge it
|
||||
grows to a fake signal
|
||||
|
||||
$$
|
||||
\langle I\rangle_\mathrm{bias} \;\approx\; n_\mathrm{pix}\,\bigl(\mathbb{E}[B]-\operatorname{median}(B)\bigr),
|
||||
R_1^2 = 2\,\langle \epsilon_t^2\rangle,
|
||||
\qquad
|
||||
\Bigl(\tfrac{I}{\sigma}\Bigr)_\mathrm{merged,\,bias} \propto \sqrt{m}.
|
||||
\langle \epsilon_t^2\rangle = \frac{\sum_p (I_p-B)\,\epsilon_{t,p}^2}{\sum_p (I_p-B)} .
|
||||
$$
|
||||
|
||||
It is worst where the real signal is weakest (high resolution), because there the offset is
|
||||
all that remains — which is exactly the observed signature. *Nothing leaks into the
|
||||
high-resolution shells; the background is simply under-estimated.*
|
||||
|
||||
**Fix.** Use the **mean** of the ring (outliers excluded by the existing spot-core mask and
|
||||
saturation sentinels). $\langle I/\sigma\rangle$ then collapses to ~0 wherever
|
||||
$\mathrm{CC}\approx0$, tracks CC down the shells, and the honest resolution limit becomes
|
||||
visible. This was the single largest contributor to untrustworthy σ — a one-line change in
|
||||
each integrator, *not* a variance-model problem.
|
||||
We bin the strong spots ($\mathrm{signif}\ge 5$) by resolution ($1/d^2$, 6 bins) and take
|
||||
the **median** $\langle\epsilon_t^2\rangle$ per bin, so each reflection integrates with the
|
||||
$R_1$ of its resolution shell (low-res spots are tight; high-res anisotropic streaks are
|
||||
wider). Weak spots fall back to the global $R_1$. Measuring the width rather than fitting
|
||||
it is what makes profile-width refinement stable — and it is a selling point: the mask
|
||||
adapts to the data per shell.
|
||||
|
||||
---
|
||||
|
||||
## 7. Error model (global $a, b$; XDS form)
|
||||
## 4. Term 1 — the intensity / scaling residual
|
||||
|
||||
The per-pixel least squares is replaced by **one residual per reflection**: the profile-fit
|
||||
amplitude $J$ (using the Term-2 $R_1$) should equal the scaled reference,
|
||||
|
||||
$$
|
||||
r_h = \frac{J_h - G\,B_\mathrm{DW}\,p_h\,\mathrm{pol}_h\,I^\mathrm{ref}_h}{\sigma_{J,h}},
|
||||
\qquad
|
||||
L = \sum_h r_h^2 + \text{(scale prior)} .
|
||||
$$
|
||||
|
||||
Only the per-image scale $G$ and Debye–Waller $B$ are optimised; geometry and $R$ are
|
||||
fixed. Three consequences:
|
||||
|
||||
* **Integration and scaling become one objective.** $J$ *is* the integrated intensity and
|
||||
the residual *is* the scaling residual.
|
||||
* **The empty-pixel problem disappears by construction.** Empty pixels enter only through
|
||||
$J$ (with ~zero profile weight); they make no residual of their own and cannot dominate.
|
||||
* **Fisher weighting puts the reference at maximum leverage.** $\sigma_J$ uses the
|
||||
*model-expected* variance $v_p = B + \max(J,0)\,P_p$ (background plus expected signal from
|
||||
$I^\mathrm{ref}$), not the observed counts — so a strong *expected* reflection observed
|
||||
absent is penalised, and a noise spike with low $I^\mathrm{ref}$ gets no weight.
|
||||
|
||||
The scale is regularised towards 1 with a data-scaled weight $w_G=\sqrt{N_\mathrm{refl}/
|
||||
\sigma_G}$ (mirrors `ScaleOnTheFly`) so weakly-measured images cannot drift and scramble
|
||||
the merge.
|
||||
|
||||
**Geometry is not refined here.** PixelRefine's earlier per-image geometry refinement
|
||||
(regularised orientation LSQ, signal-weighting, a global orientation/cell sweep) was
|
||||
removed: on true stills the predictions are already good (radial centroid error ≈ 0,
|
||||
tangential ≈ a 0.4 px sampling floor that is *not* a recoverable misprediction), and per-image
|
||||
geometry refinement only overfit the sparse signal. Geometry is the job of `XtalOptimizer`.
|
||||
|
||||
---
|
||||
|
||||
## 5. Background estimator — mean, not median (both integrators)
|
||||
|
||||
**Symptom.** Pushed past the true resolution limit, the no-signal shells reported
|
||||
$\langle I/\sigma\rangle\approx4\text{–}6$ at $\mathrm{CC}_{1/2}\approx0$ — confident "data"
|
||||
where there is none. Present in *both* PixelRefine and the classical `BraggIntegrate2D`.
|
||||
|
||||
**Cause.** Both used the **median** of the background ring. For a right-skewed (Poisson)
|
||||
background $\operatorname{median}(B) < \mathbb{E}[B]$, so subtraction under-subtracts by a
|
||||
tiny but *coherent* per-pixel offset that grows over an $n_\mathrm{pix}$ peak and a
|
||||
multiplicity-$m$ merge into a fake $\langle I/\sigma\rangle\propto\sqrt{m}$, worst where the
|
||||
real signal is weakest.
|
||||
|
||||
**Fix.** Use the **mean** of the ring (spot cores and saturation sentinels already
|
||||
excluded). $\langle I/\sigma\rangle$ then collapses to ~0 wherever $\mathrm{CC}\approx0$ and
|
||||
the honest resolution limit becomes visible. This was the single largest contributor to
|
||||
untrustworthy σ — a one-line change in each integrator.
|
||||
|
||||
---
|
||||
|
||||
## 6. Error model (global $a,b$; XDS form)
|
||||
|
||||
Counting statistics under-estimate the variance of strong reflections, which carry
|
||||
systematic errors (scaling, partiality, detector) proportional to intensity, not to
|
||||
$\sqrt{I}$. After merging, this leaves CC$_{1/2}$ and $\langle I/\sigma\rangle$
|
||||
inconsistent. The standard correction (XDS / DIALS / AIMLESS) inflates the variance with a
|
||||
**global** two-parameter model:
|
||||
systematic errors proportional to $I$, not $\sqrt{I}$. The standard correction inflates the
|
||||
variance with a **global** two-parameter model, applied at the merge level so both
|
||||
integrators benefit:
|
||||
|
||||
$$
|
||||
\sigma'^{\,2} \;=\; a\,\sigma^{2} + \bigl(b\,\langle I\rangle\bigr)^{2},
|
||||
\sigma'^{\,2} = a\,\sigma^{2} + (b\,\langle I\rangle)^{2},
|
||||
\qquad
|
||||
\mathrm{ISa} \;=\; \frac{1}{b}\;=\;\lim_{I\to\infty}\frac{I}{\sigma'} .
|
||||
\mathrm{ISa} = \frac{1}{b} = \lim_{I\to\infty}\frac{I}{\sigma'} .
|
||||
$$
|
||||
|
||||
Two points matter for an unbiased fit:
|
||||
|
||||
* The $I^2$ term uses the reflection **mean** $\langle I\rangle$ (constant over its
|
||||
observations), **not** the per-observation $I_i$. Using $I_i$ gives a down-fluctuated
|
||||
point a small $\sigma'$ and hence a large $1/\sigma'^2$ weight, biasing the merged mean —
|
||||
which collapses CC. (This was the decisive bug in the first attempt.)
|
||||
* $a$ and $b$ are fit from the spread of symmetry equivalents: for an observation in a
|
||||
group of $n$, $\mathbb{E}[(I_i-\langle I\rangle)^2] = \sigma_i^2(1-h_i)$ with leverage
|
||||
$h_i = w_i/\sum w$. Binning by intensity and regressing the bin medians of
|
||||
$(I_i-\langle I\rangle)^2/(1-h_i)$ on $(\sigma^2, \langle I\rangle^2)$ — *weighted by
|
||||
$1/\mathrm{dev}^4$ so the fit is relative* — gives $(a, b^2)$; the relative weight stops
|
||||
the strong bins (which fix $b$) from swamping the weak bins (which fix $a$).
|
||||
|
||||
It is applied at the merge level (`MergeOnTheFly`), so **both** integrators benefit, and
|
||||
`jfjoch_process` prints the model and ISa. Earlier per-resolution-shell variants were
|
||||
dropped: the standard tools use a single global $a, b$, and the per-shell version was
|
||||
partly masking the background bias of Section 6.
|
||||
The $I^2$ term uses the reflection **mean** $\langle I\rangle$ (not the per-observation
|
||||
$I_i$, which would bias the merged mean and collapse CC); $a,b$ are fit from the spread of
|
||||
symmetry equivalents with a relative ($1/\mathrm{dev}^4$) weight so the strong bins (which
|
||||
fix $b$) do not swamp the weak bins (which fix $a$). `jfjoch_process` prints the model and
|
||||
ISa.
|
||||
|
||||
---
|
||||
|
||||
## Results (lysozyme jet, 1.8 Å, identical input)
|
||||
## Results (lysozyme rotation crystal `fixed_master.h5`, treated as stills, 1.7 Å)
|
||||
|
||||
| Configuration | N_obs | Compl. | CC$_{1/2}$ | CC$_\mathrm{ref}$ |
|
||||
| Configuration | N_obs | $\langle I/\sigma\rangle$ | CC$_{1/2}$ | CC$_\mathrm{ref}$ |
|
||||
|---|---:|---:|---:|---:|
|
||||
| Classical integrator (box-sum) | 421 k | 100 % | 81.4 % | 60.9 % |
|
||||
| PixelRefine — original | 61 k | 95 % | 0.0 % | −0.4 % |
|
||||
| PixelRefine — consolidated (this work) | 951 k | 100 % | 80.0 % | **62.6 %** |
|
||||
| Baseline per-pixel loss | 799 k | 7.2 | erratic, →0 at 1.7 Å | erratic |
|
||||
| **Factored Terms 1+2 (this model)** | 1.22 M | 10.7 | **84–92 % flat** | **77–92 % flat** |
|
||||
|
||||
The consolidated integrator beats the classical one on the accuracy metric
|
||||
(CC$_\mathrm{ref}$) with > 2× the multiplicity, and turns a previously unusable result
|
||||
(CC$_{1/2}=0$) into a competitive one.
|
||||
The factored objective turns the erratic, high-res-collapsing per-pixel result into flat
|
||||
~90 % CC$_{1/2}$/CC$_\mathrm{ref}$ to 1.7 Å — matching the proper rotation integration path
|
||||
from the *stills* path. (See `FINDINGS-2026-06.md`.)
|
||||
|
||||
---
|
||||
|
||||
## Default recipe
|
||||
|
||||
Sections 1–5 are `PixelRefine`-specific; Sections 6–7 act at the integration/merge level
|
||||
and apply to the classical route too.
|
||||
§§1–4 are PixelRefine-specific; §§5–6 act at the integration/merge level and apply to the
|
||||
classical route too.
|
||||
|
||||
| Field / behaviour | Default | Section |
|
||||
|---|---|---|
|
||||
| fit/extraction variance | local background $B$ | 1 |
|
||||
| `ewald_dist_cutoff` | $2\times10^{-3}\,\text{Å}^{-1}$ | 2 |
|
||||
| `scale_reg_sigma` | 2.0 | 3 |
|
||||
| `orient_reg_sigma_deg` | 1.0 | 3 |
|
||||
| `refine_R` | `false` | 3 |
|
||||
| `fit_signal_sigma_pix` | 1.5 | 4 |
|
||||
| `sweep_orientation` | `false` (available) | 5 |
|
||||
| local background estimator | **mean** of the ring | 6 |
|
||||
| merge error model | global $a,b$ (ISa printed) | 7 |
|
||||
| tangential width $R_1$ | measured per resolution shell | 3 |
|
||||
| objective | per-reflection intensity residual, Fisher-weighted | 4 |
|
||||
| refined parameters | per-image $G$ (and $B$); geometry fixed | 4 |
|
||||
| `scale_reg_sigma` | 2.0 | 4 |
|
||||
| local background estimator | **mean** of the ring | 5 |
|
||||
| merge error model | global $a,b$ (ISa printed) | 6 |
|
||||
|
||||
`orient_reg_sigma_deg` (accuracy vs. precision) and `ewald_dist_cutoff` (multiplicity vs.
|
||||
cost) are the two knobs worth tuning per dataset.
|
||||
`ewald_dist_cutoff` (multiplicity vs. cost) and `bandwidth` (Si vs. DMM) are the two knobs
|
||||
worth setting per dataset.
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@@ -3,247 +3,100 @@
|
||||
|
||||
#pragma once
|
||||
|
||||
#include <optional>
|
||||
|
||||
#include "../bragg_prediction/BraggPrediction.h"
|
||||
#include "../common/DiffractionExperiment.h"
|
||||
#include "../scale_merge/HKLKey.h"
|
||||
|
||||
// =============================================================================
|
||||
// PixelRefine — one optimization to rule geometry, integration and scaling
|
||||
// PixelRefine — reference-driven profile-fit integration + scaling for stills
|
||||
// =============================================================================
|
||||
//
|
||||
// Intent
|
||||
// ------
|
||||
// Classical crystallographic data processing is a one-way pipeline:
|
||||
// PixelRefine is the still-image integrator: given a reference set of merged
|
||||
// intensities I_ref (the current best hypothesis for each reflection's full
|
||||
// intensity), it integrates one image and returns already-scaled intensities. It
|
||||
// is an *intensity-wise* operation - the detector geometry is taken as fixed (it
|
||||
// was refined upstream by XtalOptimizer in IndexAndRefine::RefineGeometryIfNeeded);
|
||||
// PixelRefine does not touch orientation, cell or detector parameters.
|
||||
//
|
||||
// spot finding -> indexing -> geometry refinement -> integration -> scaling -> merging
|
||||
// The objective is the factored per-reflection likelihood of FACTORED_MODEL.md,
|
||||
// Terms 1 + 2:
|
||||
//
|
||||
// Each stage consumes the previous stage's output and never talks back. The
|
||||
// integrator trusts the refined geometry; the scaler trusts the integrated
|
||||
// intensities; nothing downstream is ever allowed to correct an upstream
|
||||
// parameter. Post-refinement and profile fitting were the field's partial
|
||||
// answers to this: post-refinement lets merged intensities nudge per-image
|
||||
// orientation/cell/mosaicity, and profile fitting lets a learned spot shape
|
||||
// improve weak-reflection intensities. But both are narrow back-channels bolted
|
||||
// onto a feed-forward pipeline — there is no end-to-end gradient that flows from
|
||||
// the raw detector pixels all the way back to every parameter at once.
|
||||
// Term 2 (shape) — for each resolution shell, the tangential profile width R1 is
|
||||
// *measured* from the intensity-weighted second moment of the strong spots:
|
||||
// R1 = sqrt(2*<eps_t^2>). A second moment is normalised by the total intensity,
|
||||
// so it is decoupled from the per-image scale - which is why measuring R1 is
|
||||
// stable where *fitting* it (degenerate with G) is not.
|
||||
//
|
||||
// PixelRefine is an experiment in doing the whole thing as a *single* least
|
||||
// squares problem. We write down, for every pixel in a reflection's shoebox, the
|
||||
// expected counts as an explicit forward model
|
||||
// Term 1 (intensity / scaling) — one residual per reflection: the profile-fit
|
||||
// amplitude J (using the Term-2 R1) should equal the scaled reference
|
||||
// J_model = G * exp(-B/4d^2) * partiality * pol * I_ref,
|
||||
// weighted by the model-expected (Fisher) sigma_J. Only the per-image scale G
|
||||
// and Debye-Waller B are optimised. Integration and scaling become one objective;
|
||||
// the many empty shoebox pixels enter only through J (with ~zero profile weight)
|
||||
// instead of dominating a per-pixel loss.
|
||||
//
|
||||
// I_pred(pixel) = G * I_true * B_term * P_radial * P_tangential * pol + I_bkg
|
||||
// I_ref is NOT refined here - it is a fixed hypothesis for the pass. The intended
|
||||
// outer loop is EM-like: run PixelRefine on every image against the current I_ref,
|
||||
// re-merge to a new I_ref, repeat.
|
||||
//
|
||||
// in raw detector counts (pol = per-reflection polarization correction, I_bkg =
|
||||
// local per-shoebox background read from the image).
|
||||
// Forward model per pixel (raw detector counts, no per-pixel solid-angle/Lorentz
|
||||
// weighting - same units as the classical integrator):
|
||||
// signal = G * I_ref * B_term * P_radial * P_tangential * pol , + I_bkg
|
||||
// B_term = exp(-B |q|^2 / 4) (Debye-Waller)
|
||||
// P_radial = exp(-eps_r^2 / R0_eff^2) (still partiality, <= 1)
|
||||
// P_tangential = exp(-eps_t^2 / R1^2) / (pi R1^2) (area-normalized profile)
|
||||
// where eps_r / eps_t are the radial / tangential deviations of the pixel from the
|
||||
// predicted node, and pol is the per-reflection polarization correction.
|
||||
//
|
||||
// and let Ceres autodiff back-propagate the per-pixel residuals into ALL of:
|
||||
// * detector geometry (beam centre, distance, tilt)
|
||||
// * crystal orientation + unit cell
|
||||
// * overall scale G and Debye-Waller B
|
||||
// * the reciprocal-space spot widths R = (radial, tangential)
|
||||
// simultaneously. Geometry refinement, profile-fitted integration and scaling
|
||||
// then stop being separate stages: they are different parameters of one model,
|
||||
// coupled through the same pixels, with full backpropagation between them. Once
|
||||
// the model is differentiable end-to-end, things that used to need bespoke code
|
||||
// — mosaicity refinement, profile fitting, partiality — fall out "for free" as
|
||||
// extra parameters of the same forward model.
|
||||
//
|
||||
// How I_true enters
|
||||
// -----------------
|
||||
// I_true is NOT refined here. It is a *fixed hypothesis* for the duration of a
|
||||
// pass: the current best merged estimate of each reflection's full intensity.
|
||||
// The intended outer loop is iterative, like EM / self-consistent field:
|
||||
//
|
||||
// repeat over the whole dataset:
|
||||
// run PixelRefine on every image with the current I_true reference
|
||||
// re-merge the resulting intensities -> new I_true
|
||||
// until the reference stops changing
|
||||
//
|
||||
// So a single PixelRefine call answers "given this intensity hypothesis, what
|
||||
// geometry/scale/profile best explains these pixels, and what intensities do I
|
||||
// read back out?", and the dataset-level loop refines the hypothesis itself.
|
||||
//
|
||||
// Inner predict<->refine loop
|
||||
// ---------------------------
|
||||
// Within one image we also iterate (max_iterations): Bragg prediction places the
|
||||
// shoeboxes, we refine, the refined geometry/cell feed the next prediction, etc.
|
||||
// Initially the model geometry equals the experiment's DiffractionGeometry, but
|
||||
// as refinement proceeds it diverges, so later predictions must use the *refined*
|
||||
// data.geom rather than the static experiment geometry.
|
||||
//
|
||||
// On shoeboxes, tails, and gatekeeping
|
||||
// ------------------------------------
|
||||
// We deliberately do NOT chase the full spot. A shoebox only needs to cover
|
||||
// enough of a reflection for the fit to be meaningful; clipped tails are fine,
|
||||
// because the partiality term already downweights whatever falls outside the
|
||||
// well-modelled core. The premise - especially for serial crystallography - is
|
||||
// that the problem is rarely *missing* information; it is *failing to gate out*
|
||||
// information that is not meaningful. As long as the model knows a piece is
|
||||
// missing (low partiality), it is safe to leave it missing. That flips the usual
|
||||
// trade-off: rather than shrinking boxes to avoid contamination, we can grow them
|
||||
// and let partiality decide, per pixel and per reflection, what actually carries
|
||||
// signal. (For downstream integration this pixel-level gating is the point - keep
|
||||
// only meaningful pixels, instead of a fixed geometric mask.) The bandwidth term
|
||||
// below is part of the same idea: it tells the model where the radial tails are
|
||||
// *expected* to be, so it can weight rather than blindly include them.
|
||||
//
|
||||
// X-ray bandwidth (optional)
|
||||
// --------------------------
|
||||
// A finite bandwidth thickens the Ewald shell radially and smears spots along the
|
||||
// radial direction, growing like 1/d^2 (the pink-beam/DMM signature). It enters
|
||||
// as a fixed, resolution-dependent addition to the radial width R0 (see
|
||||
// PixelRefineData::bandwidth and PixelRefine.cpp). It is OFF by default
|
||||
// (bandwidth = 0, monochromatic); set it for DMM-type data, leave it for Si.
|
||||
//
|
||||
// Status: experimental prototype. The forward model (esp. the still-image
|
||||
// partiality normalization) is deliberately simple and expected to evolve: it
|
||||
// works in raw detector counts with a local per-shoebox background and a
|
||||
// per-reflection polarization correction (no per-pixel solid-angle/Lorentz
|
||||
// weighting), matching the classical integrator. See PixelRefine.cpp for the
|
||||
// physics conventions and known caveats.
|
||||
// X-ray bandwidth (optional): a finite bandwidth thickens the Ewald shell radially,
|
||||
// adding a fixed, resolution-dependent term to R0 that grows like 1/d^2 (the
|
||||
// pink-beam/DMM signature): R0_eff^2 = R0^2 + (b*lambda)^2/(2 d^4). b = 0 (the
|
||||
// default) is a monochromatic no-op; set it for DMM-type data, leave it for Si.
|
||||
// =============================================================================
|
||||
|
||||
struct PixelRefineData {
|
||||
// --- model state (input as initial guess, output as refined result) ---
|
||||
DiffractionGeometry geom;
|
||||
CrystalLattice latt;
|
||||
DiffractionGeometry geom; // fixed (refined upstream by XtalOptimizer)
|
||||
CrystalLattice latt; // fixed
|
||||
gemmi::CrystalSystem crystal_system = gemmi::CrystalSystem::Triclinic;
|
||||
char centering = 'P';
|
||||
|
||||
double B_factor = 0.0; // Debye-Waller B (A^2)
|
||||
double scale_factor = 1.0; // overall scale G
|
||||
double R[2] = {0.005, 0.005}; // R[0] = radial (partiality) width, R[1] = tangential (profile) width (A^-1)
|
||||
double B_factor = 0.0; // Debye-Waller B (A^2), refined
|
||||
double scale_factor = 1.0; // overall per-image scale G, refined
|
||||
double R[2] = {0.005, 0.005}; // R[0] = radial (partiality) width; R[1] = fallback
|
||||
// tangential profile width before Term 2 measures it (A^-1)
|
||||
bool refine_B = true; // refine the per-image B-factor along with G
|
||||
|
||||
// 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).
|
||||
// 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).
|
||||
double bandwidth = 0.0;
|
||||
|
||||
// --- what to refine ---
|
||||
bool refine_orientation = true; // crystal orientation (p0)
|
||||
bool refine_unit_cell = false; // cell lengths + angles
|
||||
bool refine_beam_center = false;
|
||||
bool refine_distance = false;
|
||||
bool refine_detector_angles = false;
|
||||
bool refine_scale = true;
|
||||
bool refine_B = false;
|
||||
bool refine_R = false; // per-image R refinement is unstable on sparse
|
||||
// stills; R is held at its nominal value
|
||||
|
||||
// Orientation refinement is anchored to the pre-refinement (spot-centroid)
|
||||
// orientation with weight sqrt(n_refl)/sigma, so the pixel fit can only nudge
|
||||
// the orientation by ~this many degrees before the prior pushes back. This is
|
||||
// what turns the (otherwise overfitting) 3-DOF orientation refinement into a
|
||||
// small, signal-supported sub-spot correction. Larger => freer; very large
|
||||
// approaches the unregularized (collapsing) fit. ~1 deg gives the best CCref here;
|
||||
// beyond ~2 deg the per-image fit overfits and the merge collapses.
|
||||
double orient_reg_sigma_deg = 1.0;
|
||||
|
||||
// Signal-weighting of the *fit* residuals: each pixel's weight is multiplied by
|
||||
// a detector-space Gaussian exp(-r^2/2 sigma^2) centred on the predicted spot, so
|
||||
// the many empty shoebox-corner pixels stop diluting (and destabilising) the fit
|
||||
// and the signal-bearing core drives the refined scale/R. <= 0 disables (uniform).
|
||||
double fit_signal_sigma_pix = 1.5;
|
||||
|
||||
// Per-image scale G is regularized towards 1 with weight sqrt(n_refl/scale_reg_sigma)
|
||||
// (mirrors ScaleOnTheFly). Without this the unconstrained G wanders on weakly
|
||||
// measured images and 1/G scrambles the cross-image merge. <= 0 disables.
|
||||
double scale_reg_sigma = 2.0;
|
||||
|
||||
// Radial Ewald-sphere acceptance band for prediction (A^-1): a reflection is
|
||||
// given a shoebox when ||S|-1/lambda| <= this. The narrow default predicts only
|
||||
// reflections already on the Ewald sphere; widening it (towards the integrator's
|
||||
// ~2-3x profile radius) lets in the slightly-misaligned high-resolution
|
||||
// reflections - more multiplicity, and something for orientation refinement to
|
||||
// actually centre. Safe to widen only with the per-image fit kept well-behaved
|
||||
// (de-biased variance + signal-weighting + regularization, all default here).
|
||||
// Radial Ewald-sphere acceptance band for prediction (A^-1): a reflection is given
|
||||
// a shoebox when ||S|-1/lambda| <= this. Widened from the on-sphere default towards
|
||||
// the integrator's profile radius so slightly-misaligned high-resolution reflections
|
||||
// are still integrated (multiplicity), while the partiality downweights their tails.
|
||||
double ewald_dist_cutoff = 0.0020;
|
||||
|
||||
// Pre-LSQ global orientation+cell sweep (maximises CC vs reference over the
|
||||
// strongest reflections). Bounds are derived from the detector geometry: the
|
||||
// step moves the highest-resolution spot by 1 px, the range moves the lowest-
|
||||
// resolution spot by ~2 px (rotation) / ~1 px (cell scale), so the low-res
|
||||
// XtalOptimizer solution is preserved while high-res spots are recentred.
|
||||
bool sweep_orientation = false;
|
||||
// Sweep half-range as the *orientation uncertainty* (degrees) and cell-scale
|
||||
// uncertainty (fraction). These set how far the highest-resolution spot may move
|
||||
// (a few px); low-res spots barely move and stay anchored. Keep small - a large
|
||||
// range lets the per-image CC overfit and degrades the merge.
|
||||
double sweep_max_deg = 0.15;
|
||||
double sweep_max_cell_frac = 0.003;
|
||||
|
||||
double max_time_s = 5.0;
|
||||
int shoebox_radius = 3; // half-size of the per-reflection signal box (peak region that enters the fit)
|
||||
// Half-size of the local-background sampling box. Background is estimated from
|
||||
// the ring shoebox_radius < |dx|,|dy| <= bkg_outer_radius around each spot
|
||||
// (excluding pixels belonging to any predicted spot core), mirroring the local
|
||||
// shoebox background of BraggIntegrate2D. Must be > shoebox_radius.
|
||||
int shoebox_radius = 3; // half-size of the per-reflection signal box
|
||||
// Half-size of the local-background sampling box. Background is the MEAN of the ring
|
||||
// shoebox_radius < |dx|,|dy| <= bkg_outer_radius (excluding spot cores), like
|
||||
// BraggIntegrate2D. Must be > shoebox_radius.
|
||||
int bkg_outer_radius = 6;
|
||||
int max_iterations = 3; // inner predict<->refine cycles (re-predict with refined geom/latt)
|
||||
|
||||
// Diagnostic: compute the parameter correlation matrix (Pearson) from the final
|
||||
// per-image solve, to expose degeneracies (e.g. G<->B, G<->R1) that let the fit
|
||||
// lower its chi-square along directions that do not generalise across images.
|
||||
// Requires G, B and R all refined; results in corr_* below.
|
||||
bool compute_covariance = false;
|
||||
bool fix_R0 = false; // diagnostic: with refine_R, hold R0 constant and refine R1 only
|
||||
|
||||
// Factored-likelihood Term 1 (FACTORED_MODEL.md): replace the per-pixel fit with one
|
||||
// per-reflection *intensity* residual J vs G*B_term*partiality*pol*I_ref, Fisher-
|
||||
// weighted. Geometry & R fixed; only G (and B if refine_B) are fit.
|
||||
bool intensity_residual = false;
|
||||
|
||||
// Factored-likelihood Term 2: set the tangential profile width R1 from the *measured*
|
||||
// per-resolution second moment of the strong spots (a shape statistic, decoupled from
|
||||
// the scale) and feed it to the profile template used by Term 1 and the extraction.
|
||||
bool shape_R1 = false;
|
||||
double shape_R1_lores = NAN, shape_R1_hires = NAN; // measured R1 in the lowest/highest-res bin (diag)
|
||||
|
||||
// Adaptive integration mask: set R1 (tangential profile width) from the *measured*
|
||||
// tangential second moment of the strong spots, instead of fitting it (which is
|
||||
// degenerate with the per-image scale). A shape statistic, independent of scale.
|
||||
bool adaptive_R1 = false;
|
||||
|
||||
// Diagnostic: residual centering error after refinement. For the strong spots,
|
||||
// the offset between observed intensity centroid and predicted position - bias is
|
||||
// the systematic (mean) part, rms the total. If rms is comparable to the spot size,
|
||||
// a tight profile mask lands off the spot and box-summing wins.
|
||||
bool measure_centroid = false;
|
||||
double centroid_bias_px = NAN;
|
||||
double centroid_rms_px = NAN;
|
||||
// Offset split by spot significance (lo/hi about the median) and measured two ways:
|
||||
// the intensity *centroid* (mean) and the sub-pixel *peak* (mode, parabolic fit).
|
||||
// - centroid offset shrinking lo->hi => noise floor (offset ~ 1/sqrt(counts));
|
||||
// flat with significance => a systematic position error to model.
|
||||
// - peak offset << centroid offset => the spot is asymmetric (non-Gaussian /
|
||||
// parallax): the prediction sits on the peak, only the centroid is pulled, so
|
||||
// recentring on the centroid would be wrong - the shape is what to model.
|
||||
// - radial centroid offset ~ 0 => no distance/parallax error (parallax is radial).
|
||||
double centroid_lo_signif = NAN, centroid_hi_signif = NAN; // mean significance per bin
|
||||
double centroid_lo_tang_c = NAN, centroid_hi_tang_c = NAN; // mean |tangential centroid offset|
|
||||
double centroid_lo_tang_p = NAN, centroid_hi_tang_p = NAN; // mean |tangential peak offset|
|
||||
double centroid_lo_rad_c = NAN, centroid_hi_rad_c = NAN; // mean signed radial centroid offset
|
||||
|
||||
// Test: recentre the extraction profile on each spot's observed centroid (instead of
|
||||
// the predicted position) so a tight mask lands on the real spot - but only for spots
|
||||
// whose in-shoebox significance exceeds recenter_min_signif (recentring on a noise
|
||||
// centroid would bias weak reflections positive).
|
||||
bool recenter_profile = false;
|
||||
double recenter_min_signif = 5.0;
|
||||
|
||||
// --- output ---
|
||||
std::vector<Reflection> reflections; // profile-fitted integration result
|
||||
std::vector<Reflection> reflections; // profile-fitted, scaled integration result
|
||||
bool solved = false;
|
||||
double final_cost = NAN;
|
||||
size_t residual_count = 0;
|
||||
double cc = NAN; // per-image CC of scaled intensities vs reference
|
||||
int64_t cc_n = 0; // number of reflections in the CC
|
||||
|
||||
bool covariance_valid = false;
|
||||
double corr_GB = NAN, corr_GR0 = NAN, corr_GR1 = NAN;
|
||||
double corr_BR0 = NAN, corr_BR1 = NAN, corr_R0R1 = NAN;
|
||||
};
|
||||
|
||||
class PixelRefine {
|
||||
@@ -253,67 +106,24 @@ class PixelRefine {
|
||||
const HKLKeyGenerator hkl_key_generator;
|
||||
std::map<HKLKey, double> reference_data;
|
||||
|
||||
// Fills the Ceres parameter blocks (geometry + symmetry-aware lattice
|
||||
// parametrization) from the current model state. Shared by Run and
|
||||
// PredictImage so both walk identical geometry/lattice code.
|
||||
// Fills the fixed geometry + symmetry-aware lattice parametrization (beam,
|
||||
// distance, detector tilt, and the Rodrigues orientation / cell-length / angle
|
||||
// vectors) from the current model state, for the per-pixel geometry evaluation.
|
||||
void BuildParameterBlocks(const PixelRefineData &data,
|
||||
double beam[2], double &dist_mm,
|
||||
double detector_rot[2],
|
||||
double latt_vec0[3], double latt_vec1[3], double latt_vec2[3]) const;
|
||||
|
||||
// Global orientation + uniform cell-scale sweep run before the LSQ. Re-projects
|
||||
// the strongest reference reflections through candidate lattices and keeps the
|
||||
// one maximising CC vs the reference intensities (coordinate descent over the 3
|
||||
// Rodrigues axes + cell scale, within geometry-derived pixel bounds). Writes the
|
||||
// refined orientation/cell back into data.latt. See PixelRefineData::sweep_*.
|
||||
template<class T>
|
||||
void SweepOrientationCell(const T *image, BraggPrediction &prediction,
|
||||
PixelRefineData &data) const;
|
||||
public:
|
||||
PixelRefine(const DiffractionExperiment &experiment,
|
||||
const std::vector<MergedReflection> &reference);
|
||||
|
||||
// The BraggPrediction is supplied per call (it is mutated): this keeps a
|
||||
// single PixelRefine instance usable from several threads, each passing its
|
||||
// own prediction buffer. Only `data` is written; PixelRefine state is const.
|
||||
// The image is in raw detector counts (masked/saturated pixels carry the type
|
||||
// sentinel); background is estimated locally per shoebox from the image itself.
|
||||
// The BraggPrediction is supplied per call (it is mutated): this keeps a single
|
||||
// PixelRefine instance usable from several threads, each passing its own prediction
|
||||
// buffer. Only `data` is written; PixelRefine state is const. The image is in raw
|
||||
// detector counts (masked/saturated pixels carry the type sentinel); background is
|
||||
// estimated locally per shoebox from the image itself.
|
||||
template<class T>
|
||||
void Run(const T *image,
|
||||
BraggPrediction &prediction,
|
||||
PixelRefineData &data);
|
||||
|
||||
// Render the forward model as a full detector image (raw detector counts, so
|
||||
// it overlays directly on the original image). Uses the *same* per-pixel
|
||||
// model path (PixelResidual::Model) as the optimizer, evaluated in double
|
||||
// precision - slow but exact. For each reference reflection it adds the Bragg
|
||||
// signal over its shoebox; with include_background it also lays down the local
|
||||
// per-shoebox background read from the supplied image - the same background the
|
||||
// fit uses. Diagnostic tool, not on the hot path.
|
||||
template<class T>
|
||||
std::vector<float> PredictImage(const T *image,
|
||||
BraggPrediction &prediction,
|
||||
const PixelRefineData &data,
|
||||
bool include_background = true) const;
|
||||
|
||||
// Render the per-pixel chi-square (cost density) that the optimizer actually
|
||||
// minimizes: for every shoebox pixel that enters the fit it stores the squared
|
||||
// weighted residual ((I_pred - I_obs)/sigma)^2 in raw counts - identical to the
|
||||
// Ceres residual_i^2 - accumulating where shoeboxes overlap. Pixels that are not
|
||||
// part of any shoebox stay 0; masked/saturated pixels (skipped by the fit) also
|
||||
// stay 0. Summing the image gives ~2*final_cost. Diagnostic tool.
|
||||
template<class T>
|
||||
std::vector<float> ChiSquaredImage(const T *image,
|
||||
BraggPrediction &prediction,
|
||||
const PixelRefineData &data) const;
|
||||
|
||||
// Reference (merged) intensity used as the fixed hypothesis for a reflection,
|
||||
// or nullopt if this hkl is not in the reference. Lets callers show the fitted
|
||||
// estimate next to the reference it was scaled against.
|
||||
std::optional<double> ReferenceIntensity(const Reflection &r) const {
|
||||
const auto it = reference_data.find(hkl_key_generator(r));
|
||||
if (it == reference_data.end())
|
||||
return std::nullopt;
|
||||
return it->second;
|
||||
}
|
||||
};
|
||||
|
||||
@@ -86,13 +86,8 @@ ADD_EXECUTABLE(jfjoch_viewer jfjoch_viewer.cpp JFJochViewerWindow.cpp JFJochView
|
||||
${APP_RESOURCES}
|
||||
windows/JFJochViewerReciprocalSpaceWindow.cpp
|
||||
windows/JFJochViewerReciprocalSpaceWindow.h
|
||||
windows/JFJochPixelRefineWindow.cpp
|
||||
windows/JFJochPixelRefineWindow.h
|
||||
windows/JFJochPixelRefineTableWindow.cpp
|
||||
windows/JFJochPixelRefineTableWindow.h
|
||||
windows/JFJochMagnifierWindow.cpp
|
||||
windows/JFJochMagnifierWindow.h
|
||||
windows/PixelRefineParams.h
|
||||
)
|
||||
|
||||
TARGET_LINK_LIBRARIES(jfjoch_viewer Qt6::Core Qt6::Gui Qt6::Widgets Qt6::Charts Qt6::DBus Qt6::Concurrent
|
||||
|
||||
@@ -9,8 +9,6 @@
|
||||
|
||||
#include "JFJochImageReadingWorker.h"
|
||||
#include "../reader/JFJochReaderImage.h" // JFJochReaderImage + GAP/ERROR/SATURATED sentinels
|
||||
#include "../image_analysis/LoadFCalcFromMtz.h"
|
||||
#include "../image_analysis/bragg_prediction/BraggPredictionFactory.h"
|
||||
#include "../image_analysis/geom_refinement/AssignSpotsToRings.h"
|
||||
#include "../image_analysis/spot_finding/StrongPixelSet.h"
|
||||
#include "../image_analysis/spot_finding/SpotUtils.h"
|
||||
@@ -69,8 +67,6 @@ JFJochImageReadingWorker::JFJochImageReadingWorker(const SpotFindingSettings &se
|
||||
: QObject(parent),
|
||||
indexing_settings(experiment.GetIndexingSettings()),
|
||||
azint_settings(experiment.GetAzimuthalIntegrationSettings()) {
|
||||
qRegisterMetaType<PixelRefineParams>("PixelRefineParams");
|
||||
qRegisterMetaType<PixelRefineReport>("PixelRefineReport");
|
||||
qRegisterMetaType<QVector<QRect>>("QVector<QRect>");
|
||||
spot_finding_settings = settings;;
|
||||
|
||||
@@ -307,13 +303,7 @@ void JFJochImageReadingWorker::UpdateAzint_i(const JFJochReaderDataset *dataset)
|
||||
image_analysis = std::make_unique<MXAnalysisWithoutFPGA>(curr_experiment, *azint_mapping, dataset->pixel_mask,
|
||||
*index_and_refine.get());
|
||||
|
||||
// PixelRefine state is tied to the experiment/mapping; rebuild lazily.
|
||||
pixel_refine_.reset();
|
||||
pixel_pred_.reset();
|
||||
last_profile_.reset();
|
||||
// Keep scale-on-the-fly alive across dataset reloads.
|
||||
if (!pixel_reference_.empty())
|
||||
index_and_refine->ReferenceIntensities(pixel_reference_);
|
||||
}
|
||||
}
|
||||
|
||||
@@ -722,287 +712,3 @@ void JFJochImageReadingWorker::LoadSpots(int64_t start_image, int64_t end_image,
|
||||
emit spotsLoaded(result);
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
// Experimental PixelRefine
|
||||
// ---------------------------------------------------------------------------
|
||||
void JFJochImageReadingWorker::EnsurePixelRefine_i() {
|
||||
if (!pixel_refine_ && !pixel_reference_.empty())
|
||||
pixel_refine_ = std::make_unique<PixelRefine>(curr_experiment, pixel_reference_);
|
||||
if (!pixel_pred_)
|
||||
pixel_pred_ = CreateBraggPrediction(curr_experiment.IsRotationIndexing());
|
||||
}
|
||||
|
||||
bool JFJochImageReadingWorker::BuildPixelSeed_i(PixelRefineData &d, const PixelRefineParams &p, QString &reason) const {
|
||||
if (!current_image_ptr || !index_and_refine || !current_image.has_value()) {
|
||||
reason = "No image loaded";
|
||||
return false;
|
||||
}
|
||||
|
||||
const auto &outcomes = index_and_refine->GetIntegrationOutcome();
|
||||
const int64_t n = current_image.value();
|
||||
if (n < 0 || n >= static_cast<int64_t>(outcomes.size())) {
|
||||
reason = "No analysis result for current image";
|
||||
return false;
|
||||
}
|
||||
|
||||
const auto &io = outcomes[n];
|
||||
if (io.reflections.empty()) {
|
||||
// The "indexed" badge in the viewer comes from indexing (lattice_type);
|
||||
// PixelRefine instead seeds from the *integration* outcome, which is only
|
||||
// populated when Quick integration is enabled and succeeds for this image.
|
||||
// Distinguish the two so the user knows what to turn on.
|
||||
if (current_image_ptr->ImageData().lattice_type.has_value())
|
||||
reason = "Image is indexed but not integrated - enable Quick integration";
|
||||
else
|
||||
reason = "Current image is not indexed";
|
||||
return false;
|
||||
}
|
||||
|
||||
d.geom = io.geom;
|
||||
d.latt = io.latt;
|
||||
|
||||
const auto < = current_image_ptr->ImageData().lattice_type;
|
||||
if (lt) {
|
||||
d.crystal_system = lt->crystal_system;
|
||||
d.centering = lt->centering;
|
||||
}
|
||||
if (d.crystal_system == gemmi::CrystalSystem::Trigonal)
|
||||
d.crystal_system = gemmi::CrystalSystem::Hexagonal;
|
||||
|
||||
d.R[0] = p.R0;
|
||||
d.R[1] = p.R1;
|
||||
d.bandwidth = p.bandwidth_fwhm / 2.3548; // FWHM -> sigma
|
||||
d.scale_factor = p.scale_factor;
|
||||
d.B_factor = p.B_factor;
|
||||
if (std::isfinite(p.beam_x) && std::isfinite(p.beam_y))
|
||||
d.geom.BeamX_pxl(static_cast<float>(p.beam_x)).BeamY_pxl(static_cast<float>(p.beam_y));
|
||||
|
||||
d.refine_orientation = p.refine_orientation;
|
||||
d.refine_unit_cell = p.refine_unit_cell;
|
||||
d.refine_beam_center = p.refine_beam_center;
|
||||
d.refine_scale = p.refine_scale;
|
||||
d.refine_B = p.refine_B;
|
||||
d.refine_R = p.refine_R;
|
||||
d.max_iterations = p.max_iterations;
|
||||
return true;
|
||||
}
|
||||
|
||||
std::shared_ptr<SimpleImage> JFJochImageReadingWorker::WrapFloatImage_i(const std::vector<float> &img) const {
|
||||
auto si = std::make_shared<SimpleImage>();
|
||||
// CompressedImage is a non-owning view over its data pointer. predictedImageReady
|
||||
// is a queued cross-thread connection, so the source float vector must outlive the
|
||||
// emit: copy it into SimpleImage::buffer (which the shared_ptr keeps alive) instead
|
||||
// of aliasing the caller's temporary, otherwise loadImageInternal() reads freed
|
||||
// memory in the GUI thread (SIGSEGV).
|
||||
si->buffer.resize(img.size() * sizeof(float));
|
||||
std::memcpy(si->buffer.data(), img.data(), si->buffer.size());
|
||||
si->image = CompressedImage(si->buffer, curr_experiment.GetXPixelsNum(), curr_experiment.GetYPixelsNum(),
|
||||
CompressedImageMode::Float32);
|
||||
return si;
|
||||
}
|
||||
|
||||
void JFJochImageReadingWorker::SquaredResidualWithImage_i(std::vector<float> &pred) const {
|
||||
// PredictImage() returns raw detector units (same as the measured counts), so
|
||||
// pred - measured is the per-pixel residual the model fails to explain. We plot
|
||||
// |pred - measured|^2: sign-free, so it needs no diverging colour scale and just
|
||||
// highlights where the model disagrees most. Masked / saturated pixels carry
|
||||
// sentinels rather than counts, so no comparison is possible -> NaN (gap).
|
||||
if (!current_image_ptr)
|
||||
return;
|
||||
const auto &img = current_image_ptr->Image();
|
||||
const size_t n = std::min(pred.size(), img.size());
|
||||
for (size_t i = 0; i < n; ++i) {
|
||||
const int32_t v = img[i];
|
||||
if (v == GAP_PXL_VALUE || v == ERROR_PXL_VALUE || v == SATURATED_PXL_VALUE) {
|
||||
pred[i] = NAN;
|
||||
} else {
|
||||
const float diff = pred[i] - static_cast<float>(v);
|
||||
pred[i] = diff * diff;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void JFJochImageReadingWorker::MaskMeasuredSentinels_i(std::vector<float> &img) const {
|
||||
// The chi^2 image is 0 outside shoeboxes; show masked/saturated pixels as a gap
|
||||
// (NaN) instead, so they read as "not comparable" rather than "zero cost".
|
||||
if (!current_image_ptr)
|
||||
return;
|
||||
const auto &measured = current_image_ptr->Image();
|
||||
const size_t n = std::min(img.size(), measured.size());
|
||||
for (size_t i = 0; i < n; ++i) {
|
||||
const int32_t v = measured[i];
|
||||
if (v == GAP_PXL_VALUE || v == ERROR_PXL_VALUE || v == SATURATED_PXL_VALUE)
|
||||
img[i] = NAN;
|
||||
}
|
||||
}
|
||||
|
||||
QVector<QRect> JFJochImageReadingWorker::BuildShoeboxes_i(const PixelRefineData &data) const {
|
||||
// One rectangle per fitted reflection: the shoebox the optimizer summed over,
|
||||
// centred on the predicted position with half-size data.shoebox_radius.
|
||||
QVector<QRect> boxes;
|
||||
boxes.reserve(static_cast<int>(data.reflections.size()));
|
||||
const int r = data.shoebox_radius;
|
||||
const int side = 2 * r + 1;
|
||||
for (const auto &refl : data.reflections) {
|
||||
if (!std::isfinite(refl.predicted_x) || !std::isfinite(refl.predicted_y))
|
||||
continue;
|
||||
const int cx = static_cast<int>(std::lround(refl.predicted_x));
|
||||
const int cy = static_cast<int>(std::lround(refl.predicted_y));
|
||||
boxes.push_back(QRect(cx - r, cy - r, side, side));
|
||||
}
|
||||
return boxes;
|
||||
}
|
||||
|
||||
PixelRefineReport JFJochImageReadingWorker::BuildReport_i(const PixelRefineData &data) const {
|
||||
PixelRefineReport report;
|
||||
report.pr_G = data.scale_factor;
|
||||
report.pr_B = data.B_factor;
|
||||
report.pr_cc = data.cc;
|
||||
report.pr_cc_n = data.cc_n;
|
||||
|
||||
// Standard ScaleOnTheFly pipeline result for the same image, as a baseline.
|
||||
if (current_image_ptr) {
|
||||
const auto &d = current_image_ptr->ImageData();
|
||||
if (d.image_scale_factor) report.pipe_G = d.image_scale_factor.value();
|
||||
if (d.image_scale_b_factor) report.pipe_B = d.image_scale_b_factor.value();
|
||||
if (d.image_scale_cc) report.pipe_cc = d.image_scale_cc.value();
|
||||
}
|
||||
|
||||
report.rows.reserve(data.reflections.size());
|
||||
for (const auto &r : data.reflections) {
|
||||
if (!r.observed)
|
||||
continue;
|
||||
PixelRefineReport::Row row;
|
||||
row.h = r.h; row.k = r.k; row.l = r.l;
|
||||
row.d = r.d;
|
||||
row.completeness = r.completeness;
|
||||
row.partiality = r.partiality;
|
||||
row.I = r.I;
|
||||
row.sigma = r.sigma;
|
||||
if (std::isfinite(r.image_scale_corr))
|
||||
row.I_true_est = static_cast<double>(r.I) * static_cast<double>(r.image_scale_corr);
|
||||
if (pixel_refine_)
|
||||
row.I_true_ref = pixel_refine_->ReferenceIntensity(r).value_or(NAN);
|
||||
report.rows.push_back(row);
|
||||
}
|
||||
return report;
|
||||
}
|
||||
|
||||
std::vector<float> JFJochImageReadingWorker::BuildDisplayImage_i(const PixelRefineData &data,
|
||||
int display_mode) const {
|
||||
const auto &img32 = current_image_ptr->Image();
|
||||
if (display_mode == PixelRefineParams::ChiSquared) {
|
||||
// The cost density the optimizer actually minimizes (weighted residual^2).
|
||||
auto chi2 = pixel_refine_->ChiSquaredImage<int32_t>(img32.data(), *pixel_pred_, data);
|
||||
MaskMeasuredSentinels_i(chi2);
|
||||
return chi2;
|
||||
}
|
||||
|
||||
auto pred = pixel_refine_->PredictImage<int32_t>(img32.data(), *pixel_pred_, data, true);
|
||||
if (display_mode == PixelRefineParams::SquaredDifference)
|
||||
SquaredResidualWithImage_i(pred);
|
||||
return pred;
|
||||
}
|
||||
|
||||
void JFJochImageReadingWorker::LoadReference(QString path) {
|
||||
QMutexLocker ul(&m);
|
||||
try {
|
||||
pixel_reference_ = LoadFCalcFromMtz(path.toStdString());
|
||||
if (index_and_refine)
|
||||
index_and_refine->ReferenceIntensities(pixel_reference_); // enables scale-on-the-fly too
|
||||
pixel_refine_.reset(); // rebuild with new reference
|
||||
emit pixelRefineStatus(QString("Loaded %1 reference reflections").arg(pixel_reference_.size()));
|
||||
} catch (const std::exception &e) {
|
||||
emit pixelRefineStatus(QString("Failed to load reference: %1").arg(e.what()));
|
||||
}
|
||||
}
|
||||
|
||||
void JFJochImageReadingWorker::PixelRefinePreview(PixelRefineParams params) {
|
||||
QMutexLocker ul(&m);
|
||||
if (!last_profile_) {
|
||||
emit pixelRefineStatus("Analyze an image first");
|
||||
return;
|
||||
}
|
||||
EnsurePixelRefine_i();
|
||||
if (!pixel_refine_) {
|
||||
emit pixelRefineStatus("Load reference data first");
|
||||
return;
|
||||
}
|
||||
|
||||
PixelRefineData d;
|
||||
QString seed_reason;
|
||||
if (!BuildPixelSeed_i(d, params, seed_reason)) {
|
||||
emit pixelRefineStatus(seed_reason);
|
||||
return;
|
||||
}
|
||||
|
||||
// Preview = evaluate only: do not move any parameter.
|
||||
d.refine_orientation = d.refine_unit_cell = d.refine_beam_center = false;
|
||||
d.refine_scale = d.refine_B = d.refine_R = false;
|
||||
d.max_iterations = 0;
|
||||
|
||||
try {
|
||||
const auto &img32 = current_image_ptr->Image();
|
||||
pixel_refine_->Run<int32_t>(img32.data(), *pixel_pred_, d);
|
||||
emit pixelRefineResidual(d.final_cost, d.cc, static_cast<int64_t>(d.reflections.size()));
|
||||
emit pixelRefineReport(BuildReport_i(d));
|
||||
|
||||
auto display = BuildDisplayImage_i(d, params.display_mode);
|
||||
emit predictedImageReady(WrapFloatImage_i(display));
|
||||
emit predictedShoeboxes(BuildShoeboxes_i(d));
|
||||
} catch (const std::exception &e) {
|
||||
emit pixelRefineStatus(QString("PixelRefine preview failed: %1").arg(e.what()));
|
||||
}
|
||||
}
|
||||
|
||||
void JFJochImageReadingWorker::PixelRefineRun(PixelRefineParams params) {
|
||||
QMutexLocker ul(&m);
|
||||
if (!last_profile_) {
|
||||
emit pixelRefineStatus("Analyze an image first");
|
||||
return;
|
||||
}
|
||||
EnsurePixelRefine_i();
|
||||
if (!pixel_refine_) {
|
||||
emit pixelRefineStatus("Load reference data first");
|
||||
return;
|
||||
}
|
||||
|
||||
PixelRefineData d;
|
||||
QString seed_reason;
|
||||
if (!BuildPixelSeed_i(d, params, seed_reason)) {
|
||||
emit pixelRefineStatus(seed_reason);
|
||||
return;
|
||||
}
|
||||
if (d.max_iterations <= 0)
|
||||
d.max_iterations = 3;
|
||||
|
||||
try {
|
||||
const auto &img32 = current_image_ptr->Image();
|
||||
pixel_refine_->Run<int32_t>(img32.data(), *pixel_pred_, d);
|
||||
|
||||
// Push refined values back so the sliders follow the optimizer.
|
||||
PixelRefineParams out = params;
|
||||
out.R0 = d.R[0];
|
||||
out.R1 = d.R[1];
|
||||
out.bandwidth_fwhm = d.bandwidth * 2.3548; // sigma -> FWHM
|
||||
out.scale_factor = d.scale_factor;
|
||||
out.B_factor = d.B_factor;
|
||||
out.beam_x = d.geom.GetBeamX_pxl();
|
||||
out.beam_y = d.geom.GetBeamY_pxl();
|
||||
emit pixelRefineParamsRefined(out);
|
||||
emit pixelRefineResidual(d.final_cost, d.cc, static_cast<int64_t>(d.reflections.size()));
|
||||
emit pixelRefineReport(BuildReport_i(d));
|
||||
|
||||
auto display = BuildDisplayImage_i(d, params.display_mode);
|
||||
emit predictedImageReady(WrapFloatImage_i(display));
|
||||
emit predictedShoeboxes(BuildShoeboxes_i(d));
|
||||
|
||||
// Show the refined predictions on the main image too.
|
||||
auto new_image = std::make_shared<JFJochReaderImage>(*current_image_ptr);
|
||||
new_image->ImageData().reflections = d.reflections;
|
||||
current_image_ptr = new_image;
|
||||
emit imageLoaded(current_image_ptr);
|
||||
} catch (const std::exception &e) {
|
||||
emit pixelRefineStatus(QString("PixelRefine failed: %1").arg(e.what()));
|
||||
}
|
||||
}
|
||||
|
||||
@@ -15,10 +15,8 @@
|
||||
#include "../common/Logger.h"
|
||||
#include "../reader/JFJochHttpReader.h"
|
||||
#include "../image_analysis/MXAnalysisWithoutFPGA.h"
|
||||
#include "../image_analysis/pixel_refinement/PixelRefine.h"
|
||||
#include "../image_analysis/bragg_prediction/BraggPrediction.h"
|
||||
#include "SimpleImage.h"
|
||||
#include "windows/PixelRefineParams.h"
|
||||
#include "../common/MovingAverage.h"
|
||||
|
||||
Q_DECLARE_METATYPE(std::shared_ptr<const JFJochReaderDataset>)
|
||||
@@ -60,30 +58,8 @@ private:
|
||||
std::unique_ptr<MXAnalysisWithoutFPGA> image_analysis;
|
||||
std::unique_ptr<IndexAndRefine> index_and_refine;
|
||||
|
||||
// Experimental PixelRefine support. last_profile_ keeps the azimuthal profile
|
||||
// of the most recently analyzed image (PixelRefine needs it); the engine and a
|
||||
// dedicated prediction buffer are built lazily once a reference is loaded.
|
||||
// Azimuthal profile buffer of the most recently analyzed image (filled by Analyze).
|
||||
std::unique_ptr<AzimuthalIntegrationProfile> last_profile_;
|
||||
std::vector<MergedReflection> pixel_reference_;
|
||||
std::unique_ptr<PixelRefine> pixel_refine_;
|
||||
std::unique_ptr<BraggPrediction> pixel_pred_;
|
||||
|
||||
void EnsurePixelRefine_i();
|
||||
bool BuildPixelSeed_i(PixelRefineData &d, const PixelRefineParams &p, QString &reason) const;
|
||||
std::shared_ptr<SimpleImage> WrapFloatImage_i(const std::vector<float> &img) const;
|
||||
// Turn a predicted image into the squared residual |predicted - measured|^2 in
|
||||
// place. Masked/saturated pixels become NaN (rendered as a gap: no comparison
|
||||
// possible), not 0.
|
||||
void SquaredResidualWithImage_i(std::vector<float> &pred) const;
|
||||
// Mark masked/saturated pixels of the current image as NaN (gap) in a float
|
||||
// image, leaving the rest untouched (used for the chi^2 view).
|
||||
void MaskMeasuredSentinels_i(std::vector<float> &img) const;
|
||||
// Build the per-reflection shoebox rectangles for the last refine/preview.
|
||||
QVector<QRect> BuildShoeboxes_i(const PixelRefineData &data) const;
|
||||
// Assemble the per-reflection table + per-image summary for the table window.
|
||||
PixelRefineReport BuildReport_i(const PixelRefineData &data) const;
|
||||
// Build the float image to display for the given PixelRefineParams::DisplayMode.
|
||||
std::vector<float> BuildDisplayImage_i(const PixelRefineData &data, int display_mode) const;
|
||||
|
||||
std::unique_ptr<ROIElement> roi;
|
||||
|
||||
@@ -145,14 +121,6 @@ signals:
|
||||
void fileLoadError(QString title, QString message);
|
||||
void fileLoadRetryStatus(bool active, QString message);
|
||||
|
||||
// PixelRefine (experimental)
|
||||
void predictedImageReady(std::shared_ptr<const SimpleImage> image);
|
||||
void predictedShoeboxes(QVector<QRect> boxes); // per-reflection optimization windows
|
||||
void pixelRefineResidual(double cost, double cc, int64_t n_reflections);
|
||||
void pixelRefineParamsRefined(PixelRefineParams params);
|
||||
void pixelRefineReport(PixelRefineReport report); // per-reflection table + summary
|
||||
void pixelRefineStatus(QString message);
|
||||
|
||||
public:
|
||||
JFJochImageReadingWorker(const SpotFindingSettings &settings, const DiffractionExperiment& experiment, QObject *parent = nullptr);
|
||||
~JFJochImageReadingWorker() override = default;
|
||||
@@ -189,9 +157,4 @@ public slots:
|
||||
void LoadCalibration(QString dataset);
|
||||
void setAutoLoadMode(AutoloadMode mode);
|
||||
void setAutoLoadJump(int64_t val);
|
||||
|
||||
// PixelRefine (experimental)
|
||||
void LoadReference(QString path);
|
||||
void PixelRefinePreview(PixelRefineParams params);
|
||||
void PixelRefineRun(PixelRefineParams params);
|
||||
};
|
||||
|
||||
@@ -25,8 +25,6 @@
|
||||
#include "toolbar/JFJochViewerToolbarImage.h"
|
||||
#include "windows/JFJoch2DAzintImageWindow.h"
|
||||
#include "windows/JFJochAzIntWindow.h"
|
||||
#include "windows/JFJochPixelRefineWindow.h"
|
||||
#include "windows/JFJochPixelRefineTableWindow.h"
|
||||
#include "windows/JFJochMagnifierWindow.h"
|
||||
#include "image_viewer/JFJochImage.h"
|
||||
#include "image_viewer/JFJochSimpleImage.h"
|
||||
@@ -108,8 +106,6 @@ JFJochViewerWindow::JFJochViewerWindow(QWidget *parent, bool dbus, const QString
|
||||
|
||||
auto azintWindow = new JFJochAzIntWindow(experiment.GetAzimuthalIntegrationSettings(), this);
|
||||
auto azintImageWindow = new JFJoch2DAzintImageWindow(this);
|
||||
auto pixelRefineWindow = new JFJochPixelRefineWindow(this);
|
||||
auto pixelRefineTableWindow = new JFJochPixelRefineTableWindow(this);
|
||||
auto magnifierWindow = new JFJochMagnifierWindow(this);
|
||||
|
||||
menuBar->AddWindowEntry(tableWindow, "Image list");
|
||||
@@ -121,8 +117,6 @@ JFJochViewerWindow::JFJochViewerWindow(QWidget *parent, bool dbus, const QString
|
||||
menuBar->AddWindowEntry(reciprocalWindow, "Reciprocal space viewer");
|
||||
menuBar->AddWindowEntry(azintWindow, "Azimuthal integration settings");
|
||||
menuBar->AddWindowEntry(azintImageWindow, "Azimuthal integration 2D image");
|
||||
menuBar->AddWindowEntry(pixelRefineWindow, "PixelRefine (experimental)");
|
||||
menuBar->AddWindowEntry(pixelRefineTableWindow, "PixelRefine reflections");
|
||||
menuBar->AddWindowEntry(magnifierWindow, "Magnifier");
|
||||
|
||||
if (dbus) {
|
||||
@@ -344,38 +338,6 @@ JFJochViewerWindow::JFJochViewerWindow(QWidget *parent, bool dbus, const QString
|
||||
connect(azintImageWindow, &JFJoch2DAzintImageWindow::zoomOnBin,
|
||||
viewer, &JFJochDiffractionImage::centerOnSpot);
|
||||
|
||||
// --- PixelRefine (experimental) ---
|
||||
connect(reading_worker, &JFJochImageReadingWorker::imageLoaded,
|
||||
pixelRefineWindow, &JFJochHelperWindow::imageLoaded);
|
||||
connect(pixelRefineWindow, &JFJochPixelRefineWindow::paramsChanged,
|
||||
reading_worker, &JFJochImageReadingWorker::PixelRefinePreview);
|
||||
connect(pixelRefineWindow, &JFJochPixelRefineWindow::refineRequested,
|
||||
reading_worker, &JFJochImageReadingWorker::PixelRefineRun);
|
||||
connect(pixelRefineWindow, &JFJochPixelRefineWindow::loadReferenceRequested,
|
||||
reading_worker, &JFJochImageReadingWorker::LoadReference);
|
||||
connect(reading_worker, &JFJochImageReadingWorker::predictedImageReady,
|
||||
pixelRefineWindow, &JFJochPixelRefineWindow::setPredictedImage);
|
||||
connect(reading_worker, &JFJochImageReadingWorker::predictedShoeboxes,
|
||||
pixelRefineWindow->imageView(), &JFJochSimpleImage::setShoeboxes);
|
||||
connect(reading_worker, &JFJochImageReadingWorker::pixelRefineResidual,
|
||||
pixelRefineWindow, &JFJochPixelRefineWindow::setResidual);
|
||||
connect(reading_worker, &JFJochImageReadingWorker::pixelRefineParamsRefined,
|
||||
pixelRefineWindow, &JFJochPixelRefineWindow::setRefinedParams);
|
||||
connect(reading_worker, &JFJochImageReadingWorker::pixelRefineStatus,
|
||||
pixelRefineWindow, &JFJochPixelRefineWindow::setStatus);
|
||||
|
||||
// Reflection-table window: refreshed on every preview/refine, raised by button.
|
||||
connect(reading_worker, &JFJochImageReadingWorker::pixelRefineReport,
|
||||
pixelRefineTableWindow, &JFJochPixelRefineTableWindow::setReport);
|
||||
connect(pixelRefineWindow, &JFJochPixelRefineWindow::showTableRequested,
|
||||
pixelRefineTableWindow, &JFJochHelperWindow::open);
|
||||
|
||||
// Lock the predicted-image viewport to the original image (both directions).
|
||||
connect(viewer, &JFJochImage::viewportChanged,
|
||||
pixelRefineWindow->imageView(), &JFJochImage::applyViewport);
|
||||
connect(pixelRefineWindow->imageView(), &JFJochImage::viewportChanged,
|
||||
viewer, &JFJochImage::applyViewport);
|
||||
|
||||
// --- Magnifier ---
|
||||
connect(reading_worker, &JFJochImageReadingWorker::imageLoaded,
|
||||
magnifierWindow, &JFJochHelperWindow::imageLoaded);
|
||||
|
||||
@@ -55,9 +55,8 @@ protected:
|
||||
|
||||
QColor feature_color = Qt::magenta;
|
||||
|
||||
// Decimal places for non-integer per-pixel value labels. Float images (e.g. the
|
||||
// PixelRefine prediction) are unreadable with many decimals, so subclasses can
|
||||
// lower this.
|
||||
// Decimal places for non-integer per-pixel value labels. Float images are
|
||||
// unreadable with many decimals, so subclasses can lower this.
|
||||
int label_decimals_ = 3;
|
||||
|
||||
float foreground = 10.0;
|
||||
|
||||
@@ -40,30 +40,6 @@ void JFJochSimpleImage::setImage(std::shared_ptr<const SimpleImage> img) {
|
||||
}
|
||||
}
|
||||
|
||||
void JFJochSimpleImage::setShoeboxes(QVector<QRect> boxes) {
|
||||
shoeboxes_ = std::move(boxes);
|
||||
// Redraw overlays on the current image (no-op if no image yet).
|
||||
updateOverlay();
|
||||
}
|
||||
|
||||
void JFJochSimpleImage::addCustomOverlay() {
|
||||
if (shoeboxes_.isEmpty() || !scene())
|
||||
return;
|
||||
|
||||
// Cosmetic 1-px outline so the box edges stay thin at any zoom; only draw the
|
||||
// ones currently in view (there can be hundreds of reflections).
|
||||
const QRectF visibleRect = mapToScene(viewport()->geometry()).boundingRect();
|
||||
QPen pen(QColor(0, 220, 255), 0); // cyan, distinct from the prediction colours
|
||||
pen.setCosmetic(true);
|
||||
|
||||
for (const QRect &b : shoeboxes_) {
|
||||
const QRectF r(b.x(), b.y(), b.width(), b.height());
|
||||
if (!visibleRect.intersects(r))
|
||||
continue;
|
||||
auto *item = scene()->addRect(r, pen);
|
||||
addOverlayItem(item);
|
||||
}
|
||||
}
|
||||
|
||||
void JFJochSimpleImage::mouseHover(QMouseEvent *event) {
|
||||
if (image_) {
|
||||
|
||||
@@ -20,20 +20,14 @@ class JFJochSimpleImage : public JFJochImage {
|
||||
|
||||
std::shared_ptr<const SimpleImage> image_;
|
||||
|
||||
// Per-reflection shoebox rectangles (pixel coordinates) to overlay: the pixels
|
||||
// PixelRefine actually summed over. Empty = nothing drawn.
|
||||
QVector<QRect> shoeboxes_;
|
||||
|
||||
// Prepare image
|
||||
template<class T>
|
||||
void loadImageInternal(const uint8_t *input);
|
||||
void loadImageInternal();
|
||||
|
||||
void mouseHover(QMouseEvent *event) override;
|
||||
void addCustomOverlay() override;
|
||||
public:
|
||||
explicit JFJochSimpleImage(QWidget *parent = nullptr);
|
||||
public slots:
|
||||
void setImage(std::shared_ptr<const SimpleImage> img);
|
||||
void setShoeboxes(QVector<QRect> boxes);
|
||||
};
|
||||
|
||||
@@ -1,101 +0,0 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include "JFJochPixelRefineTableWindow.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <QVBoxLayout>
|
||||
#include <QHeaderView>
|
||||
|
||||
JFJochPixelRefineTableWindow::JFJochPixelRefineTableWindow(QWidget *parent)
|
||||
: JFJochHelperWindow(parent) {
|
||||
setWindowTitle("PixelRefine reflections");
|
||||
resize(950, 600);
|
||||
|
||||
auto central = new QWidget(this);
|
||||
setCentralWidget(central);
|
||||
auto layout = new QVBoxLayout(central);
|
||||
|
||||
m_summary = new QLabel(tr("No PixelRefine result yet."), this);
|
||||
m_summary->setTextFormat(Qt::RichText);
|
||||
m_summary->setWordWrap(true);
|
||||
layout->addWidget(m_summary);
|
||||
|
||||
m_table = new QTableView(this);
|
||||
m_model = new QStandardItemModel(this);
|
||||
setupModel();
|
||||
|
||||
m_proxy = new QSortFilterProxyModel(this);
|
||||
m_proxy->setSourceModel(m_model);
|
||||
m_proxy->setSortRole(Qt::UserRole); // numeric sort on the underlying values
|
||||
|
||||
m_table->setModel(m_proxy);
|
||||
m_table->setEditTriggers(QAbstractItemView::NoEditTriggers);
|
||||
m_table->setSortingEnabled(true);
|
||||
m_table->sortByColumn(6, Qt::DescendingOrder); // default: I desc
|
||||
m_table->verticalHeader()->setVisible(false);
|
||||
m_table->horizontalHeader()->setSortIndicatorShown(true);
|
||||
m_table->horizontalHeader()->setSectionResizeMode(QHeaderView::Stretch);
|
||||
m_table->setStyleSheet("background-color: white;");
|
||||
layout->addWidget(m_table);
|
||||
}
|
||||
|
||||
void JFJochPixelRefineTableWindow::setupModel() {
|
||||
const QStringList headers = {
|
||||
"h", "k", "l", "d [Å]", "Compl.", "Part.",
|
||||
"I", "Sigma", "Est. I_true", "Ref. I_true"
|
||||
};
|
||||
m_model->setColumnCount(headers.size());
|
||||
for (int i = 0; i < headers.size(); ++i)
|
||||
m_model->setHeaderData(i, Qt::Horizontal, headers[i]);
|
||||
}
|
||||
|
||||
static QStandardItem *numItem(double value, int decimals) {
|
||||
auto *it = new QStandardItem();
|
||||
const QString text = std::isfinite(value) ? QString::number(value, 'f', decimals)
|
||||
: QStringLiteral("—");
|
||||
it->setData(text, Qt::DisplayRole);
|
||||
it->setData(value, Qt::UserRole);
|
||||
it->setTextAlignment(Qt::AlignRight | Qt::AlignVCenter);
|
||||
return it;
|
||||
}
|
||||
|
||||
static QStandardItem *intItem(int value) {
|
||||
auto *it = new QStandardItem();
|
||||
it->setData(static_cast<qlonglong>(value), Qt::DisplayRole);
|
||||
it->setData(static_cast<qlonglong>(value), Qt::UserRole);
|
||||
it->setTextAlignment(Qt::AlignRight | Qt::AlignVCenter);
|
||||
return it;
|
||||
}
|
||||
|
||||
void JFJochPixelRefineTableWindow::setReport(PixelRefineReport report) {
|
||||
// --- per-image summary: PixelRefine vs the standard pipeline ---------------
|
||||
auto fmt = [](double v, int dec) {
|
||||
return std::isfinite(v) ? QString::number(v, 'f', dec) : QStringLiteral("—");
|
||||
};
|
||||
auto pct = [](double v) {
|
||||
return std::isfinite(v) ? QString::number(v * 100.0, 'f', 1) + "%" : QStringLiteral("—");
|
||||
};
|
||||
m_summary->setText(
|
||||
tr("<b>PixelRefine:</b> scale G = %1, B = %2 Ų, CC = %3 (%4 refl)"
|
||||
" | "
|
||||
"<b>Pipeline:</b> scale G = %5, B = %6 Ų, CC = %7")
|
||||
.arg(fmt(report.pr_G, 4), fmt(report.pr_B, 1), pct(report.pr_cc))
|
||||
.arg(report.pr_cc_n)
|
||||
.arg(fmt(report.pipe_G, 4), fmt(report.pipe_B, 1), pct(report.pipe_cc)));
|
||||
|
||||
// --- per-reflection rows ---------------------------------------------------
|
||||
m_model->removeRows(0, m_model->rowCount());
|
||||
for (const auto &r : report.rows) {
|
||||
QList<QStandardItem *> row;
|
||||
row << intItem(r.h) << intItem(r.k) << intItem(r.l);
|
||||
row << numItem(r.d, 3);
|
||||
row << numItem(r.completeness, 2);
|
||||
row << numItem(r.partiality, 2);
|
||||
row << numItem(r.I, 1);
|
||||
row << numItem(r.sigma, 1);
|
||||
row << numItem(r.I_true_est, 1);
|
||||
row << numItem(r.I_true_ref, 1);
|
||||
m_model->appendRow(row);
|
||||
}
|
||||
}
|
||||
@@ -1,34 +0,0 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#pragma once
|
||||
|
||||
#include "JFJochHelperWindow.h"
|
||||
#include "PixelRefineParams.h"
|
||||
|
||||
#include <QLabel>
|
||||
#include <QTableView>
|
||||
#include <QStandardItemModel>
|
||||
#include <QSortFilterProxyModel>
|
||||
|
||||
// Reflection table for the experimental PixelRefine path. Opened from a button on
|
||||
// the PixelRefine window and refreshed on every preview/refine. Each row is one
|
||||
// matched reflection (the two partiality-style fractions, the fitted intensity,
|
||||
// and the estimated vs reference full intensity); the header line compares
|
||||
// PixelRefine's per-image scale/B/CC with the standard ScaleOnTheFly pipeline.
|
||||
class JFJochPixelRefineTableWindow : public JFJochHelperWindow {
|
||||
Q_OBJECT
|
||||
|
||||
QLabel *m_summary = nullptr;
|
||||
QTableView *m_table = nullptr;
|
||||
QStandardItemModel *m_model = nullptr;
|
||||
QSortFilterProxyModel *m_proxy = nullptr;
|
||||
|
||||
void setupModel();
|
||||
|
||||
public:
|
||||
explicit JFJochPixelRefineTableWindow(QWidget *parent = nullptr);
|
||||
|
||||
public slots:
|
||||
void setReport(PixelRefineReport report);
|
||||
};
|
||||
@@ -1,233 +0,0 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include "JFJochPixelRefineWindow.h"
|
||||
#include "../image_viewer/JFJochSimpleImage.h"
|
||||
|
||||
#include <QWidget>
|
||||
#include <QHBoxLayout>
|
||||
#include <QVBoxLayout>
|
||||
#include <QFormLayout>
|
||||
#include <QGroupBox>
|
||||
#include <QFileDialog>
|
||||
#include <cmath>
|
||||
|
||||
JFJochPixelRefineWindow::JFJochPixelRefineWindow(QWidget *parent)
|
||||
: JFJochHelperWindow(parent) {
|
||||
setWindowTitle("PixelRefine (experimental)");
|
||||
|
||||
auto central = new QWidget(this);
|
||||
setCentralWidget(central);
|
||||
auto layout = new QHBoxLayout(central);
|
||||
|
||||
// --- predicted image (left, expanding) ---------------------------------
|
||||
m_image = new JFJochSimpleImage(this);
|
||||
m_image->setSizePolicy(QSizePolicy::Expanding, QSizePolicy::Expanding);
|
||||
layout->addWidget(m_image, 1);
|
||||
|
||||
// --- control panel (right) ---------------------------------------------
|
||||
auto controls = new QWidget(this);
|
||||
controls->setMinimumWidth(320);
|
||||
auto controlsLayout = new QVBoxLayout(controls);
|
||||
layout->addWidget(controls, 0);
|
||||
|
||||
// --- what the left image shows ------------------------------------------
|
||||
m_displayMode = new QComboBox(this);
|
||||
m_displayMode->addItem(tr("Prediction"));
|
||||
m_displayMode->addItem(tr("Squared difference |pred - image|²"));
|
||||
m_displayMode->addItem(tr("χ² (weighted residual = LSQ cost)"));
|
||||
auto displayForm = new QFormLayout();
|
||||
displayForm->addRow(tr("Display:"), m_displayMode);
|
||||
controlsLayout->addLayout(displayForm);
|
||||
|
||||
auto paramBox = new QGroupBox(tr("Model parameters"), this);
|
||||
auto form = new QFormLayout(paramBox);
|
||||
|
||||
m_R0 = new SliderPlusBox(1e-4, 0.05, 1e-4, 4, this); m_R0->setValue(0.005);
|
||||
m_R1 = new SliderPlusBox(1e-4, 0.05, 1e-4, 4, this); m_R1->setValue(0.005);
|
||||
m_bw = new SliderPlusBox(0.0, 0.05, 1e-4, 4, this); m_bw->setValue(0.0);
|
||||
m_scale = new SliderPlusBox(1e-3, 1e4, 1e-3, 3, this, SliderPlusBox::Logarithmic); m_scale->setValue(1.0);
|
||||
m_B = new SliderPlusBox(0.0, 200.0, 0.1, 1, this); m_B->setValue(0.0);
|
||||
m_beamx = new SliderPlusBox(0.0, 5000.0, 0.5, 1, this); m_beamx->setValue(0.0);
|
||||
m_beamy = new SliderPlusBox(0.0, 5000.0, 0.5, 1, this); m_beamy->setValue(0.0);
|
||||
|
||||
form->addRow(tr("R0 radial [Å⁻¹]:"), m_R0);
|
||||
form->addRow(tr("R1 tangential [Å⁻¹]:"), m_R1);
|
||||
form->addRow(tr("Bandwidth FWHM (Δλ/λ):"), m_bw);
|
||||
form->addRow(tr("Scale G:"), m_scale);
|
||||
form->addRow(tr("B-factor [Ų]:"), m_B);
|
||||
|
||||
m_overrideBeam = new QCheckBox(tr("Override beam centre"), this);
|
||||
form->addRow(QString(), m_overrideBeam);
|
||||
form->addRow(tr("Beam X [px]:"), m_beamx);
|
||||
form->addRow(tr("Beam Y [px]:"), m_beamy);
|
||||
m_beamx->setEnabled(false);
|
||||
m_beamy->setEnabled(false);
|
||||
|
||||
controlsLayout->addWidget(paramBox);
|
||||
|
||||
// --- what "Refine" is allowed to move ----------------------------------
|
||||
auto refBox = new QGroupBox(tr("Refine (Ceres)"), this);
|
||||
auto refLayout = new QVBoxLayout(refBox);
|
||||
m_refOrientation = new QCheckBox(tr("Orientation"), this); m_refOrientation->setChecked(true);
|
||||
m_refCell = new QCheckBox(tr("Unit cell"), this);
|
||||
m_refBeam = new QCheckBox(tr("Beam centre"), this);
|
||||
m_refScale = new QCheckBox(tr("Scale G"), this); m_refScale->setChecked(true);
|
||||
m_refB = new QCheckBox(tr("B-factor"), this);
|
||||
m_refR = new QCheckBox(tr("Widths R0/R1"), this); m_refR->setChecked(true);
|
||||
for (auto *cb : {m_refOrientation, m_refCell, m_refBeam, m_refScale, m_refB, m_refR})
|
||||
refLayout->addWidget(cb);
|
||||
controlsLayout->addWidget(refBox);
|
||||
|
||||
// --- buttons + readouts -------------------------------------------------
|
||||
m_loadRef = new QPushButton(tr("Load reference MTZ…"), this);
|
||||
m_refine = new QPushButton(tr("Refine"), this);
|
||||
m_showTable = new QPushButton(tr("Reflection table…"), this);
|
||||
controlsLayout->addWidget(m_loadRef);
|
||||
controlsLayout->addWidget(m_refine);
|
||||
controlsLayout->addWidget(m_showTable);
|
||||
|
||||
m_residual = new QLabel(tr("Residual: —"), this);
|
||||
m_pipelineCC = new QLabel(tr("Pipeline CC (ref): —"), this);
|
||||
m_status = new QLabel(QString(), this);
|
||||
m_status->setWordWrap(true);
|
||||
m_status->setStyleSheet("color: rgb(80, 80, 80);");
|
||||
controlsLayout->addWidget(m_residual);
|
||||
controlsLayout->addWidget(m_pipelineCC);
|
||||
controlsLayout->addWidget(m_status);
|
||||
controlsLayout->addStretch(1);
|
||||
|
||||
// --- debounce timer for live preview -----------------------------------
|
||||
m_debounce = new QTimer(this);
|
||||
m_debounce->setSingleShot(true);
|
||||
m_debounce->setInterval(150);
|
||||
connect(m_debounce, &QTimer::timeout, this, [this] {
|
||||
emit paramsChanged(currentParams());
|
||||
});
|
||||
|
||||
for (auto *s : {m_R0, m_R1, m_bw, m_scale, m_B, m_beamx, m_beamy})
|
||||
connect(s, &SliderPlusBox::valueChanged, this, [this](double) { onControlChanged(); });
|
||||
|
||||
connect(m_displayMode, &QComboBox::currentIndexChanged, this, [this](int) { onControlChanged(); });
|
||||
|
||||
connect(m_overrideBeam, &QCheckBox::toggled, this, [this](bool on) {
|
||||
m_beamx->setEnabled(on);
|
||||
m_beamy->setEnabled(on);
|
||||
onControlChanged();
|
||||
});
|
||||
|
||||
connect(m_loadRef, &QPushButton::clicked, this, [this] {
|
||||
const QString path = QFileDialog::getOpenFileName(
|
||||
this, tr("Load reference MTZ"), QString(), tr("MTZ files (*.mtz);;All files (*)"));
|
||||
if (!path.isEmpty())
|
||||
emit loadReferenceRequested(path);
|
||||
});
|
||||
|
||||
connect(m_showTable, &QPushButton::clicked, this, [this] {
|
||||
emit showTableRequested();
|
||||
});
|
||||
|
||||
connect(m_refine, &QPushButton::clicked, this, [this] {
|
||||
// Cancel any pending live-preview: otherwise a debounce armed by a slider
|
||||
// move just before this click fires after the refine and overwrites the
|
||||
// refined residual/preview with the stale pre-refine slider values.
|
||||
m_debounce->stop();
|
||||
PixelRefineParams p = currentParams();
|
||||
p.max_iterations = 5;
|
||||
emit refineRequested(p);
|
||||
});
|
||||
}
|
||||
|
||||
PixelRefineParams JFJochPixelRefineWindow::currentParams() const {
|
||||
PixelRefineParams p;
|
||||
p.R0 = m_R0->value();
|
||||
p.R1 = m_R1->value();
|
||||
p.bandwidth_fwhm = m_bw->value();
|
||||
p.scale_factor = m_scale->value();
|
||||
p.B_factor = m_B->value();
|
||||
if (m_overrideBeam->isChecked()) {
|
||||
p.beam_x = m_beamx->value();
|
||||
p.beam_y = m_beamy->value();
|
||||
} else {
|
||||
p.beam_x = NAN;
|
||||
p.beam_y = NAN;
|
||||
}
|
||||
p.refine_orientation = m_refOrientation->isChecked();
|
||||
p.refine_unit_cell = m_refCell->isChecked();
|
||||
p.refine_beam_center = m_refBeam->isChecked();
|
||||
p.refine_scale = m_refScale->isChecked();
|
||||
p.refine_B = m_refB->isChecked();
|
||||
p.refine_R = m_refR->isChecked();
|
||||
p.display_mode = m_displayMode->currentIndex();
|
||||
return p;
|
||||
}
|
||||
|
||||
void JFJochPixelRefineWindow::onControlChanged() {
|
||||
if (m_suppress)
|
||||
return;
|
||||
m_debounce->start();
|
||||
}
|
||||
|
||||
void JFJochPixelRefineWindow::imageLoaded(std::shared_ptr<const JFJochReaderImage> image) {
|
||||
if (!image)
|
||||
return;
|
||||
|
||||
// Initialise the beam-centre sliders from the geometry once.
|
||||
if (!m_beamInit) {
|
||||
const auto geom = image->Dataset().experiment.GetDiffractionGeometry();
|
||||
m_beamx->setMax(static_cast<double>(image->Dataset().experiment.GetXPixelsNum()));
|
||||
m_beamy->setMax(static_cast<double>(image->Dataset().experiment.GetYPixelsNum()));
|
||||
m_suppress = true;
|
||||
m_beamx->setValue(geom.GetBeamX_pxl());
|
||||
m_beamy->setValue(geom.GetBeamY_pxl());
|
||||
m_suppress = false;
|
||||
m_beamInit = true;
|
||||
}
|
||||
|
||||
// Show the standard ScaleOnTheFly pipeline's per-image CC vs the reference (set
|
||||
// on the message during analysis when a reference is loaded), as a baseline to
|
||||
// compare PixelRefine against.
|
||||
const auto pipeline_cc = image->ImageData().image_scale_cc;
|
||||
if (pipeline_cc.has_value() && std::isfinite(pipeline_cc.value()))
|
||||
m_pipelineCC->setText(tr("Pipeline CC (ref): %1%")
|
||||
.arg(pipeline_cc.value() * 100.0, 0, 'f', 1));
|
||||
else
|
||||
m_pipelineCC->setText(tr("Pipeline CC (ref): —"));
|
||||
|
||||
// Request a predicted-image preview for the (re)loaded image. Without this the
|
||||
// preview only refreshed on a slider change, so opening/reanalyzing an image
|
||||
// left the predicted view empty. The worker no-ops if there is no reference or
|
||||
// the image is not integrated yet.
|
||||
onControlChanged();
|
||||
}
|
||||
|
||||
void JFJochPixelRefineWindow::setPredictedImage(std::shared_ptr<const SimpleImage> image) {
|
||||
m_image->setImage(std::move(image));
|
||||
}
|
||||
|
||||
void JFJochPixelRefineWindow::setResidual(double cost, double cc, int64_t n_reflections) {
|
||||
const QString cc_str = std::isfinite(cc) ? QString::number(cc * 100.0, 'f', 1) + "%"
|
||||
: QStringLiteral("—");
|
||||
m_residual->setText(tr("Residual: %1 CC: %2 (%3 reflections)")
|
||||
.arg(cost, 0, 'g', 6)
|
||||
.arg(cc_str)
|
||||
.arg(n_reflections));
|
||||
}
|
||||
|
||||
void JFJochPixelRefineWindow::setRefinedParams(PixelRefineParams params) {
|
||||
m_suppress = true;
|
||||
m_R0->setValue(params.R0);
|
||||
m_R1->setValue(params.R1);
|
||||
m_bw->setValue(params.bandwidth_fwhm);
|
||||
m_scale->setValue(params.scale_factor);
|
||||
m_B->setValue(params.B_factor);
|
||||
if (std::isfinite(params.beam_x) && std::isfinite(params.beam_y)) {
|
||||
m_beamx->setValue(params.beam_x);
|
||||
m_beamy->setValue(params.beam_y);
|
||||
}
|
||||
m_suppress = false;
|
||||
}
|
||||
|
||||
void JFJochPixelRefineWindow::setStatus(QString message) {
|
||||
m_status->setText(message);
|
||||
}
|
||||
@@ -1,79 +0,0 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#pragma once
|
||||
|
||||
#include "JFJochHelperWindow.h"
|
||||
#include "PixelRefineParams.h"
|
||||
#include "../SimpleImage.h"
|
||||
#include "../widgets/SliderPlusBox.h"
|
||||
|
||||
#include <memory>
|
||||
#include <QCheckBox>
|
||||
#include <QComboBox>
|
||||
#include <QLabel>
|
||||
#include <QPushButton>
|
||||
#include <QTimer>
|
||||
|
||||
class JFJochSimpleImage;
|
||||
|
||||
// Experimental PixelRefine control window: sliders for the forward-model
|
||||
// parameters, a live (debounced) predicted-image preview, a residual readout,
|
||||
// "Load reference" and "Refine" buttons. The predicted-image view is exposed so
|
||||
// the main window can lock its viewport to the original image.
|
||||
class JFJochPixelRefineWindow : public JFJochHelperWindow {
|
||||
Q_OBJECT
|
||||
|
||||
JFJochSimpleImage *m_image;
|
||||
|
||||
SliderPlusBox *m_R0;
|
||||
SliderPlusBox *m_R1;
|
||||
SliderPlusBox *m_bw;
|
||||
SliderPlusBox *m_scale;
|
||||
SliderPlusBox *m_B;
|
||||
SliderPlusBox *m_beamx;
|
||||
SliderPlusBox *m_beamy;
|
||||
|
||||
QComboBox *m_displayMode; // Prediction vs. Difference (prediction - image)
|
||||
|
||||
QCheckBox *m_overrideBeam;
|
||||
QCheckBox *m_refOrientation;
|
||||
QCheckBox *m_refCell;
|
||||
QCheckBox *m_refBeam;
|
||||
QCheckBox *m_refScale;
|
||||
QCheckBox *m_refB;
|
||||
QCheckBox *m_refR;
|
||||
|
||||
QLabel *m_residual;
|
||||
QLabel *m_pipelineCC; // first-image CC vs reference from the standard ScaleOnTheFly pipeline
|
||||
QLabel *m_status;
|
||||
QPushButton *m_loadRef;
|
||||
QPushButton *m_refine;
|
||||
QPushButton *m_showTable;
|
||||
|
||||
QTimer *m_debounce;
|
||||
bool m_suppress = false; // guard while pushing refined params into sliders
|
||||
bool m_beamInit = false; // beam-centre sliders initialised from geometry
|
||||
|
||||
PixelRefineParams currentParams() const;
|
||||
void onControlChanged();
|
||||
|
||||
public:
|
||||
explicit JFJochPixelRefineWindow(QWidget *parent = nullptr);
|
||||
|
||||
JFJochSimpleImage *imageView() const { return m_image; }
|
||||
|
||||
void imageLoaded(std::shared_ptr<const JFJochReaderImage> image) override;
|
||||
|
||||
signals:
|
||||
void paramsChanged(PixelRefineParams params); // debounced live preview
|
||||
void refineRequested(PixelRefineParams params); // "Refine" button
|
||||
void loadReferenceRequested(QString path); // "Load reference" button
|
||||
void showTableRequested(); // "Reflection table" button
|
||||
|
||||
public slots:
|
||||
void setPredictedImage(std::shared_ptr<const SimpleImage> image);
|
||||
void setResidual(double cost, double cc, int64_t n_reflections);
|
||||
void setRefinedParams(PixelRefineParams params);
|
||||
void setStatus(QString message);
|
||||
};
|
||||
@@ -1,66 +0,0 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#pragma once
|
||||
|
||||
#include <cmath>
|
||||
#include <cstdint>
|
||||
#include <vector>
|
||||
#include <QMetaType>
|
||||
|
||||
// Parameters exchanged between the PixelRefine window (sliders/buttons) and the
|
||||
// reading worker. bandwidth here is the *FWHM* of dlambda/lambda (user-facing);
|
||||
// the worker converts it to the sigma that PixelRefineData expects. beam_x/beam_y
|
||||
// are NaN to mean "keep the current refined geometry".
|
||||
struct PixelRefineParams {
|
||||
double R0 = 0.005; // radial / partiality width (A^-1)
|
||||
double R1 = 0.005; // tangential / profile width (A^-1)
|
||||
double bandwidth_fwhm = 0.0; // relative bandwidth FWHM (dlambda/lambda)
|
||||
double scale_factor = 1.0; // overall scale G
|
||||
double B_factor = 0.0; // Debye-Waller B (A^2)
|
||||
double beam_x = NAN; // detector beam centre X (px); NaN = keep current
|
||||
double beam_y = NAN; // detector beam centre Y (px); NaN = keep current
|
||||
|
||||
// What the "Refine" button is allowed to move (ignored by the preview path).
|
||||
bool refine_orientation = true;
|
||||
bool refine_unit_cell = false;
|
||||
bool refine_beam_center = false;
|
||||
bool refine_scale = true;
|
||||
bool refine_B = false;
|
||||
bool refine_R = true;
|
||||
|
||||
int max_iterations = 3; // <=0 means evaluate-only (preview / residual)
|
||||
|
||||
// Display only (no effect on the fit): what the preview/refine image shows.
|
||||
enum DisplayMode : int {
|
||||
Prediction = 0, // forward-model image
|
||||
SquaredDifference = 1, // |prediction - measured|^2 (raw, unweighted)
|
||||
ChiSquared = 2 // ((prediction - measured)/sigma)^2 = the LSQ cost density
|
||||
};
|
||||
int display_mode = Prediction;
|
||||
};
|
||||
|
||||
Q_DECLARE_METATYPE(PixelRefineParams)
|
||||
|
||||
// One PixelRefine result, shipped from the worker to the reflection-table window:
|
||||
// a per-image summary that puts PixelRefine's scale/B/CC next to the standard
|
||||
// ScaleOnTheFly pipeline's, plus one row per matched reflection.
|
||||
struct PixelRefineReport {
|
||||
// Per-image summary (NaN = not available).
|
||||
double pr_G = NAN, pr_B = NAN, pr_cc = NAN; // PixelRefine
|
||||
int64_t pr_cc_n = 0;
|
||||
double pipe_G = NAN, pipe_B = NAN, pipe_cc = NAN; // ScaleOnTheFly pipeline baseline
|
||||
|
||||
struct Row {
|
||||
int h = 0, k = 0, l = 0;
|
||||
double d = 0.0;
|
||||
double completeness = 1.0; // spot footprint on live pixels (1 = not clipped)
|
||||
double partiality = 1.0; // recorded rocking fraction
|
||||
double I = 0.0, sigma = 0.0;
|
||||
double I_true_est = NAN; // r.I * image_scale_corr (this image's estimate)
|
||||
double I_true_ref = NAN; // reference (merged) intensity
|
||||
};
|
||||
std::vector<Row> rows;
|
||||
};
|
||||
|
||||
Q_DECLARE_METATYPE(PixelRefineReport)
|
||||
Reference in New Issue
Block a user