From 6f2033db004899c37530bf2431f54e8c1a16054b Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 13 Jun 2026 21:34:33 +0200 Subject: [PATCH] PixelRefine: checkpoint before cleanup (factored model + all diagnostic levers) Snapshot of the messy state: factored likelihood Terms 1+2+3 behind PR_* env flags (PR_INTENSITY/PR_SHAPE/PR_RECENTER) alongside the old per-pixel ShoeboxResidual, plus diagnostic scaffolding (PR_R0/R1/COV/FIX_R0/FIX_R/ ADAPT_R1/CENTROID/RECENTER) and the FACTORED_MODEL.md spec. Next commit makes Terms 1+2 the model and strips all of this. Co-Authored-By: Claude Opus 4.8 --- image_analysis/IndexAndRefine.cpp | 34 ++ .../pixel_refinement/FACTORED_MODEL.md | 140 +++++++ .../pixel_refinement/PixelRefine.cpp | 381 +++++++++++++++++- image_analysis/pixel_refinement/PixelRefine.h | 54 +++ 4 files changed, 601 insertions(+), 8 deletions(-) create mode 100644 image_analysis/pixel_refinement/FACTORED_MODEL.md diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index f7cb25c1..ccf084fe 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -446,6 +446,28 @@ 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()) { @@ -465,6 +487,18 @@ 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 new file mode 100644 index 00000000..944514a2 --- /dev/null +++ b/image_analysis/pixel_refinement/FACTORED_MODEL.md @@ -0,0 +1,140 @@ +# A factored likelihood for joint integration + scaling + geometry + +**Status: design spec (not implemented).** 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 +thing that dissolves the empty-pixel and parameter-degeneracy problems by construction +rather than by patching. + +## 0. Notation + +Per image, parameters `θ`: scale `G`, Debye-Waller `B`, orientation + cell (geometry), +profile width `R1` (tangential, possibly a 2×2 tensor), partiality width `R0` +(radial/mosaicity; `R0_eff² = R0² + R_bw²`, `R_bw² = (bλ)²/2d⁴` *known* from bandwidth). +Per reflection `h`: reference intensity `I_ref` (the hypothesis), resolution `d`, +predicted centre `c_pred`, partiality `p = exp(−ε_r²/R0_eff²)`, polarisation `pol`, +`B_term = exp(−B/4d²)`, shoebox pixels `{I_p}` with mean local background `Bg`, and the +area-normalised tangential profile template `P_p = P_tang(ε_t,p; R1)`. + +## 1. The factorisation principle + +A reflection's shoebox carries three (to first order) **orthogonal** pieces of +information — the 0th, 1st and 2nd moments of its intensity distribution: + +| moment | statistic | constrains | +|---|---|---| +| 0th — total | profile-fit amplitude `J` | scale chain `G, B` (and `p`) | +| 1st — position | centroid `c_obs` | geometry (orientation; radial→distance/cell) | +| 2nd — shape | second moment `M₂` | profile width `R1` (and anisotropy) | + +The current per-pixel residual mixes all three into one objective over shared pixels — +*that* is what couples the parameters (measured G–R0 ≈ −0.46, G–R1 ≈ +0.51) and lets the +many empty pixels dominate. Residual-ing each **moment** against its model instead gives +a block-diagonal Jacobian: the couplings vanish because each statistic carries one +parameter block's information. + +## 2. The three residual terms + +### 2.1 Intensity / scaling residual (one scalar per reflection) + +Optimal (Diamond) profile-fit amplitude and its model: +``` +J = Σ_p w_p P_p (I_p − Bg) / Σ_p w_p P_p² w_p = 1/v_p +J_model = G · B_term · p · pol · I_ref +r¹_h = (J − J_model) / σ_J +``` +`J` is ~invariant to `R1` (a well-sampled spot integrates to the same total whatever +width is assumed) → **R1 leaves this residual**. Empty pixels make no residual; they +enter only through `J` with ~zero profile weight → **the empty-pixel problem is gone by +construction.** This residual *is* the scaling residual — integration and scaling are now +one objective. + +### 2.2 Shape residual (constrains R1; decoupled from scale) + +``` +M₂_obs = Σ_p (I_p − Bg) ε_t,p² / Σ_p (I_p − Bg) (intensity-weighted variance, Å⁻²) +M₂_model = R1² / 2 (variance of exp(−ε_t²/R1²)) +r²_h = (M₂_obs − M₂_model) / σ_M2 +``` +A moment is normalised by the total → **scale-invariant → `∂r²/∂G = 0`**. The G↔R1 +degeneracy disappears. Anisotropic extension: use the 2×2 moment tensor +`Σ(I−Bg)(ε_t⊗ε_t)/Σ(I−Bg)` vs `diag(R1a²/2, R1b²/2)` → elliptical R1 (the DMM streak). +Weak spots have huge `σ_M2` → contribute ~nothing → R1 is set by strong spots +automatically (and may be made `R1(d)` per resolution). + +### 2.3 Position residual (constrains geometry; decoupled from scale and shape) + +``` +c_obs = Σ_p (I_p − Bg)(x_p, y_p) / Σ_p (I_p − Bg) +r³_h = (c_obs − c_pred(geometry)) / σ_c (2-vector; split radial / tangential) +``` +Centroid is scale- and width-invariant → `∂r³/∂G = ∂r³/∂R1 ≈ 0`. The **radial** component +constrains distance/cell, the **tangential** constrains orientation — exactly the split +the diagnostic measured (radial≈0 = no distance error; tangential∝radius = orientation). + +## 3. Fisher / expected-variance weighting (makes it a likelihood) + +Every `σ` uses the **model-expected** variance, never observed counts — this is what +makes strong *expected* reflections carry the information and makes the model "feel pain +when something that should be there is not": +``` +v_p = Bg + J_model · P_p (background + expected signal from I_ref, not I_obs) +σ_J² = 1 / Σ_p (P_p² / v_p) +σ_M2 ≈ M₂ · √(2 / N_eff), σ_c ≈ R1 / √(N_eff), N_eff = (Σ(I−Bg))² / Σ v_p +``` +Fisher information about `G` from term 1 is `∝ (B_term·p·pol·I_ref)² / σ_J²` — driven by +`I_ref`, so a noise spike (high counts, low `I_ref`) gets *no* weight while a strong +expected reflection observed absent (`J≈0`, large residual, moderate `σ_J`) gets a large +penalty. The reference enters at maximum leverage: it sets both the target and the weight. + +## 4. Joint objective and priors + +``` +L(θ) = Σ_h [ (r¹_h)² + (r²_h)² + |r³_h|² ] + priors +``` +No free λ if the σ's are correct — the relative weighting *is* the Fisher information. +Priors are the Bayesian hooks and the principled degeneracy breaks: +- **R0 (partiality/mosaicity) is GLOBAL + prior.** R0 multiplies `J_model` (`p`), so it is + still degenerate with the per-image `G` *within term 1* — the one degeneracy the + factorisation does **not** remove. Resolve it physically, not with a directional G prior + (which would bias every output intensity): `R0 ~ N(mosaicity, σ)`, `R_bw` fixed from the + known bandwidth, and `R0` fit **globally** (one per crystal, from many reflections' + partiality distribution) so per-image G can't trade against it. +- orientation `~ N(spot-centroid, σ)`; `G ~ N(1, σ_G)` or tied to the beam monitor; + distance `~ N(nominal, σ_L)` (loose, since serial/jet alignment is poorly constrained). +- Optional Bayesian intensities: treat `I_true` as a parameter with the reference as its + prior → posterior over intensities, not point estimates. + +## 5. Why the degeneracies vanish (Jacobian structure) + +`Jᵀ W J` is approximately block-diagonal in `(G,B,p | R1 | geometry)`: +``` +∂r¹/∂{G,B,p} ≠ 0 ; ∂r¹/∂R1 ≈ 0 ; ∂r¹/∂geom ≈ 0 +∂r²/∂R1 ≠ 0 ; ∂r²/∂G = 0 ; ∂r²/∂geom ≈ 0 +∂r³/∂geom ≠ 0 ; ∂r³/∂G = 0 ; ∂r³/∂R1 ≈ 0 +``` +So G↔R1 (+0.51) and all the cross-couplings drop to ~0 by construction. Only **G↔R0** +survives (R0 is a scale-multiplier, not a shape), handled by the global+physical prior of +§4. The degeneracies we measured were artifacts of projecting all information onto a +single per-pixel residual. + +## 6. Implementation notes + +- Per reflection: 1 (intensity) + 1 (shape) + 2 (position) residuals = 4, vs ~49 per-pixel + residuals → **cheaper**, and Ceres autodiffs the moment formulas through the pixels. +- The per-pixel forward model still *defines* `P_tang`, `p`, etc.; the **loss** moves to + the moments. +- Geometry (term 3) can run as the global sweep we have (it already maximises a + position/CC objective); terms 1–2 are the per-image photometry. Or solve all three + jointly per image with the global R0/mosaicity shared across images (two-level fit). +- Drop-in path: keep the current extraction, add the three residuals as a new objective + behind a flag, compare against the per-pixel loss on both test crystals. + +## 7. Why this serves the goal + +It is one differentiable likelihood, factored along the physics, that (a) maximises use of +the reference (target + Fisher weight), (b) is the substrate for priors / posteriors over +intensities (Bayesian), and (c) lets any term — profile `P`, partiality `p`, corrections — +be replaced by a learned function trained through the same likelihood. That is the +qualitative move XDS-style empirical profile fitting cannot make. diff --git a/image_analysis/pixel_refinement/PixelRefine.cpp b/image_analysis/pixel_refinement/PixelRefine.cpp index bef4392b..42e8772f 100644 --- a/image_analysis/pixel_refinement/PixelRefine.cpp +++ b/image_analysis/pixel_refinement/PixelRefine.cpp @@ -10,6 +10,7 @@ #include #include +#include #include #include "../geom_refinement/LatticeReduction.h" @@ -35,6 +36,8 @@ 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) std::vector pixels; }; @@ -752,6 +755,30 @@ void PixelRefine::SweepOrientationCell(const T *image, BraggPrediction &predicti 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, @@ -929,12 +956,146 @@ void PixelRefine::Run(const T *image, 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) + continue; + if (sw / std::sqrt(std::max(box_px * g.Ibkg, 1.0)) < data.recenter_min_signif) + 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; + } + } + + // ---- 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 { for (const auto &g : groups) { auto *cost = new ceres::DynamicAutoDiffCostFunction( new ShoeboxResidual(g, lambda, pixel_size, data.crystal_system)); @@ -959,8 +1120,23 @@ void PixelRefine::Run(const T *image, residual_pixels += g.pixels.size(); } 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) { @@ -1029,11 +1205,17 @@ void PixelRefine::Run(const T *image, problem.SetParameterBlockConstant(&data.B_factor); if (data.refine_R) { - problem.SetParameterLowerBound(data.R, 0, 1e-5); - problem.SetParameterLowerBound(data.R, 1, 1e-5); + 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 @@ -1057,6 +1239,50 @@ void PixelRefine::Run(const T *image, 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 --------------- @@ -1096,6 +1322,139 @@ void PixelRefine::Run(const T *image, } } // 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; + 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)) + 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 (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; } + } + 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; + } + } + } + // ---- Extract integrated reflections --------------------------------------- // Profile fitting gives the recorded amplitude (fitting the tangential profile // P_t against the background-subtracted pixels): @@ -1138,26 +1497,32 @@ void PixelRefine::Run(const T *image, // 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 (valid even off the detector). - PixelObs probe{static_cast(x), static_cast(y), 0.0, g.Ibkg, 1.0}; + // 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)) continue; - if (!(data.R[0] > 0.0) || !(data.R[1] > 0.0)) + if (!(data.R[0] > 0.0) || !(g.R1_eff > 0.0)) continue; - // Tangential profile shape (area-normalized) -> the fit template. - const double P_t = std::exp(-eps_t_sq / (data.R[1] * data.R[1])) - / (M_PI * data.R[1] * data.R[1]); + // Tangential profile shape (area-normalized) -> the fit template. Uses the + // per-reflection R1_eff (Term 2), falling back to the global R1 by default. + 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 // Only real, unmasked detector pixels carry signal. diff --git a/image_analysis/pixel_refinement/PixelRefine.h b/image_analysis/pixel_refinement/PixelRefine.h index fc4acabd..f8db193f 100644 --- a/image_analysis/pixel_refinement/PixelRefine.h +++ b/image_analysis/pixel_refinement/PixelRefine.h @@ -183,6 +183,56 @@ struct PixelRefineData { 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 bool solved = false; @@ -190,6 +240,10 @@ struct PixelRefineData { 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 {