diff --git a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md index 7528a8f6c..ae4bdfd32 100644 --- a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md +++ b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md @@ -151,9 +151,9 @@ $ $ 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 non-signal variance all pick it up. The variance carried to the merge, $\mathrm{var}_\mathrm{bkg} = g\,(n_S\hat b + n_S^2\,\mathrm{var}(\hat b))$, moves with the corrected $\hat b$ in its first term — the Poisson share of the background *under the disk* — while $\mathrm{var}(\hat b)$, the noise of the ring estimate, and $\sigma$, taken from the counts in the disk, do not depend on a smooth curve and stay. -The curve $\bar B(r)$ is the adaptive spot finder's peak-excluded background per resolution ring of the **whole frame** (§3.2), one value per pixel of radius: every pixel of the frame, against the few dozen to few thousand reflections whose own rings would otherwise sample it. Where that curve is in use the ring's **high-side clip follows it**: a ring pixel at radial offset $u$ is clipped against the threshold plus the curve's rise $\bar B(r_0+u)-\bar B(r_0)$, in whole counts. On a background falling steeply away from the beam a flat threshold cuts the beam-side half of the ring as outliers and the background reads low — up to 150× too strong a low-resolution reflection on a binned CCD near the beam stop. Without a frame curve (fixed-threshold spot finding) the curve is accumulated from the clipped annulus pixels of the predicted reflections themselves, and the clip stays flat. +The curve $\bar B(r)$ is the adaptive spot finder's peak-excluded background per resolution ring of the **whole frame** (§3.2), one value per pixel of radius: every pixel of the frame, against the few dozen to few thousand reflections whose own rings would otherwise sample it. Where that curve is in use the ring's **high-side clip follows it**: a ring pixel at radial offset $u$ is clipped against the threshold plus the curve's rise $\bar B(r_0+u)-\bar B(r_0)$, in whole counts. On a background falling steeply away from the beam a flat threshold cuts the beam-side half of the ring as outliers and the background reads low — up to 150× too strong a low-resolution reflection on a binned CCD near the beam stop. With the frame's curve the correction is not the kernel above but the curve itself averaged over the pixels the reflection actually reads: its mean over the ring pixels the clip kept minus its mean over the readable disk pixels, summed while the box is read anyway. The kernel is an average over a whole annulus, so it counts the curve wherever the ring has lost pixels — to the clip, to a neighbour's signal, to a masked or shadowed pixel — and next to the beam-stop shadow, where the curve falls steeply into the shadow and the ring's beam-side part is masked, that turned a correction of a few counts into tens and made low-resolution reflections near the stop read strongly negative. The curve is read at the $2\theta$ of each reflection's position — the radius $F\tan 2\theta$ it is tabulated against, read off the geometry as the curve's rings are — and not at its distance from the PONI, which on a tilted detector is a few pixels off exactly where the curve is steepest. Without a frame curve (fixed-threshold spot finding) the curve is accumulated from the clipped annulus pixels of the predicted reflections themselves, the kernel applies it, and the clip stays flat. -**When it is applied.** `--background-radial` takes `on`, `off` or `auto`; in Rugnux it is **`auto` by default** (the broker keeps it off), and auto is decided **per run, in the pre-scan, on empty apertures** — the same measurement that, on stills, also chooses how far the ring reaches (the $r_3$, up to the one the run starts from, whose ring predicts the empty disk best). On a grid of apertures away from every found spot, the ring's prediction of the empty disk is taken twice: plain, and with the correction applied exactly as the integrator applies it (the clip following the curve, the kernel's correction). With nothing in the aperture each difference is the error the ring makes in a reflection's background, noise plus bias. Rotation measures every stencil the spot width can give it (and the one it starts from, which the geometry pre-pass and the starvation guard return to) and reads the one it gets. +**When it is applied.** `--background-radial` takes `on`, `off` or `auto`; in Rugnux it is **`auto` by default** (the broker keeps it off), and auto is decided **per run, in the pre-scan, on empty apertures** — the same measurement that, on stills, also chooses how far the ring reaches (the $r_3$, up to the one the run starts from, whose ring predicts the empty disk best). On a grid of apertures away from every found spot, the ring's prediction of the empty disk is taken twice: plain, and with the correction applied exactly as the integrator applies it (the clip following the curve, and the curve's mean over the ring pixels kept minus its mean over the disk). With nothing in the aperture each difference is the error the ring makes in a reflection's background, noise plus bias. Rotation measures every stencil the spot width can give it (and the one it starts from, which the geometry pre-pass and the starvation guard return to) and reads the one it gets. - **Stills** have no other correction, so it is applied to every image of a run where it lowers the apertures' mean square error at all. - **Rotation** has one: under auto, each image whose peak-excluded ice score (§3.3) reaches `--ice-min-score` is corrected with the curve of its reflections' rings. The frame's curve on every image replaces that only where it takes **at least 10 %** of the apertures' mean square error away. On photon-counting sweeps the apertures gain 0.1–5 % (statistically significant, $z\approx3$–4), and correcting every image there was measured to cost ISa and anomalous signal; where it pays — the rings of home-source CCD frames, a smooth ice ring — they gain 13–31 %. A run below the margin keeps the per-image ice correction, bit for bit. diff --git a/image_analysis/bragg_integration/BackgroundAperture.cpp b/image_analysis/bragg_integration/BackgroundAperture.cpp index 04d107f5f..e731d64fe 100644 --- a/image_analysis/bragg_integration/BackgroundAperture.cpp +++ b/image_analysis/bragg_integration/BackgroundAperture.cpp @@ -101,17 +101,17 @@ void background_aperture::MeasureFrame(ImagePreprocessorBuffer &image, int width // Every candidate ring as the integrator reads it: the mean, then one high-side clip pass. With the // radial correction, as the integrator applies it: the clip follows the curve (each pixel's // threshold raised by the curve's rise from the centre to it, in whole counts - - // RadialClipOffset), and the correction is taken over the ring's geometry, as the kernel is. + // RadialClipOffset), and the correction is the curve's mean over the ring pixels the clip kept + // minus its mean over the disk (BraggCurveCorrection). std::vector d_plain(r3.size()), d_radial(r3.size()); bool measured = true; for (size_t k = 0; k < r3.size() && measured; ++k) { const float r3_sq = r3[k] * r3[k]; - double sum = 0.0, sum_b = 0.0; + double sum = 0.0; int n = 0; for (const auto &p : ring) if (p.d2 < r3_sq) { sum += p.v; - sum_b += p.b; ++n; } if (n == 0) { @@ -120,7 +120,7 @@ void background_aperture::MeasureFrame(ImagePreprocessorBuffer &image, int width } const double mean = sum / n; const double thr = mean + clip_nsigma * std::sqrt(gain * std::max(mean, static_cast(gain))); - double kept = 0.0, kept_radial = 0.0; + double kept = 0.0, kept_radial = 0.0, kept_radial_b = 0.0; int n_kept = 0, n_kept_radial = 0; for (const auto &p : ring) { if (p.d2 >= r3_sq) @@ -131,13 +131,14 @@ void background_aperture::MeasureFrame(ImagePreprocessorBuffer &image, int width } if (clip_nsigma <= 0.0f || p.v <= thr + static_cast(std::lround(p.b - centre_b))) { kept_radial += p.v; + kept_radial_b += p.b; ++n_kept_radial; } } measured = n_kept > 0 && n_kept_radial > 0; if (measured) { d_plain[k] = disk_v - kept / n_kept; - d_radial[k] = disk_v - kept_radial / n_kept_radial - (disk_b - sum_b / n); + d_radial[k] = disk_v - kept_radial / n_kept_radial - (disk_b - kept_radial_b / n_kept_radial); } } if (!measured) diff --git a/image_analysis/bragg_integration/BackgroundAperture.h b/image_analysis/bragg_integration/BackgroundAperture.h index d73ea4897..55db18292 100644 --- a/image_analysis/bragg_integration/BackgroundAperture.h +++ b/image_analysis/bragg_integration/BackgroundAperture.h @@ -34,7 +34,8 @@ // The same apertures also say whether the radial background-curvature correction helps: each difference // is taken a second time with the correction applied as the integrator applies it - the ring's clip // following the frame's own ring background (RingBackground), and that background averaged over the -// ring's pixels minus over the disk's, which is what the integrator's kernel takes off a reflection. +// ring pixels the clip kept minus over the disk's, which is what the integrator takes off a reflection +// (BraggCurveCorrection). // Where the background really is a function of the radius the corrected ring predicts the disk better; // where it is not - a strongly polarised beam, panels of different response - it predicts it worse, and // the run measures which. On stills the ring is chosen on the plain differences first, so its reach does diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index 3a7d84f9b..4b4430808 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -302,18 +302,34 @@ public: // a few dozen reflections, against every pixel of the frame here. An empty profile goes back to // accumulating the curve from the reflections. void RadialBackgroundFromRings(const std::vector &ring_bkg, int32_t q_bins, float low_q, float dq); + // Where a reflection sits on that curve. The curve is tabulated against F_px tan(2 theta), so this is + // the 2 theta of the reflection's position, read off the geometry as the curve's rings and the empty + // apertures read every pixel's; a pixel of its stencil is at that radius plus its radial offset. Not + // the distance from the PONI (BraggStencil::r0): on a tilted detector the beam lands elsewhere, and + // near the beam stop, where the curve is steepest, those few pixels put a reflection on the wrong + // part of it. + [[nodiscard]] float CurveRadius(const BraggPanel &pn) const { + const double s = geom.GetWavelength_A() / (2.0 * geom.PxlToRes(pn.plane_x, pn.plane_y)); + return static_cast(F_px * std::tan(2.0 * std::asin(s))); + } // Where that curve is in use (the correction on), the ring's high-side clip follows it as well: a ring - // pixel at radial offset `rad` from a reflection at radius r0 is clipped against the threshold plus - // the curve's rise from r0 to it, rounded to whole counts so that every engine keeps the same pixels. - // On a background falling steeply away from the beam the flat threshold otherwise cuts the beam-side - // half of the ring as outliers, and the background reads low. 0 where no curve is in use. - [[nodiscard]] int32_t RadialClipOffset(float r0, float rad) const { + // pixel at radial offset `rad` from a reflection at curve radius rc (CurveRadius) is clipped against + // the threshold plus the curve's rise from rc to it, rounded to whole counts so that every engine + // keeps the same pixels. On a background falling steeply away from the beam the flat threshold + // otherwise cuts the beam-side half of the ring as outliers, and the background reads low. 0 where no + // curve is in use. + [[nodiscard]] int32_t RadialClipOffset(float rc, float rad) const { if (!bkg_radial || radial_curve.empty()) return 0; + return static_cast(std::lround(CurveAt(rc + rad) - CurveAt(rc))); + } + // The curve at a radius, to the nearest bin, and the same in whole 1/BRAGG_CURVE_FIXED counts. + [[nodiscard]] float CurveAt(float r) const { const int n = static_cast(radial_curve.size()); - const int b = std::clamp(static_cast(std::lround(r0 + rad)), 0, n - 1); - const int b0 = std::clamp(static_cast(std::lround(r0)), 0, n - 1); - return static_cast(std::lround(radial_curve[b] - radial_curve[b0])); + return radial_curve[std::clamp(static_cast(std::lround(r)), 0, n - 1)]; + } + [[nodiscard]] int64_t CurveFixed(float r) const { + return static_cast(CurveAt(r) * BRAGG_CURVE_FIXED); } [[nodiscard]] const std::vector &GetRadialBackground() const { return radial_curve; } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 205621eb8..a1b6583fe 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -259,6 +259,12 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, BraggBoxSums s{}; bkg_vals.clear(); clip_ring.clear(); + // Where the frame's curve is in use: the reflection's place on it, and the curve summed over the + // readable disk pixels and over the ring pixels the background is taken from + // (BraggCurveCorrection). + const bool curve = bkg_radial && !radial_curve.empty(); + const float rc = curve ? CurveRadius(pn) : st.r0; + int64_t curve_disk = 0, curve_ring = 0, curve_kept = 0; for (int y = y0; y <= y1; ++y) for (int x = x0; x <= x1; ++x) { const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y); @@ -292,6 +298,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, s.x_sum += x; s.y_sum += y; ++s.n_inner_valid; + if (curve) curve_disk += CurveFixed(rc + d.rad); } else if (d.inner >= r2_sq && d.outer < r3_sq) { if (!valid(px)) continue; ++s.n_ring; @@ -301,6 +308,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, s.bkg_sum += px; if (bkg_trim > 0.0f) bkg_vals.push_back(px); ++s.n_bkg; + if (curve) curve_ring += CurveFixed(rc + d.rad); } } @@ -310,9 +318,10 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, if (ClipsRing(s)) { const int32_t thr = ClipThreshold(s); for (const auto &[px, rad] : clip_ring) - if (static_cast(px) <= static_cast(thr) + RadialClipOffset(st.r0, rad)) { + if (static_cast(px) <= static_cast(thr) + RadialClipOffset(rc, rad)) { clip.sum += px; ++clip.n; + if (curve) curve_kept += CurveFixed(rc + rad); if (bkg_radial && radial_curve.empty()) { // The radial curve is binned on the TRUE detector radius, so the // offset here is the unshrunk radial projection. @@ -323,6 +332,11 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, } } FinishPassA(a, i, r, st, s, clip, bkg_vals); + if (curve) { + const bool kept = ClipsRing(s) && clip.n > 5; // the pixels FinishPassA took the background from + a.rough[i].curve_corr = BraggCurveCorrection(kept ? curve_kept : curve_ring, kept ? clip.n : s.n_bkg, + curve_disk, static_cast(s.n_inner_valid)); + } } return ProfilePasses(img, owner, predicted, npredicted, panels, a, image_number); @@ -609,24 +623,23 @@ std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler & // --- Radial background curvature correction. The annulus mean is blind to the curvature of the // radial background (a linear background cancels between the concentric disk and annulus), so - // correct it by mean_annulus(B) - mean_disk(B) taken from the curve just accumulated. Reads no - // pixels: one short dot product per reflection. --- + // correct it by mean_ring(B) - mean_disk(B). With the frame's curve that was taken in pass A over + // the pixels themselves (BraggCurveCorrection); with the curve just accumulated from the + // reflections it is the kernel, one short dot product per reflection. --- if (bkg_radial) { for (size_t i = 0; i < npredicted; ++i) { auto &rh = rough[i]; 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(rh.k_bin) * k_len; - double corr = 0.0; - for (int k = 0; k < k_len; ++k) { - const int b = std::clamp(rh.r_bin + k - k_off, 0, a.n_rad - 1); - double v; - if (!radial_curve.empty()) - v = radial_curve[b]; - else - v = a.rad_cnt[b] > 0 ? static_cast(a.rad_sum[b]) / a.rad_cnt[b] : rh.bkg; - corr += static_cast(kern[k]) * v; + double corr = rh.curve_corr; + if (radial_curve.empty()) { + const float *kern = k_diff.data() + static_cast(rh.k_bin) * k_len; + for (int k = 0; k < k_len; ++k) { + const int b = std::clamp(rh.r_bin + k - k_off, 0, a.n_rad - 1); + const double v = a.rad_cnt[b] > 0 ? static_cast(a.rad_sum[b]) / a.rad_cnt[b] : rh.bkg; + corr += static_cast(kern[k]) * v; + } } rh.bkg -= corr; // annulus mean -> mean over the signal disk rh.I = static_cast(rh.I_sum) - static_cast(rh.n_inner) * rh.bkg; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h index 5e121a800..9b47b390f 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h @@ -105,6 +105,7 @@ protected: int n_disk = 0, n_own = 0; // signal-disk pixels, and how many of them are this reflection's 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 + double curve_corr = 0.0; // the correction with the frame's curve (BraggCurveCorrection) 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 diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 78a85d108..34e242fab 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -41,6 +41,12 @@ struct BraggGpuParams { __device__ inline bool valid(int32_t v) { return v != INT32_MIN && v != INT32_MAX; } +// The frame's curve at a radius, to the nearest bin, in whole 1/BRAGG_CURVE_FIXED counts +// (BraggIntegrationEngine::CurveFixed). +__device__ inline long long curve_fixed(const float *curve, int n, float r) { + return (long long) (curve[min(max((int) lroundf(r), 0), n - 1)] * BRAGG_CURVE_FIXED); +} + // 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; @@ -169,7 +175,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p uint8_t *hasobs_o, uint8_t *overloaded_o, int *saturated_o, unsigned long long *invd2mm, float *isum_o, int *ninner_o, int *rbin_o, int *kbin_o, unsigned long long *rad_sum, int *rad_cnt, int n_rad, int accumulate_rad, - const float *rad_curve, unsigned long long *counts) { + const float *rad_curve, const float *curve_r, float *curvecorr_o, + unsigned long long *counts) { const int i = blockIdx.x; if (i >= n) return; @@ -184,6 +191,11 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p __shared__ long long s_thr_i; __shared__ int s_clipn; __shared__ float s_r0; + // Where the frame's curve is in use: the reflection's place on it (CurveRadius), and the curve + // summed over the readable disk pixels, over the ring pixels behind the background, and over those + // the clip keeps - in whole 1/BRAGG_CURVE_FIXED counts, so the sums are exact (BraggCurveCorrection). + __shared__ float s_rc; + __shared__ unsigned long long s_cdisk, s_cring, s_ckept; // 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. @@ -200,6 +212,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p if (threadIdx.x == 0) { const float rx = pn.plane_x - p.beam_x, ry = pn.plane_y - p.beam_y; s_r0 = sqrtf(rx * rx + ry * ry); + s_rc = rad_curve ? curve_r[i] : s_r0; + s_cdisk = 0; s_cring = 0; s_ckept = 0; s_Isum = 0; s_Ix = 0; s_Iy = 0; s_x = 0; s_y = 0; s_ninner = 0; s_ninner_valid = 0; s_nbkg = 0; s_bkgsum = 0; s_ndisk = 0; s_nown = 0; s_nring = 0; s_nclean = 0; s_nsat = 0; s_nmax = 0; @@ -221,6 +235,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p long long l_x = 0, l_y = 0; // positions behind l_Isum, for the background-free centroid int l_ni = 0, l_niv = 0, l_nb = 0, l_nd = 0, l_no = 0, l_nring = 0, l_nclean = 0, l_nsat = 0, l_nmax = 0; long long l_bkg = 0; + long long l_cdisk = 0, l_cring = 0; for (int t = threadIdx.x; t < area; t += blockDim.x) { const int x = x0 + t % bw, y = y0 + t / bw; const BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy); @@ -242,7 +257,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p if (px == INT32_MAX || (px == INT32_MIN && on && !static_mask[y * p.W + x])) ++l_nsat; if (px == INT32_MAX) ++l_nmax; if (valid(px)) { l_Isum += px; l_Ix += (long long) x * px; l_Iy += (long long) y * px; - l_x += x; l_y += y; ++l_niv; } + l_x += x; l_y += y; ++l_niv; + if (rad_curve) l_cdisk += curve_fixed(rad_curve, n_rad, s_rc + d.rad); } } else if (d.inner >= p.r2_sq && d.outer < p.r3_sq) { if (!BraggPanelHolds(pn, x, y)) continue; const int32_t px = img[y * p.W + x]; @@ -252,6 +268,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p if (m == 0) ++l_nclean; if (m == MASK_FLUX) continue; l_bkg += px; ++l_nb; + if (rad_curve) l_cring += curve_fixed(rad_curve, n_rad, s_rc + d.rad); } } WARP_ATOMIC_ADD(s_Isum, (unsigned long long) l_Isum); @@ -260,6 +277,10 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p WARP_ATOMIC_ADD(s_x, (unsigned long long) l_x); WARP_ATOMIC_ADD(s_y, (unsigned long long) l_y); WARP_ATOMIC_ADD(s_ninner, l_ni); + if (rad_curve) { + WARP_ATOMIC_ADD(s_cdisk, (unsigned long long) l_cdisk); + WARP_ATOMIC_ADD(s_cring, (unsigned long long) l_cring); + } WARP_ATOMIC_ADD(s_ninner_valid, l_niv); WARP_ATOMIC_ADD(s_nbkg, l_nb); WARP_ATOMIC_ADD(s_bkgsum, (unsigned long long) l_bkg); @@ -352,7 +373,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p // Second ring pass for the high-side sigma-clip (re-reads the annulus; avoids storing bkg values). if (s_accept && p.bkg_clip_nsigma > 0.0f && !do_trim) { - long long c_l = 0; int cn_l = 0; + long long c_l = 0, cc_l = 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 BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy); @@ -363,12 +384,13 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p // Where the frame's radial curve is in use the clip follows it (RadialClipOffset). long long lim = s_thr_i; if (rad_curve) { - const int b = min(max((int) lroundf(s_r0 + d.rad), 0), n_rad - 1); - const int b0 = min(max((int) lroundf(s_r0), 0), n_rad - 1); + const int b = min(max((int) lroundf(s_rc + d.rad), 0), n_rad - 1); + const int b0 = min(max((int) lroundf(s_rc), 0), n_rad - 1); lim += lroundf(rad_curve[b] - rad_curve[b0]); } if ((long long) px <= lim) { c_l += px; ++cn_l; + if (rad_curve) cc_l += curve_fixed(rad_curve, n_rad, s_rc + d.rad); if (accumulate_rad) { // The radial curve is binned on the TRUE detector radius, so this is the // unshrunk radial projection - no per-pixel sqrt. @@ -381,6 +403,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p } } WARP_ATOMIC_ADD(s_clipsum, (unsigned long long) c_l); WARP_ATOMIC_ADD(s_clipn, cn_l); + if (rad_curve) WARP_ATOMIC_ADD(s_ckept, (unsigned long long) cc_l); } __syncthreads(); @@ -445,6 +468,12 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p isum_o[i] = (float) Isum; ninner_o[i] = s_ninner_valid; rbin_o[i] = min(max((int) lroundf(s_r0), 0), n_rad > 0 ? n_rad - 1 : 0); + if (rad_curve) { + // From the ring pixels the background above was taken from (see the CPU engine). + const bool kept = p.bkg_clip_nsigma > 0.0f && s_clipn > 5; + curvecorr_o[i] = (float) BraggCurveCorrection((long long) (kept ? s_ckept : s_cring), kept ? s_clipn : s_nbkg, + (long long) s_cdisk, s_ninner_valid); + } kbin_o[i] = BraggStencilKernelIndex(st, p.n_kern); obsx_o[i] = (float) ox; obsy_o[i] = (float) oy; hasobs_o[i] = hasobs; overloaded_o[i] = s_nsat > 0 ? 1 : 0; @@ -610,9 +639,11 @@ __global__ void build_profiles(const unsigned long long *shell_grid, const unsig // 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. The kernel depends on how far this -// reflection's ring was elongated, so k_diff is a table and kbin_a picks the row. +// reflection's ring was elongated, so k_diff is a table and kbin_a picks the row. With the +// frame's curve handed over, boxsum already took the correction over the pixels themselves +// (curve_corr, BraggCurveCorrection) and only its application is left here. // Mirrors BraggIntegrationEngineCPU. --- -__global__ void radial_correct(const unsigned long long *rad_sum, const int *rad_cnt, const float *curve, int n_rad, +__global__ void radial_correct(const unsigned long long *rad_sum, const int *rad_cnt, const float *curve_corr, 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, float pixel_gain, @@ -620,14 +651,18 @@ __global__ void radial_correct(const unsigned long long *rad_sum, const int *rad 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; + // With the frame's curve the correction was taken over the pixels themselves in boxsum. 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 = curve ? curve[b] - : rad_cnt[b] > 0 ? (float) (long long) rad_sum[b] / (float) rad_cnt[b] : bkg; - corr += kern[k] * v; + if (curve_corr) + corr = curve_corr[i]; + else { + const float *kern = k_diff + (size_t) kbin_a[i] * k_len; + 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 ? (float) (long long) rad_sum[b] / (float) rad_cnt[b] : bkg; + corr += kern[k] * v; + } } const float bkg_new = bkg - corr; bkg_o[i] = bkg_new; @@ -952,6 +987,8 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) { d_isum = CudaDevicePtr(new_capacity); d_ninner = CudaDevicePtr(new_capacity); d_rbin = CudaDevicePtr(new_capacity); + d_curve_corr = CudaDevicePtr(new_capacity); + d_curve_r = CudaDevicePtr(new_capacity); d_kbin = CudaDevicePtr(new_capacity); d_obs_x = CudaDevicePtr(new_capacity); d_obs_y = CudaDevicePtr(new_capacity); @@ -964,6 +1001,7 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) { h_px_x = CudaHostPtr(new_capacity); h_px_y = CudaHostPtr(new_capacity); h_d = CudaHostPtr(new_capacity); + h_curve_r = CudaHostPtr(new_capacity); h_panel = CudaHostPtr(new_capacity); h_mark = CudaHostPtr(new_capacity); h_I = CudaHostPtr(new_capacity); @@ -1007,6 +1045,12 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu cuda_err(cudaMemcpyAsync(d_d, h_d.get(), sizeof(float) * npredicted, cudaMemcpyHostToDevice, *stream)); cuda_err(cudaMemcpyAsync(d_panel, h_panel.get(), sizeof(BraggPanel) * npredicted, cudaMemcpyHostToDevice, *stream)); cuda_err(cudaMemcpyAsync(d_mark, h_mark.get(), sizeof(uint8_t) * npredicted, cudaMemcpyHostToDevice, *stream)); + // Where the frame's curve is in use, where each reflection sits on it (CurveRadius). + if (bkg_radial && !radial_curve.empty()) { + for (size_t i = 0; i < npredicted; ++i) + h_curve_r[i] = CurveRadius(panels[i]); + cuda_err(cudaMemcpyAsync(d_curve_r, h_curve_r.get(), sizeof(float) * npredicted, cudaMemcpyHostToDevice, *stream)); + } BraggGpuParams p{ .W = static_cast(xpixel), .H = static_cast(ypixel), @@ -1074,14 +1118,14 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu d_ok, d_strong, d_has_obs, d_overloaded, d_saturated, d_invd2, d_isum, d_ninner, d_rbin, d_kbin, d_rad_sum, d_rad_cnt, rad_n, rad_n > 0 && !radial_given, - radial_given ? d_rad_curve.get() : nullptr, d_counts); + radial_given ? d_rad_curve.get() : nullptr, d_curve_r, d_curve_corr, d_counts); cuda_err(cudaGetLastError()); // Correct the flat annulus background for the curvature of the radial background before anything // downstream (profile fit, variance) reads it. if (radial_on) { radial_correct<<<(n + threads - 1) / threads, threads, 0, *stream>>>( - d_rad_sum, d_rad_cnt, radial_given ? d_rad_curve.get() : nullptr, n_rad, d_k_diff, k_len, k_off, + d_rad_sum, d_rad_cnt, radial_given ? d_curve_corr.get() : nullptr, n_rad, d_k_diff, k_len, k_off, d_isum, d_ninner, d_rbin, d_kbin, static_cast(pixel_gain), d_ok, d_bkg, d_I, d_var_bkg, n); cuda_err(cudaGetLastError()); } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h index 502978d72..49d357cd7 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h @@ -56,6 +56,8 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { CudaDevicePtr d_k_diff; CudaDevicePtr d_rad_cnt; CudaDevicePtr d_rad_curve; // the curve handed over (RadialBackgroundFromRings) + CudaDevicePtr d_curve_corr; // its correction per reflection (BraggCurveCorrection) + CudaDevicePtr d_curve_r; // where each reflection sits on it (CurveRadius) // --- fixed-size device arrays --- // The learning/fit math is single precision: FP64 is heavily throttled on consumer GPUs and the @@ -81,6 +83,7 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { // driver has staged it through a bounce buffer, so eight of them in a row are eight serialised // round-trips rather than eight queued transfers. CudaHostPtr h_px_x, h_px_y, h_d; + CudaHostPtr h_curve_r; // where each reflection sits on the frame's curve (CurveRadius) CudaHostPtr h_panel; CudaHostPtr h_mark; CudaHostPtr h_I, h_sigma, h_bkg, h_var_bkg, h_obs_x, h_obs_y; diff --git a/image_analysis/bragg_integration/BraggStencil.h b/image_analysis/bragg_integration/BraggStencil.h index e30477268..fcfc48072 100644 --- a/image_analysis/bragg_integration/BraggStencil.h +++ b/image_analysis/bragg_integration/BraggStencil.h @@ -290,6 +290,25 @@ struct BraggStencilDist { float signal, inner, outer, rad; }; +// --- The frame's radial background curve (BraggIntegrationEngine::RadialBackgroundFromRings) ------- +// The curve is summed over a reflection's pixels in whole 1/BRAGG_CURVE_FIXED counts: integer sums, +// so they do not depend on the order the pixels are visited in, as the engines' other sums do not. +BRAGG_STENCIL_CONSTANT constexpr float BRAGG_CURVE_FIXED = 256.0f; + +#ifndef __METAL_VERSION__ +// The radial background correction with the frame's curve: the curve's mean over the ring pixels the +// background was taken from, minus its mean over the disk pixels the intensity was summed over, from +// their fixed-point sums. It is what the ring mean got wrong about the disk if the background is the +// curve - taken over the pixels actually read, so it holds whatever the ring lost (the clip, a +// neighbour's signal, a masked or shadowed pixel), where a kernel averaged over a whole annulus does +// not. Nothing to correct on a disk with no readable pixel. +BRAGG_STENCIL_HD double BraggCurveCorrection(long long ring_sum, int ring_n, long long disk_sum, int disk_n) { + if (ring_n <= 0 || disk_n <= 0) + return 0.0; + return ((double) ring_sum / ring_n - (double) disk_sum / disk_n) / BRAGG_CURVE_FIXED; +} +#endif + // --- Ownership of a shared signal region --------------------------------------------------------- // Two reflections whose centres are closer than 2*r1 have signal disks that intersect, and a pixel in // the intersection carries both. It belongs to the NEARER centre. The union mask above cannot say diff --git a/tests/BraggIntegrationEngineCPUTest.cpp b/tests/BraggIntegrationEngineCPUTest.cpp index a8b3e9ac6..c89320728 100644 --- a/tests/BraggIntegrationEngineCPUTest.cpp +++ b/tests/BraggIntegrationEngineCPUTest.cpp @@ -181,7 +181,7 @@ TEST_CASE("BraggIntegrationEngineCPU_RadialCurveFromRings", "[Integration][porta r.h = 1; r.k = 0; r.l = 0; r.predicted_x = bx + 300.0f; r.predicted_y = by; - r.d = 2.0f; + r.d = static_cast(lambda / (2.0 * std::sin(0.5 * std::atan(300.0 / F_px)))); // where the curve has it r.prescaling_corr = 1.0f; r.partiality = 1.0f; const std::vector predicted{r}; @@ -235,7 +235,7 @@ TEST_CASE("BraggIntegrationEngineCPU_ClipFollowsRadialCurve", "[Integration][por r.h = 1; r.k = 0; r.l = 0; r.predicted_x = bx + 200.0f; r.predicted_y = by; - r.d = 2.0f; + r.d = static_cast(lambda / (2.0 * std::sin(0.5 * std::atan(200.0 / F_px)))); r.prescaling_corr = 1.0f; r.partiality = 1.0f; const std::vector predicted{r}; @@ -254,6 +254,77 @@ TEST_CASE("BraggIntegrationEngineCPU_ClipFollowsRadialCurve", "[Integration][por CHECK(std::fabs(background_read(true) - truth) < 20.0); } +// An empty aperture next to the edge of the beam-stop shadow: the background rises towards the beam until +// the shadow cuts it off, and the shadowed pixels are masked, so the ring has lost its beam-side part. The +// correction is taken over the ring pixels actually read (BraggCurveCorrection) and recovers the +// background under the disk (51.7 against 51.5; the plain ring reads 54.3). An average of the curve over +// the whole annulus counts the curve inside the shadow, where the ring has no pixels, and on the same +// edge without the tilt moved the background the wrong way (57 against 50.4). The detector is tilted, so +// the curve has to be read at the reflection's own 2 theta (CurveRadius): the PONI is 17 px from where +// the beam lands, and read at the distance from it the reflection sits on the flat part of the curve and +// nothing is corrected (54.3). +TEST_CASE("BraggIntegrationEngineCPU_RadialCurveOverThePixelsRead", "[Integration][portable]") { + DiffractionExperiment experiment(DetJF(2)); + experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(400.0f).BeamY_pxl(400.0f) + .PoniRot1_rad(0.01f); + const DiffractionGeometry geom = experiment.GetDiffractionGeometry(); + const double F_px = 100.0 / 0.075, lambda = 1.0; + // The radius the curve is tabulated against, F tan(2 theta), of a pixel. + const auto radius = [&](double x, double y) { + return F_px * std::tan(2.0 * std::asin(lambda / (2.0 * geom.PxlToRes(static_cast(x), static_cast(y))))); + }; + const double edge = 296.0; // the shadow's edge + const auto background = [](double r) { return 50.0 + 4.0 * std::max(0.0, 302.0 - r); }; + const size_t width = experiment.GetXPixelsNum(), npixel = experiment.GetPixelsNum(); + ImagePreprocessorBuffer image(npixel); + for (size_t i = 0; i < npixel; ++i) { + const double r = radius(static_cast(i % width), static_cast(i / width)); + image[i] = r < edge ? INT32_MIN : static_cast(std::lround(background(r))); + } + // The rings inside the shadow hold what is left there, little; fine enough that the edge is sharp. + const float low_q = 0.05f, dq = 0.0002f; + std::vector rings(10000); + for (size_t i = 0; i < rings.size(); ++i) { + const double q = low_q + (i + 0.5) * dq; + const double r = F_px * std::tan(2.0 * std::asin(q * lambda / (4.0 * PI))); + rings[i] = static_cast(r < edge ? 10.0 : background(r)); + } + // 304 px out from where the beam lands, on the side away from the PONI. + const auto [dbx, dby] = geom.GetDirectBeam_pxl(); + REQUIRE(std::fabs(dbx - 400.0f) > 10.0f); + Reflection r{}; + r.h = 1; r.k = 0; r.l = 0; + r.predicted_x = std::round(dbx + (dbx > 400.0f ? 304.0f : -304.0f)); + r.predicted_y = std::round(dby); + r.d = geom.PxlToRes(r.predicted_x, r.predicted_y); + r.prescaling_corr = 1.0f; + r.partiality = 1.0f; + const std::vector predicted{r}; + auto integrate = [&](bool correct) { + experiment.ImportBraggIntegrationSettings( + BraggIntegrationSettings().Integrator(IntegratorMode::BoxSum).BackgroundRadialCorrection(correct)); + BraggIntegrationEngineCPU engine(experiment, PixelMask(experiment)); + engine.RadialBackgroundFromRings(rings, static_cast(rings.size()), low_q, dq); + const auto out = engine.Run(image, predicted, predicted.size(), 0); + REQUIRE(out.size() == 1); + return out[0]; + }; + // The background under the disk, from the image itself. + const float r1 = experiment.GetBraggIntegrationSettings().GetR1(); + double truth = 0.0; + int n = 0; + for (int dy = -8; dy <= 8; ++dy) + for (int dx = -8; dx <= 8; ++dx) + if (dx * dx + dy * dy < r1 * r1) { + truth += image[static_cast(r.predicted_y + dy) * width + static_cast(r.predicted_x + dx)]; + ++n; + } + truth /= n; + const Reflection plain = integrate(false), corrected = integrate(true); + CHECK(plain.bkg > truth + 2.0); + CHECK(corrected.bkg == Catch::Approx(truth).margin(0.5)); +} + // A still has no rocking event to drop, so an overloaded observation is kept as a lower bound: the box sum // with every saturated pixel at the saturation limit. A spot with its core clipped at the limit, the // clipped pixels marked saturated as the preprocessor marks them: the still gives exactly the clipped diff --git a/tests/BraggIntegrationEngineGPUTest.cpp b/tests/BraggIntegrationEngineGPUTest.cpp index 889360085..b979522d6 100644 --- a/tests/BraggIntegrationEngineGPUTest.cpp +++ b/tests/BraggIntegrationEngineGPUTest.cpp @@ -496,8 +496,8 @@ TEST_CASE("BraggIntegrationEngineGPU_RadialCurveFromRingsMatchesCPU", "[Integrat ImagePreprocessorBuffer cpu_image(npixel); ImagePreprocessorBufferGPU gpu_image(npixel); for (size_t i = 0; i < npixel; ++i) { - cpu_image[i] = static_cast(std::lround( - background(std::hypot(static_cast(i % width) - bx, static_cast(i / width) - by)))); + const double r = std::hypot(static_cast(i % width) - bx, static_cast(i / width) - by); + cpu_image[i] = r < 296.0 ? INT32_MIN : static_cast(std::lround(background(r))); gpu_image[i] = cpu_image[i]; } auto stream = std::make_shared(); @@ -512,9 +512,14 @@ TEST_CASE("BraggIntegrationEngineGPU_RadialCurveFromRingsMatchesCPU", "[Integrat rings[i] = static_cast(background(F_px * std::tan(2.0 * std::asin(q * lambda / (4.0 * PI))))); } + // Half of them next to a masked disc around the beam - a beam-stop shadow - whose edge their rings + // reach, so the correction is taken over a ring that lost part of itself (BraggCurveCorrection). std::vector predicted; - for (int i = 0; i < 8; ++i) - predicted.push_back(MakeReflection(bx + 300.0f * std::cos(0.7f * i), by + 300.0f * std::sin(0.7f * i), 2.0f, i + 1)); + for (int i = 0; i < 8; ++i) { + const float rr = i % 2 ? 306.0f : 300.0f; + const float d = static_cast(lambda / (2.0 * std::sin(0.5 * std::atan(rr / F_px)))); + predicted.push_back(MakeReflection(bx + rr * std::cos(0.7f * i), by + rr * std::sin(0.7f * i), d, i + 1)); + } for (const auto mode : {IntegratorMode::BoxSum, IntegratorMode::ProfileGaussian}) { experiment.ImportBraggIntegrationSettings(BraggIntegrationSettings().Integrator(mode).BackgroundRadialCorrection(true));