diff --git a/common/BraggIntegrationSettings.h b/common/BraggIntegrationSettings.h index 1ac24669..8f33c87d 100644 --- a/common/BraggIntegrationSettings.h +++ b/common/BraggIntegrationSettings.h @@ -116,11 +116,17 @@ class BraggIntegrationSettings { // Offline that is what is wanted. ONLINE it is not: the broker bootstraps a concrete value // (BRAGG_ONLINE_DEFAULT_MAX_HKL) so per-image cost stays predictable across samples. std::optional max_hkl; - // Overlap treatment and, for OverlapMode::Reject, the least fraction of a reflection's expected - // profile that must be cleanly its own for the reflection to be kept (XDS calls it MINPK). - // Excluding the shared pixels is the default: over the rotation battery it costs 1.1% of the wall - // clock (23% on a genuinely crowded crystal, nothing where no two predictions touch) and buys ISa - // on 15 crystals against 5, cutting the summed shortfall against XDS by a third. + // Overlap treatment and the MINPK threshold: the least fraction of a reflection's expected profile + // that must be usable for the reflection to be kept. Excluding the shared pixels is the default: + // over the rotation battery it costs 1.1% of the wall clock (23% on a genuinely crowded crystal, + // nothing where no two predictions touch) and buys ISa on 15 crystals against 5, cutting the summed + // shortfall against XDS by a third. + // + // As in XDS, one threshold governs both ways a reflection can lose part of its profile. Under + // OverlapMode::Reject it is the fraction that must be cleanly the reflection's own rather than a + // neighbour's. In every profile mode it is also the fraction that must be READABLE - not masked, + // untrusted, in a detector gap or overloaded - because the profile fit renormalises to the pixels + // it can read (dials calls that valid_foreground_threshold, and defaults it to the same 0.75). OverlapMode overlap_mode = OverlapMode::Exclude; float overlap_min_peak = 0.75f; diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 74ea92a3..e0094f00 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -9,6 +9,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as * Bragg integration: the profile fit's `background_variance` now takes the fitted intensity itself out of the fit variance instead of `max(0, I)`, so a reflection that fluctuated below zero no longer reports a background variance two to three times too small and is no longer weighted up for it. * Scaling: the rotation merge weights each combined full by its variance rebuilt at the reflection's mean intensity rather than by the full's own sigma, as the stills merge already did. * Scaling: a rotation frame whose fitted scale collapses - it recorded no diffraction from the indexed lattice - is now **dropped from the merge** instead of being merged unscaled, which had asserted a scale of 1 for a frame demonstrably nowhere near it. +* Bragg integration: a reflection whose signal disk is cut by a mask, an untrusted region, a detector gap or an overload is now profile-fitted over the pixels that remain instead of being discarded, as long as at least `--overlap-minpk` of its expected profile is readable (XDS's MINPK); `--integrator boxsum` still discards it. * rugnux: Rotation data are integrated on **every frame whose spots the sweep's lattice explains**, instead of only on frames that would also index on their own; the reported indexing rate still counts the latter. * Scaling: a rotation frame too sparse to fit a rocking width of its own now takes the run's median instead of a fixed default. * rugnux: The error-model **a** and **b** are reported in XDS's convention, and `_reflns.jfjoch_diffrn_ISa` now carries the whole-range `1/sqrt(a*b)` that XDS's ISa denotes; the strong-reflection asymptote moves to `_reflns.jfjoch_diffrn_ISa_asymptotic`. **A file written by an earlier version carries the asymptote under the old name.** diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 2c00b97e..90e55694 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -589,7 +589,7 @@ One caveat on "exactly": the *classification* is bit-identical at $k=0$, and so Growing the ring also grows the neighbour exclusion, so on a crowded pattern fewer background pixels survive; a reflection left with too few is rejected outright. On the data this was measured on the loss is under 0.1% of reflections, but it is not structurally zero. -Only the ring moves. The signal disk $r_1$ stays circular, deliberately: it sets $n_S$, it sets $\mathrm{var}(\hat b)$, and it drives the all-or-nothing “every signal pixel valid” acceptance gate (§9.2), so growing it would reject any reflection carrying a single bad pixel anywhere along a long streak. What a circular $r_1$ loses is flux, and that loss is a function of resolution alone, which the per-shell scale absorbs. +Only the ring moves. The signal disk $r_1$ stays circular, deliberately: it sets $n_S$ and it sets $\mathrm{var}(\hat b)$, and with `--integrator boxsum` it drives the all-or-nothing “every signal pixel valid” acceptance gate (§9.2), so growing it would reject any box sum carrying a single bad pixel anywhere along a long streak. What a circular $r_1$ loses is flux, and that loss is a function of resolution alone, which the per-shell scale absorbs. ### 9.2 Box summation (seed and fallback) @@ -604,7 +604,7 @@ $ \hat{b} = \frac{B}{n_B},\qquad \hat{I} = S - n_S \hat{b}, $ -with a Poisson-like uncertainty $\sigma(\hat{I})=\max\!\big(1,\ r_\sigma\hat{I},\ \sqrt{S + n_S^2\,\mathrm{var}(\hat{b})}\big)$, i.e. $\sqrt{S}$ floored both at 1 and at a small fraction $r_\sigma$ of the intensity. The second term under the root is the **uncertainty of the background estimate itself**: $\hat b$ is measured from a finite number of ring pixels, $\mathrm{var}(\hat b)=\hat b/n_B$, and it is subtracted $n_S$ times over, so it enters squared. Omitting it understates $\sigma$ by $\sqrt{1+n_S/n_B}$ — 1.109 with the shipped circular stencil; with an elongated ring $n_B$ grows with resolution, so the factor is no longer one number for a run — uniformly, on every reflection of every dataset. The same term is carried into the profile fit (§9.3), where it adds $(\sum wP/\sum P^2/v)^2\,\mathrm{var}(\hat b)$; $n_B$ is the count of pixels behind the *final* background value, so a clip or trim that discards ring pixels raises it. A reflection is accepted as “observed” only if all signal pixels were valid and $n_B$ exceeds a minimum. This box sum is the classical estimator; it is used directly with `--integrator boxsum`, and otherwise seeds the profile fit below. +with a Poisson-like uncertainty $\sigma(\hat{I})=\max\!\big(1,\ r_\sigma\hat{I},\ \sqrt{S + n_S^2\,\mathrm{var}(\hat{b})}\big)$, i.e. $\sqrt{S}$ floored both at 1 and at a small fraction $r_\sigma$ of the intensity. The second term under the root is the **uncertainty of the background estimate itself**: $\hat b$ is measured from a finite number of ring pixels, $\mathrm{var}(\hat b)=\hat b/n_B$, and it is subtracted $n_S$ times over, so it enters squared. Omitting it understates $\sigma$ by $\sqrt{1+n_S/n_B}$ — 1.109 with the shipped circular stencil; with an elongated ring $n_B$ grows with resolution, so the factor is no longer one number for a run — uniformly, on every reflection of every dataset. The same term is carried into the profile fit (§9.3), where it adds $(\sum wP/\sum P^2/v)^2\,\mathrm{var}(\hat b)$; $n_B$ is the count of pixels behind the *final* background value, so a clip or trim that discards ring pixels raises it. A box sum is accepted as “observed” only if all signal pixels were valid and $n_B$ exceeds a minimum — it measures what is in the disk with no model of what should be there, so it cannot renormalise a disk it has lost pixels out of. The profile modes can, and do (§9.3). This box sum is the classical estimator; it is used directly with `--integrator boxsum`, and otherwise seeds the profile fit below, where $S$ and $n_S$ then count only the pixels that were actually read. **High-side clipped background (default on).** Because $\hat{I}=S-n_S\hat{b}$ is a small difference of large numbers for weak reflections, a per-pixel background bias $\delta\hat{b}$ becomes a *fractional* intensity bias $\approx n_S\,\delta\hat{b}/\hat{I}$ that grows as $\hat{I}$ shrinks — worst at the resolution edge. A plain ring mean reads high there, because neighbour-spot wings that survive the signal-disk mask, tails and zingers are one-sided (positive) contaminants. The ring mean is therefore made robust: pixels above $\hat{b}+n\sqrt{\hat{b}}$ are rejected and the mean recomputed, with $n=4$ (`--background-clip`; $n=0$ disables), lowered by `rugnux` to $n=3$ on broadband (non-zero bandwidth: pink-beam / DMM) data, where a bandwidth-streaked high-resolution spot leaks into the ring more readily. That is only a default — the flag sets $n$ whatever the bandwidth is. A clean Poisson ring is essentially unchanged by the cut (measured false-rejection rate 0.04–0.39 % at $4\sigma$), while a 40-pixel neighbour core at $+100$ counts shifts the estimate by $+0.009$ ct/px. @@ -638,6 +638,8 @@ v = \max\!\left(B + I\,P,\ \tfrac{1}{2}B\right), $ where $c$ is the pixel value and the de-biased variance $v$ (background plus model signal, rather than the down-fluctuating observed count) is iterated (a few passes). The plug-in $I$ enters **as it is**: half-wave rectifying it, $v=B+\max(I,0)P$, lets $v$ — and with it the reported $1/\sum P^2/v$ — respond only to *upward* fluctuations of a noisy estimate, which adds $\approx0.4\,\sigma\sum P^3/(\sum P^2)^2$ to every $\sigma$ whatever the count rate. That offset is invisible on strong reflections and a large fractional inflation on weak ones; the $\tfrac12 B$ clamp keeps $v$ positive without reintroducing it. As a guard, if the profile intensity runs away from the box-sum seed (by more than ~10 box-sum $\sigma$) it falls back to the seed, and the background term is floored at $0.01$ ct/px — enough to keep $P^2/v$ finite when the ring mean reads exactly zero, which a ring of $n_B$ pixels cannot distinguish from any background below $\approx1/n_B$. The rotation/excitation partiality is carried exactly as in the box-sum path. +**Pixels the fit cannot use (MINPK).** A profile fit is the amplitude of a *normalised* profile, so a pixel left out of the sum renormalises the estimator by construction: it costs information — $\sum P^2/v$ shrinks and $\sigma$ grows — but biases nothing. That is what keeps a reflection whose signal disk is cut by a mask, an untrusted region, a detector gap or an overload: those pixels are simply not read, and the fit is taken over the rest, exactly as the shared pixels of a crowded reflection are (`--overlap exclude`). The reflection is kept only while enough of the expected profile survives — at least `--overlap-minpk` of the profile mass that falls on the detector at all, default 0.75, which is XDS's `MINPK` and dials' `valid_foreground_threshold`. The complete reflections alone teach the profile, its resolution shells and their widths. `--integrator boxsum` has no profile to renormalise with and keeps the all-or-nothing rule of §9.2. + The integrator is selected by `--integrator boxsum|gaussian|empirical` (default `gaussian`). ### 9.4 Lorentz–polarization factor handling diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 84c5bb36..18dc85ae 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -128,6 +128,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, int k_bin = 0; // which radial-background kernel this reflection's stencil needs int cx = 0, cy = 0, shell = -1; bool ok = false, strong = false, has_obs = false; + bool full = false; // every pixel of the signal disk was readable }; std::vector rough(npredicted); double inv_d2_min = std::numeric_limits::max(), inv_d2_max = 0.0; @@ -196,7 +197,15 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, } int n_bkg_used = n_bkg; // pixels behind the FINAL background value (trim/clip shrink it) - if (n_inner_valid == n_inner && n_bkg > 5) { + // A masked, untrusted, gapped or overloaded pixel inside the signal disk used to discard the + // reflection outright. A profile fit does not need it to: the fit is the amplitude of a + // NORMALISED profile, so leaving pixels out renormalises the estimator by construction and + // widens sigma by the information they carried - the same argument Exclude makes for a + // neighbour's pixels. Pass B keeps the reflection only while enough of the expected profile + // survived to constrain the amplitude (XDS's MINPK, dials' valid_foreground_threshold). A box + // sum has no profile to renormalise with, so there it stays all or nothing. + const bool full = n_inner_valid == n_inner; + if ((full || mode != IntegratorMode::BoxSum) && n_bkg > 5) { out.bkg = bkg_sum / n_bkg; if (bkg_trim_frac > 0.0 && bkg_vals.size() > 5 && bkg_vals.size() <= static_cast(bragg_engine::BKG_TRIM_MAX)) { @@ -239,31 +248,39 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, } if (n > 5) { out.bkg = s / n; n_bkg_used = n; } } - out.I = static_cast(I_sum) - static_cast(n_inner) * out.bkg; + // The sum is over the pixels actually READ, so that is the count the background is + // subtracted with; with nothing missing it is the whole disk, exactly as before. + out.I = static_cast(I_sum) - static_cast(n_inner_valid) * out.bkg; // I = I_sum - n_inner*bkg, and bkg is itself estimated from n_bkg_used pixels, so its // error enters n_inner times over: var(I) = I_sum + n_inner^2 * bkg/n_bkg_used. Leaving // the second term out understates sigma by sqrt(1 + n_inner/n_bkg) - 1.109x at the // default r1=4/r2=6/r3=10 stencil, on every reflection of every dataset. out.bkg_var = out.bkg / n_bkg_used; - out.var_bkg = static_cast(n_inner) * out.bkg - + static_cast(n_inner) * n_inner * out.bkg_var; + out.var_bkg = static_cast(n_inner_valid) * out.bkg + + static_cast(n_inner_valid) * n_inner_valid * out.bkg_var; out.I_sum = I_sum; - out.n_inner = static_cast(n_inner); + out.n_inner = static_cast(n_inner_valid); out.n_disk = n_disk; out.n_own = n_own; - const double var_bkg_term = static_cast(n_inner) * n_inner * out.bkg_var; + const double var_bkg_term = static_cast(n_inner_valid) * n_inner_valid * out.bkg_var; out.sigma = 1.0; if (I_sum > 0) { out.sigma = std::max(out.sigma, std::sqrt(static_cast(I_sum) + var_bkg_term)); out.obs_x = static_cast(I_sum_x) / static_cast(I_sum); out.obs_y = static_cast(I_sum_y) / static_cast(I_sum); - out.has_obs = true; + // A disk with a hole in it gives a centroid pulled away from the hole, and the hole + // sits at a fixed place on the detector - post-refinement would read that as geometry. + out.has_obs = full; } out.cx = static_cast(std::lround(r.predicted_x)); out.cy = static_cast(std::lround(r.predicted_y)); out.ok = true; - out.strong = out.sigma > 0.0 && out.I / out.sigma >= STRONG_I_OVER_SIGMA; - if (r.d > 0.0f) { + out.full = full; + // The profile, its resolution shells and their widths are learned from COMPLETE + // reflections: a partial one has no intensity to normalise its grid by and a hole to + // stack into it. So the learning sees exactly what it saw before this rescue existed. + out.strong = full && out.sigma > 0.0 && out.I / out.sigma >= STRONG_I_OVER_SIGMA; + if (full && r.d > 0.0f) { const double inv_d2 = 1.0 / (static_cast(r.d) * r.d); inv_d2_min = std::min(inv_d2_min, inv_d2); inv_d2_max = std::max(inv_d2_max, inv_d2); @@ -454,29 +471,30 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, const int Gf = 2 * Rf + 1; - // --- How much of the expected profile is cleanly this reflection's own. p_own is the whole - // grid's clean mass, i.e. XDS's MINPK quantity, and it is what Reject cuts on. m_own / - // m_all is the same fraction over the r1 disk alone, which is what the summation seed the - // runaway guard below compares against actually saw; with nothing excluded it is 1 and - // the guard is untouched. --- - double p_own = 1.0, m_all = 0.0, m_own = 0.0; - if (overlap != OverlapMode::Off) { - p_own = 0.0; - for (int dy = -Rf; dy <= Rf; ++dy) - for (int dx = -Rf; dx <= Rf; ++dx) { - const double Pp = (*Pvec)[(dy + Rf) * Gf + (dx + Rf)]; - if (Pp <= 0.0) continue; - const int x = rh.cx + dx, y = rh.cy + dy; - if (x < 0 || y < 0 || x >= W || y >= H) continue; - if (!valid(img[y * W + x])) continue; - const bool own = clean(x, y, i); - if (own) p_own += Pp; - if (dx * dx + dy * dy < r1_sq) { - m_all += Pp; - if (own) m_own += Pp; - } - } - } + // --- How much of the expected profile the fit can actually see. p_valid is the readable + // fraction of it - XDS's MINPK, dials' valid_foreground_threshold - measured against the + // mass that falls on the detector at all, so a reflection is judged on the pixels that + // exist and not on where the sensor ends. p_own is the same quantity over neighbour-owned + // pixels, what Reject cuts on. m_read / m_all is the fraction of the r1 disk the + // summation seed the runaway guard compares against actually summed; with nothing missing + // and nothing excluded it is 1 and the guard is untouched. --- + double p_grid = 0.0, p_valid = 0.0, p_own = 0.0, m_all = 0.0, m_read = 0.0; + for (int dy = -Rf; dy <= Rf; ++dy) + for (int dx = -Rf; dx <= Rf; ++dx) { + const double Pp = (*Pvec)[(dy + Rf) * Gf + (dx + Rf)]; + if (Pp <= 0.0) continue; + const int x = rh.cx + dx, y = rh.cy + dy; + if (x < 0 || y < 0 || x >= W || y >= H) continue; + const bool in_disk = dx * dx + dy * dy < r1_sq; + p_grid += Pp; + if (in_disk) m_all += Pp; + if (!valid(img[y * W + x])) continue; + p_valid += Pp; + const bool own = clean(x, y, i); + if (own) p_own += Pp; + if (in_disk && (own || !exclude)) m_read += Pp; + } + if (p_valid < overlap_min_peak * p_grid) continue; if (overlap == OverlapMode::Reject && p_own < overlap_min_peak) continue; const double B = std::max(rh.bkg, PIXEL_VARIANCE_FLOOR); @@ -519,10 +537,10 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, // fluctuation, and 1/den - I is what does not (1/den grows with I one for one). double var_bkg = std::max(0.0, 1.0 / den - I + (wsum / den) * (wsum / den) * rh.bkg_var); - // The seed is a sum over the disk the box sum actually read, so when Exclude has taken pixels - // out of both, the fit's full-profile intensity has to be scaled down to that same disk - // before the two are comparable. Nothing excluded gives exactly 1. - const double guard_scale = exclude && m_all > 0.0 ? m_own / m_all : 1.0; + // The seed is a sum over the disk the box sum actually read, so when pixels are missing from + // both - excluded to a neighbour, or unreadable - the fit's full-profile intensity has to be + // scaled down to that same disk before the two are comparable. Nothing dropped gives 1. + const double guard_scale = m_all > 0.0 ? m_read / m_all : 1.0; if (std::abs(I * guard_scale - rh.I) > PROFILE_SUMMATION_MAX_NSIGMA * rh.sigma) { I = rh.I; sigma = rh.sigma; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 986c272b..178475c9 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -32,7 +32,8 @@ struct BraggGpuParams { int boxsum_reject; // BoxSum mode under Reject: drop on the disk-AREA fraction (see the CPU engine) int exclude_px; // drop a neighbour's pixels from the disk: Exclude, and not a box sum float claim_sq, inv_claim; // how far a reflection claims pixels in the owner map - float minpk; // Reject: least clean profile fraction that is kept + float minpk; // least readable/clean profile fraction that is kept (XDS MINPK) + int partial_ok; // keep a signal disk with unreadable pixels: anything but a box sum }; __device__ inline bool valid(int32_t v) { return v != INT32_MIN && v != INT32_MAX; } @@ -82,7 +83,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, __shared__ unsigned long long s_Isum, s_Ix, s_Iy; __shared__ int s_ninner, s_ninner_valid, s_nbkg, s_ndisk, s_nown; __shared__ double s_bkgsum; - __shared__ int s_accept; + __shared__ int s_accept, s_full; __shared__ double s_bkg, s_thr, s_clipsum; __shared__ int s_clipn; __shared__ float s_r0; @@ -149,7 +150,11 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, for (int t = threadIdx.x; t < p.rad_w; t += blockDim.x) { s_radv[t] = 0.0f; s_radn[t] = 0; } if (threadIdx.x == 0) s_radbase = (int) lroundf(s_r0) - p.rad_w / 2; if (threadIdx.x == 0) { - s_accept = (s_ninner_valid == s_ninner && s_nbkg > 5) ? 1 : 0; + // A hole in the signal disk no longer discards the reflection in the profile modes - the fit + // renormalises to the pixels it can read and Pass B cuts on how much of the profile survived + // (XDS's MINPK). A box sum has no profile to renormalise with. See the CPU engine. + s_full = (s_ninner_valid == s_ninner) ? 1 : 0; + s_accept = ((s_full || p.partial_ok) && s_nbkg > 5) ? 1 : 0; s_bkg = s_accept ? (s_bkgsum / (double) s_nbkg) : 0.0; s_thr = s_bkg + (double) p.bkg_clip_nsigma * sqrt(fmax(s_bkg, 1.0)); s_clipsum = 0.0; s_clipn = 0; @@ -258,34 +263,39 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, } if (p.bkg_clip_nsigma > 0.0f && s_clipn > 5) { bkg = s_clipsum / (double) s_clipn; n_bkg_used = s_clipn; } const long long Isum = (long long) s_Isum; - const double I = (double) Isum - (double) s_ninner * bkg; + // The sum is over the pixels actually READ, so that is the count the background is subtracted + // with; with nothing missing it is the whole disk, exactly as before. + const double I = (double) Isum - (double) s_ninner_valid * bkg; // See the CPU engine: bkg is estimated from n_bkg_used pixels and subtracted n_inner times, so // var(I) = Isum + n_inner^2 * bkg/n_bkg_used. Both engines must agree. const double bkg_var = bkg / (double) n_bkg_used; - const double var_bkg_term = (double) s_ninner * (double) s_ninner * bkg_var; + const double var_bkg_term = (double) s_ninner_valid * (double) s_ninner_valid * bkg_var; double sigma = 1.0; uint8_t hasobs = 0; double ox = 0.0, oy = 0.0; if (Isum > 0) { sigma = fmax(sigma, sqrt((double) Isum + var_bkg_term)); ox = (double) (long long) s_Ix / (double) Isum; oy = (double) (long long) s_Iy / (double) Isum; - hasobs = 1; + // A disk with a hole gives a centroid pulled away from it; see the CPU engine. + hasobs = s_full ? 1 : 0; } cx_o[i] = (int) lroundf(cx); cy_o[i] = (int) lroundf(cy); I_o[i] = (float) I; sigma_o[i] = (float) sigma; bkg_o[i] = (float) bkg; bkgvar_o[i] = (float) bkg_var; - varbkg_o[i] = (float) ((double) s_ninner * bkg + var_bkg_term); + varbkg_o[i] = (float) ((double) s_ninner_valid * bkg + var_bkg_term); isum_o[i] = (float) Isum; - ninner_o[i] = s_ninner; + ninner_o[i] = s_ninner_valid; rbin_o[i] = min(max((int) lroundf(s_r0), 0), n_rad > 0 ? n_rad - 1 : 0); kbin_o[i] = BraggStencilKernelIndex(st, p.n_kern); obsx_o[i] = (float) ox; obsy_o[i] = (float) oy; hasobs_o[i] = hasobs; ok_o[i] = 1; - strong_o[i] = (sigma > 0.0 && I / sigma >= STRONG_I_OVER_SIGMA) ? 1 : 0; + // The profile and its resolution shells are learned from COMPLETE reflections only; see the CPU + // engine. + strong_o[i] = (s_full && sigma > 0.0 && I / sigma >= STRONG_I_OVER_SIGMA) ? 1 : 0; const float d = dd[i]; - if (d > 0.0f) { + if (s_full && d > 0.0f) { // Positive doubles keep IEEE bit-pattern ordering, so atomicMin/Max on the ull view works. const unsigned long long b = (unsigned long long) __double_as_longlong(1.0 / ((double) d * d)); atomicMin(&invd2mm[0], b); @@ -456,7 +466,7 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x if (i >= n) return; extern __shared__ float Pbuf[]; __shared__ float s_gs, s_num, s_den, s_I, s_wsum; - __shared__ float s_pown, s_mall, s_mown; + __shared__ float s_pgrid, s_pvalid, s_pown, s_mall, s_mread; __shared__ int s_Rf, s_Gf; if (!ok_a[i]) { if (threadIdx.x == 0) ok_o[i] = 0; return; } @@ -509,37 +519,45 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x const int Rf = s_Rf, Gf = s_Gf, GfGf = Gf * Gf; const float B = fmaxf(bkg, (float) PIXEL_VARIANCE_FLOOR); - // How much of the expected profile is cleanly this reflection's own: s_pown over the whole grid - // (XDS's MINPK quantity, what Reject cuts on) and s_mown / s_mall over the r1 disk alone, which - // is what the summation seed the runaway guard compares against actually saw. See the CPU engine. - if (p.overlap) { - if (threadIdx.x == 0) { s_pown = 0.0f; s_mall = 0.0f; s_mown = 0.0f; } - __syncthreads(); - float l_pown = 0.0f, l_mall = 0.0f, l_mown = 0.0f; - for (int k = threadIdx.x; k < GfGf; k += blockDim.x) { - const float Pp = Pbuf[k]; - if (Pp <= 0.0f) continue; - const int dx = k % Gf - Rf, dy = k / Gf - Rf; - const int x = cx + dx, y = cy + dy; - if (x < 0 || y < 0 || x >= p.W || y >= p.H) continue; - if (!valid(img[y * p.W + x])) continue; - const bool own = BraggOwnedBy(owner[y * p.W + x], i); - // Zeroing the profile here is how Exclude drops the pixel: the fit skips any cell with - // P <= 0 already, and P is not renormalised, so the fitted amplitude comes out on the - // scale of the WHOLE profile - the renormalisation is the estimator's own doing. - if (!own && p.overlap == 2) Pbuf[k] = 0.0f; - if (own) l_pown += Pp; - if ((float) (dx * dx + dy * dy) < p.r1_sq) { - l_mall += Pp; - if (own) l_mown += Pp; - } - } - atomicAdd(&s_pown, l_pown); atomicAdd(&s_mall, l_mall); atomicAdd(&s_mown, l_mown); - __syncthreads(); - if (p.overlap == 1 && s_pown < p.minpk) { - if (threadIdx.x == 0) ok_o[i] = 0; - return; - } + // How much of the expected profile the fit can actually see: s_pvalid over the whole grid against + // s_pgrid, the mass that falls on the detector at all (XDS's MINPK), s_pown the same over + // neighbour-owned pixels (what Reject cuts on), and s_mread / s_mall over the r1 disk alone, + // which is what the summation seed the runaway guard compares against actually summed. See the + // CPU engine. + if (threadIdx.x == 0) { + s_pgrid = 0.0f; s_pvalid = 0.0f; s_pown = 0.0f; s_mall = 0.0f; s_mread = 0.0f; + } + __syncthreads(); + float l_pgrid = 0.0f, l_pvalid = 0.0f, l_pown = 0.0f, l_mall = 0.0f, l_mread = 0.0f; + for (int k = threadIdx.x; k < GfGf; k += blockDim.x) { + const float Pp = Pbuf[k]; + if (Pp <= 0.0f) continue; + const int dx = k % Gf - Rf, dy = k / Gf - Rf; + const int x = cx + dx, y = cy + dy; + if (x < 0 || y < 0 || x >= p.W || y >= p.H) continue; + const bool in_disk = (float) (dx * dx + dy * dy) < p.r1_sq; + l_pgrid += Pp; + if (in_disk) l_mall += Pp; + if (!valid(img[y * p.W + x])) continue; + l_pvalid += Pp; + const bool own = !p.overlap || BraggOwnedBy(owner[y * p.W + x], i); + // Zeroing the profile here is how Exclude drops the pixel: the fit skips any cell with + // P <= 0 already, and P is not renormalised, so the fitted amplitude comes out on the + // scale of the WHOLE profile - the renormalisation is the estimator's own doing. + if (!own && p.overlap == 2) Pbuf[k] = 0.0f; + if (own) l_pown += Pp; + if (in_disk && (own || p.overlap != 2)) l_mread += Pp; + } + atomicAdd(&s_pgrid, l_pgrid); atomicAdd(&s_pvalid, l_pvalid); atomicAdd(&s_pown, l_pown); + atomicAdd(&s_mall, l_mall); atomicAdd(&s_mread, l_mread); + __syncthreads(); + if (s_pvalid < p.minpk * s_pgrid) { + if (threadIdx.x == 0) ok_o[i] = 0; + return; + } + if (p.overlap == 1 && s_pown < p.minpk) { + if (threadIdx.x == 0) ok_o[i] = 0; + return; } if (threadIdx.x == 0) s_I = I_seed[i]; @@ -577,10 +595,10 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x // The signal part to remove is I itself, not max(0, I) - see the CPU engine. float var_bkg = fmaxf(0.0f, 1.0f / s_den - I + wr * wr * bkgvar_a[i]); // Guard against profile-fit runaways (see the CPU engine): fall back to the summation seed - // when the profile result diverges from it. Exclude has taken pixels out of both, so the - // full-profile intensity is scaled back to the disk the seed read; nothing excluded gives - // exactly 1. - const float gs = (p.overlap == 2 && s_mall > 0.0f) ? s_mown / s_mall : 1.0f; + // when the profile result diverges from it. Pixels missing from both - excluded to a + // neighbour, or unreadable - scale the full-profile intensity back to the disk the seed + // read; nothing dropped gives exactly 1. + const float gs = s_mall > 0.0f ? s_mread / s_mall : 1.0f; if (fabsf(I * gs - I_seed[i]) > (float) PROFILE_SUMMATION_MAX_NSIGMA * sigma_seed[i]) { I = I_seed[i]; sigma = sigma_seed[i]; @@ -726,6 +744,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu .boxsum_reject = (overlap == OverlapMode::Reject && mode == IntegratorMode::BoxSum) ? 1 : 0, .exclude_px = (overlap == OverlapMode::Exclude && mode != IntegratorMode::BoxSum) ? 1 : 0, .claim_sq = claim * claim, .inv_claim = inv_claim, .minpk = overlap_min_peak, + .partial_ok = mode != IntegratorMode::BoxSum ? 1 : 0, }; // Whether the radial correction runs for THIS image. n_rad only says the buffers exist - under diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 253da9a6..9f3b1323 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -169,7 +169,7 @@ void print_usage() { std::cout << " --background-radial[=on|off|auto] Correct the background ring for the CURVATURE of the radial background (default off). The signal disk and the background ring are concentric, so a background linear in position cancels between them and only curvature survives - which on a smooth ice ring reaches +26 counts on a single reflection. =auto applies it per image where that image's ice score shows a smooth powder ring, which is where a radius-only background model holds; on ice made of discrete crystallite spots there is no such ring and the correction makes the bias worse. Costs one short dot product per reflection and no extra pixel reads" << std::endl; std::cout << " --background-trim Use the old symmetric trimmed mean for the background ring instead of the clip (0<=f<0.5; 0.10 was the former default). Switches --background-clip off. A symmetric trim is biased low on Poisson data and adds ~5 counts to every partial, so this is for back compatibility only; 0 = plain ring mean. Applies whatever --bandwidth is set to" << std::endl; std::cout << " --overlap What to do where two predicted reflections share signal pixels: off|reject|exclude (default exclude). A pixel inside two signal disks belongs to the nearer centre; before this, nothing kept a neighbour's flux out of a reflection's own disk, so on a dense pattern a crowded reflection read high. exclude drops the shared PIXELS from the profile fit, which renormalises itself, so the reflection is kept; reject instead drops the whole reflection when less than --overlap-minpk of its expected profile is cleanly its own (what XDS calls MINPK). --integrator boxsum has no profile to renormalise with, so exclude does nothing there and only reject acts" << std::endl; - std::cout << " --overlap-minpk Least fraction of a reflection's expected profile that must be cleanly its own for --overlap reject to keep it (default 0.75, XDS MINPK). With --integrator boxsum the fraction is by disk AREA instead, which cuts harder" << std::endl; + std::cout << " --overlap-minpk Least fraction of a reflection's expected profile that must be usable for the reflection to be kept (default 0.75, XDS MINPK). Governs both ways part of a profile is lost: the fraction that must be READABLE - not masked, untrusted, in a detector gap or overloaded - in every profile mode, and, under --overlap reject, the fraction that must be cleanly the reflection's own. --integrator boxsum has no profile to renormalise with, so there a signal disk with any unreadable pixel is still discarded outright and the reject fraction is by disk AREA, which cuts harder" << std::endl; std::cout << " --integrator Spot integrator boxsum|gaussian|empirical (default: gaussian profile-fit; boxsum is the classical fallback)" << std::endl; std::cout << " --simple-stills stills: treat every reflection as a full (p=1, single-pass scale/merge); disables the default physical partiality post-refinement" << std::endl; std::cout << " -q, --azim-q-spacing Azimuthal-integration Q bin spacing (1/A) (default: 0.01)" << std::endl;