diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index ccf084fe..f7cb25c1 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -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 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); diff --git a/image_analysis/pixel_refinement/FACTORED_MODEL.md b/image_analysis/pixel_refinement/FACTORED_MODEL.md index 944514a2..cb398977 100644 --- a/image_analysis/pixel_refinement/FACTORED_MODEL.md +++ b/image_analysis/pixel_refinement/FACTORED_MODEL.md @@ -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 diff --git a/image_analysis/pixel_refinement/FINDINGS-2026-06.md b/image_analysis/pixel_refinement/FINDINGS-2026-06.md new file mode 100644 index 00000000..a7345de4 --- /dev/null +++ b/image_analysis/pixel_refinement/FINDINGS-2026-06.md @@ -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. diff --git a/image_analysis/pixel_refinement/METHODS.md b/image_analysis/pixel_refinement/METHODS.md index 8604ffc7..1c1fd9e1 100644 --- a/image_analysis/pixel_refinement/METHODS.md +++ b/image_analysis/pixel_refinement/METHODS.md @@ -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. diff --git a/image_analysis/pixel_refinement/PixelRefine.cpp b/image_analysis/pixel_refinement/PixelRefine.cpp index 42e8772f..75facde4 100644 --- a/image_analysis/pixel_refinement/PixelRefine.cpp +++ b/image_analysis/pixel_refinement/PixelRefine.cpp @@ -6,11 +6,11 @@ #include #include #include +#include #include #include #include -#include #include #include "../geom_refinement/LatticeReduction.h" @@ -24,7 +24,6 @@ struct PixelObs { double x, y; // detector pixel coordinate double Iobs; // raw pixel value (signal + background) double Ibkg; // local background estimate (per-shoebox level, raw counts) - double weight; // 1 / sigma_pixel }; // One reflection together with the pixels of its shoebox. @@ -36,17 +35,10 @@ struct ReflGroup { double pol; // per-reflection polarization correction (raw = true * pol) double Ibkg; // local flat background (raw counts, constant over the shoebox) double predicted_x, predicted_y; - double R1_eff = 0.0; // tangential profile width to use (Term 2; 0 => fall back to data.R[1]) - double dcx = 0.0, dcy = 0.0; // Term 3: profile recentre shift (observed centroid - predicted) + double R1_eff = 0.0; // tangential profile width to use (Term 2) std::vector pixels; }; -double SafeInv(double x, double fallback) { - if (!std::isfinite(x) || std::fabs(x) < 1e-30) - return fallback; - return 1.0 / x; -} - // Median of a vector (in place, partially reorders it). double MedianInPlace(std::vector &v) { if (v.empty()) @@ -90,8 +82,7 @@ std::vector BuildSpotMask(const std::vector &predicted, int // Square shoebox bounds (inclusive) around a predicted spot, clamped to the // detector. The centre is rounded to the nearest pixel with std::lround so the // signal box is centred identically to the spot-core mask (BuildSpotMask) and -// the local-background ring (EstimateLocalBackground), which also lround. Used by -// Run and the diagnostic renderers so all three share one shoebox definition. +// the local-background ring (EstimateLocalBackground), which also lround. struct ShoeboxBox { int min_x, max_x, min_y, max_y; }; ShoeboxBox ShoeboxBounds(double px, double py, int radius, size_t xpixel, size_t ypixel) { const int cx = static_cast(std::lround(px)); @@ -284,6 +275,30 @@ bool PredictedNode(const T *p0, const T *p1, const T *p2, return true; } +// Geometry terms for one shoebox pixel under a FIXED geometry (PixelRefine no +// longer refines geometry, so this is a plain double evaluation, not a Ceres cost): +// q_sq = |g_hkl|^2 (predicted node, for the B-factor) +// eps_radial = deviation along the Ewald normal (the partiality direction) +// eps_tang_sq = squared deviation in the tangent plane (the profile direction) +bool GeometryProbe(double obs_x, double obs_y, double lambda, double pixel_size, + int h, int k, int l, gemmi::CrystalSystem symmetry, + const double beam[2], double dist_mm, const double detector_rot[2], + const double p0[3], const double p1[3], const double p2[3], + double &q_sq, double &eps_radial, double &eps_tang_sq) { + const double inv_lambda = 1.0 / lambda; + Eigen::Vector3d e_obs; + ObservedRecip(beam, &dist_mm, detector_rot, obs_x, obs_y, pixel_size, inv_lambda, e_obs); + + Eigen::Vector3d e_pred, n_radial; + if (!PredictedNode(p0, p1, p2, h, k, l, symmetry, inv_lambda, e_pred, n_radial, q_sq)) + return false; + + const Eigen::Vector3d delta_q = e_obs - e_pred; + eps_radial = delta_q.dot(n_radial); + eps_tang_sq = (delta_q - eps_radial * n_radial).squaredNorm(); + return true; +} + // Pulls a scalar parameter towards an expected value with a fixed weight (the // data-scaled prior). Identical in spirit to ScaleOnTheFly's regularizer: it is what // keeps the per-image scale G from wandering on weakly-constrained images and @@ -299,250 +314,30 @@ struct ScalarRegularizer { double expected; }; -// Anchors the orientation (angle-axis vector) to its pre-refinement value with a -// data-scaled weight. Without it the three orientation DOF chase the sparse signal -// (and the few noisy background pixels) and the per-image intensities collapse; -// with it the fit can only make a small, signal-supported sub-spot correction - the -// push that brings slightly-misaligned high-resolution reflections onto their -// shoeboxes. Mirrors the G/B regularizers in ScaleOnTheFly. -struct OrientationRegularizer { - OrientationRegularizer(double weight, const double prior[3]) : weight(weight) { - for (int i = 0; i < 3; ++i) - prior_[i] = prior[i]; - } +// Term 1 of the factored likelihood (FACTORED_MODEL.md): the per-reflection +// *intensity* (0th-moment) residual. The profile-fit amplitude J should equal the +// scaled reference J_model = G * exp(-B/4d^2) * partiality * pol * I_ref. One scalar +// residual per reflection, weighted by the model-expected (Fisher) sigma_J. This is +// the scaling residual - integration and scaling become one objective, and the empty +// pixels (which make no residual of their own) stop dominating the fit. Geometry is +// fixed, so J, partiality and sigma_J are constants and only G and B are free. +struct IntensityResidual { + IntensityResidual(double J, double sigma_J, double partiality, double pol, + double I_ref, double inv_4d2) + : J(J), inv_sigma(1.0 / sigma_J), partiality(partiality), pol(pol), + I_ref(I_ref), inv_4d2(inv_4d2) {} template - bool operator()(const T *p0, T *residual) const { - for (int i = 0; i < 3; ++i) - residual[i] = T(weight) * (p0[i] - T(prior_[i])); + bool operator()(const T *const G, const T *const B, T *residual) const { + const T B_term = ceres::exp(-B[0] * T(inv_4d2)); + const T J_model = G[0] * B_term * T(partiality) * T(pol) * T(I_ref); + residual[0] = (J_model - T(J)) * T(inv_sigma); return true; } - double weight; - double prior_[3]; + double J, inv_sigma, partiality, pol, I_ref, inv_4d2; }; } // namespace -// --------------------------------------------------------------------------- -// Cost functor -// -// I_pred(pixel) = G * Itrue * 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) (partiality: fraction of -// the mosaic blob on the -// Ewald sphere; <= 1) -// P_tangential = exp(-eps_t^2/R1^2) / (pi R1^2) (Gaussian spatial profile -// in the Ewald tangent plane) -// pol = per-reflection polarization correction (raw = true * pol), -// evaluated once at the predicted spot position (as in -// BraggIntegrate2D). 1 if polarization is disabled. -// -// Everything is in *raw* detector counts: there is no per-pixel solid-angle or -// area (Lorentz/Jacobian) weighting - each pixel counts equally, like the normal -// integrator. The tangential factor is what makes this "profile fitting"; the -// 1/(pi R1^2) normalization keeps the profile width R1 from soaking up the -// overall scale G. -// -// X-ray bandwidth: a spread in lambda is a spread in the Ewald-sphere radius, -// i.e. a purely *radial* thickening of the shell. It adds (in quadrature) a -// resolution-dependent term to the radial width: -// R0_eff^2 = R0^2 + R_bw^2 , R_bw^2 = (b*lambda)^2 / (2 d^4) -// where b = relative bandwidth (sigma of dlambda/lambda). R_bw grows like 1/d^2, -// so bandwidth leaves low-resolution spots sharp and smears high-resolution ones -// radially - the pink-beam/DMM signature. R_bw_sq is a fixed per-reflection -// constant (b is known), so R0 keeps meaning "intrinsic" width (mosaic + -// divergence + beam). b = 0 makes R_bw = 0: a monochromatic no-op. -// --------------------------------------------------------------------------- -struct PixelResidual { - PixelResidual(const PixelObs &obs, double Itrue, - double lambda, double pixel_size, - double exp_h, double exp_k, double exp_l, - double R_bw_sq, double pol, - gemmi::CrystalSystem symmetry) - : Itrue(Itrue), Iobs(obs.Iobs), Ibkg(obs.Ibkg), weight(obs.weight), - obs_x(obs.x), obs_y(obs.y), - inv_lambda(1.0 / lambda), pixel_size(pixel_size), - exp_h(exp_h), exp_k(exp_k), exp_l(exp_l), - R_bw_sq(R_bw_sq), pol(pol), symmetry(symmetry) { - if (std::fabs(lambda) < 1e-6) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "Lambda cannot be close to zero"); - } - - // Maps a detector pixel through the current geometry + lattice into the - // reference reciprocal frame and returns: - // q_sq = |g_hkl|^2 (predicted node, for B-factor) - // eps_radial = deviation along Ewald normal (partiality direction) - // eps_tang_sq = squared deviation in the detector-tangential plane (profile) - template - bool GeometryTerms(const T *const beam, - const T *const distance_mm, - const T *const detector_rot, - const T *const p0, - const T *const p1, - const T *const p2, - T &q_sq, T &eps_radial, T &eps_tang_sq) const { - Eigen::Matrix e_obs_recip; - ObservedRecip(beam, distance_mm, detector_rot, - obs_x, obs_y, pixel_size, inv_lambda, - e_obs_recip); - - Eigen::Matrix e_pred_recip, n_radial; - if (!PredictedNode(p0, p1, p2, exp_h, exp_k, exp_l, symmetry, inv_lambda, - e_pred_recip, n_radial, q_sq)) - return false; - - const Eigen::Matrix delta_q = e_obs_recip - e_pred_recip; - eps_radial = delta_q.dot(n_radial); - eps_tang_sq = (delta_q - eps_radial * n_radial).squaredNorm(); - return true; - } - - // Assembles the full model intensity for the pixel from the geometry terms. - template - bool Model(const T *const beam, const T *const distance_mm, - const T *const detector_rot, - const T *const p0, const T *const p1, const T *const p2, - const T *const scale_factor, const T *const B, const T *const R, - T &Ipred) const { - T q_sq, eps_radial, eps_tang_sq; - if (!GeometryTerms(beam, distance_mm, detector_rot, - p0, p1, p2, q_sq, eps_radial, eps_tang_sq)) - return false; - - if (R[0] < T(1e-10) || R[1] < T(1e-10)) - return false; - - const T B_term = ceres::exp(-B[0] * q_sq / T(4.0)); - - // Separable Gaussian spot model: - // radial P_r(e) = exp(-e^2/R0_eff^2) (peak-normalized, in (0,1]) - // tangent g_t(e) = exp(-|e|^2/R1^2) / (pi R1^2) [1/A^-2] - // Every pixel counts equally (no area/Lorentz weighting); the radial factor - // is the still-image partiality (how far the reflection sits from the Ewald - // sphere); the overall scale is carried by the free G. - // - // IMPORTANT: the radial factor MUST use the same convention here as the - // extraction's `partiality` (peak-normalized), otherwise image_scale_corr - // = 1/(partiality*G*B) does not invert the model and a leftover, R0_eff- - // dependent (hence resolution-dependent) factor biases the intensities. - // R0_eff folds in the energy-bandwidth broadening via R_bw_sq. - const T R0_eff_sq = R[0] * R[0] + T(R_bw_sq); - const T P_radial = ceres::exp(-eps_radial * eps_radial / R0_eff_sq); - const T P_tang = ceres::exp(-eps_tang_sq / (R[1] * R[1])) - / (T(M_PI) * R[1] * R[1]); - - const T signal = scale_factor[0] * T(Itrue) * B_term * P_radial * P_tang * T(pol); - Ipred = signal + T(Ibkg); - return true; - } - - template - bool operator()(const T *const beam, - const T *const distance_mm, - const T *const detector_rot, - const T *const p0, - const T *const p1, - const T *const p2, - const T *const scale_factor, - const T *const B, - const T *const R, - T *residual) const { - T Ipred; - if (!Model(beam, distance_mm, detector_rot, p0, p1, p2, scale_factor, B, R, Ipred)) - return false; - - residual[0] = (Ipred - T(Iobs)) * T(weight); - return true; - } - - const double Itrue, Iobs, Ibkg, weight; - const double obs_x, obs_y; - const double inv_lambda; - const double pixel_size; - const double exp_h, exp_k, exp_l; - const double R_bw_sq; // bandwidth radial-width^2 contribution (0 = monochromatic) - const double pol; // per-reflection polarization correction - gemmi::CrystalSystem symmetry; -}; - -// --------------------------------------------------------------------------- -// Per-shoebox cost functor -// -// One residual block per reflection emitting N residuals (one per shoebox pixel). -// The expensive per-reflection geometry (PredictedNode: symmetry-aware B matrix, -// three rotations, cross products) is computed ONCE; only the cheap per-pixel -// ObservedRecip + Gaussian profile run in the pixel loop. This is identical in -// value to the old one-block-per-pixel formulation but ~(pixels-per-shoebox)x -// fewer evaluations of the costly node computation. Uses the same shared helpers -// (and hence the same conventions) as PixelResidual. -// --------------------------------------------------------------------------- -struct ShoeboxResidual { - ShoeboxResidual(const ReflGroup &g, double lambda, double pixel_size, - gemmi::CrystalSystem symmetry) - : pixels(g.pixels), Itrue(g.Itrue), R_bw_sq(g.R_bw_sq), pol(g.pol), - exp_h(g.h), exp_k(g.k), exp_l(g.l), - inv_lambda(1.0 / lambda), pixel_size(pixel_size), - symmetry(symmetry) {} - - template - bool operator()(const T *const *params, T *residual) const { - // Parameter blocks (order matches AddParameterBlock in Run): - // 0 beam[2] 1 distance[1] 2 detector_rot[2] - // 3 p0[3] 4 p1[3] 5 p2[3] 6 scale[1] 7 B[1] 8 R[2] - const T *beam = params[0]; - const T *distance_mm = params[1]; - const T *detector_rot = params[2]; - const T *p0 = params[3]; - const T *p1 = params[4]; - const T *p2 = params[5]; - const T *scale_factor = params[6]; - const T *B = params[7]; - const T *R = params[8]; - - if (R[0] < T(1e-10) || R[1] < T(1e-10)) - return false; - - // --- per-reflection: computed once --------------------------------- - Eigen::Matrix e_pred_recip, n_radial; - T q_sq; - if (!PredictedNode(p0, p1, p2, exp_h, exp_k, exp_l, symmetry, inv_lambda, - e_pred_recip, n_radial, q_sq)) - return false; - - const T B_term = ceres::exp(-B[0] * q_sq / T(4.0)); - const T R0_eff_sq = R[0] * R[0] + T(R_bw_sq); - - // --- per-pixel loop ------------------------------------------------- - for (size_t i = 0; i < pixels.size(); ++i) { - const PixelObs &obs = pixels[i]; - - Eigen::Matrix e_obs_recip; - ObservedRecip(beam, distance_mm, detector_rot, - obs.x, obs.y, pixel_size, inv_lambda, e_obs_recip); - - const Eigen::Matrix delta_q = e_obs_recip - e_pred_recip; - const T eps_radial = delta_q.dot(n_radial); - const T eps_tang_sq = (delta_q - eps_radial * n_radial).squaredNorm(); - - const T P_radial = ceres::exp(-eps_radial * eps_radial / R0_eff_sq); - const T P_tang = ceres::exp(-eps_tang_sq / (R[1] * R[1])) - / (T(M_PI) * R[1] * R[1]); - - const T signal = scale_factor[0] * T(Itrue) * B_term * P_radial * P_tang * T(pol); - const T Ipred = signal + T(obs.Ibkg); - residual[i] = (Ipred - T(obs.Iobs)) * T(obs.weight); - } - return true; - } - - std::vector pixels; - const double Itrue, R_bw_sq, pol; - const double exp_h, exp_k, exp_l; - const double inv_lambda, pixel_size; - gemmi::CrystalSystem symmetry; -}; - PixelRefine::PixelRefine(const DiffractionExperiment &experiment, const std::vector &reference) : xpixel(experiment.GetXPixelsNum()), @@ -598,187 +393,6 @@ void PixelRefine::BuildParameterBlocks(const PixelRefineData &data, } } -template -void PixelRefine::SweepOrientationCell(const T *image, BraggPrediction &prediction, - PixelRefineData &data) const { - const int radius = data.shoebox_radius; - const double beam_x = data.geom.GetBeamX_pxl(); - const double beam_y = data.geom.GetBeamY_pxl(); - const auto qnan = std::numeric_limits::quiet_NaN(); - - // Box-sum minus local (perimeter) background, raw counts. NaN if the box runs - // off the detector or hits a masked/saturated pixel. - auto integrate = [&](double px, double py) -> double { - const int cx = static_cast(std::lround(px)); - const int cy = static_cast(std::lround(py)); - const int outer = radius + 1; - if (cx - outer < 0 || cy - outer < 0 || - cx + outer >= static_cast(xpixel) || cy + outer >= static_cast(ypixel)) - return qnan; - double sig = 0.0; - int nsig = 0; - std::vector ring; - ring.reserve((2 * outer + 1) * (2 * outer + 1)); - for (int y = cy - outer; y <= cy + outer; ++y) { - for (int x = cx - outer; x <= cx + outer; ++x) { - const T raw = image[static_cast(xpixel) * y + x]; - if (raw == std::numeric_limits::max()) - return qnan; - if (std::is_signed_v && raw == std::numeric_limits::min()) - return qnan; - const double v = static_cast(raw); - if (std::abs(x - cx) <= radius && std::abs(y - cy) <= radius) { - sig += v; - ++nsig; - } else { - ring.push_back(v); - } - } - } - if (ring.size() < 5) - return qnan; - return sig - nsig * MedianInPlace(ring); - }; - - // Predict (wide band) and collect every reflection that has a reference value, - // with its detector radius. The full set is scored - the strong low-res spots - // anchor the CC, the weak high-res spots are what "appear" at the right cell. - DiffractionExperiment exp_iter = experiment; - exp_iter.BeamX_pxl(beam_x).BeamY_pxl(beam_y) - .DetectorDistance_mm(data.geom.GetDetectorDistance_mm()) - .PoniRot1_rad(data.geom.GetPoniRot1_rad()) - .PoniRot2_rad(data.geom.GetPoniRot2_rad()); - const BraggPredictionSettings settings{ - .high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(), - .ewald_dist_cutoff = static_cast(data.ewald_dist_cutoff), - .max_hkl = 100, - .centering = data.centering, - .bandwidth_sigma = static_cast(data.bandwidth) - }; - const int nrefl = prediction.Calc(exp_iter, data.latt, settings); - const auto &predicted = prediction.GetReflections(); - - struct Matched { int h, k, l; double refI; }; - std::vector matched; - double r_max = 0.0, r_min = std::numeric_limits::max(); - for (int i = 0; i < nrefl; ++i) { - const auto &r = predicted[i]; - const auto it = reference_data.find(hkl_key_generator(r)); - if (it == reference_data.end()) - continue; - matched.push_back({r.h, r.k, r.l, it->second}); - const double dx = r.predicted_x - beam_x; - const double dy = r.predicted_y - beam_y; - const double rad = std::sqrt(dx * dx + dy * dy); - r_max = std::max(r_max, rad); - r_min = std::min(r_min, rad); - } - if (matched.size() < 20 || r_min <= 1.0 || r_max <= r_min) - return; // too little to anchor a meaningful sweep - - // CC of the box-summed intensities against the reference, over all matched hkls. - auto score = [&](const CrystalLattice &L) -> double { - const Coord A = L.Astar(), B = L.Bstar(), C = L.Cstar(); - double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0; - int n = 0; - for (const auto &m : matched) { - const Coord g = A * static_cast(m.h) + B * static_cast(m.k) - + C * static_cast(m.l); - const auto [x, y] = data.geom.RecipToDetector(g); - if (!std::isfinite(x) || !std::isfinite(y)) - continue; - const double I = integrate(x, y); - if (!std::isfinite(I)) - continue; - const double yv = m.refI; - sx += I; sy += yv; sxx += I * I; syy += yv * yv; sxy += I * yv; ++n; - } - if (n < 10) - return -2.0; - const double nd = n; - const double cov = sxy - sx * sy / nd; - const double vx = sxx - sx * sx / nd; - const double vy = syy - sy * sy / nd; - if (!(vx > 0.0 && vy > 0.0)) - return -2.0; - return cov / std::sqrt(vx * vy); - }; - - // Step = 1 px at the highest resolution. Range = assumed orientation/cell-scale - // uncertainty (a few px at high res), NOT the low-res 2 px cap: the latter is - // ~2*r_max/r_min px at high res - far too permissive, and lets the per-image CC - // overfit. Here the low-res spots barely move (stay anchored). - const double step = 1.0 / r_max; - const int n_rot = std::clamp( - static_cast(std::lround(data.sweep_max_deg * M_PI / 180.0 * r_max)), 1, 25); - const int n_scale = std::clamp( - static_cast(std::lround(data.sweep_max_cell_frac * r_max)), 1, 25); - const Coord axes[3] = {Coord(1, 0, 0), Coord(0, 1, 0), Coord(0, 0, 1)}; - - CrystalLattice best = data.latt; - double best_cc = score(best); - - for (int round = 0; round < 2; ++round) { - for (const auto &axis : axes) { - CrystalLattice axis_best = best; - double axis_cc = best_cc; - for (int i = -n_rot; i <= n_rot; ++i) { - if (i == 0) - continue; - CrystalLattice cand = best.Multiply(RotMatrix(static_cast(i * step), axis)); - const double cc = score(cand); - if (cc > axis_cc) { - axis_cc = cc; - axis_best = cand; - } - } - best = axis_best; - best_cc = axis_cc; - } - CrystalLattice scale_best = best; - double scale_cc = best_cc; - for (int i = -n_scale; i <= n_scale; ++i) { - if (i == 0) - continue; - const double s = 1.0 / (1.0 + i * step); // cell scale (1+eps) -> recip * 1/(1+eps) - CrystalLattice cand = best.Multiply(gemmi::Mat33(s, 0, 0, 0, s, 0, 0, 0, s)); - const double cc = score(cand); - if (cc > scale_cc) { - scale_cc = cc; - scale_best = cand; - } - } - best = scale_best; - best_cc = scale_cc; - } - - data.latt = best; -} - -// --------------------------------------------------------------------------- -// Term 1 of the factored likelihood (FACTORED_MODEL.md): the per-reflection -// *intensity* (0th-moment) residual. The profile-fit amplitude J should equal the -// scaled reference J_model = G * exp(-B/4d^2) * partiality * pol * I_ref. One scalar -// residual per reflection, weighted by the model-expected (Fisher) sigma_J. This is -// the scaling residual - integration and scaling become one objective, and the empty -// pixels (which make no residual of their own) stop dominating the fit. With geometry -// and R held fixed, J, partiality and sigma_J are constants, so only G and B are free. -// --------------------------------------------------------------------------- -struct IntensityResidual { - IntensityResidual(double J, double sigma_J, double partiality, double pol, - double I_ref, double inv_4d2) - : J(J), inv_sigma(1.0 / sigma_J), partiality(partiality), pol(pol), - I_ref(I_ref), inv_4d2(inv_4d2) {} - template - bool operator()(const T *const G, const T *const B, T *residual) const { - const T B_term = ceres::exp(-B[0] * T(inv_4d2)); - const T J_model = G[0] * B_term * T(partiality) * T(pol) * T(I_ref); - residual[0] = (J_model - T(J)) * T(inv_sigma); - return true; - } - double J, inv_sigma, partiality, pol, I_ref, inv_4d2; -}; - template void PixelRefine::Run(const T *image, BraggPrediction &prediction, @@ -786,11 +400,6 @@ void PixelRefine::Run(const T *image, data.solved = false; data.reflections.clear(); - // Global orientation + cell-scale sweep before the local LSQ, to recentre the - // high-resolution shoeboxes onto signal that small misalignments hide. - if (data.sweep_orientation) - SweepOrientationCell(image, prediction, data); - const double lambda = data.geom.GetWavelength_A(); const double pixel_size = data.geom.GetPixelSize_mm(); @@ -826,701 +435,244 @@ void PixelRefine::Run(const T *image, return bl * bl / (2.0 * d * d * d * d); }; - // Mutable experiment whose geometry is re-synced from the refined data.geom - // before each prediction, so shoeboxes track the refined geometry/cell. + // Geometry is FIXED here: orientation/cell/detector were already refined upstream + // by XtalOptimizer (IndexAndRefine::RefineGeometryIfNeeded). PixelRefine is an + // intensity-only operation - it predicts shoeboxes with this geometry, measures the + // tangential profile width, and fits the per-image scale G (and B) to the reference. + double beam[2], dist_mm, detector_rot[2]; + double latt_vec0[3], latt_vec1[3], latt_vec2[3]; + BuildParameterBlocks(data, beam, dist_mm, detector_rot, latt_vec0, latt_vec1, latt_vec2); + + // ---- 1. Predict shoeboxes for the current geometry ------------------------ DiffractionExperiment exp_iter = experiment; + exp_iter.BeamX_pxl(data.geom.GetBeamX_pxl()) + .BeamY_pxl(data.geom.GetBeamY_pxl()) + .DetectorDistance_mm(data.geom.GetDetectorDistance_mm()) + .PoniRot1_rad(data.geom.GetPoniRot1_rad()) + .PoniRot2_rad(data.geom.GetPoniRot2_rad()); + const int nrefl = prediction.Calc(exp_iter, data.latt, settings_prediction); + + // ---- 2. Collect per-reflection shoebox pixels + local background ---------- + // GetReflections() returns the full pre-sized buffer; only the first nrefl + // entries are valid for this image. A spot-core mask over ALL predictions keeps + // each reflection's background ring from picking up a neighbour's signal. + const auto &predicted = prediction.GetReflections(); + const auto spot_mask = BuildSpotMask(predicted, nrefl, xpixel, ypixel, radius); - // State retained after the loop for the final reflection extraction. std::vector groups; - double beam[2] = {0, 0}; - double dist_mm = data.geom.GetDetectorDistance_mm(); - double detector_rot[2] = {0, 0}; - double latt_vec0[3] = {0, 0, 0}; // orientation (Rodrigues) - double latt_vec1[3] = {0, 0, 0}; // lengths - double latt_vec2[3] = {0, 0, 0}; // angles (rad) - double orient_prior[3] = {0, 0, 0}; // pre-refinement orientation (regularization anchor) + for (int ri = 0; ri < nrefl; ++ri) { + const auto &refl = predicted[ri]; + const auto hkl = hkl_key_generator(refl); + if (!reference_data.contains(hkl)) + continue; - const bool eval_only = (data.max_iterations <= 0); - const int n_iter = std::max(1, data.max_iterations); - for (int iter = 0; iter < n_iter; ++iter) { - // ---- 1. Re-sync prediction geometry from the (refined) model ---------- - exp_iter.BeamX_pxl(data.geom.GetBeamX_pxl()) - .BeamY_pxl(data.geom.GetBeamY_pxl()) - .DetectorDistance_mm(data.geom.GetDetectorDistance_mm()) - .PoniRot1_rad(data.geom.GetPoniRot1_rad()) - .PoniRot2_rad(data.geom.GetPoniRot2_rad()); + // Local flat background from the ring around the shoebox (raw counts). If we + // cannot estimate a clean local background the reflection is dropped, exactly + // as BraggIntegrate2D marks it unobserved when too few background pixels survive. + double Ibkg = 0.0; + if (!EstimateLocalBackground(image, spot_mask, xpixel, ypixel, + refl.predicted_x, refl.predicted_y, + radius, bkg_outer_radius, Ibkg)) + continue; - const int nrefl = prediction.Calc(exp_iter, data.latt, settings_prediction); + ReflGroup g; + g.h = refl.h; + g.k = refl.k; + g.l = refl.l; + g.d = refl.d; + g.Itrue = reference_data[hkl]; + g.R_bw_sq = bandwidth_radial_sq(refl.d); + g.pol = polarization(refl.predicted_x, refl.predicted_y); + g.Ibkg = Ibkg; + g.predicted_x = refl.predicted_x; + g.predicted_y = refl.predicted_y; - // ---- 2. Collect per-reflection shoebox pixels ------------------------- - // GetReflections() returns the full pre-sized buffer; only the first - // nrefl entries are valid for this image (the rest are stale/zeroed). - groups.clear(); - const auto &predicted = prediction.GetReflections(); - - // Spot-core mask over ALL predicted reflections, so each reflection's - // local background ignores pixels that belong to a neighbouring spot. - const auto spot_mask = BuildSpotMask(predicted, nrefl, xpixel, ypixel, radius); - - for (int ri = 0; ri < nrefl; ++ri) { - const auto &refl = predicted[ri]; - const auto hkl = hkl_key_generator(refl); - if (!reference_data.contains(hkl)) - continue; - - // Local flat background from the ring around the shoebox (raw counts). - // No azimuthal fallback: if we cannot estimate a clean local background - // the reflection is dropped, exactly as BraggIntegrate2D marks it - // unobserved when fewer than a handful of background pixels survive. - double Ibkg = 0.0; - if (!EstimateLocalBackground(image, spot_mask, xpixel, ypixel, - refl.predicted_x, refl.predicted_y, - radius, bkg_outer_radius, Ibkg)) - continue; - - ReflGroup g; - g.h = refl.h; - g.k = refl.k; - g.l = refl.l; - g.d = refl.d; - g.Itrue = reference_data[hkl]; - g.R_bw_sq = bandwidth_radial_sq(refl.d); - g.pol = polarization(refl.predicted_x, refl.predicted_y); - g.Ibkg = Ibkg; - g.predicted_x = refl.predicted_x; - g.predicted_y = refl.predicted_y; - - const auto box = ShoeboxBounds(refl.predicted_x, refl.predicted_y, radius, xpixel, ypixel); - - for (int y = box.min_y; y <= box.max_y; ++y) { - for (int x = box.min_x; x <= box.max_x; ++x) { - const size_t npixel = xpixel * y + x; - - // Skip sentinel (masked / saturated) pixels. We assume the pixel - // mask is already applied upstream (encoded as the sentinel). - if (image[npixel] == std::numeric_limits::max()) - continue; - if (std::is_signed_v && (image[npixel] == std::numeric_limits::min())) - continue; - - const double Iobs = static_cast(image[npixel]); // raw counts - - // Variance for the fit weight. Weighting by the observed count - // (var = Iobs) lets down-fluctuated background pixels carry the - // largest 1/sqrt(var) weight, which biases the fit towards "no - // signal" and drove the per-image scale G to 0 on weak images - // (collapsing the merge). Use the local background as the - // (background-limited) variance, constant over the shoebox - the - // same de-biasing applied to the extraction. - double var = std::max(Ibkg, 1.0); - double weight = 1.0 / std::sqrt(var); - - // Signal-weighting: down-weight pixels far from the predicted spot - // centre so the empty shoebox corners cannot dilute or destabilise - // the fit; the signal-bearing core drives the refined parameters. - if (data.fit_signal_sigma_pix > 0.0) { - const double dx = x - g.predicted_x; - const double dy = y - g.predicted_y; - const double s2 = data.fit_signal_sigma_pix * data.fit_signal_sigma_pix; - weight *= std::exp(-0.5 * (dx * dx + dy * dy) / s2); - } - - PixelObs obs{ - .x = static_cast(x), - .y = static_cast(y), - .Iobs = Iobs, - .Ibkg = Ibkg, - .weight = weight - }; - g.pixels.push_back(obs); - } - } - - if (!g.pixels.empty()) - groups.push_back(std::move(g)); - } - - if (groups.empty()) - return; - - // ---- 3. Set up parameter blocks (geometry part mirrors XtalOptimizer) - - BuildParameterBlocks(data, beam, dist_mm, detector_rot, - latt_vec0, latt_vec1, latt_vec2); - - // Anchor for orientation regularization = the orientation the LSQ starts from - // (captured before the predict<->refine iterations move it). When the global - // sweep ran first this is the swept orientation, not the original spot-centroid - // one - which is intended: the regularizer keeps the LSQ near its own starting - // point, it is not meant to pull a deliberate sweep back. - if (iter == 0) - for (int i = 0; i < 3; ++i) - orient_prior[i] = latt_vec0[i]; - - // ---- Term 3: per-reflection recentre on the observed centroid ---------------- - // The geometry predicts the spot to ~0.4 px (per-reflection scatter a global fit - // cannot remove); a tight Term-2 template centred on the prediction then sits off - // the real spot. For confident spots, shift the profile centre to the observed - // centroid (used consistently by Term 2, Term 1 and the extraction below). Weak - // spots keep the prediction (recentring on a noise centroid would bias them). - for (auto &g : groups) { - g.dcx = 0.0; - g.dcy = 0.0; - } - if (data.recenter_profile && !groups.empty()) { - const double box_px = (2.0 * radius + 1.0) * (2.0 * radius + 1.0); - for (auto &g : groups) { - double sw = 0.0, swx = 0.0, swy = 0.0; - for (const auto &px : g.pixels) { - const double w = std::max(px.Iobs - g.Ibkg, 0.0); - sw += w; swx += w * px.x; swy += w * px.y; - } - if (sw <= 0.0) + const auto box = ShoeboxBounds(refl.predicted_x, refl.predicted_y, radius, xpixel, ypixel); + for (int y = box.min_y; y <= box.max_y; ++y) { + for (int x = box.min_x; x <= box.max_x; ++x) { + const size_t npixel = xpixel * y + x; + // Skip sentinel (masked / saturated) pixels. + if (image[npixel] == std::numeric_limits::max()) continue; - if (sw / std::sqrt(std::max(box_px * g.Ibkg, 1.0)) < data.recenter_min_signif) + if (std::is_signed_v && (image[npixel] == std::numeric_limits::min())) continue; - double dcx = swx / sw - g.predicted_x, dcy = swy / sw - g.predicted_y; - const double dl = std::sqrt(dcx * dcx + dcy * dcy); - if (dl > 2.0) { dcx *= 2.0 / dl; dcy *= 2.0 / dl; } // clamp runaway centroids - g.dcx = dcx; - g.dcy = dcy; + g.pixels.push_back({static_cast(x), static_cast(y), + static_cast(image[npixel]), Ibkg}); } } + if (!g.pixels.empty()) + groups.push_back(std::move(g)); + } + if (groups.empty()) + return; - // ---- Term 2: per-resolution tangential profile width R1 from spot shapes ------ - // Default: every reflection uses the global R1; with shape_R1 on, override it with - // R1 = sqrt(2*) from the intensity-weighted second moment of the strong - // spots, binned by resolution (low res small spots, high res larger). A shape - // statistic - normalised by the total, so decoupled from the per-image scale. - for (auto &g : groups) - g.R1_eff = data.R[1]; - if (data.shape_R1 && !groups.empty()) { - constexpr int n_bins = 6; - const double box_px = (2.0 * radius + 1.0) * (2.0 * radius + 1.0); - double s2min = 1e30, s2max = 0.0; - for (const auto &g : groups) { - const double s2 = 1.0 / (g.d * g.d); - s2min = std::min(s2min, s2); - s2max = std::max(s2max, s2); - } - const double span = std::max(s2max - s2min, 1e-12); - auto bin_of = [&](double d) { - return std::clamp(static_cast((1.0 / (d * d) - s2min) / span * n_bins), 0, n_bins - 1); - }; - std::vector> bin_M2(n_bins); - for (const auto &g : groups) { - double sw = 0.0, sw_et2 = 0.0; - for (const auto &px : g.pixels) { - PixelObs probe{px.x - g.dcx, px.y - g.dcy, 0.0, g.Ibkg, 1.0}; - PixelResidual pr(probe, 1.0, lambda, pixel_size, g.h, g.k, g.l, - g.R_bw_sq, g.pol, data.crystal_system); - double q_sq, eps_r, eps_t_sq; - if (!pr.GeometryTerms(beam, &dist_mm, detector_rot, - latt_vec0, latt_vec1, latt_vec2, q_sq, eps_r, eps_t_sq)) - continue; - const double w = std::max(px.Iobs - g.Ibkg, 0.0); - sw += w; - sw_et2 += w * eps_t_sq; - } - if (sw <= 0.0) - continue; - const double signif = sw / std::sqrt(std::max(box_px * g.Ibkg, 1.0)); - if (signif >= 5.0) // only well-measured spots define the shape - bin_M2[bin_of(g.d)].push_back(sw_et2 / sw); - } - std::vector bin_R1(n_bins, data.R[1]); - for (int b = 0; b < n_bins; ++b) - if (bin_M2[b].size() >= 5) { - const double r1 = std::sqrt(2.0 * std::max(MedianInPlace(bin_M2[b]), 0.0)); - if (std::isfinite(r1) && r1 > 1e-4) - bin_R1[b] = std::clamp(r1, 1e-4, 0.05); - } - for (auto &g : groups) - g.R1_eff = bin_R1[bin_of(g.d)]; - data.shape_R1_lores = bin_R1[0]; // lowest resolution bin - data.shape_R1_hires = bin_R1[n_bins - 1]; // highest resolution bin - } - - // ---- 4. Build the problem --------------------------------------------- - // One residual block per shoebox (N residuals), so the expensive - // per-reflection node geometry is evaluated once per reflection instead - // of once per pixel. - ceres::Problem problem; - size_t residual_pixels = 0; - if (data.intensity_residual) { - // Term-1 path: one per-reflection intensity residual. Geometry & R fixed, so - // J / partiality / sigma_J are computed here as constants and only G, B vary. - const double R0 = data.R[0]; - for (const auto &g : groups) { - const double R1 = g.R1_eff; // Term 2: per-resolution profile width - double num = 0.0, den = 0.0, rad = 0.0; - std::vector> pt_sig; // (P_t, Iobs-Bg) for Fisher pass - pt_sig.reserve(g.pixels.size()); - for (const auto &px : g.pixels) { - PixelObs probe{px.x - g.dcx, px.y - g.dcy, 0.0, g.Ibkg, 1.0}; // Term 3 recentre - PixelResidual pr(probe, 1.0, lambda, pixel_size, g.h, g.k, g.l, - g.R_bw_sq, g.pol, data.crystal_system); - double q_sq, eps_r, eps_t_sq; - if (!pr.GeometryTerms(beam, &dist_mm, detector_rot, - latt_vec0, latt_vec1, latt_vec2, q_sq, eps_r, eps_t_sq)) - continue; - if (!(R1 > 0.0) || !(R0 > 0.0)) - continue; - const double P_t = std::exp(-eps_t_sq / (R1 * R1)) / (M_PI * R1 * R1); - const double R0_eff_sq = R0 * R0 + g.R_bw_sq; - const double P_rad = std::exp(-eps_r * eps_r / R0_eff_sq); - const double v = std::max(g.Ibkg, 1.0); - const double sig = px.Iobs - g.Ibkg; - num += P_t * sig / v; - den += P_t * P_t / v; - rad += P_rad * P_t * P_t / v; - pt_sig.emplace_back(P_t, sig); - } - if (!(den > 0.0)) - continue; - const double J = num / den; - const double partiality = rad / den; - // Model-expected (Fisher) variance: v_p = background + expected signal J*P_t, - // not the per-pixel observed counts (which down-bias) - so the weight tracks - // information, and an expected-strong reflection that is absent hurts. - double den_f = 0.0; - for (const auto &[P_t, sig] : pt_sig) { - const double v_f = std::max(g.Ibkg + std::max(J, 0.0) * P_t, 1.0); - den_f += P_t * P_t / v_f; - } - const double sigma_J = std::sqrt(1.0 / std::max(den_f, 1e-30)); - const double inv_4d2 = (g.d > 0.0) ? 1.0 / (4.0 * g.d * g.d) : 0.0; - auto *cost = new ceres::AutoDiffCostFunction( - new IntensityResidual(J, sigma_J, partiality, g.pol, g.Itrue, inv_4d2)); - problem.AddResidualBlock(cost, nullptr, &data.scale_factor, &data.B_factor); - ++residual_pixels; - } - data.residual_count = residual_pixels; - } else { + // ---- 3. Term 2: per-resolution tangential profile width R1 ---------------- + // R1 = sqrt(2*) from the intensity-weighted tangential second moment of + // the strong spots, binned by resolution (low res small spots, high res larger). + // A *shape* statistic, normalised by the total intensity, so it is decoupled from + // the per-image scale - which is what makes measuring it (rather than fitting it, + // where it is degenerate with G) stable. Weak spots fall back to the global R[1]. + for (auto &g : groups) + g.R1_eff = data.R[1]; + { + constexpr int n_bins = 6; + const double box_px = (2.0 * radius + 1.0) * (2.0 * radius + 1.0); + double s2min = 1e30, s2max = 0.0; for (const auto &g : groups) { - auto *cost = new ceres::DynamicAutoDiffCostFunction( - new ShoeboxResidual(g, lambda, pixel_size, data.crystal_system)); - cost->AddParameterBlock(2); // beam - cost->AddParameterBlock(1); // distance - cost->AddParameterBlock(2); // detector_rot - cost->AddParameterBlock(3); // p0 (orientation) - cost->AddParameterBlock(3); // p1 (lengths) - cost->AddParameterBlock(3); // p2 (angles) - cost->AddParameterBlock(1); // scale G - cost->AddParameterBlock(1); // B - cost->AddParameterBlock(2); // R - cost->SetNumResiduals(static_cast(g.pixels.size())); - // No robust loss here: a per-block (whole-shoebox) Huber would act on - // the sum of ~N squared residuals and mis-scale, unlike the previous - // per-pixel Huber. Per-pixel sigma weighting is retained; per-pixel - // outlier rejection (zingers) is a TODO if needed. - problem.AddResidualBlock(cost, nullptr, - beam, &dist_mm, detector_rot, - latt_vec0, latt_vec1, latt_vec2, - &data.scale_factor, &data.B_factor, data.R); - residual_pixels += g.pixels.size(); + const double s2 = 1.0 / (g.d * g.d); + s2min = std::min(s2min, s2); + s2max = std::max(s2max, s2); } - data.residual_count = residual_pixels; - } - - // ---- 5. Constrain / bound parameter blocks ---------------------------- - if (data.intensity_residual) { - // Only G and B are in this problem; geometry/R are not parameters here. - problem.SetParameterLowerBound(&data.scale_factor, 0, 0.0); - if (!data.refine_B) - problem.SetParameterBlockConstant(&data.B_factor); - // Regularize G->1, weight sqrt(n_refl/sigma): commensurate because the data - // term is now one residual per reflection (unlike the per-pixel path). - if (data.scale_reg_sigma > 0.0 && !groups.empty()) { - const double w = std::sqrt(static_cast(groups.size()) / data.scale_reg_sigma); - auto *reg = new ceres::AutoDiffCostFunction( - new ScalarRegularizer(w, 1.0)); - problem.AddResidualBlock(reg, nullptr, &data.scale_factor); - } - } else { - if (!data.refine_orientation) { - problem.SetParameterBlockConstant(latt_vec0); - } else if (data.orient_reg_sigma_deg > 0.0) { - // Anchor orientation to its spot-centroid prior. The weight is scaled to - // the *pixel* data term (sqrt(n_pixels)/sigma_rad), not the reflection - // count - the data has one residual per shoebox pixel, so a reflection- - // scaled prior (~50x too weak) was simply not felt. At a misorientation of - // orient_reg_sigma_deg the prior matches the data, so the fit only moves - // further when the pixels strongly agree it should. - const double sigma_rad = std::max(data.orient_reg_sigma_deg * M_PI / 180.0, 1e-9); - const double w = std::sqrt(static_cast(residual_pixels)) / sigma_rad; - auto *reg = new ceres::AutoDiffCostFunction( - new OrientationRegularizer(w, orient_prior)); - problem.AddResidualBlock(reg, nullptr, latt_vec0); - } - - if (!data.refine_unit_cell) { - problem.SetParameterBlockConstant(latt_vec1); - problem.SetParameterBlockConstant(latt_vec2); - } else { - for (int i = 0; i < 3; ++i) { - problem.SetParameterLowerBound(latt_vec1, i, 5.0); - problem.SetParameterUpperBound(latt_vec1, i, 1000.0); - } - if (data.crystal_system != gemmi::CrystalSystem::Monoclinic && - data.crystal_system != gemmi::CrystalSystem::Triclinic) - problem.SetParameterBlockConstant(latt_vec2); - } - - if (!data.refine_beam_center) - problem.SetParameterBlockConstant(beam); - - if (!data.refine_distance) { - problem.SetParameterBlockConstant(&dist_mm); - } else { - problem.SetParameterLowerBound(&dist_mm, 0, dist_mm * 0.9); - problem.SetParameterUpperBound(&dist_mm, 0, dist_mm * 1.1); - } - - if (!data.refine_detector_angles) { - problem.SetParameterBlockConstant(detector_rot); - } else { - const double rng = 3.0 / 180.0 * M_PI; - for (int i = 0; i < 2; ++i) { - problem.SetParameterLowerBound(detector_rot, i, detector_rot[i] - rng); - problem.SetParameterUpperBound(detector_rot, i, detector_rot[i] + rng); - } - } - - if (data.refine_scale) { - problem.SetParameterLowerBound(&data.scale_factor, 0, 0.0); - // Regularize G towards 1 so weakly-constrained images cannot wander - // (an unconstrained 1/G is what collapsed the cross-image merge). Weight - // scaled to the pixel data term (n_pixels), not the reflection count. - if (data.scale_reg_sigma > 0.0) { - const double w = std::sqrt(static_cast(residual_pixels) / data.scale_reg_sigma); - auto *reg = new ceres::AutoDiffCostFunction( - new ScalarRegularizer(w, 1.0)); - problem.AddResidualBlock(reg, nullptr, &data.scale_factor); - } - } else { - problem.SetParameterBlockConstant(&data.scale_factor); - } - - if (!data.refine_B) - problem.SetParameterBlockConstant(&data.B_factor); - - if (data.refine_R) { - if (data.fix_R0) { - // Diagnostic: hold R0 constant, refine R1 only. - problem.SetManifold(data.R, new ceres::SubsetManifold(2, {0})); - } else { - problem.SetParameterLowerBound(data.R, 0, 1e-5); - problem.SetParameterLowerBound(data.R, 1, 1e-5); - } - } else { - problem.SetParameterBlockConstant(data.R); - } - } // end per-pixel (non-intensity_residual) constraints - - // ---- 6. Solve (or, for max_iterations<=0, just evaluate the cost) ----- - // Evaluate-only is the live-residual path: it reports the current cost - // without moving any parameter, so a UI can show how good the present - // R0/R1/bandwidth/geometry are as the user drags sliders. - if (eval_only) { - double cost = 0.0; - problem.Evaluate(ceres::Problem::EvaluateOptions(), &cost, nullptr, nullptr, nullptr); - data.final_cost = cost; - data.solved = true; - } else { - ceres::Solver::Options options; - options.linear_solver_type = ceres::DENSE_QR; - options.minimizer_progress_to_stdout = false; - options.logging_type = ceres::LoggingType::SILENT; - options.max_solver_time_in_seconds = data.max_time_s; - options.num_threads = 1; - - ceres::Solver::Summary summary; - ceres::Solve(options, &problem, &summary); - - data.final_cost = summary.final_cost; - data.solved = summary.IsSolutionUsable(); - - // Diagnostic: Pearson correlations on the final solve. Always needs G and B - // refined for G-B; the R correlations are added only when R is also refined. - if (data.compute_covariance && data.solved && iter == n_iter - 1 && - data.refine_scale && data.refine_B) { - ceres::Covariance::Options copt; - copt.algorithm_type = ceres::DENSE_SVD; - copt.null_space_rank = -1; // tolerate (and reveal) degenerate directions - ceres::Covariance cov(copt); - std::vector> blocks = { - {&data.scale_factor, &data.scale_factor}, {&data.B_factor, &data.B_factor}, - {&data.scale_factor, &data.B_factor}}; - if (data.refine_R) { - blocks.push_back({data.R, data.R}); - blocks.push_back({&data.scale_factor, data.R}); - blocks.push_back({&data.B_factor, data.R}); - } - if (cov.Compute(blocks, &problem)) { - double cGG, cBB, cGB; - cov.GetCovarianceBlock(&data.scale_factor, &data.scale_factor, &cGG); - cov.GetCovarianceBlock(&data.B_factor, &data.B_factor, &cBB); - cov.GetCovarianceBlock(&data.scale_factor, &data.B_factor, &cGB); - const double sG = std::sqrt(std::max(cGG, 0.0)); - const double sB = std::sqrt(std::max(cBB, 0.0)); - auto rho = [](double c, double s1, double s2) { - return (s1 > 0.0 && s2 > 0.0) ? c / (s1 * s2) : NAN; - }; - data.corr_GB = rho(cGB, sG, sB); - if (data.refine_R) { - double cRR[4], cGR[2], cBR[2]; - cov.GetCovarianceBlock(data.R, data.R, cRR); - cov.GetCovarianceBlock(&data.scale_factor, data.R, cGR); - cov.GetCovarianceBlock(&data.B_factor, data.R, cBR); - const double sR0 = std::sqrt(std::max(cRR[0], 0.0)); - const double sR1 = std::sqrt(std::max(cRR[3], 0.0)); - data.corr_GR0 = rho(cGR[0], sG, sR0); - data.corr_GR1 = rho(cGR[1], sG, sR1); - data.corr_BR0 = rho(cBR[0], sB, sR0); - data.corr_BR1 = rho(cBR[1], sB, sR1); - data.corr_R0R1 = rho(cRR[1], sR0, sR1); - } - data.covariance_valid = true; - } - } - } - - // ---- 7. Write refined geometry + lattice back into data --------------- - if (data.refine_beam_center) - data.geom.BeamX_pxl(beam[0]).BeamY_pxl(beam[1]); - if (data.refine_distance) - data.geom.DetectorDistance_mm(dist_mm); - if (data.refine_detector_angles) - data.geom.PoniRot1_rad(detector_rot[0]).PoniRot2_rad(detector_rot[1]); - - if (data.refine_orientation || data.refine_unit_cell) { - switch (data.crystal_system) { - case gemmi::CrystalSystem::Orthorhombic: - data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI/2, M_PI/2, M_PI/2); - break; - case gemmi::CrystalSystem::Tetragonal: - latt_vec1[1] = latt_vec1[0]; - data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI/2, M_PI/2, M_PI/2); - break; - case gemmi::CrystalSystem::Cubic: - latt_vec1[1] = latt_vec1[0]; - latt_vec1[2] = latt_vec1[0]; - data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI/2, M_PI/2, M_PI/2); - break; - case gemmi::CrystalSystem::Hexagonal: - latt_vec1[1] = latt_vec1[0]; - data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI/2, M_PI/2, 2.0*M_PI/3.0); - break; - case gemmi::CrystalSystem::Monoclinic: - data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI/2, latt_vec2[0], M_PI/2); - break; - default: - data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, - latt_vec2[0], latt_vec2[1], latt_vec2[2]); - break; - } - } - } // predict<->refine iterations - - // ---- Adaptive integration mask -------------------------------------------- - // Measure R1 (tangential profile width) from the intensity-weighted tangential - // second moment of the strong spots, rather than fitting it. R1 is a *shape* - // statistic: sigma_t^2 = sum_p (I_p-B) eps_t,p^2 / sum_p (I_p-B), normalised by the - // total so it is independent of the per-image scale - which is exactly what breaks - // the R1<->G degeneracy (a measured width cannot be traded against G). One value per - // image here (from the strong, mostly low-res spots); a per-resolution version is the - // natural next step for the high-res / DMM-streak shapes. - if (data.adaptive_R1 && !groups.empty()) { - std::vector itrue; - itrue.reserve(groups.size()); - for (const auto &g : groups) - itrue.push_back(g.Itrue); - const size_t cut_idx = itrue.size() * 7 / 10; // keep the strongest ~30% - std::nth_element(itrue.begin(), itrue.begin() + cut_idx, itrue.end()); - const double itrue_cut = itrue[cut_idx]; - - std::vector sigma_t2; + const double span = std::max(s2max - s2min, 1e-12); + auto bin_of = [&](double d) { + return std::clamp(static_cast((1.0 / (d * d) - s2min) / span * n_bins), 0, n_bins - 1); + }; + std::vector> bin_M2(n_bins); for (const auto &g : groups) { - if (g.Itrue < itrue_cut) - continue; double sw = 0.0, sw_et2 = 0.0; for (const auto &px : g.pixels) { - PixelObs probe{px.x, px.y, 0.0, g.Ibkg, 1.0}; - PixelResidual pr(probe, 1.0, lambda, pixel_size, g.h, g.k, g.l, - g.R_bw_sq, g.pol, data.crystal_system); double q_sq, eps_r, eps_t_sq; - if (!pr.GeometryTerms(beam, &dist_mm, detector_rot, - latt_vec0, latt_vec1, latt_vec2, q_sq, eps_r, eps_t_sq)) + if (!GeometryProbe(px.x, px.y, lambda, pixel_size, g.h, g.k, g.l, data.crystal_system, + beam, dist_mm, detector_rot, latt_vec0, latt_vec1, latt_vec2, + q_sq, eps_r, eps_t_sq)) continue; const double w = std::max(px.Iobs - g.Ibkg, 0.0); sw += w; sw_et2 += w * eps_t_sq; } - if (sw > 0.0) - sigma_t2.push_back(sw_et2 / sw); - } - if (sigma_t2.size() >= 5) { - const double r1 = std::sqrt(2.0 * MedianInPlace(sigma_t2)); // R1^2 = 2 sigma_t^2 - if (std::isfinite(r1) && r1 > 1e-5) - data.R[1] = r1; - } - } - - // ---- Centering diagnostic -------------------------------------------------- - // Observed-centroid vs predicted-position offset for the strong spots, after all - // refinement. Large rms (relative to the spot size) means a tight profile mask - // sits off the real spot - which is why a generous box can beat profile fitting. - if (data.measure_centroid && !groups.empty()) { - const double beam_x = data.geom.GetBeamX_pxl(); - const double beam_y = data.geom.GetBeamY_pxl(); - const double box_px = (2.0 * radius + 1.0) * (2.0 * radius + 1.0); - // Raw pixel value (sentinel/bounds safe) for the parabolic peak fit. A constant - // background cancels in the parabola, so no need to subtract it here. - auto pix_val = [&](int x, int y) -> double { - if (x < 0 || x >= static_cast(xpixel) || y < 0 || y >= static_cast(ypixel)) - return std::numeric_limits::quiet_NaN(); - const size_t np = static_cast(xpixel) * y + x; - if (image[np] == std::numeric_limits::max()) - return std::numeric_limits::quiet_NaN(); - if (std::is_signed_v && image[np] == std::numeric_limits::min()) - return std::numeric_limits::quiet_NaN(); - return static_cast(image[np]); - }; - struct Off { double signif, tang_c, tang_p, rad_c; }; - std::vector offs; - double sdx = 0.0, sdy = 0.0, sd2 = 0.0; - size_t nc = 0; - for (const auto &g : groups) { - double sw = 0.0, swx = 0.0, swy = 0.0, bmax = -1e30; - int bx = 0, by = 0; - for (const auto &px : g.pixels) { - const double s = px.Iobs - g.Ibkg; - const double w = std::max(s, 0.0); - sw += w; swx += w * px.x; swy += w * px.y; - if (s > bmax) { bmax = s; bx = static_cast(std::lround(px.x)); by = static_cast(std::lround(px.y)); } - } if (sw <= 0.0) continue; const double signif = sw / std::sqrt(std::max(box_px * g.Ibkg, 1.0)); - if (signif < 5.0) - continue; // a measurable spot at any resolution (not just the strong low-res ones) - // Sub-pixel peak (mode): parabola through the brightest pixel and its two - // neighbours per axis. The mode tracks the prediction even when an asymmetric - // tail drags the centroid (mean) sideways - so peak vs centroid separates a - // shape asymmetry from a true position error. - double peak_x = bx, peak_y = by; - { const double l = pix_val(bx - 1, by), c = pix_val(bx, by), r = pix_val(bx + 1, by); - const double den = l - 2.0 * c + r; - if (std::isfinite(den) && den < -1e-9) - peak_x = bx + std::clamp(0.5 * (l - r) / den, -1.0, 1.0); } - { const double l = pix_val(bx, by - 1), c = pix_val(bx, by), r = pix_val(bx, by + 1); - const double den = l - 2.0 * c + r; - if (std::isfinite(den) && den < -1e-9) - peak_y = by + std::clamp(0.5 * (l - r) / den, -1.0, 1.0); } - const double dcx = swx / sw - g.predicted_x, dcy = swy / sw - g.predicted_y; - const double dpx = peak_x - g.predicted_x, dpy = peak_y - g.predicted_y; - sdx += dcx; sdy += dcy; sd2 += dcx * dcx + dcy * dcy; ++nc; - const double rx = g.predicted_x - beam_x, ry = g.predicted_y - beam_y; - const double rr = std::sqrt(rx * rx + ry * ry); - if (rr < 1.0) - continue; - const double rad_c = (dcx * rx + dcy * ry) / rr; // signed radial (outward +) - const double tang_c = (dcx * -ry + dcy * rx) / rr; // signed tangential - const double tang_p = (dpx * -ry + dpy * rx) / rr; - offs.push_back({signif, std::fabs(tang_c), std::fabs(tang_p), rad_c}); + if (signif >= 5.0) // only well-measured spots define the shape + bin_M2[bin_of(g.d)].push_back(sw_et2 / sw); } - if (nc >= 5) { - data.centroid_bias_px = std::sqrt((sdx / nc) * (sdx / nc) + (sdy / nc) * (sdy / nc)); - data.centroid_rms_px = std::sqrt(sd2 / nc); - } - if (offs.size() >= 10) { - std::vector sig; - sig.reserve(offs.size()); - for (const auto &o : offs) - sig.push_back(o.signif); - std::nth_element(sig.begin(), sig.begin() + sig.size() / 2, sig.end()); - const double smed = sig[sig.size() / 2]; - double slo = 0, shi = 0, tclo = 0, tchi = 0, tplo = 0, tphi = 0, rclo = 0, rchi = 0; - int nlo = 0, nhi = 0; - for (const auto &o : offs) { - if (o.signif < smed) { slo += o.signif; tclo += o.tang_c; tplo += o.tang_p; rclo += o.rad_c; ++nlo; } - else { shi += o.signif; tchi += o.tang_c; tphi += o.tang_p; rchi += o.rad_c; ++nhi; } + std::vector bin_R1(n_bins, data.R[1]); + for (int b = 0; b < n_bins; ++b) + if (bin_M2[b].size() >= 5) { + const double r1 = std::sqrt(2.0 * std::max(MedianInPlace(bin_M2[b]), 0.0)); + if (std::isfinite(r1) && r1 > 1e-4) + bin_R1[b] = std::clamp(r1, 1e-4, 0.05); } - if (nlo > 0 && nhi > 0) { - data.centroid_lo_signif = slo / nlo; data.centroid_hi_signif = shi / nhi; - data.centroid_lo_tang_c = tclo / nlo; data.centroid_hi_tang_c = tchi / nhi; - data.centroid_lo_tang_p = tplo / nlo; data.centroid_hi_tang_p = tphi / nhi; - data.centroid_lo_rad_c = rclo / nlo; data.centroid_hi_rad_c = rchi / nhi; - } - } + for (auto &g : groups) + g.R1_eff = bin_R1[bin_of(g.d)]; } - // ---- Extract integrated reflections --------------------------------------- - // Profile fitting gives the recorded amplitude (fitting the tangential profile - // P_t against the background-subtracted pixels): + // ---- 4. Term 1: one intensity residual per reflection; fit G (and B) ------ + // J / partiality / sigma_J are computed here as constants (geometry & R fixed), + // and only the per-image scale G and Debye-Waller B are optimised. + ceres::Problem problem; + size_t n_blocks = 0; + const double R0 = data.R[0]; + for (const auto &g : groups) { + const double R1 = g.R1_eff; // Term 2: per-resolution profile width + if (!(R1 > 0.0) || !(R0 > 0.0)) + continue; + double num = 0.0, den = 0.0, rad = 0.0; + std::vector> pt_sig; // (P_t, Iobs-Bg) for the Fisher pass + pt_sig.reserve(g.pixels.size()); + for (const auto &px : g.pixels) { + double q_sq, eps_r, eps_t_sq; + if (!GeometryProbe(px.x, px.y, lambda, pixel_size, g.h, g.k, g.l, data.crystal_system, + beam, dist_mm, detector_rot, latt_vec0, latt_vec1, latt_vec2, + q_sq, eps_r, eps_t_sq)) + continue; + const double P_t = std::exp(-eps_t_sq / (R1 * R1)) / (M_PI * R1 * R1); + const double R0_eff_sq = R0 * R0 + g.R_bw_sq; + const double P_rad = std::exp(-eps_r * eps_r / R0_eff_sq); + const double v = std::max(g.Ibkg, 1.0); + const double sig = px.Iobs - g.Ibkg; + num += P_t * sig / v; + den += P_t * P_t / v; + rad += P_rad * P_t * P_t / v; + pt_sig.emplace_back(P_t, sig); + } + if (!(den > 0.0)) + continue; + const double J = num / den; + const double partiality = rad / den; + // Model-expected (Fisher) variance: v_p = background + expected signal J*P_t, + // not the per-pixel observed counts (which down-bias) - so the weight tracks + // information, and an expected-strong reflection that is absent hurts. + double den_f = 0.0; + for (const auto &[P_t, sig] : pt_sig) { + const double v_f = std::max(g.Ibkg + std::max(J, 0.0) * P_t, 1.0); + den_f += P_t * P_t / v_f; + } + const double sigma_J = std::sqrt(1.0 / std::max(den_f, 1e-30)); + const double inv_4d2 = (g.d > 0.0) ? 1.0 / (4.0 * g.d * g.d) : 0.0; + auto *cost = new ceres::AutoDiffCostFunction( + new IntensityResidual(J, sigma_J, partiality, g.pol, g.Itrue, inv_4d2)); + problem.AddResidualBlock(cost, nullptr, &data.scale_factor, &data.B_factor); + ++n_blocks; + } + data.residual_count = n_blocks; + if (n_blocks == 0) + return; + + // G >= 0; B fixed unless requested; G regularized -> 1 with weight sqrt(n/sigma) + // (mirrors ScaleOnTheFly) so weakly-measured images cannot drift and scramble the merge. + problem.SetParameterLowerBound(&data.scale_factor, 0, 0.0); + if (!data.refine_B) + problem.SetParameterBlockConstant(&data.B_factor); + if (data.scale_reg_sigma > 0.0) { + const double w = std::sqrt(static_cast(groups.size()) / data.scale_reg_sigma); + auto *reg = new ceres::AutoDiffCostFunction( + new ScalarRegularizer(w, 1.0)); + problem.AddResidualBlock(reg, nullptr, &data.scale_factor); + } + + ceres::Solver::Options options; + options.linear_solver_type = ceres::DENSE_QR; + options.minimizer_progress_to_stdout = false; + options.logging_type = ceres::LoggingType::SILENT; + options.max_solver_time_in_seconds = data.max_time_s; + options.num_threads = 1; + ceres::Solver::Summary summary; + ceres::Solve(options, &problem, &summary); + data.final_cost = summary.final_cost; + data.solved = summary.IsSolutionUsable(); + + // ---- 5. Extract integrated reflections ------------------------------------ + // Profile fitting gives the recorded amplitude (fitting the tangential profile P_t + // against the background-subtracted pixels): // J = sum_p[ P_t,p (Iobs_p - Ibkg)/v_p ] / sum_p[ P_t,p^2 / v_p ] // ~ G * Itrue * B_term * partiality * pol (recorded raw counts) // var(J) = 1 / sum_p[ P_t,p^2 / v_p ] // - // Two SEPARATE fractions reduce the full intensity to what these pixels record: - // - // partiality - the radial / rocking dimension that a still does NOT sample. - // Only the slice of the reflection that crosses the Ewald - // sphere on this shot is recorded; <= 1. We DIVIDE it out to - // recover the full intensity. = profile-weighted P_radial. - // - // completeness - the fraction of the spot's detector footprint that landed on - // live pixels (= profile captured by live pixels / profile over - // the whole shoebox). 1.0 when the spot sits fully on the - // detector; < 1.0 only when a detector edge, gap or mask clips - // it. Profile fitting already extrapolates over the missing - // pixels, so this is NOT applied to r.I - it is a quality flag. - // // Output split (Merge multiplies r.I * image_scale_corr and weights by - // 1/(sigma*image_scale_corr)^2 - see Merge.cpp): + // 1/(sigma*image_scale_corr)^2): // r.I = J / (B_term * partiality * pol) = G * Itrue // r.sigma = sqrt(var(J)) / (B_term * partiality * pol) // r.partiality = profile-weighted P_radial in (0,1] (the rocking fraction) // r.completeness = live/total tangential profile in (0,1] (detector clipping) - // r.image_scale_corr = 1/G (per-image scale ONLY) - // so r.I * image_scale_corr = Itrue. B, partiality and polarization live on the - // intensity, G lives on image_scale_corr - one clean meaning per field. - // - // We walk the full (unclamped) shoebox once: every grid point feeds the total - // tangential profile (completeness denominator); points that are real, live - // detector pixels also feed the profile fit and the captured profile. + // r.image_scale_corr = 1/G (per-image scale ONLY) + // so r.I * image_scale_corr = Itrue: B, partiality and polarization live on the + // intensity, G lives on image_scale_corr - one clean meaning per field. We walk the + // full (unclamped) shoebox once: every grid point feeds the total tangential profile + // (completeness denominator); points that are real, live pixels also feed the fit. data.reflections.reserve(groups.size()); for (const auto &g : groups) { const int cx = static_cast(std::lround(g.predicted_x)); const int cy = static_cast(std::lround(g.predicted_y)); - - // Debye-Waller factor for this reflection (constant over its shoebox). const double B_term = std::exp(-data.B_factor / (4.0 * g.d * g.d)); - // Term 3 recentre shift (precomputed in the pre-pass; 0 if off or the spot was not - // confident). Profile evaluated at (x-dcx, y-dcy), data summed at (x,y). - const double dcx = g.dcx, dcy = g.dcy; - double num = 0.0, den = 0.0, bkg_sum = 0.0, radial_sum = 0.0; double prof_live = 0.0, prof_full = 0.0; // tangential profile: captured / total size_t n = 0; for (int y = cy - radius; y <= cy + radius; ++y) { for (int x = cx - radius; x <= cx + radius; ++x) { - // Geometry/profile for this grid point (profile recentred by (dcx,dcy)). - PixelObs probe{static_cast(x) - dcx, static_cast(y) - dcy, - 0.0, g.Ibkg, 1.0}; - PixelResidual pr(probe, 1.0, lambda, pixel_size, g.h, g.k, g.l, - g.R_bw_sq, g.pol, data.crystal_system); double q_sq, eps_r, eps_t_sq; - if (!pr.GeometryTerms(beam, &dist_mm, detector_rot, - latt_vec0, latt_vec1, latt_vec2, q_sq, eps_r, eps_t_sq)) + if (!GeometryProbe(static_cast(x), static_cast(y), + lambda, pixel_size, g.h, g.k, g.l, data.crystal_system, + beam, dist_mm, detector_rot, latt_vec0, latt_vec1, latt_vec2, + q_sq, eps_r, eps_t_sq)) continue; if (!(data.R[0] > 0.0) || !(g.R1_eff > 0.0)) continue; - // Tangential profile shape (area-normalized) -> the fit template. Uses the - // per-reflection R1_eff (Term 2), falling back to the global R1 by default. + // Tangential profile shape (area-normalized) -> the fit template, using the + // per-reflection R1_eff (Term 2). const double R1 = g.R1_eff; const double P_t = std::exp(-eps_t_sq / (R1 * R1)) / (M_PI * R1 * R1); prof_full += P_t; // whole shoebox, on- or off-detector @@ -1535,30 +687,21 @@ void PixelRefine::Run(const T *image, continue; const double Iobs = static_cast(image[np]); // raw counts - // Variance for the profile-fit weights. Weighting by the *observed* - // per-pixel count (v = Iobs) biases the amplitude negative: a - // down-fluctuated background pixel gets the smallest v and hence the - // largest 1/v weight, so num = sum P_t (Iobs - Ibkg)/v is pulled below - // zero - worst where the signal is weakest, i.e. the high-resolution - // shells (the negative we see there). For background-limited - // reflections the variance is the local background, constant over the - // shoebox, so use that instead of the observed count. - double v = std::max(g.Ibkg, 1.0); + // Background-limited variance (constant over the shoebox): weighting by + // the observed count biases the amplitude negative where signal is weakest. + const double v = std::max(g.Ibkg, 1.0); - // Peak-normalized radial factor (the partiality), in (0,1]. The - // bandwidth-broadened radial width matches the model in Model(). + // Peak-normalized radial factor (the partiality), in (0,1]. MUST use the + // same P_t^2/v weights as the amplitude, else an R0_eff-dependent (hence + // resolution-dependent) factor is left behind in r.I. const double R0_eff_sq = data.R[0] * data.R[0] + g.R_bw_sq; const double P_radial = std::exp(-eps_r * eps_r / R0_eff_sq); - // Profile-fit accumulators. The amplitude estimator weights pixels by - // P_t^2/v, so the partiality (which de-scales that amplitude) MUST use - // the SAME weights - otherwise an R0_eff-dependent (resolution- - // dependent) factor is left behind in r.I. const double w = P_t * P_t / v; num += P_t * (Iobs - g.Ibkg) / v; den += w; - radial_sum += P_radial * w; // partiality weighted exactly like num/den - prof_live += P_t; // captured tangential profile + radial_sum += P_radial * w; + prof_live += P_t; bkg_sum += g.Ibkg; ++n; } @@ -1580,7 +723,7 @@ void PixelRefine::Run(const T *image, if (den > 0.0 && n > 0) { const double I_amp = num / den; // ~ G*Itrue*B_term*partiality*pol const double sigma_amp = std::sqrt(1.0 / den); - const double corr = static_cast(r.partiality) * B_term * g.pol; // B, partiality & pol + const double corr = static_cast(r.partiality) * B_term * g.pol; r.bkg = static_cast(bkg_sum / static_cast(n)); r.observed = true; @@ -1602,9 +745,9 @@ void PixelRefine::Run(const T *image, data.reflections.push_back(r); } - // ---- Per-image CC vs reference (the half/ref correlation diagnostic) ------- - // Pearson CC of the scaled estimate (r.I * image_scale_corr = Itrue_est) - // against the reference intensities, over the matched reflections. + // ---- 6. Per-image CC vs reference (diagnostic) ---------------------------- + // Pearson CC of the scaled estimate (r.I * image_scale_corr = Itrue_est) against + // the reference intensities, over the matched reflections. { double sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0; size_t cn = 0; @@ -1633,231 +776,9 @@ void PixelRefine::Run(const T *image, } } -template -std::vector PixelRefine::PredictImage(const T *image, - BraggPrediction &prediction, - const PixelRefineData &data, - bool include_background) const { - std::vector img(xpixel * ypixel, 0.0f); - - const double lambda = data.geom.GetWavelength_A(); - const double pixel_size = data.geom.GetPixelSize_mm(); - const int radius = data.shoebox_radius; - const int bkg_outer_radius = std::max(radius + 1, data.bkg_outer_radius); - const double bw = data.bandwidth; - - const auto pol_factor = experiment.GetPolarizationFactor(); - auto polarization = [&](double x, double y) -> double { - if (!pol_factor) - return 1.0; - return data.geom.CalcAzIntPolarizationCorr(static_cast(x), static_cast(y), - pol_factor.value()); - }; - auto bandwidth_radial_sq = [&](double d) -> double { - if (bw <= 0.0 || d <= 0.0) - return 0.0; - const double bl = bw * lambda; - return bl * bl / (2.0 * d * d * d * d); - }; - - double beam[2], dist_mm, detector_rot[2]; - double latt_vec0[3], latt_vec1[3], latt_vec2[3]; - BuildParameterBlocks(data, beam, dist_mm, detector_rot, latt_vec0, latt_vec1, latt_vec2); - - DiffractionExperiment exp_iter = experiment; - exp_iter.BeamX_pxl(data.geom.GetBeamX_pxl()) - .BeamY_pxl(data.geom.GetBeamY_pxl()) - .DetectorDistance_mm(data.geom.GetDetectorDistance_mm()) - .PoniRot1_rad(data.geom.GetPoniRot1_rad()) - .PoniRot2_rad(data.geom.GetPoniRot2_rad()); - - const BraggPredictionSettings settings_prediction{ - .high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(), - .ewald_dist_cutoff = static_cast(data.ewald_dist_cutoff), - .max_hkl = 100, - .centering = data.centering, - .bandwidth_sigma = static_cast(data.bandwidth) // relative Δλ/λ sigma - }; - const int nrefl = prediction.Calc(exp_iter, data.latt, settings_prediction); - const auto &predicted = prediction.GetReflections(); - const auto spot_mask = BuildSpotMask(predicted, nrefl, xpixel, ypixel, radius); - - for (int ri = 0; ri < nrefl; ++ri) { - const auto &refl = predicted[ri]; - const auto it = reference_data.find(hkl_key_generator(refl)); - if (it == reference_data.end()) - continue; - - const double Itrue = it->second; - const double R_bw_sq = bandwidth_radial_sq(refl.d); - const double pol = polarization(refl.predicted_x, refl.predicted_y); - - // Local background straight from the actual image (flat per shoebox), laid - // into the box so the prediction overlays the real frame - the same model - // path Run() fits, now reproduced faithfully because we have the image. - double Ibkg = 0.0; - const bool have_bkg = include_background && - EstimateLocalBackground(image, spot_mask, xpixel, ypixel, - refl.predicted_x, refl.predicted_y, - radius, bkg_outer_radius, Ibkg); - - const auto box = ShoeboxBounds(refl.predicted_x, refl.predicted_y, radius, xpixel, ypixel); - - for (int y = box.min_y; y <= box.max_y; ++y) { - for (int x = box.min_x; x <= box.max_x; ++x) { - const size_t npixel = xpixel * y + x; - - PixelObs obs{ - .x = static_cast(x), - .y = static_cast(y), - .Iobs = 0.0, - .Ibkg = have_bkg ? Ibkg : 0.0, - .weight = 1.0 - }; - PixelResidual pr(obs, Itrue, lambda, pixel_size, - refl.h, refl.k, refl.l, R_bw_sq, pol, data.crystal_system); - - double Ipred = 0.0; // raw counts: signal (+ local background) - if (pr.Model(beam, &dist_mm, detector_rot, - latt_vec0, latt_vec1, latt_vec2, - &data.scale_factor, &data.B_factor, data.R, Ipred)) - img[npixel] += static_cast(Ipred); - } - } - } - - return img; -} - -template -std::vector PixelRefine::ChiSquaredImage(const T *image, - BraggPrediction &prediction, - const PixelRefineData &data) const { - std::vector img(xpixel * ypixel, 0.0f); - - const double lambda = data.geom.GetWavelength_A(); - const double pixel_size = data.geom.GetPixelSize_mm(); - const int radius = data.shoebox_radius; - const int bkg_outer_radius = std::max(radius + 1, data.bkg_outer_radius); - const double bw = data.bandwidth; - - const auto pol_factor = experiment.GetPolarizationFactor(); - auto polarization = [&](double x, double y) -> double { - if (!pol_factor) - return 1.0; - return data.geom.CalcAzIntPolarizationCorr(static_cast(x), static_cast(y), - pol_factor.value()); - }; - auto bandwidth_radial_sq = [&](double d) -> double { - if (bw <= 0.0 || d <= 0.0) - return 0.0; - const double bl = bw * lambda; - return bl * bl / (2.0 * d * d * d * d); - }; - - double beam[2], dist_mm, detector_rot[2]; - double latt_vec0[3], latt_vec1[3], latt_vec2[3]; - BuildParameterBlocks(data, beam, dist_mm, detector_rot, latt_vec0, latt_vec1, latt_vec2); - - DiffractionExperiment exp_iter = experiment; - exp_iter.BeamX_pxl(data.geom.GetBeamX_pxl()) - .BeamY_pxl(data.geom.GetBeamY_pxl()) - .DetectorDistance_mm(data.geom.GetDetectorDistance_mm()) - .PoniRot1_rad(data.geom.GetPoniRot1_rad()) - .PoniRot2_rad(data.geom.GetPoniRot2_rad()); - - const BraggPredictionSettings settings_prediction{ - .high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(), - .ewald_dist_cutoff = static_cast(data.ewald_dist_cutoff), - .max_hkl = 100, - .centering = data.centering, - .bandwidth_sigma = static_cast(data.bandwidth) - }; - const int nrefl = prediction.Calc(exp_iter, data.latt, settings_prediction); - const auto &predicted = prediction.GetReflections(); - const auto spot_mask = BuildSpotMask(predicted, nrefl, xpixel, ypixel, radius); - - for (int ri = 0; ri < nrefl; ++ri) { - const auto &refl = predicted[ri]; - const auto it = reference_data.find(hkl_key_generator(refl)); - if (it == reference_data.end()) - continue; - - const double Itrue = it->second; - const double R_bw_sq = bandwidth_radial_sq(refl.d); - const double pol = polarization(refl.predicted_x, refl.predicted_y); - - // Local flat background, identical to Run(); skip the reflection if it - // cannot be estimated (matches Run() dropping the reflection). - double Ibkg = 0.0; - if (!EstimateLocalBackground(image, spot_mask, xpixel, ypixel, - refl.predicted_x, refl.predicted_y, - radius, bkg_outer_radius, Ibkg)) - continue; - - const auto box = ShoeboxBounds(refl.predicted_x, refl.predicted_y, radius, xpixel, ypixel); - - for (int y = box.min_y; y <= box.max_y; ++y) { - for (int x = box.min_x; x <= box.max_x; ++x) { - const size_t npixel = xpixel * y + x; - - // Same gating as Run(): only pixels that actually enter the fit. - if (image[npixel] == std::numeric_limits::max()) - continue; - if (std::is_signed_v && (image[npixel] == std::numeric_limits::min())) - continue; - - const double Iobs = static_cast(image[npixel]); // raw counts - - double var = std::max(Iobs, 0.0); - if (!(var > 1.0)) - var = 1.0; - const double weight = 1.0 / std::sqrt(var); - - PixelObs obs{ - .x = static_cast(x), - .y = static_cast(y), - .Iobs = Iobs, - .Ibkg = Ibkg, - .weight = weight - }; - PixelResidual pr(obs, Itrue, lambda, pixel_size, - refl.h, refl.k, refl.l, R_bw_sq, pol, data.crystal_system); - - double Ipred = 0.0; - if (pr.Model(beam, &dist_mm, detector_rot, - latt_vec0, latt_vec1, latt_vec2, - &data.scale_factor, &data.B_factor, data.R, Ipred)) { - // residual_i = (I_pred - I_obs) * weight (== Ceres residual); - // its square is this pixel's contribution to the cost. - const double rw = (Ipred - Iobs) * weight; - img[npixel] += static_cast(rw * rw); - } - } - } - } - - return img; -} - -// Explicit instantiations for the supported (uncompressed) image pixel types. template void PixelRefine::Run(const int8_t *, BraggPrediction &, PixelRefineData &); template void PixelRefine::Run(const int16_t *, BraggPrediction &, PixelRefineData &); template void PixelRefine::Run(const int32_t *, BraggPrediction &, PixelRefineData &); template void PixelRefine::Run(const uint8_t *, BraggPrediction &, PixelRefineData &); template void PixelRefine::Run(const uint16_t *, BraggPrediction &, PixelRefineData &); template void PixelRefine::Run(const uint32_t *, BraggPrediction &, PixelRefineData &); - -template std::vector PixelRefine::PredictImage(const int8_t *, BraggPrediction &, const PixelRefineData &, bool) const; -template std::vector PixelRefine::PredictImage(const int16_t *, BraggPrediction &, const PixelRefineData &, bool) const; -template std::vector PixelRefine::PredictImage(const int32_t *, BraggPrediction &, const PixelRefineData &, bool) const; -template std::vector PixelRefine::PredictImage(const uint8_t *, BraggPrediction &, const PixelRefineData &, bool) const; -template std::vector PixelRefine::PredictImage(const uint16_t *, BraggPrediction &, const PixelRefineData &, bool) const; -template std::vector PixelRefine::PredictImage(const uint32_t *, BraggPrediction &, const PixelRefineData &, bool) const; - -template std::vector PixelRefine::ChiSquaredImage(const int8_t *, BraggPrediction &, const PixelRefineData &) const; -template std::vector PixelRefine::ChiSquaredImage(const int16_t *, BraggPrediction &, const PixelRefineData &) const; -template std::vector PixelRefine::ChiSquaredImage(const int32_t *, BraggPrediction &, const PixelRefineData &) const; -template std::vector PixelRefine::ChiSquaredImage(const uint8_t *, BraggPrediction &, const PixelRefineData &) const; -template std::vector PixelRefine::ChiSquaredImage(const uint16_t *, BraggPrediction &, const PixelRefineData &) const; -template std::vector PixelRefine::ChiSquaredImage(const uint32_t *, BraggPrediction &, const PixelRefineData &) const; diff --git a/image_analysis/pixel_refinement/PixelRefine.h b/image_analysis/pixel_refinement/PixelRefine.h index f8db193f..dfc560c5 100644 --- a/image_analysis/pixel_refinement/PixelRefine.h +++ b/image_analysis/pixel_refinement/PixelRefine.h @@ -3,247 +3,100 @@ #pragma once -#include - #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*). 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 reflections; // profile-fitted integration result + std::vector 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 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 - void SweepOrientationCell(const T *image, BraggPrediction &prediction, - PixelRefineData &data) const; public: PixelRefine(const DiffractionExperiment &experiment, const std::vector &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 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 - std::vector 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 - std::vector 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 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; - } }; diff --git a/viewer/CMakeLists.txt b/viewer/CMakeLists.txt index 6bfbd37e..d7f0875c 100644 --- a/viewer/CMakeLists.txt +++ b/viewer/CMakeLists.txt @@ -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 diff --git a/viewer/JFJochImageReadingWorker.cpp b/viewer/JFJochImageReadingWorker.cpp index 1e92cfed..a10ded30 100644 --- a/viewer/JFJochImageReadingWorker.cpp +++ b/viewer/JFJochImageReadingWorker.cpp @@ -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"); - qRegisterMetaType("PixelRefineReport"); qRegisterMetaType>("QVector"); spot_finding_settings = settings;; @@ -307,13 +303,7 @@ void JFJochImageReadingWorker::UpdateAzint_i(const JFJochReaderDataset *dataset) image_analysis = std::make_unique(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(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(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(p.beam_x)).BeamY_pxl(static_cast(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 JFJochImageReadingWorker::WrapFloatImage_i(const std::vector &img) const { - auto si = std::make_shared(); - // 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 &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(v); - pred[i] = diff * diff; - } - } -} - -void JFJochImageReadingWorker::MaskMeasuredSentinels_i(std::vector &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 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 boxes; - boxes.reserve(static_cast(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(std::lround(refl.predicted_x)); - const int cy = static_cast(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(r.I) * static_cast(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 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(img32.data(), *pixel_pred_, data); - MaskMeasuredSentinels_i(chi2); - return chi2; - } - - auto pred = pixel_refine_->PredictImage(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(img32.data(), *pixel_pred_, d); - emit pixelRefineResidual(d.final_cost, d.cc, static_cast(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(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(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(*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())); - } -} diff --git a/viewer/JFJochImageReadingWorker.h b/viewer/JFJochImageReadingWorker.h index 50cd77ca..fcf74d5d 100644 --- a/viewer/JFJochImageReadingWorker.h +++ b/viewer/JFJochImageReadingWorker.h @@ -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) @@ -60,30 +58,8 @@ private: std::unique_ptr image_analysis; std::unique_ptr 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 last_profile_; - std::vector pixel_reference_; - std::unique_ptr pixel_refine_; - std::unique_ptr pixel_pred_; - - void EnsurePixelRefine_i(); - bool BuildPixelSeed_i(PixelRefineData &d, const PixelRefineParams &p, QString &reason) const; - std::shared_ptr WrapFloatImage_i(const std::vector &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 &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 &img) const; - // Build the per-reflection shoebox rectangles for the last refine/preview. - QVector 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 BuildDisplayImage_i(const PixelRefineData &data, int display_mode) const; std::unique_ptr 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 image); - void predictedShoeboxes(QVector 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); }; diff --git a/viewer/JFJochViewerWindow.cpp b/viewer/JFJochViewerWindow.cpp index 1eb2c3c7..903b7528 100644 --- a/viewer/JFJochViewerWindow.cpp +++ b/viewer/JFJochViewerWindow.cpp @@ -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); diff --git a/viewer/image_viewer/JFJochImage.h b/viewer/image_viewer/JFJochImage.h index 5b3816ea..5e98e599 100644 --- a/viewer/image_viewer/JFJochImage.h +++ b/viewer/image_viewer/JFJochImage.h @@ -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; diff --git a/viewer/image_viewer/JFJochSimpleImage.cpp b/viewer/image_viewer/JFJochSimpleImage.cpp index d3494a8f..bfa029e7 100644 --- a/viewer/image_viewer/JFJochSimpleImage.cpp +++ b/viewer/image_viewer/JFJochSimpleImage.cpp @@ -40,30 +40,6 @@ void JFJochSimpleImage::setImage(std::shared_ptr img) { } } -void JFJochSimpleImage::setShoeboxes(QVector 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_) { diff --git a/viewer/image_viewer/JFJochSimpleImage.h b/viewer/image_viewer/JFJochSimpleImage.h index 40abac06..3232bc7c 100644 --- a/viewer/image_viewer/JFJochSimpleImage.h +++ b/viewer/image_viewer/JFJochSimpleImage.h @@ -20,20 +20,14 @@ class JFJochSimpleImage : public JFJochImage { std::shared_ptr image_; - // Per-reflection shoebox rectangles (pixel coordinates) to overlay: the pixels - // PixelRefine actually summed over. Empty = nothing drawn. - QVector shoeboxes_; - // Prepare image template 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 img); - void setShoeboxes(QVector boxes); }; diff --git a/viewer/windows/JFJochPixelRefineTableWindow.cpp b/viewer/windows/JFJochPixelRefineTableWindow.cpp deleted file mode 100644 index f656164a..00000000 --- a/viewer/windows/JFJochPixelRefineTableWindow.cpp +++ /dev/null @@ -1,101 +0,0 @@ -// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#include "JFJochPixelRefineTableWindow.h" - -#include -#include -#include - -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(value), Qt::DisplayRole); - it->setData(static_cast(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("PixelRefine: scale G = %1, B = %2 Ų, CC = %3 (%4 refl)" - "  |  " - "Pipeline: 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 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); - } -} diff --git a/viewer/windows/JFJochPixelRefineTableWindow.h b/viewer/windows/JFJochPixelRefineTableWindow.h deleted file mode 100644 index 5187b546..00000000 --- a/viewer/windows/JFJochPixelRefineTableWindow.h +++ /dev/null @@ -1,34 +0,0 @@ -// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#pragma once - -#include "JFJochHelperWindow.h" -#include "PixelRefineParams.h" - -#include -#include -#include -#include - -// 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); -}; diff --git a/viewer/windows/JFJochPixelRefineWindow.cpp b/viewer/windows/JFJochPixelRefineWindow.cpp deleted file mode 100644 index 0e8b4a16..00000000 --- a/viewer/windows/JFJochPixelRefineWindow.cpp +++ /dev/null @@ -1,233 +0,0 @@ -// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#include "JFJochPixelRefineWindow.h" -#include "../image_viewer/JFJochSimpleImage.h" - -#include -#include -#include -#include -#include -#include -#include - -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 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(image->Dataset().experiment.GetXPixelsNum())); - m_beamy->setMax(static_cast(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 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); -} diff --git a/viewer/windows/JFJochPixelRefineWindow.h b/viewer/windows/JFJochPixelRefineWindow.h deleted file mode 100644 index 9398389d..00000000 --- a/viewer/windows/JFJochPixelRefineWindow.h +++ /dev/null @@ -1,79 +0,0 @@ -// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#pragma once - -#include "JFJochHelperWindow.h" -#include "PixelRefineParams.h" -#include "../SimpleImage.h" -#include "../widgets/SliderPlusBox.h" - -#include -#include -#include -#include -#include -#include - -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 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 image); - void setResidual(double cost, double cc, int64_t n_reflections); - void setRefinedParams(PixelRefineParams params); - void setStatus(QString message); -}; diff --git a/viewer/windows/PixelRefineParams.h b/viewer/windows/PixelRefineParams.h deleted file mode 100644 index e180847e..00000000 --- a/viewer/windows/PixelRefineParams.h +++ /dev/null @@ -1,66 +0,0 @@ -// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#pragma once - -#include -#include -#include -#include - -// 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 rows; -}; - -Q_DECLARE_METATYPE(PixelRefineReport)