mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 67 Wiki Activity

Bragg integration: fit a reflection over the pixels it has, not only over all of them

Browse Source 
A predicted reflection was discarded outright if ANY pixel of its signal disk was
unreadable - masked, untrusted, in a detector gap, or overloaded. On a battery
crystal that is 11.1% of all predictions, thrown away for a defect in one pixel
of fifty, and the pixels concerned sit at fixed places on the detector, so the
loss is systematic in reciprocal space rather than random.

Neither XDS nor dials does that. Both estimate the missing part from the profile
instead and keep the reflection while enough of it was seen: XDS's MINPK (default
75%, "the missing intensity is estimated from the learned profiles"), dials'
integration.profile.valid_foreground_threshold (default 0.75). MOSFLM is the one
program that rejects by default, and even it relaxes to 50% with PROFILE EDGE.

We already had the argument and the machinery: a profile fit is the amplitude of
a NORMALISED profile, so leaving pixels out renormalises the estimator by
construction - it costs information, which sum P^2/v duly loses and sigma duly
gains, and biases nothing. That is exactly why --overlap exclude drops a
neighbour's pixels from the fit rather than the reflection. Unreadable pixels are
the same case with a different reason, so they take the same treatment, cut on
the same threshold, in the same place: the readable fraction of the expected
profile, measured against the profile mass that lands on the detector at all so a
reflection is judged on the pixels that exist. A box sum has no profile to
renormalise with and keeps the all-or-nothing rule.

Two consequences handled. The summation seed and its variance now count the
pixels actually read, and the runaway guard scales the fit back to that same disk
before comparing - both exactly as before wherever nothing is missing. (Its
fallback then hands back that partial sum unrescaled, which would read low; the
guard fires on 8 of 96 688 recovered reflections, and on none at all on a weak
crystal, so it is not worth a branch.) And the profile, its resolution shells and
their widths are learned from COMPLETE reflections only, as is the box-sum
centroid post-refinement reads as an observed position: a disk with a hole gives
a centroid pulled away from the hole, and the hole does not move between frames.

That sigma gains what the missing pixels carried is the claim the whole change
rests on, and it is measurable. Force the conventional CENTRED cell of a
body-centred crystal in P1: the predictor then enumerates every lattice point,
and the reflections the centring makes systematically absent have a true
intensity of exactly zero, so their scatter about zero must equal their reported
sigma. Over 7.1 M such observations, matched by resolution shell, the trimmed
std(I)/rms(sigma) of the recovered reflections is 0.99 / 1.20 / 2.33 / 1.04 /
1.69 against 0.98 / 1.22 / 2.29 / 1.03 / 1.56 for the reflections that were
complete - the same calibration to a few percent. The lever there is small,
because the typical recovered reflection is missing only 5% of its disk. Lowering
the threshold to 0.50 admits a band missing 25-50%, which is a real lever: there
sigma comes out 8-43% larger than a complete reflection's in the same shell, and
the scatter about zero tracks it, 0.97 / 1.09 / 1.92 / 0.99 / 1.37, at or below
the complete population. Sigma grows, and by the amount it should.

The threshold stays at XDS's and dials' 0.75, on that evidence and on quality.
Below it the estimator starts to run out: on those same zero-intensity
reflections the recovered ones read +0.8 counts high at 0.75 and +1.9 counts high
in the 0.50-0.75 band, against a sigma of 12-17, and at 0.25 the fit degenerates
outright, single reflections carrying sigma in the thousands. Above it there is
nothing to buy: 0.90 leaves a fifth of the recoverable observations behind and
measures no better for them. On the high-multiplicity control, R_rim over
as-shipped / 0.90 / 0.75 / 0.50 runs 4.49% / 4.51% / 4.56% / 4.78% while
<I/sigma> runs 33.47 / 34.02 / 33.89 / 33.43 - 0.50 is where the recovered
observations stop paying for themselves.

