Bragg integration: elongate the background ring per reflection
The signal disk and the r2..r3 background ring were fixed pixel circles, identical for every reflection at every resolution. A reflection is not round: a finite bandwidth streaks it radially by bw_sigma*Rpx, so at high resolution the ring sits within 1.3-2.2 sigma of the reflection's own profile and measures its tails as background. --integration-stencil <k> makes the RING an ellipse, elongated along the beam->reflection direction by k times that streak, capped at 2*r3. The tangential half-widths stay r2 and r3, and the r1 signal disk stays a circle: r1 drives the all-or-nothing n_inner_valid == n_inner gate, so growing it rejects any reflection carrying one bad pixel along a long streak, and the flux a circular r1 loses is a function of resolution alone, which the per-shell scale absorbs. The geometry lives in one shared header compiled by both the host compiler and nvcc, so the seven pixel-classification sites - the CPU mask/main/clip loops and the GPU mark_mask/main/trim/clip kernels - cannot drift apart. Rather than evaluate an ellipse, each pixel's squared distance has its radial part scaled down, d2 - q*rad^2 against r2^2/r3^2 with q = 1 - (r/(r+grow))^2, so grow = 0 gives q = 0 and both tests collapse onto d2 exactly in floating point. The width is the bandwidth streak alone, not the profile's full radial variance, which also carries the sensor parallax and weak-spot capture terms. Deriving the growth from those was implemented first and measured on the rotation battery: at k=1 it took Thau_9's high-shell CC1/2 from 75.8 to 27.9 and Benas_3's from 14.1 to 6.0, against cytC_10 +1.2 and lyso_ref flat. On a monochromatic beam they are the only terms there are, and C_CAPTURE is 64% of them. Keeping only the streak also makes the option exactly inert without a bandwidth, rather than merely small. Default 0. Measured on broadband rotation data with the bandwidth set to its spectroscopic value, matched resolution limits: high-shell CC1/2 30.6 -> 46.4 at k=4, and better in EVERY shell in both CC1/2 and R_meas (top shell R_meas 194.7% -> 138.7%), with completeness, multiplicity and space group unchanged and 28 of 98833 unique reflections lost. Anomalous peak height over 18 sites +0.107 +- 0.039 sigma (p = 0.013). The full 38-crystal rotation battery is unchanged to every reported digit, base against k=3. Two consequences of an elongated ring are handled rather than inherited. The neighbour exclusion marks the inner ELLIPSE in each neighbour's own frame, or an elongated neighbour leaks its tails into this reflection's ring. And the radial-background curvature kernel becomes a small table indexed by the growth, because its azimuthal average makes one kernel serve every reflection only while their stencils are identical; the GPU's radial window, previously a fixed 32 bins, is now sized on the host from the widest ring on the detector. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
Claude Opus 5
@@ -44,6 +44,14 @@ BraggIntegrationSettings &BraggIntegrationSettings::R3(float input) {
|
||||
return *this;
|
||||
}
|
||||
|
||||
BraggIntegrationSettings &BraggIntegrationSettings::StencilKSigma(float input) {
|
||||
check_finite("Integration stencil elongation", input);
|
||||
check_min("Integration stencil elongation", input, 0.0);
|
||||
check_max("Integration stencil elongation", input, 10.0);
|
||||
stencil_k_sigma = input;
|
||||
return *this;
|
||||
}
|
||||
|
||||
BraggIntegrationSettings &BraggIntegrationSettings::DMinLimit_A(std::optional<float> input) {
|
||||
if (input) {
|
||||
check_finite("Minimum d-spacing", *input);
|
||||
@@ -90,6 +98,10 @@ float BraggIntegrationSettings::GetR3() const {
|
||||
return r_3;
|
||||
}
|
||||
|
||||
float BraggIntegrationSettings::GetStencilKSigma() const {
|
||||
return stencil_k_sigma;
|
||||
}
|
||||
|
||||
std::optional<float> BraggIntegrationSettings::GetDMinLimit_A() const {
|
||||
return d_min_limit_A;
|
||||
}
|
||||
|
||||
@@ -21,6 +21,14 @@ class BraggIntegrationSettings {
|
||||
float r_1 = 4;
|
||||
float r_2 = 6;
|
||||
float r_3 = 10;
|
||||
// How many times the beam's radial streak to push the r2..r3 background ring out by, per
|
||||
// reflection. A bandwidth streaks a spot radially by bw_sigma*Rpx, and against a fixed pixel ring
|
||||
// that puts the background annulus on the reflection's own tails at high resolution, where it
|
||||
// measures signal as background. The ring's radial semi-axes become r2 + this*bw_sigma*Rpx and
|
||||
// r3 + this*bw_sigma*Rpx, the tangential ones stay r2 and r3, and the r1 signal disk stays a
|
||||
// circle (growing it trips the all-or-nothing n_inner_valid gate). 0 reproduces the fixed
|
||||
// circular stencil exactly, and so does any monochromatic beam, where the streak is zero.
|
||||
float stencil_k_sigma = 0.0f;
|
||||
// Integration/prediction resolution limit. Unset means "as far as the detector reaches", resolved
|
||||
// from the geometry where it is used. The predictor independently rejects any reflection that misses
|
||||
// the detector, so this is a bound on how far the lattice walk goes rather than a second opinion on
|
||||
@@ -92,6 +100,7 @@ public:
|
||||
BraggIntegrationSettings& R1(float input);
|
||||
BraggIntegrationSettings& R2(float input);
|
||||
BraggIntegrationSettings& R3(float input);
|
||||
BraggIntegrationSettings& StencilKSigma(float input);
|
||||
BraggIntegrationSettings& DMinLimit_A(std::optional<float> input);
|
||||
BraggIntegrationSettings& FixedProfileRadius_recipA(std::optional<float> input);
|
||||
BraggIntegrationSettings& Integrator(IntegratorMode input);
|
||||
@@ -105,6 +114,7 @@ public:
|
||||
[[nodiscard]] float GetR1() const;
|
||||
[[nodiscard]] float GetR2() const;
|
||||
[[nodiscard]] float GetR3() const;
|
||||
[[nodiscard]] float GetStencilKSigma() const;
|
||||
[[nodiscard]] std::optional<float> GetFixedProfileRadius_recipA() const;
|
||||
[[nodiscard]] std::optional<float> GetDMinLimit_A() const;
|
||||
|
||||
|
||||
@@ -3,6 +3,7 @@
|
||||
### 1.0.0-rc.161
|
||||
This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use.
|
||||
|
||||
* rugnux: The background ring can be **elongated radially per reflection** for broadband data (`--integration-stencil <k>`, default 0 = the fixed circular ring), by `k` times the beam's radial streak; the `r1` signal box stays circular.
|
||||
* rugnux: New **beam-stop shadow detection**, **on by default** (`--detect-beam-stop[=N|off]`), finds the beam stop and its holder in a projection of N images (default 60) and adds them to the pixel mask as bit 9, which is cleared at the start of every run.
|
||||
* Viewer: the detected beam-stop shadow is drawn in coral, with a "Show beam stop" switch in the side panel.
|
||||
* rugnux: Scaling and merging now apply a **low-resolution limit** of 50 Å (`--scaling-low-resolution <num>`, 0 removes it).
|
||||
|
||||
@@ -565,7 +565,24 @@ Pixels are classified by their squared distance $r^2=(x-x_p)^2+(y-y_p)^2$:
|
||||
- **signal region:** $r^2 < r_1^2$,
|
||||
- **background annulus:** $r_2^2 \le r^2 < r_3^2$.
|
||||
|
||||
Invalid pixels (masked/bad/saturated) are excluded from both sums. In addition, pixels lying inside the signal disk ($r<r_2$) of any *other* predicted reflection are removed from this reflection's background annulus, so a neighbouring spot cannot leak into the background estimate.
|
||||
Invalid pixels (masked/bad/saturated) are excluded from both sums. In addition, pixels lying inside the signal disk ($r<r_2$) of any *other* predicted reflection are removed from this reflection's background annulus, so a neighbouring spot cannot leak into the background estimate. (Both the annulus and that exclusion become ellipses when the option below is used; with it off, which is the default, they are the circles just described.)
|
||||
|
||||
**Radially elongated background ring (opt-in, `--integration-stencil <k>`, default 0).** The three radii above are fixed pixel counts, identical for every reflection at every resolution. A reflection is not round, though: a finite bandwidth streaks it radially by $\sigma_\mathrm{bw}=\text{bandwidth}\cdot R_\mathrm{px}$. On a radially smeared spot the fixed $6\ldots10$ px ring therefore sits only $\approx1.3$–$2.2$ radial $\sigma$ from the centre — on the reflection's own tails, which it then measures as background.
|
||||
|
||||
With $k>0$ the **background ring becomes an ellipse**, elongated along the beam→reflection direction by $k\sigma_\mathrm{bw}$. The **radial** semi-axes become $r_2+k\sigma_\mathrm{bw}$ and $r_3+k\sigma_\mathrm{bw}$; the **tangential** half-widths stay $r_2$ and $r_3$; and the growth is capped at $2r_3$, which bounds what a mis-declared bandwidth can do to the bounding box. Pixels are then classified as
|
||||
|
||||
- **signal region:** $r^2 < r_1^2$ — a circle, unchanged,
|
||||
- **background ring:** $r^2-q_\mathrm{in}\rho^2 \ge r_2^2$ **and** $r^2-q_\mathrm{out}\rho^2 < r_3^2$,
|
||||
|
||||
where $\rho$ is the pixel's radial offset (its projection on the beam→reflection direction), $g=\min(k\sigma_\mathrm{bw},\,2r_3)$ is the capped growth, and $q=1-\big(r/(r+g)\big)^2$ for the boundary concerned. Written this way $k=0$ gives $q=0$ and both tests collapse onto $r^2$ **exactly in floating point**, so the default classifies every pixel exactly as the circular stencil did. The neighbour exclusion above follows: each neighbour's **inner ellipse**, taken in that neighbour's own radial frame, is what is masked out of this reflection's ring.
|
||||
|
||||
The width is the bandwidth streak alone, and deliberately **not** the profile's full radial variance of §9.3, which also carries the sensor parallax and weak-spot capture terms. Those two are the only terms there are on a monochromatic beam, and widening the ring by them was measured on the rotation battery: it neither helped the crystals with clean high-resolution shells nor left the weak ones alone. The bandwidth streak, by contrast, is a measured elongation of the recorded spot — principal axis along the radius to within a couple of degrees, and azimuth-independent. Keeping only it also makes the option exactly inert on a monochromatic beam, where $\sigma_\mathrm{bw}$ is zero.
|
||||
|
||||
One caveat on "exactly": the *classification* is bit-identical at $k=0$, and so is every intensity on the default path. The radial-background curve of §9.2 is a separate matter — its per-pixel radial offset used to be computed in double precision on the CPU and is now computed in single, as it always was on the GPU, so with `--background-radial` on a pixel can change radial bin. That brings the two engines together rather than apart, but it does mean the CPU is not bit-compatible with earlier builds on that one path.
|
||||
|
||||
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.
|
||||
|
||||
### 9.2 Box summation (seed and fallback)
|
||||
|
||||
@@ -580,7 +597,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 stencil — 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 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.
|
||||
|
||||
**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.
|
||||
|
||||
@@ -588,11 +605,11 @@ The clip cuts only the high tail, which matters: the **symmetric** trimmed mean
|
||||
|
||||
Both estimators are computed in the shared background pass, but only the trim reaches plain box summation: the high-side clip is skipped for `--integrator boxsum`, which therefore uses the plain ring mean unless `--background-trim` is given.
|
||||
|
||||
**Radial background correction (opt-in).** A ring mean estimates the background *under* the signal disk correctly only if the background is flat there. The signal disk and the ring are concentric, so for a background that is **linear** in position $\langle B\rangle_\mathrm{ring}=\langle B\rangle_\mathrm{disk}$ identically — a plane or gradient fit buys exactly nothing. The leading error is the **curvature** of the radial background, which is negligible on a smooth background but reaches tens of counts on a single reflection sitting on a sharp powder ring. Because every reflection uses the same stencil, that error is a fixed kernel over radial offset,
|
||||
**Radial background correction (opt-in).** A ring mean estimates the background *under* the signal disk correctly only if the background is flat there. The signal disk and the ring are concentric, so for a background that is **linear** in position $\langle B\rangle_\mathrm{ring}=\langle B\rangle_\mathrm{disk}$ identically — a plane or gradient fit buys exactly nothing. The leading error is the **curvature** of the radial background, which is negligible on a smooth background but reaches tens of counts on a single reflection sitting on a sharp powder ring. That error is a kernel over radial offset,
|
||||
$
|
||||
\delta \hat b \;=\; \textstyle\sum_k \kappa_k\, \bar B(r_0+k),
|
||||
$
|
||||
with $\kappa$ obtained once by azimuth-averaging the stencil and $\bar B(r)$ the image's own radial background curve. Applying it costs one short dot product per reflection and no extra pixel reads; correcting the background *scalar* means the box sum, the profile fit and the variance all pick it up. The curve is accumulated from the same annulus pixels the background pass already reads (a pixel's radius is the reflection's radius plus the pixel's projection on the beam→reflection direction, so no per-pixel square root is needed) and specifically from the **clipped** pixels, or it would carry neighbour tails and zingers — which is why the correction is inert under `--integrator boxsum`, that path having no clip pass.
|
||||
with $\kappa$ the annulus-minus-disk histogram of the stencil over radial offset, averaged over azimuth, and $\bar B(r)$ the image's own radial background curve. With the fixed circular stencil ($k=0$, §9.1) that single kernel serves every reflection. An elongated ring does not: its radial-offset histogram depends on how far that particular reflection's ring was grown, so $\kappa$ becomes a small table of kernels, indexed by the growth rounded to whole pixels. The azimuthal average survives the change unaltered, because the stencil is rebuilt in the reflection's own frame at each azimuth and so stays radially aligned: what is averaged over is the sub-pixel phase of the detector grid against the radius, which is what genuinely differs between reflections. Applying it costs one short dot product per reflection and no extra pixel reads; correcting the background *scalar* means the box sum, the profile fit and the variance all pick it up. The curve is accumulated from the same annulus pixels the background pass already reads (a pixel's radius is the reflection's radius plus the pixel's projection on the beam→reflection direction, so no per-pixel square root is needed) and specifically from the **clipped** pixels, or it would carry neighbour tails and zingers — which is why the correction is inert under `--integrator boxsum`, that path having no clip pass.
|
||||
|
||||
The model is a function of **radius alone**, so it is applied only where that is true of the background. `--background-radial` takes `on`, `off` or `auto`. It is **off by default**; under `auto` each image's peak-excluded ice score (§3.3) is taken after spot detection and before integration, and the correction is applied to that image when the score reaches the same `--ice-min-score` gate. Smooth powder ice *is* a radial feature and is corrected; ice made of discrete crystallite spots — which the profile channel is blind to and the spot channel catches — leaves no smooth ring to model, and correcting it makes matters worse. Measured against a fixed atomic model, comparing ice bands with resolution-matched decoy bands carrying no ice: on a crystal with pure smooth ice the correction removes **43 % of the bands' excess amplitude**, and the improvement is **7× larger inside the bands than outside**, which is its stated mechanism; on a crystal whose ice is textured the same correction *increased* the excess amplitude by half; on a clean crystal it is inert to four decimal places. Auto engages only where a peak-excluded score exists (adaptive spot finding, §3.2) — a plain azimuthal profile carries the Bragg peaks and cannot support an absolute threshold, so without one auto leaves the correction off.
|
||||
|
||||
|
||||
+2
-1
@@ -345,8 +345,9 @@ Integration:
|
||||
| --- | --- |
|
||||
| `--integrator <txt>` | Spot integrator: `gaussian` (profile-fit, default) \| `empirical` \| `boxsum` (classical fallback) |
|
||||
| `--integration-radius <r>` | Signal-box radius `r1`, or `r1,r2,r3` (px). One value ⇒ `r2=r1+2`, `r3=r1+4` |
|
||||
| `--integration-stencil <k>` | Push the `r2..r3` background ring out by `k` times the beam's radial streak `bandwidth·R_px`, per reflection (default `0` = the fixed circular ring). A fixed ring otherwise ends up on a streaked reflection's own tails at high resolution and measures them as background. Only the ring moves, and only radially — the `r1` signal box stays a circle — and the growth is capped at `2·r3`. The neighbour exclusion grows with it, so on a crowded pattern a few reflections can be left with too little background and dropped. Needs `--bandwidth`: on a monochromatic beam the streak is zero and this does nothing |
|
||||
| `--background-clip <n>` | Monochromatic (rotation + still): high-side clip of the background ring at `mean + n·√mean` (default 4; 0 = off). The default background estimator — it rejects neighbour cores and zingers without the symmetric trim's Poisson skew bias. Broadband data always clip, at 3σ; ignored by `--integrator boxsum` |
|
||||
| `--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 |
|
||||
| `--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. Rings holding more than 512 pixels fall back to the plain mean (the GPU sorts the ring in shared memory and the CPU now matches it), which the default radii never reach but wide ones do |
|
||||
| `--background-radial[=on\|off\|auto]` | Correct the background ring for the **curvature** of the radial background (default **off**). Disk and 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, since the model is a function of radius alone: on ice made of discrete crystallite spots there is no smooth ring and the correction makes the bias worse. Ignored by `--integrator boxsum` (no clip pass to take the curve from) |
|
||||
| `--integration-high-resolution <num>` | High-resolution limit for prediction and integration. Omitted (or 0) means integration extends as far as the detector reaches — which is what the predictor can place on the detector anyway, since it rejects reflections that miss it. Set a value to integrate less than the detector offers |
|
||||
| `--max-hkl <n>` | Predict reflections with \|h\|,\|k\|,\|l\| ≤ `n` (max 511). By default this is derived per crystal from the refined cell as `ceil(max(a,b,c)/d_min) + 1`, which is the exact bound: the predictor keeps only \|q\| ≤ 1/d_min and `h = a·q`, so no reflection can lie outside it and no candidate inside it is wasted on a shorter axis. Set it only to override that |
|
||||
|
||||
@@ -70,6 +70,17 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe
|
||||
beam_y = geom.GetBeamY_pxl();
|
||||
use_ellipse = !empirical;
|
||||
|
||||
// Per-reflection signal/background geometry: the ring elongated radially by k_sigma times the
|
||||
// beam's own radial streak, capped. k_sigma = 0 is the fixed circular stencil, bit for bit, and
|
||||
// so is any monochromatic beam, where the streak is zero.
|
||||
stencil.beam_x = beam_x;
|
||||
stencil.beam_y = beam_y;
|
||||
stencil.r2 = r2;
|
||||
stencil.r3 = r3;
|
||||
stencil.bw_sigma = static_cast<float>(bw_sigma);
|
||||
stencil.k_sigma = settings.GetStencilKSigma();
|
||||
stencil.max_grow = bragg_engine::MAX_STENCIL_GROW_OVER_R3 * r3;
|
||||
|
||||
// Robust background ring, one estimator or the other (see BraggIntegrationSettings): a high-side
|
||||
// sigma-clip (rugnux --background-clip, the default) or, when the clip is switched off, a
|
||||
// symmetric trimmed mean (rugnux --background-trim). The caller owns the choice - the engine no
|
||||
@@ -88,33 +99,74 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe
|
||||
const auto radial = settings.GetBackgroundRadialCorrection();
|
||||
bkg_radial_auto = !radial.has_value();
|
||||
bkg_radial = radial.value_or(false);
|
||||
k_off = static_cast<int>(std::ceil(r3)) + 1;
|
||||
k_diff.assign(2 * k_off + 1, 0.0f);
|
||||
{
|
||||
std::vector<double> hist_disk(k_diff.size(), 0.0), hist_ann(k_diff.size(), 0.0);
|
||||
constexpr int n_phi = 512;
|
||||
const int span = static_cast<int>(std::ceil(r3)) + 1;
|
||||
for (int p = 0; p < n_phi; ++p) {
|
||||
const double phi = 2.0 * PI * p / n_phi, cp = std::cos(phi), sp = std::sin(phi);
|
||||
for (int dy = -span; dy <= span; ++dy)
|
||||
for (int dx = -span; dx <= span; ++dx) {
|
||||
const double d2 = static_cast<double>(dx) * dx + static_cast<double>(dy) * dy;
|
||||
const int k = k_off + static_cast<int>(std::lround(dx * cp + dy * sp));
|
||||
if (k < 0 || k >= static_cast<int>(k_diff.size()))
|
||||
continue;
|
||||
if (d2 < r1_sq) hist_disk[k] += 1.0;
|
||||
else if (d2 >= r2_sq && d2 < r3_sq) hist_ann[k] += 1.0;
|
||||
}
|
||||
}
|
||||
const double sd = std::accumulate(hist_disk.begin(), hist_disk.end(), 0.0);
|
||||
const double sa = std::accumulate(hist_ann.begin(), hist_ann.end(), 0.0);
|
||||
for (size_t k = 0; k < k_diff.size(); ++k)
|
||||
k_diff[k] = static_cast<float>(hist_ann[k] / sa - hist_disk[k] / sd);
|
||||
}
|
||||
|
||||
// The table spans zero growth up to whatever the widest reflection on this detector reaches, one
|
||||
// kernel per pixel of growth; with nothing elongated a single kernel is all there is, which is
|
||||
// the layout and the values of every build before the stencil existed. It is built only when the
|
||||
// correction can ever run - the rows are not cheap, and nothing may read them otherwise:
|
||||
// bkg_radial is raised after construction only by the auto mode (MXAnalysisWithoutFPGA), which
|
||||
// requires bkg_radial_auto, and the GPU allocates its curve buffers under the same condition.
|
||||
// n_kern is the largest row BraggStencilKernelIndex can select, plus one.
|
||||
r_max = std::hypot(std::max<double>(beam_x, static_cast<double>(xpixel) - beam_x),
|
||||
std::max<double>(beam_y, static_cast<double>(ypixel) - beam_y));
|
||||
const float grow_max = (bkg_radial || bkg_radial_auto)
|
||||
? BraggStencilGrow_px(static_cast<float>(r_max), stencil)
|
||||
: 0.0f;
|
||||
n_kern = static_cast<int>(std::lround(grow_max)) + 1;
|
||||
// Every row must fit: the last one is built at grow = n_kern - 1, which rounding can put just
|
||||
// above grow_max.
|
||||
k_off = static_cast<int>(std::ceil(r3 + std::max<double>(grow_max, n_kern - 1))) + 1;
|
||||
k_len = 2 * k_off + 1;
|
||||
k_diff.clear();
|
||||
k_diff.reserve(static_cast<size_t>(n_kern) * k_len);
|
||||
for (int j = 0; j < n_kern; ++j)
|
||||
BuildRadialKernel(static_cast<float>(j));
|
||||
|
||||
polarization = experiment.GetPolarizationFactor();
|
||||
}
|
||||
|
||||
void BraggIntegrationEngine::BuildRadialKernel(float grow) {
|
||||
// Histogram the stencil over radial offset, averaged over azimuth so the kernel does not depend
|
||||
// on where the reflection sits. The average is over the SUB-PIXEL PHASE of the detector grid
|
||||
// against the radial direction, not over the stencil's own orientation: the stencil is built in
|
||||
// the reflection's frame at each azimuth, so an elongated one stays aligned with the radius, as
|
||||
// it is on the detector. k_diff is the annulus histogram minus the disk histogram, each
|
||||
// normalised, so dot(k_diff, B) is directly mean_annulus(B) - mean_disk(B).
|
||||
// The signal disk is a circle whatever the ring does, so its histogram is the same for every
|
||||
// kernel in the table - build it once.
|
||||
const bool first = hist_disk.empty();
|
||||
if (first)
|
||||
hist_disk.assign(k_len, 0.0);
|
||||
std::vector<double> hist_ann(k_len, 0.0);
|
||||
constexpr int n_phi = 512;
|
||||
const int span = static_cast<int>(std::ceil(r3 + grow)) + 1;
|
||||
const float si = r2 / (r2 + grow), so = r3 / (r3 + grow);
|
||||
const double q_in = 1.0 - static_cast<double>(si) * si;
|
||||
const double q_out = 1.0 - static_cast<double>(so) * so;
|
||||
for (int p = 0; p < n_phi; ++p) {
|
||||
const double phi = 2.0 * PI * p / n_phi, cp = std::cos(phi), sp = std::sin(phi);
|
||||
for (int dy = -span; dy <= span; ++dy)
|
||||
for (int dx = -span; dx <= span; ++dx) {
|
||||
const double d2 = static_cast<double>(dx) * dx + static_cast<double>(dy) * dy;
|
||||
const double rad = dx * cp + dy * sp;
|
||||
const int k = k_off + static_cast<int>(std::lround(rad));
|
||||
if (k < 0 || k >= k_len)
|
||||
continue;
|
||||
const double rad2 = rad * rad;
|
||||
if (d2 < r1_sq) {
|
||||
if (first) hist_disk[k] += 1.0;
|
||||
} else if (d2 - q_in * rad2 >= r2_sq && d2 - q_out * rad2 < r3_sq) {
|
||||
hist_ann[k] += 1.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
if (first) sum_disk = std::accumulate(hist_disk.begin(), hist_disk.end(), 0.0);
|
||||
const double sd = sum_disk;
|
||||
const double sa = std::accumulate(hist_ann.begin(), hist_ann.end(), 0.0);
|
||||
for (int k = 0; k < k_len; ++k)
|
||||
k_diff.push_back(static_cast<float>(hist_ann[k] / sa - hist_disk[k] / sd));
|
||||
}
|
||||
|
||||
std::vector<Reflection> BraggIntegrationEngine::Finalize(const std::vector<Reflection> &predicted,
|
||||
size_t npredicted,
|
||||
const std::vector<BraggFitResult> &results,
|
||||
|
||||
@@ -42,6 +42,7 @@
|
||||
#include "../../common/DiffractionGeometry.h"
|
||||
#include "../../common/Reflection.h"
|
||||
#include "../image_preprocessing/ImagePreprocessorBuffer.h"
|
||||
#include "BraggStencil.h"
|
||||
|
||||
namespace bragg_engine {
|
||||
// Shared with both engines so the CPU and GPU paths stay numerically aligned.
|
||||
@@ -64,6 +65,20 @@ constexpr double PIXEL_VARIANCE_FLOOR = 0.01;
|
||||
// estimate, which adds ~0.4*sigma*sum(P^3)/sum(P^2)^2 to every sigma whatever the count rate.
|
||||
constexpr double WEIGHT_VARIANCE_MIN_FRACTION = 0.5;
|
||||
|
||||
// Most background-ring pixels a block can hold for the GPU trimmed-mean sort. Shared with the CPU so
|
||||
// that a ring which overflows it falls back to the plain mean in BOTH engines: the CPU sorts an
|
||||
// unbounded vector and would otherwise keep trimming where the GPU had silently stopped. Note that
|
||||
// an elongated ring makes this a function of resolution rather than a property of the dataset - the
|
||||
// ring area grows with the elongation, so on a wide enough stencil the estimator changes at a fixed
|
||||
// detector radius. The growth cap below keeps the default r2=6/r3=10 ring under the bound at any
|
||||
// bandwidth; the wider stills radii can cross it, and only ever with --background-trim, which is
|
||||
// off by default and kept for back compatibility.
|
||||
constexpr int BKG_TRIM_MAX = 512;
|
||||
// Ceiling on how far the background ring may be pushed out radially, as a multiple of r3 - so the
|
||||
// outer ellipse never exceeds (1 + this) * r3. Bounds what a mis-declared bandwidth can do to the
|
||||
// per-reflection bounding box, and with it the shared memory the GPU sizes from the widest ring.
|
||||
constexpr float MAX_STENCIL_GROW_OVER_R3 = 2.0f;
|
||||
|
||||
// Guard against profile-fit runaways: on a weak / near-zero reflection the reweighted Kabsch iteration
|
||||
// has no real peak to lock onto and can manufacture intensity the box sum never sees. Fall back to the
|
||||
// summation (box-sum) intensity when the profile result disagrees with the summation seed by more than
|
||||
@@ -108,6 +123,12 @@ protected:
|
||||
double c_radial; // radial variance coefficient of tan^2(2theta): parallax + capture
|
||||
double F_px; // detector distance expressed in pixels
|
||||
float beam_x, beam_y;
|
||||
double r_max; // distance from the beam centre to the far corner [px]
|
||||
// Per-reflection signal/background geometry (BraggStencil.h). With k_sigma = 0 this is the fixed
|
||||
// r1 disk + r2..r3 ring the integrator has always used, bit for bit; above 0 the RING is
|
||||
// elongated radially, per reflection, by the analytic radial smear. Both engines build every
|
||||
// stencil through MakeBraggStencil, so the geometry has one definition.
|
||||
BraggStencilParams stencil;
|
||||
// Effective symmetric trimmed-mean background fraction (BraggIntegrationSettings), used only when
|
||||
// the high-side clip is switched off - the two are alternatives. 0 = plain ring mean. Read by both
|
||||
// the CPU and GPU engines.
|
||||
@@ -115,17 +136,37 @@ protected:
|
||||
|
||||
// --- radial background curvature correction (BraggIntegrationSettings) ---
|
||||
// The disk and the annulus are concentric, so any background LINEAR in position cancels between
|
||||
// them; what survives is the curvature of the radial background. Every reflection uses the same
|
||||
// stencil, so mean_annulus(B) - mean_disk(B) of a radial B is a FIXED kernel over radial offset:
|
||||
// them; what survives is the curvature of the radial background. mean_annulus(B) - mean_disk(B)
|
||||
// of a radial B is a kernel over radial offset:
|
||||
// bkg_error = sum_k k_diff[k] * B(r0 + k - k_off)
|
||||
// That is one short dot product per reflection and reads no pixels. Built in the constructor.
|
||||
// The kernel is built whatever the setting says, so bkg_radial can be flipped between images at
|
||||
// no cost - which is what the auto mode does, applying the correction only to the images whose
|
||||
// background really is a smooth function of radius (see BackgroundRadial below).
|
||||
// That is one short dot product per reflection and reads no pixels. Built in the constructor,
|
||||
// so bkg_radial can be flipped between images at no cost - which is what the auto mode does,
|
||||
// applying the correction only to the images whose background really is a smooth function of
|
||||
// radius (see BackgroundRadial below). It is built only when that mode, or an explicit setting,
|
||||
// could ever raise the correction; an engine that can never apply it keeps the single circular
|
||||
// kernel and never reads it.
|
||||
//
|
||||
// With a circular stencil ONE kernel serves every reflection. An elongated ring does not: its
|
||||
// radial-offset histogram depends on how far that reflection's ring was grown. So k_diff holds
|
||||
// n_kern kernels of k_len each, indexed by the growth quantized to whole pixels
|
||||
// (BraggStencilKernelIndex); n_kern is 1 when nothing is elongated, which is the old
|
||||
// single-kernel layout unchanged. The azimuthal average is kept, and still means what it did:
|
||||
// the stencil is rebuilt in the reflection's own frame at each azimuth, so it stays radially
|
||||
// aligned and what is averaged over is the sub-pixel phase of the detector grid against the
|
||||
// radius. What an elongated ring rules out is one kernel for ALL of them, not the average.
|
||||
bool bkg_radial = false;
|
||||
bool bkg_radial_auto = false; // settings left it unset: decide per image from the ice score
|
||||
int k_off = 0; // index of offset 0 in k_diff
|
||||
std::vector<float> k_diff; // annulus-minus-disk weight per integer radial offset
|
||||
int k_off = 0; // index of offset 0 within one kernel
|
||||
int k_len = 0; // entries per kernel
|
||||
int n_kern = 1; // kernels in the table (1 = circular stencil)
|
||||
std::vector<float> k_diff; // n_kern * k_len, annulus-minus-disk weight per radial offset
|
||||
|
||||
// One radial-offset kernel for a ring grown by `grow` px, appended to k_diff. Kept out of line
|
||||
// so the circular and elongated cases cannot drift apart. The signal disk does not change with
|
||||
// the growth, so its histogram is built on the first call and reused.
|
||||
void BuildRadialKernel(float grow);
|
||||
std::vector<double> hist_disk;
|
||||
double sum_disk = 0.0;
|
||||
|
||||
DiffractionGeometry geom; // kept for the per-reflection polarization correction
|
||||
std::optional<float> polarization;
|
||||
|
||||
@@ -66,19 +66,22 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
|
||||
|
||||
auto grid_idx = [this](int dx, int dy) { return (dy + R) * G + (dx + R); };
|
||||
|
||||
// --- Reflection mask: mark the r2 signal disk of every predicted reflection so a neighbour's
|
||||
// disk is excluded from this reflection's r2..r3 background ring. ---
|
||||
// --- Reflection mask: mark the r2 signal region of every predicted reflection so a neighbour's
|
||||
// signal is excluded from this reflection's r2..r3 background ring. The region is the INNER
|
||||
// stencil ellipse, taken in each neighbour's OWN frame - an elongated neighbour whose streak
|
||||
// is still marked as a disk would leak its tails into this reflection's ring. ---
|
||||
std::vector<uint8_t> refl_mask(npixel, 0);
|
||||
for (size_t i = 0; i < npredicted; ++i) {
|
||||
const auto &r = predicted[i];
|
||||
const int x0 = std::max(0, static_cast<int>(std::floor(r.predicted_x - r2 - 1.0f)));
|
||||
const int x1 = std::min(W - 1, static_cast<int>(std::ceil(r.predicted_x + r2 + 1.0f)));
|
||||
const int y0 = std::max(0, static_cast<int>(std::floor(r.predicted_y - r2 - 1.0f)));
|
||||
const int y1 = std::min(H - 1, static_cast<int>(std::ceil(r.predicted_y + r2 + 1.0f)));
|
||||
const BraggStencil st = MakeBraggStencil(r.predicted_x, r.predicted_y, stencil);
|
||||
const int x0 = std::max(0, static_cast<int>(std::floor(r.predicted_x - st.ex_in - 1.0f)));
|
||||
const int x1 = std::min(W - 1, static_cast<int>(std::ceil(r.predicted_x + st.ex_in + 1.0f)));
|
||||
const int y0 = std::max(0, static_cast<int>(std::floor(r.predicted_y - st.ey_in - 1.0f)));
|
||||
const int y1 = std::min(H - 1, static_cast<int>(std::ceil(r.predicted_y + st.ey_in + 1.0f)));
|
||||
for (int y = y0; y <= y1; ++y)
|
||||
for (int x = x0; x <= x1; ++x) {
|
||||
const double d2 = (x - r.predicted_x) * (x - r.predicted_x) + (y - r.predicted_y) * (y - r.predicted_y);
|
||||
if (d2 < r2_sq) refl_mask[y * W + x] = 1;
|
||||
const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y);
|
||||
if (d.inner < r2_sq) refl_mask[y * W + x] = 1;
|
||||
}
|
||||
}
|
||||
|
||||
@@ -90,6 +93,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
|
||||
int64_t I_sum = 0; // kept so I can be rebuilt after the radial background correction
|
||||
int n_inner = 0;
|
||||
int r_bin = 0; // rounded distance from the beam centre, indexes the radial curve
|
||||
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;
|
||||
};
|
||||
@@ -112,16 +116,18 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
|
||||
for (size_t i = 0; i < npredicted; ++i) {
|
||||
const auto &r = predicted[i];
|
||||
Rough out;
|
||||
const int x0 = std::max(0, static_cast<int>(std::floor(r.predicted_x - r3 - 1.0)));
|
||||
const int x1 = std::min(W - 1, static_cast<int>(std::ceil(r.predicted_x + r3 + 1.0)));
|
||||
const int y0 = std::max(0, static_cast<int>(std::floor(r.predicted_y - r3 - 1.0)));
|
||||
const int y1 = std::min(H - 1, static_cast<int>(std::ceil(r.predicted_y + r3 + 1.0)));
|
||||
// This reflection's stencil: the r1 signal disk, and the r2..r3 ring elongated radially by
|
||||
// the analytic smear. The bounding box spans the OUTER ellipse, tightly - taking the largest
|
||||
// semi-axis in both directions instead would read up to 60% more pixels for nothing.
|
||||
const BraggStencil st = MakeBraggStencil(r.predicted_x, r.predicted_y, stencil);
|
||||
const int x0 = std::max(0, static_cast<int>(std::floor(r.predicted_x - st.ex_out - 1.0f)));
|
||||
const int x1 = std::min(W - 1, static_cast<int>(std::ceil(r.predicted_x + st.ex_out + 1.0f)));
|
||||
const int y0 = std::max(0, static_cast<int>(std::floor(r.predicted_y - st.ey_out - 1.0f)));
|
||||
const int y1 = std::min(H - 1, static_cast<int>(std::ceil(r.predicted_y + st.ey_out + 1.0f)));
|
||||
|
||||
// Unit vector beam -> reflection: a stencil pixel's radial offset is its projection on it.
|
||||
const double rx = r.predicted_x - beam_x, ry = r.predicted_y - beam_y;
|
||||
const double r0 = std::hypot(rx, ry);
|
||||
const double ux = r0 > 1e-6 ? rx / r0 : 1.0, uy = r0 > 1e-6 ? ry / r0 : 0.0;
|
||||
out.r_bin = std::clamp(static_cast<int>(std::lround(r0)), 0, n_rad - 1);
|
||||
// Both from the stencil's own radius, so nothing downstream is derived from a second one.
|
||||
out.r_bin = std::clamp(static_cast<int>(std::lround(st.r0)), 0, n_rad - 1);
|
||||
out.k_bin = BraggStencilKernelIndex(st, n_kern);
|
||||
|
||||
int64_t I_sum = 0, I_sum_x = 0, I_sum_y = 0, n_inner = 0, n_inner_valid = 0;
|
||||
double bkg_sum = 0.0;
|
||||
@@ -129,16 +135,16 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
|
||||
bkg_vals.clear();
|
||||
for (int y = y0; y <= y1; ++y)
|
||||
for (int x = x0; x <= x1; ++x) {
|
||||
const double d2 = (x - r.predicted_x) * (x - r.predicted_x) + (y - r.predicted_y) * (y - r.predicted_y);
|
||||
const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y);
|
||||
const int32_t px = img[y * W + x];
|
||||
if (d2 < r1_sq) {
|
||||
if (d.signal < r1_sq) {
|
||||
++n_inner;
|
||||
if (!valid(px)) continue;
|
||||
I_sum += px;
|
||||
I_sum_x += static_cast<int64_t>(x) * px;
|
||||
I_sum_y += static_cast<int64_t>(y) * px;
|
||||
++n_inner_valid;
|
||||
} else if (d2 >= r2_sq && d2 < r3_sq) {
|
||||
} else if (d.inner >= r2_sq && d.outer < r3_sq) {
|
||||
if (refl_mask[y * W + x]) continue;
|
||||
if (!valid(px)) continue;
|
||||
bkg_sum += static_cast<double>(px);
|
||||
@@ -150,7 +156,8 @@ 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) {
|
||||
out.bkg = bkg_sum / n_bkg;
|
||||
if (bkg_trim_frac > 0.0 && bkg_vals.size() > 5) {
|
||||
if (bkg_trim_frac > 0.0 && bkg_vals.size() > 5
|
||||
&& bkg_vals.size() <= static_cast<size_t>(bragg_engine::BKG_TRIM_MAX)) {
|
||||
// Symmetric trimmed mean over the background ring (idea 1): drop the lowest and highest
|
||||
// bkg_trim_frac of the pixels, average the rest. Robust to the high-side contamination
|
||||
// that biases the plain ring mean and makes it over-subtract at high resolution.
|
||||
@@ -171,8 +178,8 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
|
||||
int n = 0;
|
||||
for (int y = y0; y <= y1; ++y)
|
||||
for (int x = x0; x <= x1; ++x) {
|
||||
const double d2 = (x - r.predicted_x) * (x - r.predicted_x) + (y - r.predicted_y) * (y - r.predicted_y);
|
||||
if (!(d2 >= r2_sq && d2 < r3_sq)) continue;
|
||||
const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y);
|
||||
if (!(d.inner >= r2_sq && d.outer < r3_sq)) continue;
|
||||
if (refl_mask[y * W + x]) continue;
|
||||
const int32_t px = img[y * W + x];
|
||||
if (!valid(px)) continue;
|
||||
@@ -180,8 +187,9 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
|
||||
s += px;
|
||||
++n;
|
||||
if (bkg_radial) {
|
||||
const double off = (x - r.predicted_x) * ux + (y - r.predicted_y) * uy;
|
||||
const int b = std::clamp(static_cast<int>(std::lround(r0 + off)), 0, n_rad - 1);
|
||||
// The radial curve is binned on the TRUE detector radius, so the
|
||||
// offset here is the unshrunk radial projection.
|
||||
const int b = std::clamp(static_cast<int>(std::lround(st.r0 + d.rad)), 0, n_rad - 1);
|
||||
rad_sum[b] += static_cast<double>(px);
|
||||
++rad_cnt[b];
|
||||
}
|
||||
@@ -230,11 +238,12 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
|
||||
if (!rh.ok) continue;
|
||||
// An empty bin contributes the reflection's own background, so a fully empty
|
||||
// neighbourhood gives corr == 0 exactly (the kernel weights sum to zero).
|
||||
const float *kern = k_diff.data() + static_cast<size_t>(rh.k_bin) * k_len;
|
||||
double corr = 0.0;
|
||||
for (int k = 0; k < static_cast<int>(k_diff.size()); ++k) {
|
||||
for (int k = 0; k < k_len; ++k) {
|
||||
const int b = std::clamp(rh.r_bin + k - k_off, 0, n_rad - 1);
|
||||
const double v = rad_cnt[b] > 0 ? rad_sum[b] / rad_cnt[b] : rh.bkg;
|
||||
corr += static_cast<double>(k_diff[k]) * v;
|
||||
corr += static_cast<double>(kern[k]) * v;
|
||||
}
|
||||
rh.bkg -= corr; // annulus mean -> mean over the signal disk
|
||||
rh.I = static_cast<double>(rh.I_sum) - static_cast<double>(rh.n_inner) * rh.bkg;
|
||||
|
||||
@@ -25,31 +25,33 @@ struct BraggGpuParams {
|
||||
float F_px;
|
||||
float beam_x, beam_y;
|
||||
float bkg_trim; // idea 1: symmetric trimmed-mean background fraction (0 = plain ring mean)
|
||||
BraggStencilParams stencil; // per-reflection signal/background geometry (BraggStencil.h)
|
||||
int rad_w; // radial-background window held in shared memory, in bins of one pixel
|
||||
int n_kern; // radial-background kernels in the table
|
||||
};
|
||||
|
||||
__device__ inline bool valid(int32_t v) { return v != INT32_MIN && v != INT32_MAX; }
|
||||
|
||||
// idea 1: max annulus pixels held in shared memory for the per-block trimmed-mean sort (else flat mean).
|
||||
constexpr int BKG_TRIM_MAX = 512;
|
||||
|
||||
// Learned second moments, (sum v*rad^2, sum v*tan^2, sum v) per shell plus one global slot at N_SHELL.
|
||||
constexpr int MOM_STRIDE = 3;
|
||||
|
||||
// --- Mark the r2 signal disk of every predicted reflection (race-free: all writes are 1). ---
|
||||
// --- Mark the r2 signal region of every predicted reflection (race-free: all writes are 1). The
|
||||
// region is the INNER stencil ellipse in the neighbour's own frame; see the CPU engine. ---
|
||||
__global__ void mark_mask(const float *px_x, const float *px_y, uint8_t *mask, BraggGpuParams p, int n) {
|
||||
const int i = blockIdx.x;
|
||||
if (i >= n) return;
|
||||
const float cx = px_x[i], cy = px_y[i];
|
||||
const int x0 = max(0, (int) floorf(cx - p.r2 - 1.0f));
|
||||
const int x1 = min(p.W - 1, (int) ceilf(cx + p.r2 + 1.0f));
|
||||
const int y0 = max(0, (int) floorf(cy - p.r2 - 1.0f));
|
||||
const int y1 = min(p.H - 1, (int) ceilf(cy + p.r2 + 1.0f));
|
||||
const BraggStencil st = MakeBraggStencil(cx, cy, p.stencil);
|
||||
const int x0 = max(0, (int) floorf(cx - st.ex_in - 1.0f));
|
||||
const int x1 = min(p.W - 1, (int) ceilf(cx + st.ex_in + 1.0f));
|
||||
const int y0 = max(0, (int) floorf(cy - st.ey_in - 1.0f));
|
||||
const int y1 = min(p.H - 1, (int) ceilf(cy + st.ey_in + 1.0f));
|
||||
const int bw = x1 - x0 + 1, bh = y1 - y0 + 1;
|
||||
if (bw <= 0 || bh <= 0) return;
|
||||
for (int t = threadIdx.x; t < bw * bh; t += blockDim.x) {
|
||||
const int x = x0 + t % bw, y = y0 + t / bw;
|
||||
const float ddx = (float) x - cx, ddy = (float) y - cy;
|
||||
if (ddx * ddx + ddy * ddy < p.r2_sq) mask[y * p.W + x] = 1;
|
||||
const BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy);
|
||||
if (d.inner < p.r2_sq) mask[y * p.W + x] = 1;
|
||||
}
|
||||
}
|
||||
|
||||
@@ -59,7 +61,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
|
||||
int *cx_o, int *cy_o, float *I_o, float *sigma_o, float *bkg_o,
|
||||
float *bkgvar_o, float *varbkg_o, float *obsx_o, float *obsy_o, uint8_t *ok_o, uint8_t *strong_o,
|
||||
uint8_t *hasobs_o, unsigned long long *invd2mm,
|
||||
float *isum_o, int *ninner_o, int *rbin_o,
|
||||
float *isum_o, int *ninner_o, int *rbin_o, int *kbin_o,
|
||||
float *rad_sum, int *rad_cnt, int n_rad) {
|
||||
const int i = blockIdx.x;
|
||||
if (i >= n) return;
|
||||
@@ -70,29 +72,30 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
|
||||
__shared__ int s_accept;
|
||||
__shared__ double s_bkg, s_thr, s_clipsum;
|
||||
__shared__ int s_clipn;
|
||||
__shared__ float s_ux, s_uy, s_r0;
|
||||
// The whole stencil spans r0 +- r3, so a 32-bin window centred on the reflection always covers it.
|
||||
constexpr int RAD_W = 32;
|
||||
__shared__ float s_radv[RAD_W];
|
||||
__shared__ int s_radn[RAD_W];
|
||||
__shared__ float s_r0;
|
||||
// The stencil spans r0 +- the radial semi-axis of the outer ellipse, so the window has to be
|
||||
// sized from the widest aperture on the detector - which the host does, into p.rad_w. A fixed
|
||||
// 32-bin window silently dropped the outer bins as soon as anything was elongated.
|
||||
extern __shared__ float s_rad[];
|
||||
float *s_radv = s_rad;
|
||||
int *s_radn = (int *) (s_rad + p.rad_w);
|
||||
__shared__ int s_radbase;
|
||||
if (threadIdx.x == 0) {
|
||||
const float rx = px_x[i] - p.beam_x, ry = px_y[i] - p.beam_y;
|
||||
s_r0 = sqrtf(rx * rx + ry * ry);
|
||||
s_ux = s_r0 > 1e-6f ? rx / s_r0 : 1.0f;
|
||||
s_uy = s_r0 > 1e-6f ? ry / s_r0 : 0.0f;
|
||||
}
|
||||
if (threadIdx.x == 0) {
|
||||
s_Isum = 0; s_Ix = 0; s_Iy = 0;
|
||||
s_ninner = 0; s_ninner_valid = 0; s_nbkg = 0; s_bkgsum = 0.0;
|
||||
}
|
||||
__syncthreads();
|
||||
|
||||
const float cx = px_x[i], cy = px_y[i];
|
||||
const int x0 = max(0, (int) floorf(cx - p.r3 - 1.0f));
|
||||
const int x1 = min(p.W - 1, (int) ceilf(cx + p.r3 + 1.0f));
|
||||
const int y0 = max(0, (int) floorf(cy - p.r3 - 1.0f));
|
||||
const int y1 = min(p.H - 1, (int) ceilf(cy + p.r3 + 1.0f));
|
||||
// Every thread builds the same stencil from the same inputs; it is a few flops against a whole
|
||||
// bounding box of pixel reads, so it costs less than staging it through shared memory.
|
||||
const BraggStencil st = MakeBraggStencil(cx, cy, p.stencil);
|
||||
const int x0 = max(0, (int) floorf(cx - st.ex_out - 1.0f));
|
||||
const int x1 = min(p.W - 1, (int) ceilf(cx + st.ex_out + 1.0f));
|
||||
const int y0 = max(0, (int) floorf(cy - st.ey_out - 1.0f));
|
||||
const int y1 = min(p.H - 1, (int) ceilf(cy + st.ey_out + 1.0f));
|
||||
const int bw = x1 - x0 + 1, bh = y1 - y0 + 1;
|
||||
const int area = (bw > 0 && bh > 0) ? bw * bh : 0;
|
||||
|
||||
@@ -101,13 +104,12 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
|
||||
double l_bkg = 0.0;
|
||||
for (int t = threadIdx.x; t < area; t += blockDim.x) {
|
||||
const int x = x0 + t % bw, y = y0 + t / bw;
|
||||
const float ddx = (float) x - cx, ddy = (float) y - cy;
|
||||
const float d2 = ddx * ddx + ddy * ddy;
|
||||
const BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy);
|
||||
const int32_t px = img[y * p.W + x];
|
||||
if (d2 < p.r1_sq) {
|
||||
if (d.signal < p.r1_sq) {
|
||||
++l_ni;
|
||||
if (valid(px)) { l_Isum += px; l_Ix += (long long) x * px; l_Iy += (long long) y * px; ++l_niv; }
|
||||
} else if (d2 >= p.r2_sq && d2 < p.r3_sq) {
|
||||
} else if (d.inner >= p.r2_sq && d.outer < p.r3_sq) {
|
||||
if (mask[y * p.W + x]) continue;
|
||||
if (!valid(px)) continue;
|
||||
l_bkg += (double) px; ++l_nb;
|
||||
@@ -122,8 +124,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
|
||||
atomicAdd(&s_bkgsum, l_bkg);
|
||||
__syncthreads();
|
||||
|
||||
for (int t = threadIdx.x; t < RAD_W; t += blockDim.x) { s_radv[t] = 0.0f; s_radn[t] = 0; }
|
||||
if (threadIdx.x == 0) s_radbase = (int) lroundf(s_r0) - RAD_W / 2;
|
||||
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;
|
||||
s_bkg = s_accept ? (s_bkgsum / (double) s_nbkg) : 0.0;
|
||||
@@ -143,9 +145,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
|
||||
if (do_trim) {
|
||||
for (int t = threadIdx.x; t < area; t += blockDim.x) {
|
||||
const int x = x0 + t % bw, y = y0 + t / bw;
|
||||
const float ddx = (float) x - cx, ddy = (float) y - cy;
|
||||
const float d2 = ddx * ddx + ddy * ddy;
|
||||
if (!(d2 >= p.r2_sq && d2 < p.r3_sq)) continue;
|
||||
const BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy);
|
||||
if (!(d.inner >= p.r2_sq && d.outer < p.r3_sq)) continue;
|
||||
if (mask[y * p.W + x]) continue;
|
||||
const int32_t px = img[y * p.W + x];
|
||||
if (!valid(px)) continue;
|
||||
@@ -186,19 +187,18 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
|
||||
double c_l = 0.0; int cn_l = 0;
|
||||
for (int t = threadIdx.x; t < area; t += blockDim.x) {
|
||||
const int x = x0 + t % bw, y = y0 + t / bw;
|
||||
const float ddx = (float) x - cx, ddy = (float) y - cy;
|
||||
const float d2 = ddx * ddx + ddy * ddy;
|
||||
if (!(d2 >= p.r2_sq && d2 < p.r3_sq)) continue;
|
||||
const BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy);
|
||||
if (!(d.inner >= p.r2_sq && d.outer < p.r3_sq)) continue;
|
||||
if (mask[y * p.W + x]) continue;
|
||||
const int32_t px = img[y * p.W + x];
|
||||
if (!valid(px)) continue;
|
||||
if ((double) px <= s_thr) {
|
||||
c_l += px; ++cn_l;
|
||||
if (n_rad > 0) {
|
||||
// Radial offset = projection on the beam->reflection direction, no per-pixel sqrt.
|
||||
const float off = ddx * s_ux + ddy * s_uy;
|
||||
const int idx = (int) lroundf(s_r0 + off) - s_radbase;
|
||||
if (idx >= 0 && idx < RAD_W) {
|
||||
// The radial curve is binned on the TRUE detector radius, so this is the
|
||||
// unshrunk radial projection - no per-pixel sqrt.
|
||||
const int idx = (int) lroundf(s_r0 + d.rad) - s_radbase;
|
||||
if (idx >= 0 && idx < p.rad_w) {
|
||||
atomicAdd(&s_radv[idx], (float) px); // shared, not global
|
||||
atomicAdd(&s_radn[idx], 1);
|
||||
}
|
||||
@@ -210,7 +210,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
|
||||
__syncthreads();
|
||||
|
||||
if (n_rad > 0) {
|
||||
for (int t = threadIdx.x; t < RAD_W; t += blockDim.x) {
|
||||
for (int t = threadIdx.x; t < p.rad_w; t += blockDim.x) {
|
||||
if (s_radn[t] == 0) continue;
|
||||
const int b = min(max(s_radbase + t, 0), n_rad - 1);
|
||||
atomicAdd(&rad_sum[b], s_radv[t]);
|
||||
@@ -252,6 +252,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd,
|
||||
isum_o[i] = (float) Isum;
|
||||
ninner_o[i] = s_ninner;
|
||||
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;
|
||||
@@ -387,20 +388,24 @@ __global__ void build_profiles(const float *shell_grid, const float *global_grid
|
||||
// them; what is left is the curvature of the radial background. k_diff (annulus minus disk
|
||||
// histogram over radial offset) turns that into one short dot product. An empty radial bin
|
||||
// contributes the reflection's own background, so an empty neighbourhood gives exactly zero
|
||||
// correction because the kernel weights sum to zero. Mirrors BraggIntegrationEngineCPU. ---
|
||||
// correction because the kernel weights sum to zero. The kernel depends on how far this
|
||||
// reflection's ring was elongated, so k_diff is a table and kbin_a picks the row.
|
||||
// Mirrors BraggIntegrationEngineCPU. ---
|
||||
__global__ void radial_correct(const float *rad_sum, const int *rad_cnt, int n_rad,
|
||||
const float *k_diff, int k_len, int k_off,
|
||||
const float *isum_a, const int *ninner_a, const int *rbin_a,
|
||||
const int *kbin_a,
|
||||
const uint8_t *ok_a, float *bkg_o, float *I_o, int n) {
|
||||
const int i = blockIdx.x * blockDim.x + threadIdx.x;
|
||||
if (i >= n || !ok_a[i]) return;
|
||||
const float bkg = bkg_o[i];
|
||||
const float *kern = k_diff + (size_t) kbin_a[i] * k_len;
|
||||
float corr = 0.0f;
|
||||
for (int k = 0; k < k_len; ++k) {
|
||||
int b = rbin_a[i] + k - k_off;
|
||||
b = min(max(b, 0), n_rad - 1);
|
||||
const float v = rad_cnt[b] > 0 ? rad_sum[b] / (float) rad_cnt[b] : bkg;
|
||||
corr += k_diff[k] * v;
|
||||
corr += kern[k] * v;
|
||||
}
|
||||
const float bkg_new = bkg - corr;
|
||||
bkg_o[i] = bkg_new;
|
||||
@@ -540,6 +545,12 @@ BraggIntegrationEngineGPU::BraggIntegrationEngineGPU(const DiffractionExperiment
|
||||
const int max_Gf = 2 * max_Rf + 1;
|
||||
fit_shared_bytes = static_cast<size_t>(max_Gf) * max_Gf * sizeof(float);
|
||||
|
||||
// Shared radial window of the background curve in boxsum. The stencil spans r0 +- the radial
|
||||
// semi-axis of the outer ellipse, widest at the far corner of the detector, and the window
|
||||
// covers twice that so a reflection anywhere on the detector fits.
|
||||
rad_w = 2 * (static_cast<int>(std::ceil(r3 + BraggStencilGrow_px(static_cast<float>(r_max), stencil))) + 1) + 1;
|
||||
boxsum_shared_bytes = static_cast<size_t>(rad_w) * (sizeof(float) + sizeof(int));
|
||||
|
||||
// The current device, not device 0: workers are pinned round-robin across GPUs, so device 0's
|
||||
// shared-memory size can belong to a different card than the one these kernels launch on.
|
||||
int device = 0;
|
||||
@@ -549,14 +560,17 @@ BraggIntegrationEngineGPU::BraggIntegrationEngineGPU(const DiffractionExperiment
|
||||
if (fit_shared_bytes > prop.sharedMemPerBlock)
|
||||
throw JFJochException(JFJochExceptionCategory::GPUCUDAError,
|
||||
"BraggIntegrationEngineGPU: profile grid exceeds shared memory (r2 too large)");
|
||||
// boxsum's dynamic window sits alongside its static shared arrays (the trimmed-mean buffer
|
||||
// dominates them), so leave room for those rather than budgeting the whole block.
|
||||
if (boxsum_shared_bytes + sizeof(int) * BKG_TRIM_MAX + 256 > prop.sharedMemPerBlock)
|
||||
throw JFJochException(JFJochExceptionCategory::GPUCUDAError,
|
||||
"BraggIntegrationEngineGPU: background ring exceeds shared memory");
|
||||
|
||||
// Radial background curve: one bin per pixel of distance from the beam, out to the far corner.
|
||||
// Allocated whenever the correction COULD run, so the auto mode can turn it on for an individual
|
||||
// image; whether it runs for a given image is decided in Run() from the per-image bkg_radial.
|
||||
if (bkg_radial || bkg_radial_auto) {
|
||||
const double fx = std::max<double>(beam_x, static_cast<double>(xpixel) - beam_x);
|
||||
const double fy = std::max<double>(beam_y, static_cast<double>(ypixel) - beam_y);
|
||||
n_rad = static_cast<int>(std::ceil(std::hypot(fx, fy))) + 2;
|
||||
n_rad = static_cast<int>(std::ceil(r_max)) + 2;
|
||||
d_rad_sum = CudaDevicePtr<float>(n_rad);
|
||||
d_rad_cnt = CudaDevicePtr<int>(n_rad);
|
||||
d_k_diff = CudaDevicePtr<float>(k_diff.size());
|
||||
@@ -585,6 +599,7 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) {
|
||||
d_isum = CudaDevicePtr<float>(new_capacity);
|
||||
d_ninner = CudaDevicePtr<int>(new_capacity);
|
||||
d_rbin = CudaDevicePtr<int>(new_capacity);
|
||||
d_kbin = CudaDevicePtr<int>(new_capacity);
|
||||
d_obs_x = CudaDevicePtr<float>(new_capacity);
|
||||
d_obs_y = CudaDevicePtr<float>(new_capacity);
|
||||
d_ok = CudaDevicePtr<uint8_t>(new_capacity);
|
||||
@@ -633,6 +648,7 @@ std::vector<Reflection> BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu
|
||||
.F_px = static_cast<float>(F_px),
|
||||
.beam_x = beam_x, .beam_y = beam_y,
|
||||
.bkg_trim = bkg_trim, // effective trim fraction (0 for stills), set by the base ctor from settings
|
||||
.stencil = stencil, .rad_w = rad_w, .n_kern = n_kern,
|
||||
};
|
||||
|
||||
// Whether the radial correction runs for THIS image. n_rad only says the buffers exist - under
|
||||
@@ -647,18 +663,18 @@ std::vector<Reflection> BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu
|
||||
}
|
||||
reset<<<32, 256, 0, *stream>>>(d_shell_grid, d_global_grid, d_mom, d_shell_n, d_global_n, d_invd2, GG);
|
||||
mark_mask<<<n, threads, 0, *stream>>>(d_px_x, d_px_y, d_mask, p, n);
|
||||
boxsum<<<n, threads, 0, *stream>>>(d_px_x, d_px_y, d_d, img, d_mask, p, n,
|
||||
boxsum<<<n, threads, boxsum_shared_bytes, *stream>>>(d_px_x, d_px_y, d_d, img, d_mask, p, n,
|
||||
d_cx, d_cy, d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_obs_x, d_obs_y,
|
||||
d_ok, d_strong, d_has_obs, d_invd2,
|
||||
d_isum, d_ninner, d_rbin,
|
||||
d_isum, d_ninner, d_rbin, d_kbin,
|
||||
d_rad_sum, d_rad_cnt, rad_n);
|
||||
|
||||
// Correct the flat annulus background for the curvature of the radial background before anything
|
||||
// downstream (profile fit, variance) reads it.
|
||||
if (rad_n > 0)
|
||||
radial_correct<<<(n + threads - 1) / threads, threads, 0, *stream>>>(
|
||||
d_rad_sum, d_rad_cnt, rad_n, d_k_diff, static_cast<int>(k_diff.size()), k_off,
|
||||
d_isum, d_ninner, d_rbin, d_ok, d_bkg, d_I, n);
|
||||
d_rad_sum, d_rad_cnt, rad_n, d_k_diff, k_len, k_off,
|
||||
d_isum, d_ninner, d_rbin, d_kbin, d_ok, d_bkg, d_I, n);
|
||||
|
||||
if (mode != IntegratorMode::BoxSum) {
|
||||
// Pass B: learn (shell computed inline) -> build -> fit.
|
||||
|
||||
@@ -23,6 +23,8 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine {
|
||||
std::shared_ptr<CudaStream> stream;
|
||||
int threads;
|
||||
size_t fit_shared_bytes;
|
||||
int rad_w = 0; // radial-background window of boxsum, in bins of one pixel
|
||||
size_t boxsum_shared_bytes = 0;
|
||||
|
||||
size_t capacity = 0; // per-reflection device/host arrays hold at least this many reflections
|
||||
|
||||
@@ -31,7 +33,7 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine {
|
||||
CudaDevicePtr<int> d_cx, d_cy;
|
||||
CudaDevicePtr<float> d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_obs_x, d_obs_y;
|
||||
CudaDevicePtr<float> d_isum; // box-sum raw sum, for the radial correction
|
||||
CudaDevicePtr<int> d_ninner, d_rbin;
|
||||
CudaDevicePtr<int> d_ninner, d_rbin, d_kbin;
|
||||
CudaDevicePtr<uint8_t> d_ok, d_strong, d_has_obs;
|
||||
|
||||
// --- radial background curvature correction (see BraggIntegrationEngine) ---
|
||||
@@ -42,7 +44,7 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine {
|
||||
// --- fixed-size device arrays ---
|
||||
// The learning/fit math is single precision: FP64 is heavily throttled on consumer GPUs and the
|
||||
// extraction is Poisson-noise limited, so float reproduces the double CPU path to ~1e-4.
|
||||
CudaDevicePtr<uint8_t> d_mask; // per-pixel r2-disk reflection mask
|
||||
CudaDevicePtr<uint8_t> d_mask; // per-pixel inner-stencil reflection mask
|
||||
CudaDevicePtr<float> d_shell_grid, d_global_grid; // learned profile accumulators (N_SHELL*GG, GG)
|
||||
CudaDevicePtr<float> d_shell_P, d_global_P; // normalised profiles (empirical mode)
|
||||
CudaDevicePtr<float> d_mom; // learned 2nd moments, 3 per shell + global
|
||||
|
||||
@@ -0,0 +1,142 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#pragma once
|
||||
|
||||
// =============================================================================
|
||||
// BraggStencil - the per-reflection signal/background geometry, shared by both engines
|
||||
// =============================================================================
|
||||
//
|
||||
// The integrator classifies every pixel of a reflection's neighbourhood as signal (inside r1),
|
||||
// background (the r2..r3 ring) or neither. Those three radii used to be plain circles of a fixed
|
||||
// pixel size, which is wrong for a spot that is not round: a bandwidth streaks a reflection
|
||||
// radially by sigma_bw = bw_sigma*Rpx, so at high resolution the r2..r3 ring sits on the
|
||||
// reflection's own tails and measures them as background.
|
||||
//
|
||||
// So the RING is an ellipse, elongated along the beam->reflection direction by k_sigma times that
|
||||
// streak, and left alone tangentially. The SIGNAL disk stays a circle: r1 governs n_inner, var_bkg
|
||||
// and the all-or-nothing `n_inner_valid == n_inner` acceptance gate, so growing it rejects every
|
||||
// reflection with one bad pixel anywhere along a long streak. The flux it loses is a function of
|
||||
// resolution alone and the per-shell scale absorbs it.
|
||||
//
|
||||
// The width is the BANDWIDTH STREAK ALONE, and deliberately not the profile's full radial variance
|
||||
// (which also carries a sensor parallax term and a weak-spot capture term). Those two were measured
|
||||
// on the rotation battery: they widen the ring on monochromatic data, where they are the only terms
|
||||
// there are, and that neither helped the two crystals with clean high-resolution shells nor left the
|
||||
// two weak ones alone. The bandwidth streak, by contrast, is a measured elongation of the recorded
|
||||
// spot - principal axis along the radius to within a couple of degrees, azimuth-independent - and it
|
||||
// is the term that carries the pink-beam case. Keeping only it also makes the whole feature exactly
|
||||
// inert on a monochromatic beam, where bw_sigma is 0 and so is the growth.
|
||||
//
|
||||
// Rather than evaluate an ellipse, each pixel's squared distance has the radial part scaled down:
|
||||
// with rad the pixel's radial offset and d2 its plain squared distance,
|
||||
//
|
||||
// d2_inner = d2 - q_in *rad^2 tested against r2^2 (q_in = 1 - (r2/(r2+grow))^2)
|
||||
// d2_outer = d2 - q_out*rad^2 tested against r3^2 (q_out = 1 - (r3/(r3+grow))^2)
|
||||
//
|
||||
// which is the ellipse (radial semi-axis r+grow, tangential semi-axis r) written so that grow = 0
|
||||
// gives q = 0 and hence d2_inner == d2_outer == d2 EXACTLY, in floating point, not merely to within
|
||||
// rounding. k_sigma = 0 therefore reproduces the old circular stencil bit for bit.
|
||||
//
|
||||
// This header is compiled by both the host compiler and nvcc, so the geometry has exactly one
|
||||
// definition: the two engines classify by the same rule, built from the same inputs. They are not
|
||||
// bit-identical to each other - nvcc contracts a*b+c into an FMA and the host baseline has no FMA
|
||||
// instruction to contract into - so a pixel within a rounding of a boundary can still fall either
|
||||
// way, which is what the parity test's tolerances are for.
|
||||
// =============================================================================
|
||||
|
||||
#include <cmath>
|
||||
|
||||
#ifdef __CUDACC__
|
||||
#define BRAGG_STENCIL_HD __host__ __device__ inline
|
||||
#else
|
||||
#define BRAGG_STENCIL_HD inline
|
||||
#endif
|
||||
|
||||
// Fixed per-experiment inputs to the stencil law (mirrors BraggIntegrationEngine's members).
|
||||
struct BraggStencilParams {
|
||||
float beam_x = 0.0f, beam_y = 0.0f;
|
||||
float r2 = 6.0f, r3 = 10.0f;
|
||||
float bw_sigma = 0.0f; // radial streak per pixel of radius (bandwidth sigma, dimensionless)
|
||||
float k_sigma = 0.0f; // radial sigmas to push the ring out by; 0 = the old circular stencil
|
||||
float max_grow = 0.0f; // hard cap on the radial growth [px]. 0 disables growth, not the cap.
|
||||
};
|
||||
|
||||
// One reflection's stencil, in its own radial/tangential frame.
|
||||
struct BraggStencil {
|
||||
float ux = 1.0f, uy = 0.0f; // unit vector beam -> reflection
|
||||
float r0 = 0.0f; // distance from the beam centre [px]
|
||||
float grow = 0.0f; // radial growth of the ring [px] (0 = circular)
|
||||
float q_in = 0.0f, q_out = 0.0f; // radial shrink coefficients (0 = circular)
|
||||
float ex_in = 0.0f, ey_in = 0.0f; // axis-aligned half-extent of the inner (r2) ellipse
|
||||
float ex_out = 0.0f, ey_out = 0.0f; // axis-aligned half-extent of the outer (r3) ellipse
|
||||
};
|
||||
|
||||
// Radial growth of the ring at Rpx, in pixels: k_sigma times the bandwidth streak, capped. Shared
|
||||
// so the stencil and the radial-background kernel table agree on what "grown by this much" means.
|
||||
BRAGG_STENCIL_HD float BraggStencilGrow_px(float Rpx, const BraggStencilParams &p) {
|
||||
if (!(p.k_sigma > 0.0f)) return 0.0f;
|
||||
// The cap bounds what a mis-declared bandwidth can do to the bounding box.
|
||||
const float grow = p.k_sigma * p.bw_sigma * Rpx;
|
||||
return grow > p.max_grow ? p.max_grow : grow;
|
||||
}
|
||||
|
||||
BRAGG_STENCIL_HD BraggStencil MakeBraggStencil(float px_x, float px_y, const BraggStencilParams &p) {
|
||||
BraggStencil s;
|
||||
const float rx = px_x - p.beam_x, ry = px_y - p.beam_y;
|
||||
const float r0 = sqrtf(rx * rx + ry * ry);
|
||||
s.r0 = r0;
|
||||
if (r0 > 1e-6f) {
|
||||
s.ux = rx / r0;
|
||||
s.uy = ry / r0;
|
||||
}
|
||||
|
||||
s.grow = BraggStencilGrow_px(r0, p);
|
||||
const float grow = s.grow;
|
||||
if (grow == 0.0f) { // the common case, and the only one the online path can reach
|
||||
s.ex_in = p.r2; s.ey_in = p.r2;
|
||||
s.ex_out = p.r3; s.ey_out = p.r3;
|
||||
return s;
|
||||
}
|
||||
|
||||
const float a_in = p.r2 + grow, a_out = p.r3 + grow;
|
||||
const float si = p.r2 / a_in, so = p.r3 / a_out;
|
||||
s.q_in = 1.0f - si * si;
|
||||
s.q_out = 1.0f - so * so;
|
||||
|
||||
// Axis-aligned half-extents of an ellipse with semi-axes (a along u, b across it).
|
||||
const float ux2 = s.ux * s.ux, uy2 = s.uy * s.uy;
|
||||
s.ex_in = sqrtf(a_in * a_in * ux2 + p.r2 * p.r2 * uy2);
|
||||
s.ey_in = sqrtf(a_in * a_in * uy2 + p.r2 * p.r2 * ux2);
|
||||
s.ex_out = sqrtf(a_out * a_out * ux2 + p.r3 * p.r3 * uy2);
|
||||
s.ey_out = sqrtf(a_out * a_out * uy2 + p.r3 * p.r3 * ux2);
|
||||
return s;
|
||||
}
|
||||
|
||||
// Which of the n_kern radial-background kernels this reflection uses: one per whole pixel of growth.
|
||||
// Taken from the stencil, so the kernel and the ring it corrects are built from the same radius -
|
||||
// deriving it from a separately computed r0 let the two engines round differently and pick
|
||||
// different rows.
|
||||
BRAGG_STENCIL_HD int BraggStencilKernelIndex(const BraggStencil &s, int n_kern) {
|
||||
if (n_kern <= 1) return 0;
|
||||
int idx = (int) lroundf(s.grow);
|
||||
if (idx < 0) idx = 0;
|
||||
if (idx >= n_kern) idx = n_kern - 1;
|
||||
return idx;
|
||||
}
|
||||
|
||||
// The three squared distances a pixel at offset (ddx, ddy) is tested by, plus its radial offset
|
||||
// (which the radial-background curve bins on). `signal` is the plain circular distance.
|
||||
struct BraggStencilDist {
|
||||
float signal, inner, outer, rad;
|
||||
};
|
||||
|
||||
BRAGG_STENCIL_HD BraggStencilDist BraggStencilDistances(const BraggStencil &s, float ddx, float ddy) {
|
||||
BraggStencilDist d;
|
||||
d.rad = ddx * s.ux + ddy * s.uy;
|
||||
d.signal = ddx * ddx + ddy * ddy;
|
||||
const float rad2 = d.rad * d.rad;
|
||||
d.inner = d.signal - s.q_in * rad2;
|
||||
d.outer = d.signal - s.q_out * rad2;
|
||||
return d;
|
||||
}
|
||||
@@ -142,6 +142,8 @@ std::string RugnuxCommandLine(const ProcessConfig &config,
|
||||
std::ostringstream radii;
|
||||
radii << bragg.GetR1() << "," << bragg.GetR2() << "," << bragg.GetR3();
|
||||
add("--integration-radius", radii.str());
|
||||
if (bragg.GetStencilKSigma() > 0.0f)
|
||||
add("--integration-stencil", num(bragg.GetStencilKSigma()));
|
||||
// Background ring: the CLI defaults to the 4 sigma high-side clip, so emit a flag only when the
|
||||
// GUI chose otherwise. The trim is the alternative estimator (it clears the clip), and a clip of
|
||||
// 0 with no trim is the plain ring mean, which needs the flag to be reproduced.
|
||||
|
||||
+23
-1
@@ -158,6 +158,7 @@ void print_usage() {
|
||||
std::cout << " Integration" << std::endl;
|
||||
std::cout << " --bandwidth <num> Relative X-ray bandwidth FWHM (e.g. 0.01 for 1% DMM); default from file or 0" << std::endl;
|
||||
std::cout << " --integration-radius <r> Signal-box radius r1, or r1,r2,r3 (px). One value => r2=r1+2, r3=r1+4" << std::endl;
|
||||
std::cout << " --integration-stencil <k> Push the r2..r3 background ring out by k times the beam's radial streak (bandwidth*Rpx), per reflection (default 0 = a fixed circular ring). A fixed ring otherwise ends up on a streaked reflection's own tails at high resolution and measures them as background. Only the ring moves, and only radially - the r1 signal box stays a circle. Needs --bandwidth: on a monochromatic beam the streak is zero and this does nothing" << std::endl;
|
||||
std::cout << " --integration-high-resolution <num> High resolution limit for prediction/integration. If omitted (or 0), integration extends as far as the detector reaches" << std::endl;
|
||||
std::cout << " --max-hkl <n> Predict reflections with |h|,|k|,|l| <= n. Default: derived per crystal from the refined cell (ceil(longest axis / d_min) + 1), which is the exact bound - set it only to override that" << std::endl;
|
||||
std::cout << " --background-clip <n> High-side clip of the background ring at mean + n*sqrt(mean) (default 4, or 3 when --bandwidth is set; 0 = off). This is the default background estimator - it rejects neighbour cores and zingers without the symmetric trim's Poisson skew bias. Ignored by --integrator boxsum" << std::endl;
|
||||
@@ -211,6 +212,7 @@ enum {
|
||||
OPT_DETECT_BEAM_STOP,
|
||||
OPT_BANDWIDTH,
|
||||
OPT_INTEGRATION_RADIUS,
|
||||
OPT_INTEGRATION_STENCIL,
|
||||
OPT_BACKGROUND_TRIM,
|
||||
OPT_MAX_HKL,
|
||||
OPT_INTEGRATION_HIGH_RES,
|
||||
@@ -328,6 +330,7 @@ static option long_options[] = {
|
||||
{"resolution-shells", required_argument, nullptr, OPT_RESOLUTION_SHELLS},
|
||||
{"bandwidth", required_argument, nullptr, OPT_BANDWIDTH},
|
||||
{"integration-radius", required_argument, nullptr, OPT_INTEGRATION_RADIUS},
|
||||
{"integration-stencil", required_argument, nullptr, OPT_INTEGRATION_STENCIL},
|
||||
{"background-trim", required_argument, nullptr, OPT_BACKGROUND_TRIM},
|
||||
{"max-hkl", required_argument, nullptr, OPT_MAX_HKL},
|
||||
{"integration-high-resolution", required_argument, nullptr, OPT_INTEGRATION_HIGH_RES},
|
||||
@@ -624,7 +627,8 @@ static int RunRugnux(int argc, char **argv) {
|
||||
std::optional<ResolutionCutoffMethod> resolution_cutoff_method; // --resolution-cutoff cc-logistic|off
|
||||
std::optional<double> resolution_cc_target; // --resolution-cc-target
|
||||
std::optional<int> report_shell_count; // --resolution-shells
|
||||
std::optional<std::string> integration_radius_arg; // "r1" or "r1,r2,r3"
|
||||
std::optional<std::string> integration_radius_arg;
|
||||
std::optional<std::string> integration_stencil_arg; // --integration-stencil: ring elongation, in sigma
|
||||
std::optional<double> background_trim_arg; // --background-trim: background-ring trimmed-mean fraction
|
||||
std::optional<int64_t> max_hkl_arg; // --max-hkl: half-width of the predicted hkl box
|
||||
std::optional<double> integration_d_min_arg; // --integration-high-resolution; unset = detector reach
|
||||
@@ -971,6 +975,9 @@ static int RunRugnux(int argc, char **argv) {
|
||||
case OPT_INTEGRATION_RADIUS:
|
||||
integration_radius_arg = optarg;
|
||||
break;
|
||||
case OPT_INTEGRATION_STENCIL:
|
||||
integration_stencil_arg = optarg;
|
||||
break;
|
||||
case OPT_BACKGROUND_TRIM:
|
||||
background_trim_arg = parse_double_arg(optarg, "--background-trim", logger);
|
||||
break;
|
||||
@@ -1822,6 +1829,21 @@ static int RunRugnux(int argc, char **argv) {
|
||||
logger.Info("Stills integration radii default to r1=6.0 r2=8.0 r3=12.0 (override with --integration-radius)");
|
||||
}
|
||||
|
||||
if (integration_stencil_arg) {
|
||||
const float k = parse_number_arg<float>(integration_stencil_arg->c_str(), "--integration-stencil",
|
||||
logger, 0.0f, 10.0f);
|
||||
BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
|
||||
bis.StencilKSigma(k);
|
||||
experiment.ImportBraggIntegrationSettings(bis);
|
||||
if (!(k > 0.0f))
|
||||
logger.Info("Background ring left circular (--integration-stencil 0)");
|
||||
else if (!(experiment.GetBandwidthFWHM().value_or(0.0f) > 0.0f))
|
||||
logger.Info("--integration-stencil {:.2f} has no effect without --bandwidth: "
|
||||
"the ring is elongated by the beam's radial streak, which is zero here", k);
|
||||
else
|
||||
logger.Info("Background ring elongated by {:.2f} times the beam's radial streak, per reflection", k);
|
||||
}
|
||||
|
||||
if (integrator_mode) {
|
||||
BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
|
||||
bis.Integrator(*integrator_mode);
|
||||
|
||||
@@ -83,31 +83,39 @@ Scene BuildScene(size_t width, size_t height, int spacing = 60) {
|
||||
DiffractionExperiment MakeExperiment(IntegratorMode mode, std::optional<float> bandwidth_fwhm,
|
||||
float clip_nsigma = 4.0f,
|
||||
bool radial = false,
|
||||
const DetectorSetup &det = DetJF(2)) {
|
||||
const DetectorSetup &det = DetJF(2),
|
||||
float stencil_k = 0.0f,
|
||||
float r1 = 0.0f, float r2 = 0.0f, float r3 = 0.0f) {
|
||||
DiffractionExperiment experiment(det); // DetJF(2) (small) keeps the correctness test fast
|
||||
experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV)
|
||||
.BeamX_pxl(400.0f).BeamY_pxl(400.0f);
|
||||
experiment.BandwidthFWHM(bandwidth_fwhm);
|
||||
BraggIntegrationSettings settings;
|
||||
settings.Integrator(mode);
|
||||
if (r1 > 0.0f)
|
||||
settings.R1(r1).R2(r2).R3(r3);
|
||||
if (clip_nsigma > 0.0f)
|
||||
settings.BackgroundClipNSigma(clip_nsigma);
|
||||
else
|
||||
settings.BackgroundTrimFraction(0.10f);
|
||||
settings.BackgroundRadialCorrection(radial);
|
||||
settings.StencilKSigma(stencil_k);
|
||||
experiment.ImportBraggIntegrationSettings(settings);
|
||||
return experiment;
|
||||
}
|
||||
|
||||
void CompareCpuVsGpu(IntegratorMode mode, std::optional<float> bandwidth_fwhm,
|
||||
float clip_nsigma = 4.0f, bool radial = false) {
|
||||
const DiffractionExperiment experiment = MakeExperiment(mode, bandwidth_fwhm, clip_nsigma, radial);
|
||||
float clip_nsigma = 4.0f, bool radial = false, int spacing = 60,
|
||||
float stencil_k = 0.0f,
|
||||
float r1 = 0.0f, float r2 = 0.0f, float r3 = 0.0f) {
|
||||
const DiffractionExperiment experiment =
|
||||
MakeExperiment(mode, bandwidth_fwhm, clip_nsigma, radial, DetJF(2), stencil_k, r1, r2, r3);
|
||||
const size_t width = experiment.GetXPixelsNum();
|
||||
const size_t height = experiment.GetYPixelsNum();
|
||||
const size_t npixel = experiment.GetPixelsNum();
|
||||
REQUIRE(npixel == width * height);
|
||||
|
||||
const Scene scene = BuildScene(width, height);
|
||||
const Scene scene = BuildScene(width, height, spacing);
|
||||
REQUIRE(scene.image.size() == npixel);
|
||||
REQUIRE(scene.predicted.size() > 60);
|
||||
|
||||
@@ -154,12 +162,52 @@ TEST_CASE("BraggIntegrationEngineGPU_MatchesCPU") {
|
||||
SECTION("BoxSum") { CompareCpuVsGpu(IntegratorMode::BoxSum, std::nullopt); }
|
||||
SECTION("ProfileGaussian mono") { CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt); }
|
||||
SECTION("ProfileGaussian broadband") { CompareCpuVsGpu(IntegratorMode::ProfileGaussian, 0.03f); }
|
||||
// An elongated background ring: the classification, the bounding box, the neighbour mask and the
|
||||
// shared-memory radial window all become reflection-dependent, and the two engines have to agree
|
||||
// on every one of them. Spots spaced wider so the grown rings stay clear of the neighbours -
|
||||
// what is under test is the stencil, not the crowding.
|
||||
SECTION("ProfileGaussian stencil broadband") {
|
||||
CompareCpuVsGpu(IntegratorMode::ProfileGaussian, 0.005f, 4.0f, false, 120, 3.0f);
|
||||
}
|
||||
// A monochromatic beam has no streak, so k_sigma changes nothing - the point of the section is
|
||||
// that both engines agree that it changes nothing.
|
||||
SECTION("ProfileGaussian stencil mono") {
|
||||
CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 120, 3.0f);
|
||||
}
|
||||
// Crowded: at the default spacing the grown rings DO overlap their neighbours, so the elongated
|
||||
// neighbour mask, the shrinking background-pixel count and the n_bkg acceptance gate are all in
|
||||
// play. That is the case the feature meets at high resolution, and the wide-spacing sections
|
||||
// above deliberately avoid it.
|
||||
SECTION("ProfileGaussian stencil crowded") {
|
||||
CompareCpuVsGpu(IntegratorMode::ProfileGaussian, 0.02f, 4.0f, false, 60, 4.0f);
|
||||
}
|
||||
SECTION("BoxSum stencil") {
|
||||
CompareCpuVsGpu(IntegratorMode::BoxSum, 0.005f, 4.0f, false, 120, 3.0f);
|
||||
}
|
||||
// The trimmed-mean ring is sorted in a fixed-size shared buffer on the GPU; an elongated ring
|
||||
// holds more pixels, so both engines have to fall back to the plain mean at the same place.
|
||||
SECTION("ProfileGaussian stencil trim") {
|
||||
CompareCpuVsGpu(IntegratorMode::ProfileGaussian, 0.005f, 0.0f, false, 120, 3.0f);
|
||||
}
|
||||
// A ring wide enough to overflow the GPU's fixed trimmed-mean buffer, so the fallback to the
|
||||
// plain ring mean is exercised - and has to happen in both engines at the same reflection. The
|
||||
// growth cap keeps the default 6/10 ring under the buffer at any bandwidth, so this needs the
|
||||
// wider stills radii to be reachable at all.
|
||||
SECTION("ProfileGaussian stencil trim overflow") {
|
||||
CompareCpuVsGpu(IntegratorMode::ProfileGaussian, 0.04f, 0.0f, false, 120, 4.0f, 6.0f, 8.0f, 12.0f);
|
||||
}
|
||||
SECTION("ProfileEmpirical") { CompareCpuVsGpu(IntegratorMode::ProfileEmpirical, std::nullopt); }
|
||||
SECTION("ProfileGaussian mono trim") { CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 0.0f); }
|
||||
// The radial background curvature correction is computed independently in the two engines
|
||||
// (host loop vs radial_correct kernel), so it needs its own parity coverage.
|
||||
SECTION("BoxSum radial") { CompareCpuVsGpu(IntegratorMode::BoxSum, std::nullopt, 4.0f, true); }
|
||||
SECTION("ProfileGaussian radial") { CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, true); }
|
||||
// With an elongated ring the radial-curvature kernel is a table indexed per reflection, and the
|
||||
// shared window boxsum accumulates the curve in is sized from the widest aperture on the
|
||||
// detector. Both are computed independently in the two engines.
|
||||
SECTION("ProfileGaussian radial stencil") {
|
||||
CompareCpuVsGpu(IntegratorMode::ProfileGaussian, 0.005f, 4.0f, true, 120, 3.0f);
|
||||
}
|
||||
}
|
||||
|
||||
// Hidden ([.]) benchmark: the raison d'etre of the GPU port is < 2 ms/frame (vs ~142 ms on the CPU
|
||||
|
||||
@@ -0,0 +1,160 @@
|
||||
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||||
// SPDX-License-Identifier: GPL-3.0-only
|
||||
|
||||
#include <catch2/catch_test_macros.hpp>
|
||||
#include <catch2/matchers/catch_matchers_floating_point.hpp>
|
||||
|
||||
#include <cmath>
|
||||
|
||||
#include "../image_analysis/bragg_integration/BraggStencil.h"
|
||||
|
||||
namespace {
|
||||
|
||||
BraggStencilParams Params(float k_sigma, float bw_sigma = 0.002f) {
|
||||
BraggStencilParams p;
|
||||
p.beam_x = 400.0f;
|
||||
p.beam_y = 400.0f;
|
||||
p.r2 = 6.0f;
|
||||
p.r3 = 10.0f;
|
||||
p.bw_sigma = bw_sigma;
|
||||
p.k_sigma = k_sigma;
|
||||
p.max_grow = 2.0f * p.r3;
|
||||
return p;
|
||||
}
|
||||
|
||||
} // namespace
|
||||
|
||||
// The whole change rests on this: with no elongation asked for, the three squared distances the
|
||||
// integrator tests against must be the SAME BITS as the plain circular distance used before, so
|
||||
// that every pixel is classified exactly as it was, not merely nearly.
|
||||
TEST_CASE("BraggStencil_ZeroElongationIsExactlyCircular", "[Integration]") {
|
||||
const BraggStencilParams p = Params(0.0f); // a bandwidth, but k_sigma = 0
|
||||
|
||||
for (float py = 0.0f; py < 800.0f; py += 37.0f)
|
||||
for (float px = 0.0f; px < 800.0f; px += 41.0f) {
|
||||
const BraggStencil s = MakeBraggStencil(px, py, p);
|
||||
REQUIRE(s.q_in == 0.0f);
|
||||
REQUIRE(s.q_out == 0.0f);
|
||||
for (int dy = -12; dy <= 12; ++dy)
|
||||
for (int dx = -12; dx <= 12; ++dx) {
|
||||
const auto d = BraggStencilDistances(s, static_cast<float>(dx), static_cast<float>(dy));
|
||||
const float circular = static_cast<float>(dx) * dx + static_cast<float>(dy) * dy;
|
||||
REQUIRE(d.signal == circular);
|
||||
REQUIRE(d.inner == circular);
|
||||
REQUIRE(d.outer == circular);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// A monochromatic beam has no streak, so nothing is elongated whatever k_sigma says - which is what
|
||||
// makes the feature inert on every monochromatic dataset rather than merely small.
|
||||
TEST_CASE("BraggStencil_MonochromaticIsInert", "[Integration]") {
|
||||
const BraggStencilParams p = Params(4.0f, 0.0f);
|
||||
for (float py = 0.0f; py < 800.0f; py += 53.0f)
|
||||
for (float px = 0.0f; px < 800.0f; px += 59.0f) {
|
||||
const BraggStencil s = MakeBraggStencil(px, py, p);
|
||||
REQUIRE(s.grow == 0.0f);
|
||||
REQUIRE(s.q_in == 0.0f);
|
||||
REQUIRE(s.q_out == 0.0f);
|
||||
}
|
||||
}
|
||||
|
||||
// The elongated region really is the ellipse it claims: radial semi-axis r + grow, tangential r.
|
||||
TEST_CASE("BraggStencil_ElongatedSemiAxes", "[Integration]") {
|
||||
const BraggStencilParams p = Params(3.0f);
|
||||
|
||||
for (float py = 120.0f; py < 800.0f; py += 91.0f)
|
||||
for (float px = 120.0f; px < 800.0f; px += 97.0f) {
|
||||
const BraggStencil s = MakeBraggStencil(px, py, p);
|
||||
const float grow = s.grow;
|
||||
REQUIRE(grow > 0.0f);
|
||||
REQUIRE(grow <= p.max_grow);
|
||||
REQUIRE(s.grow == BraggStencilGrow_px(s.r0, p)); // the kernel table indexes on this
|
||||
|
||||
// On the radial axis the inner boundary sits at r2 + grow, the outer at r3 + grow.
|
||||
const auto rad_in = BraggStencilDistances(s, (p.r2 + grow) * s.ux, (p.r2 + grow) * s.uy);
|
||||
const auto rad_out = BraggStencilDistances(s, (p.r3 + grow) * s.ux, (p.r3 + grow) * s.uy);
|
||||
CHECK_THAT(rad_in.inner, Catch::Matchers::WithinRel(p.r2 * p.r2, 1e-4f));
|
||||
CHECK_THAT(rad_out.outer, Catch::Matchers::WithinRel(p.r3 * p.r3, 1e-4f));
|
||||
|
||||
// Across it, at the untouched tangential half-widths r2 and r3. Testing on the exact
|
||||
// tangential axis would be a tautology - rad is 0 there, so q never enters - so the
|
||||
// point that matters is that the SAME offset is inside the region radially and outside
|
||||
// it tangentially. That is the anisotropy, and it fails if q is built from the wrong
|
||||
// radius or from a constant.
|
||||
const float probe = p.r2 + 0.5f * grow;
|
||||
const auto radial_probe = BraggStencilDistances(s, probe * s.ux, probe * s.uy);
|
||||
const auto tangent_probe = BraggStencilDistances(s, -probe * s.uy, probe * s.ux);
|
||||
CHECK(radial_probe.inner < p.r2 * p.r2); // still signal, the ring starts further out
|
||||
CHECK(tangent_probe.inner > p.r2 * p.r2); // already background across the streak
|
||||
const auto tan_in = BraggStencilDistances(s, -p.r2 * s.uy, p.r2 * s.ux);
|
||||
const auto tan_out = BraggStencilDistances(s, -p.r3 * s.uy, p.r3 * s.ux);
|
||||
CHECK_THAT(tan_in.inner, Catch::Matchers::WithinRel(p.r2 * p.r2, 1e-4f));
|
||||
CHECK_THAT(tan_out.outer, Catch::Matchers::WithinRel(p.r3 * p.r3, 1e-4f));
|
||||
}
|
||||
}
|
||||
|
||||
// The bounding boxes the engines scan must contain the regions they classify - a box one pixel too
|
||||
// small silently drops background pixels on one side of every reflection, which no parity test
|
||||
// between two engines making the same mistake would catch.
|
||||
TEST_CASE("BraggStencil_BoundingBoxContainsRegion", "[Integration]") {
|
||||
const BraggStencilParams p = Params(3.0f);
|
||||
const float r2_sq = p.r2 * p.r2, r3_sq = p.r3 * p.r3;
|
||||
|
||||
for (float py = 0.0f; py < 800.0f; py += 53.0f)
|
||||
for (float px = 0.0f; px < 800.0f; px += 59.0f) {
|
||||
const BraggStencil s = MakeBraggStencil(px, py, p);
|
||||
const int span = static_cast<int>(std::ceil(p.r3 + p.max_grow)) + 4;
|
||||
for (int dy = -span; dy <= span; ++dy)
|
||||
for (int dx = -span; dx <= span; ++dx) {
|
||||
const auto d = BraggStencilDistances(s, static_cast<float>(dx), static_cast<float>(dy));
|
||||
const float ax = std::fabs(static_cast<float>(dx)), ay = std::fabs(static_cast<float>(dy));
|
||||
// No slack: the offsets are integers from an exactly centred stencil, so the
|
||||
// extents bound them outright. A tolerance of a pixel here would accept a box
|
||||
// one pixel too small, which is the error this exists to catch.
|
||||
if (d.inner < r2_sq) {
|
||||
INFO("inner region outside its box at " << dx << "," << dy);
|
||||
REQUIRE(ax <= s.ex_in);
|
||||
REQUIRE(ay <= s.ey_in);
|
||||
}
|
||||
if (d.inner >= r2_sq && d.outer < r3_sq) {
|
||||
INFO("ring outside its box at " << dx << "," << dy);
|
||||
REQUIRE(ax <= s.ex_out);
|
||||
REQUIRE(ay <= s.ey_out);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// The growth is capped, so a mis-declared bandwidth cannot run away with the bounding box.
|
||||
TEST_CASE("BraggStencil_GrowthIsCapped", "[Integration]") {
|
||||
BraggStencilParams p = Params(3.0f, 0.5f); // an absurdly declared bandwidth
|
||||
for (float r0 = 0.0f; r0 < 4000.0f; r0 += 17.0f)
|
||||
REQUIRE(BraggStencilGrow_px(r0, p) <= p.max_grow);
|
||||
|
||||
const BraggStencil s = MakeBraggStencil(4000.0f, 4000.0f, p);
|
||||
REQUIRE(s.ex_out <= p.r3 + p.max_grow + 1e-3f);
|
||||
REQUIRE(s.ey_out <= p.r3 + p.max_grow + 1e-3f);
|
||||
}
|
||||
|
||||
// The kernel table is indexed by the growth rounded to whole pixels, so the table has to have a row
|
||||
// for every index any reflection on the detector can produce. An off-by-one here is an out-of-range
|
||||
// read of k_diff - on the GPU, a device-side one.
|
||||
TEST_CASE("BraggStencil_KernelIndexInRange", "[Integration]") {
|
||||
for (const float k : {0.0f, 0.4f, 1.0f, 2.5f, 3.0f, 6.0f}) {
|
||||
const BraggStencilParams p = Params(k);
|
||||
const float r_max = std::hypot(800.0f - p.beam_x, 800.0f - p.beam_y);
|
||||
const int n_kern = static_cast<int>(std::lround(BraggStencilGrow_px(r_max, p))) + 1;
|
||||
REQUIRE(n_kern >= 1);
|
||||
for (float py = 0.0f; py <= 800.0f; py += 13.0f)
|
||||
for (float px = 0.0f; px <= 800.0f; px += 17.0f) {
|
||||
const BraggStencil s = MakeBraggStencil(px, py, p);
|
||||
const int idx = BraggStencilKernelIndex(s, n_kern);
|
||||
INFO("k " << k << " at " << px << "," << py << " grow " << s.grow);
|
||||
REQUIRE(idx >= 0);
|
||||
REQUIRE(idx < n_kern);
|
||||
// The clamp must never be what saves it: the table is sized so the row exists.
|
||||
REQUIRE(static_cast<int>(std::lround(s.grow)) == idx);
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -38,6 +38,7 @@ ADD_EXECUTABLE(jfjoch_test
|
||||
ROIIntegrationGPUTest.cpp
|
||||
BraggIntegrationEngineGPUTest.cpp
|
||||
BraggIntegrationEngineCompressedImageTest.cpp
|
||||
BraggStencilTest.cpp
|
||||
LossyFilterTest.cpp
|
||||
ImageBufferTest.cpp
|
||||
PixelMaskTest.cpp
|
||||
|
||||
Reference in New Issue
Block a user