Probe against the previous commit, six crystals. The high-multiplicity control
gains 4.2% more observations, 924 803 -> 963 946, which lands it on XDS's 961 379
from the same images, for <I/sigma> 33.47 -> 33.89, R_rim 4.49% -> 4.56% at 4.3%
more multiplicity, CC1/2 unchanged at 0.9998 and ISa 27.80 -> 27.12. Five weaker
crystals gain 3.3-4.8% of their observations and up to 1.0 point of completeness,
for <I/sigma> +0.4 to +3.6%, R_rim between -8.1% and +5.8% relative, CC1/2
+6.6 / +0.3 / +0.2 / -0.0 / -1.2 points, and ISa between +0.3% and -3.4%. Some of
that ISa is the point rather than the price: a reflection integrated over fewer
pixels carries less information, and the absence test above says the sigma that
reports so is honest. The GPU and CPU engines agree as before.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

Full 38-crystal rotation battery, against the same binary without it:

  observations     better 38 / worse  0,  +937 100
  unique refl      better 30 / worse  0,    +9 229
  overall <I/sig>  better 33 / worse  1,     +7.00
  CC1/2            better  5 / worse  1,      +6.2
  space groups     unchanged at 35/38

Every crystal gains observations and not one loses a unique reflection. The two
costs are small and both are understood. Low-resolution R_meas is worse on eight
crystals, by +0.8 pp at most and +3.2 pp summed - a reflection whose own peak
pixel is unreadable loses the part of the profile that carries most of the
amplitude, and that population sits at low resolution; the following commit
handles it. And ISa falls on 32 crystals, by 10.9 summed, which is what admitting
937 000 further observations does to the strong-reflection asymptote: R_meas
excluding the one crystal whose thread-count noise is 1.5 pp is flat.
This commit is contained in:
leonarski_f
2026-08-12 02:03:42 +02:00
parent e0e429ae2f
commit 12ae6c5228
6 changed files with 136 additions and 90 deletions
Download Patch File Download Diff File Expand all files Collapse all files
common/BraggIntegrationSettings.h
+11 -5
View File
@@ -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<int> 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;
docs/CHANGELOG.md
+1
View File
@@ -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.**
docs/CPU_DATA_ANALYSIS.md
+4 -2
View File
@@ -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
image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp
+54 -36
View File
@@ -128,6 +128,7 @@ std::vector<Reflection> 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> rough(npredicted);
double inv_d2_min = std::numeric_limits<double>::max(), inv_d2_max = 0.0;
@@ -196,7 +197,15 @@ std::vector<Reflection> 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<size_t>(bragg_engine::BKG_TRIM_MAX)) {
@@ -239,31 +248,39 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
}
if (n > 5) { out.bkg = s / n; n_bkg_used = n; }
}
out.I = static_cast<double>(I_sum) - static_cast<double>(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<double>(I_sum) - static_cast<double>(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<double>(n_inner) * out.bkg
+ static_cast<double>(n_inner) * n_inner * out.bkg_var;
out.var_bkg = static_cast<double>(n_inner_valid) * out.bkg
+ static_cast<double>(n_inner_valid) * n_inner_valid * out.bkg_var;
out.I_sum = I_sum;
out.n_inner = static_cast<int>(n_inner);
out.n_inner = static_cast<int>(n_inner_valid);
out.n_disk = n_disk;
out.n_own = n_own;
const double var_bkg_term = static_cast<double>(n_inner) * n_inner * out.bkg_var;
const double var_bkg_term = static_cast<double>(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<double>(I_sum) + var_bkg_term));
out.obs_x = static_cast<double>(I_sum_x) / static_cast<double>(I_sum);
out.obs_y = static_cast<double>(I_sum_y) / static_cast<double>(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<int>(std::lround(r.predicted_x));
out.cy = static_cast<int>(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<double>(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<Reflection> 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<Reflection> 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;
image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu
+65 -46
View File
@@ -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<Reflection> 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
rugnux/rugnux_cli.cpp
+1 -1
View File
@@ -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 <f> 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 <txt> 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 <f> 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 <f> 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 <txt> 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 <num> Azimuthal-integration Q bin spacing (1/A) (default: 0.01)" << std::endl;