From 0e23fd3ab9b40143631c31c9a0cb63abfd5d6adc Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 6 Aug 2026 15:44:13 +0200 Subject: [PATCH] Bragg integration: propagate the background-estimate uncertainty, add an opt-in radial background correction Two independent pieces in the same code path. The background-estimate variance was never propagated. A reflection's background comes from a finite ring of n_b pixels, so subtracting it adds var(B)/n_b per signal pixel - sqrt(1 + n_d/n_b) = 1.109 with the shipped stencil. Both engines omitted it, which is exactly the 1.11-1.19 gap measured between the off-ring scatter and the reported sigma. Three lines each; it affects every dataset, not only iced ones. The radial correction is new and OFF by default (--background-radial). The signal disk and the background ring are concentric, so for any background LINEAR in position _ann == _disk identically and a plane fit buys nothing; the leading error is the CURVATURE of the radial background, which on a sharp ice ring reaches +26 counts on a single reflection. Since every reflection uses the same stencil, that error is a fixed kernel over radial offset - one short dot product per reflection and no extra pixel reads. Validated on empty apertures before any C++: mean |bias| over 9 bands / 3 crystals 4.33 -> 0.79 counts with the scatter unchanged. Three things it cost a battery each to learn, all now in the code: - the radial curve must be accumulated from CLIPPED annulus pixels, inside the clip pass, or it carries neighbour tails and zingers (so it is inert under --integrator boxsum, which has no clip pass); - the GPU version was a 1.8x slowdown from atomicAdd contention on a small radial array - staged in shared memory per block it now costs nothing measurable; - it is battery-NEUTRAL as a default, because the reflections whose bias it fixes are the ones the ice handling already excludes. Hence off by default. CPU/GPU parity extended with two radial sections: 9002 assertions. Also fixes a latent French-Wilson quadrature collapse: j_max = I + 8 sigma on a fixed 400-point grid degenerates to a single cell once sigma >> 50 , giving F = 0.1 sqrt(sigma) with sigmaF -> 0. Harmless today, but any sigma-inflation scheme detonates it. Co-Authored-By: Claude Opus 5 (1M context) --- common/BraggIntegrationSettings.cpp | 9 ++ common/BraggIntegrationSettings.h | 12 ++ .../BraggIntegrationEngine.cpp | 33 +++++ .../BraggIntegrationEngine.h | 10 ++ .../BraggIntegrationEngineCPU.cpp | 77 +++++++++- .../BraggIntegrationEngineGPU.cu | 138 ++++++++++++++++-- .../BraggIntegrationEngineGPU.h | 9 +- image_analysis/scale_merge/FrenchWilson.cpp | 11 +- rugnux/rugnux_cli.cpp | 22 +++ tests/BraggIntegrationEngineGPUTest.cpp | 12 +- 10 files changed, 309 insertions(+), 24 deletions(-) diff --git a/common/BraggIntegrationSettings.cpp b/common/BraggIntegrationSettings.cpp index e1fad22f..7142f13c 100644 --- a/common/BraggIntegrationSettings.cpp +++ b/common/BraggIntegrationSettings.cpp @@ -139,3 +139,12 @@ BraggIntegrationSettings &BraggIntegrationSettings::BackgroundClipNSigma(float i float BraggIntegrationSettings::GetBackgroundClipNSigma() const { return bkg_clip_nsigma; } + +BraggIntegrationSettings &BraggIntegrationSettings::BackgroundRadialCorrection(bool input) { + bkg_radial_correction = input; + return *this; +} + +bool BraggIntegrationSettings::IsBackgroundRadialCorrection() const { + return bkg_radial_correction; +} diff --git a/common/BraggIntegrationSettings.h b/common/BraggIntegrationSettings.h index 074168f2..fa24c4ac 100644 --- a/common/BraggIntegrationSettings.h +++ b/common/BraggIntegrationSettings.h @@ -47,6 +47,16 @@ class BraggIntegrationSettings { // 10% symmetric trim +5.05..+6.34, 4 sigma clip +0.02..+0.54. Applied to monochromatic data // (rugnux --background-clip); broadband (non-zero bandwidth) data always clip, at their tuned 3 sigma. float bkg_clip_nsigma = 4.0f; + // Radial background curvature correction. The signal disk and the background annulus are + // concentric, so for ANY background linear in position their means are equal - a plane fit buys + // nothing and the leading error is the CURVATURE of the radial background, which the flat annulus + // mean is structurally blind to. Sitting on an ice ring that reaches +26 counts on a single + // reflection. When on, a radial background curve is accumulated per image from the annulus pixels + // that are already read, and each reflection's background is corrected by + // mean_annulus(B) - mean_disk(B), evaluated as a fixed kernel over radial offset (O(1), no extra + // pixel reads). Measured empty-aperture bias over 9 bands on 3 crystals: 4.33 -> 0.79 counts mean + // |bias|, scatter unchanged. + bool bkg_radial_correction = false; // Half-width of the hkl cube the predictor walks: every reflection with |h|,|k|,|l| <= this is // tested against the Ewald sphere, and nothing outside it can ever be predicted. An axis is // truncated once a/d_min exceeds this, and the GPU cost is the cube (2n+1)^3 of candidates, so @@ -67,6 +77,7 @@ public: BraggIntegrationSettings& Integrator(IntegratorMode input); BraggIntegrationSettings& BackgroundTrimFraction(float input); BraggIntegrationSettings& BackgroundClipNSigma(float input); + BraggIntegrationSettings& BackgroundRadialCorrection(bool input); BraggIntegrationSettings& MaxHKL(std::optional input); @@ -80,5 +91,6 @@ public: [[nodiscard]] float GetMinimumSigmaInRegardsToI() const; [[nodiscard]] float GetBackgroundTrimFraction() const; [[nodiscard]] float GetBackgroundClipNSigma() const; + [[nodiscard]] bool IsBackgroundRadialCorrection() const; [[nodiscard]] std::optional GetMaxHKL() const; }; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp index 58a4e2fb..b409549b 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp @@ -4,6 +4,8 @@ #include "BraggIntegrationEngine.h" #include +#include +#include #include namespace { @@ -75,6 +77,37 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe bkg_clip_nsigma = broadband ? 3.0f : settings.GetBackgroundClipNSigma(); bkg_trim = (broadband || bkg_clip_nsigma > 0.0f) ? 0.0f : settings.GetBackgroundTrimFraction(); + // Radial-offset kernels for the background curvature correction. A stencil pixel at (dx, dy) + // sits at radial offset dx*cos(phi) + dy*sin(phi) from the reflection, where phi is the + // reflection's azimuth; averaging over phi makes the kernels position-independent, which is + // exact to the extent the stencil is small against the reflection's radius (r3 = 10 px vs + // hundreds). 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). + bkg_radial = settings.IsBackgroundRadialCorrection(); + k_off = static_cast(std::ceil(r3)) + 1; + k_diff.assign(2 * k_off + 1, 0.0f); + { + std::vector hist_disk(k_diff.size(), 0.0), hist_ann(k_diff.size(), 0.0); + constexpr int n_phi = 512; + const int span = static_cast(std::ceil(r3)) + 1; + for (int p = 0; p < n_phi; ++p) { + const double phi = 2.0 * M_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(dx) * dx + static_cast(dy) * dy; + const int k = k_off + static_cast(std::lround(dx * cp + dy * sp)); + if (k < 0 || k >= static_cast(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(hist_ann[k] / sa - hist_disk[k] / sd); + } + polarization = experiment.GetPolarizationFactor(); } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index a2c39c9a..4db28f57 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -102,6 +102,16 @@ protected: // high-side sigma-clip). 0 = plain ring mean. Read by both the CPU and GPU engines. float bkg_trim = 0.0f; + // --- 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: + // 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. + bool bkg_radial = false; + int k_off = 0; // index of offset 0 in k_diff + std::vector k_diff; // annulus-minus-disk weight per integer radial offset + DiffractionGeometry geom; // kept for the per-reflection polarization correction std::optional polarization; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 2b2c9b14..a55544fe 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -85,6 +85,10 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, // --- Pass A: box-sum every reflection (rough I, background, centroid, strong flag). --- struct Rough { double I = 0.0, sigma = NAN, bkg = 0.0, obs_x = 0.0, obs_y = 0.0; + double bkg_var = 0.0; // variance of the background ESTIMATE itself, bkg / n_bkg + 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 cx = 0, cy = 0, shell = -1; bool ok = false, strong = false, has_obs = false; }; @@ -92,6 +96,18 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, double inv_d2_min = std::numeric_limits::max(), inv_d2_max = 0.0; std::vector bkg_vals; // reused per reflection for the trimmed-mean background (idea 1) + // Radial background curve, accumulated from the annulus pixels this 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 sqrt is needed. + const int n_rad = static_cast(std::ceil(std::hypot(std::max(beam_x, W - beam_x), + std::max(beam_y, H - beam_y)))) + 2; + std::vector rad_sum; + std::vector rad_cnt; + if (bkg_radial) { + rad_sum.assign(n_rad, 0.0); + rad_cnt.assign(n_rad, 0); + } + for (size_t i = 0; i < npredicted; ++i) { const auto &r = predicted[i]; Rough out; @@ -100,6 +116,12 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, const int y0 = std::max(0, static_cast(std::floor(r.predicted_y - r3 - 1.0))); const int y1 = std::min(H - 1, static_cast(std::ceil(r.predicted_y + r3 + 1.0))); + // 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(std::lround(r0)), 0, n_rad - 1); + 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; int n_bkg = 0; @@ -124,6 +146,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, } } + int n_bkg_used = n_bkg; // pixels behind the FINAL background value (trim/clip shrink it) if (n_inner_valid == n_inner && n_bkg > 5) { out.bkg = bkg_sum / n_bkg; if (bkg_trim_frac > 0.0 && bkg_vals.size() > 5) { @@ -137,6 +160,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, double s = 0.0; for (size_t t = lo; t < hi; ++t) s += bkg_vals[t]; out.bkg = s / static_cast(hi - lo); + n_bkg_used = static_cast(hi - lo); } } else if (do_clip) { // One high-outlier sigma-clip pass on the background ring: reject pixels above @@ -151,14 +175,31 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, if (refl_mask[y * W + x]) continue; const int32_t px = img[y * W + x]; if (!valid(px)) continue; - if (static_cast(px) <= thr) { s += px; ++n; } + if (static_cast(px) <= thr) { + 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(std::lround(r0 + off)), 0, n_rad - 1); + rad_sum[b] += static_cast(px); + ++rad_cnt[b]; + } + } } - if (n > 5) out.bkg = s / n; + if (n > 5) { out.bkg = s / n; n_bkg_used = n; } } out.I = static_cast(I_sum) - static_cast(n_inner) * out.bkg; + // I = I_sum - n_inner*bkg, and bkg is itself estimated from n_bkg_used pixels, so its + // error enters n_inner times over: var(I) = I_sum + n_inner^2 * bkg/n_bkg_used. Leaving + // the second term out understates sigma by sqrt(1 + n_inner/n_bkg) - 1.109x at the + // default r1=4/r2=6/r3=10 stencil, on every reflection of every dataset. + out.bkg_var = out.bkg / n_bkg_used; + out.I_sum = I_sum; + out.n_inner = static_cast(n_inner); + const double var_bkg_term = static_cast(n_inner) * n_inner * out.bkg_var; out.sigma = std::max(1.0, out.I * min_sigma_ratio); if (I_sum > 0) { - out.sigma = std::max(out.sigma, std::sqrt(static_cast(I_sum))); + out.sigma = std::max(out.sigma, std::sqrt(static_cast(I_sum) + var_bkg_term)); out.obs_x = static_cast(I_sum_x) / static_cast(I_sum); out.obs_y = static_cast(I_sum_y) / static_cast(I_sum); out.has_obs = true; @@ -176,6 +217,27 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, rough[i] = out; } + // --- 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. --- + 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). + double corr = 0.0; + for (int k = 0; k < static_cast(k_diff.size()); ++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(k_diff[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; + } + } + // --- BoxSum mode is BraggIntegrate2D: emit the rough result directly. --- if (mode == IntegratorMode::BoxSum) { for (size_t i = 0; i < npredicted; ++i) { @@ -317,10 +379,11 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, const int Gf = 2 * Rf + 1; const double B = std::max(rh.bkg, PIXEL_VARIANCE_FLOOR); - double I = rh.I, den = 0.0; + double I = rh.I, den = 0.0, wsum = 0.0; for (int iter = 0; iter < 4; ++iter) { double num = 0.0; den = 0.0; + wsum = 0.0; for (int dy = -Rf; dy <= Rf; ++dy) for (int dx = -Rf; dx <= Rf; ++dx) { const double Pp = (*Pvec)[(dy + Rf) * Gf + (dx + Rf)]; @@ -332,6 +395,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, const double v = B + std::max(0.0, I) * Pp; num += Pp * (static_cast(px) - rh.bkg) / v; den += Pp * Pp / v; + wsum += Pp / v; } if (den > 0.0) I = num / den; else break; } @@ -341,7 +405,10 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, // iteration has no real peak to lock onto and can manufacture intensity the box sum never sees. // Keep the profile intensity only if it agrees with the summation seed within the margin; // otherwise fall back to the summation, which is robust there. - double sigma = std::sqrt(1.0 / den); + // 1/den is the profile-fit variance with the background taken as exact. The fit is + // I = sum(P*(px-bkg)/v) / sum(P^2/v), so dI/dbkg = -wsum/den and the background estimate's + // own error adds (wsum/den)^2 * var(bkg) - the same term the box sum was missing. + double sigma = std::sqrt(1.0 / den + (wsum / den) * (wsum / den) * rh.bkg_var); if (std::abs(I - rh.I) > PROFILE_SUMMATION_MAX_NSIGMA * rh.sigma) { I = rh.I; sigma = rh.sigma; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 42d5131e..71f2ec5b 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -56,8 +56,10 @@ __global__ void mark_mask(const float *px_x, const float *px_y, uint8_t *mask, B __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, const int32_t *img, const uint8_t *mask, BraggGpuParams p, int n, int *cx_o, int *cy_o, float *I_o, float *sigma_o, float *bkg_o, - float *obsx_o, float *obsy_o, uint8_t *ok_o, uint8_t *strong_o, - uint8_t *hasobs_o, unsigned long long *invd2mm) { + float *bkgvar_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 *rad_sum, int *rad_cnt, int n_rad) { const int i = blockIdx.x; if (i >= n) return; @@ -67,6 +69,18 @@ __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__ 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; @@ -107,6 +121,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; 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; @@ -175,23 +191,54 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, 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 ((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) { + atomicAdd(&s_radv[idx], (float) px); // shared, not global + atomicAdd(&s_radn[idx], 1); + } + } + } } atomicAdd(&s_clipsum, c_l); atomicAdd(&s_clipn, cn_l); } __syncthreads(); + if (n_rad > 0) { + for (int t = threadIdx.x; t < 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]); + atomicAdd(&rad_cnt[b], s_radn[t]); + } + __syncthreads(); + } + if (threadIdx.x != 0) return; if (!s_accept) { ok_o[i] = 0; strong_o[i] = 0; hasobs_o[i] = 0; return; } double bkg = s_bkg; - if (p.bkg_clip_nsigma > 0.0f && s_clipn > 5) bkg = s_clipsum / (double) s_clipn; + int n_bkg_used = s_nbkg; // pixels behind the FINAL background value (trim/clip shrink it) + if (do_trim) { + const int nb = min(s_bn, BKG_TRIM_MAX); + const int lo = (int) (nb * p.bkg_trim), hi = nb - lo; + if (hi > lo) n_bkg_used = hi - lo; + } + if (p.bkg_clip_nsigma > 0.0f && s_clipn > 5) { bkg = s_clipsum / (double) s_clipn; n_bkg_used = s_clipn; } const long long Isum = (long long) s_Isum; const double I = (double) Isum - (double) s_ninner * bkg; + // See the CPU engine: bkg is estimated from n_bkg_used pixels and subtracted n_inner times, so + // var(I) = Isum + n_inner^2 * bkg/n_bkg_used. Both engines must agree. + const double bkg_var = bkg / (double) n_bkg_used; + const double var_bkg_term = (double) s_ninner * (double) s_ninner * bkg_var; double sigma = fmax(1.0, I * (double) p.min_sigma_ratio); uint8_t hasobs = 0; double ox = 0.0, oy = 0.0; if (Isum > 0) { - sigma = fmax(sigma, sqrt((double) Isum)); + sigma = fmax(sigma, sqrt((double) Isum + var_bkg_term)); ox = (double) (long long) s_Ix / (double) Isum; oy = (double) (long long) s_Iy / (double) Isum; hasobs = 1; @@ -199,6 +246,10 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, cx_o[i] = (int) lroundf(cx); cy_o[i] = (int) lroundf(cy); I_o[i] = (float) I; sigma_o[i] = (float) sigma; bkg_o[i] = (float) bkg; + bkgvar_o[i] = (float) bkg_var; + 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); 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; @@ -302,18 +353,44 @@ __global__ void build_profiles(const float *shell_grid, const float *global_grid } } +// --- Radial background curvature correction: one thread per reflection, no pixel reads. +// The disk and the annulus are concentric, so a background linear in position cancels between +// 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. --- +__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 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]; + 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; + } + const float bkg_new = bkg - corr; + bkg_o[i] = bkg_new; + I_o[i] = isum_a[i] - (float) ninner_a[i] * bkg_new; +} + // --- Pass B Kabsch profile fit: I = sum P(c-B)/v over sum P^2/v, v = B + max(I,0)P (iterate). // One block per reflection; the (possibly elongated) profile is built in shared memory. --- __global__ void fit(const int32_t *img, const float *px_x, const float *px_y, const int *cx_a, const int *cy_a, const float *dd, const unsigned long long *invd2mm, - const float *I_seed, const float *sigma_seed, const float *bkg_a, const uint8_t *ok_a, + const float *I_seed, const float *sigma_seed, const float *bkg_a, + const float *bkgvar_a, const uint8_t *ok_a, const float *shell_P, const float *global_P, const float *shell_sigma2, const float *global_sigma2, const int *shell_n, float *I_o, float *sigma_o, uint8_t *ok_o, BraggGpuParams p, int n) { const int i = blockIdx.x; if (i >= n) return; extern __shared__ float Pbuf[]; - __shared__ float s_gs, s_num, s_den, s_I; + __shared__ float s_gs, s_num, s_den, s_I, s_wsum; __shared__ int s_Rf, s_Gf; if (!ok_a[i]) { if (threadIdx.x == 0) ok_o[i] = 0; return; } @@ -366,10 +443,10 @@ __global__ void fit(const int32_t *img, const float *px_x, const float *px_y, __syncthreads(); for (int iter = 0; iter < 4; ++iter) { - if (threadIdx.x == 0) { s_num = 0.0f; s_den = 0.0f; } + if (threadIdx.x == 0) { s_num = 0.0f; s_den = 0.0f; s_wsum = 0.0f; } __syncthreads(); const float Ihere = s_I; - float l_num = 0.0f, l_den = 0.0f; + float l_num = 0.0f, l_den = 0.0f, l_wsum = 0.0f; for (int k = threadIdx.x; k < GfGf; k += blockDim.x) { const float Pp = Pbuf[k]; if (Pp <= 0.0f) continue; @@ -380,8 +457,9 @@ __global__ void fit(const int32_t *img, const float *px_x, const float *px_y, const float v = B + fmaxf(0.0f, Ihere) * Pp; l_num += Pp * ((float) px - bkg) / v; l_den += Pp * Pp / v; + l_wsum += Pp / v; } - atomicAdd(&s_num, l_num); atomicAdd(&s_den, l_den); + atomicAdd(&s_num, l_num); atomicAdd(&s_den, l_den); atomicAdd(&s_wsum, l_wsum); __syncthreads(); if (threadIdx.x == 0 && s_den > 0.0f) s_I = s_num / s_den; __syncthreads(); @@ -389,7 +467,10 @@ __global__ void fit(const int32_t *img, const float *px_x, const float *px_y, if (threadIdx.x == 0) { if (s_den > 0.0f) { - float I = s_I, sigma = sqrtf(1.0f / s_den); + // 1/s_den takes the background as exact; dI/dbkg = -s_wsum/s_den adds the background + // estimate's own error (see the CPU engine). + const float wr = s_wsum / s_den; + float I = s_I, sigma = sqrtf(1.0f / s_den + wr * wr * bkgvar_a[i]); // Guard against profile-fit runaways (see the CPU engine): fall back to the summation seed // when the profile result diverges from it. if (fabsf(I - I_seed[i]) > (float) PROFILE_SUMMATION_MAX_NSIGMA * sigma_seed[i]) { @@ -433,6 +514,18 @@ BraggIntegrationEngineGPU::BraggIntegrationEngineGPU(const DiffractionExperiment if (fit_shared_bytes > prop.sharedMemPerBlock) throw JFJochException(JFJochExceptionCategory::GPUCUDAError, "BraggIntegrationEngineGPU: profile grid exceeds shared memory (r2 too large)"); + + // Radial background curve: one bin per pixel of distance from the beam, out to the far corner. + if (bkg_radial) { + const double fx = std::max(beam_x, static_cast(xpixel) - beam_x); + const double fy = std::max(beam_y, static_cast(ypixel) - beam_y); + n_rad = static_cast(std::ceil(std::hypot(fx, fy))) + 2; + d_rad_sum = CudaDevicePtr(n_rad); + d_rad_cnt = CudaDevicePtr(n_rad); + d_k_diff = CudaDevicePtr(k_diff.size()); + cuda_err(cudaMemcpy(d_k_diff, k_diff.data(), sizeof(float) * k_diff.size(), + cudaMemcpyHostToDevice)); + } } void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) { @@ -450,6 +543,10 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) { d_I = CudaDevicePtr(new_capacity); d_sigma = CudaDevicePtr(new_capacity); d_bkg = CudaDevicePtr(new_capacity); + d_bkg_var = CudaDevicePtr(new_capacity); + d_isum = CudaDevicePtr(new_capacity); + d_ninner = CudaDevicePtr(new_capacity); + d_rbin = CudaDevicePtr(new_capacity); d_obs_x = CudaDevicePtr(new_capacity); d_obs_y = CudaDevicePtr(new_capacity); d_ok = CudaDevicePtr(new_capacity); @@ -503,11 +600,24 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu // Pass A: reset accumulators, mask, then box-sum. cuda_err(cudaMemsetAsync(d_mask, 0, npixel, *stream)); + if (n_rad > 0) { + cuda_err(cudaMemsetAsync(d_rad_sum, 0, sizeof(float) * n_rad, *stream)); + cuda_err(cudaMemsetAsync(d_rad_cnt, 0, sizeof(int) * n_rad, *stream)); + } reset<<<32, 256, 0, *stream>>>(d_shell_grid, d_global_grid, d_shell_n, d_global_n, d_invd2, GG); mark_mask<<>>(d_px_x, d_px_y, d_mask, p, n); boxsum<<>>(d_px_x, d_px_y, d_d, img, d_mask, p, n, - d_cx, d_cy, d_I, d_sigma, d_bkg, d_obs_x, d_obs_y, - d_ok, d_strong, d_has_obs, d_invd2); + d_cx, d_cy, d_I, d_sigma, d_bkg, d_bkg_var, d_obs_x, d_obs_y, + d_ok, d_strong, d_has_obs, d_invd2, + d_isum, d_ninner, d_rbin, + d_rad_sum, d_rad_cnt, n_rad); + + // Correct the flat annulus background for the curvature of the radial background before anything + // downstream (profile fit, variance) reads it. + if (n_rad > 0) + radial_correct<<<(n + threads - 1) / threads, threads, 0, *stream>>>( + d_rad_sum, d_rad_cnt, n_rad, d_k_diff, static_cast(k_diff.size()), k_off, + d_isum, d_ninner, d_rbin, d_ok, d_bkg, d_I, n); if (mode != IntegratorMode::BoxSum) { // Pass B: learn (shell computed inline) -> build -> fit. @@ -517,7 +627,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu d_shell_grid, d_global_grid, d_shell_n, d_global_n, d_shell_P, d_global_P, d_shell_sigma2, d_global_sigma2, p); fit<<>>(img, d_px_x, d_px_y, d_cx, d_cy, d_d, d_invd2, - d_I, d_sigma, d_bkg, d_ok, d_shell_P, d_global_P, + d_I, d_sigma, d_bkg, d_bkg_var, d_ok, d_shell_P, d_global_P, d_shell_sigma2, d_global_sigma2, d_shell_n, d_I, d_sigma, d_ok, p, n); } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h index 6e9ef467..f83d03fb 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h @@ -29,9 +29,16 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { // --- per-reflection device arrays (grown by EnsureCapacity) --- CudaDevicePtr d_px_x, d_px_y, d_d; CudaDevicePtr d_cx, d_cy; - CudaDevicePtr d_I, d_sigma, d_bkg, d_obs_x, d_obs_y; + CudaDevicePtr d_I, d_sigma, d_bkg, d_bkg_var, d_obs_x, d_obs_y; + CudaDevicePtr d_isum; // box-sum raw sum, for the radial correction + CudaDevicePtr d_ninner, d_rbin; CudaDevicePtr d_ok, d_strong, d_has_obs; + // --- radial background curvature correction (see BraggIntegrationEngine) --- + int n_rad = 0; // radial bins, 0 when the correction is off + CudaDevicePtr d_rad_sum, d_k_diff; + CudaDevicePtr d_rad_cnt; + // --- 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. diff --git a/image_analysis/scale_merge/FrenchWilson.cpp b/image_analysis/scale_merge/FrenchWilson.cpp index 362e42b3..35d22264 100644 --- a/image_analysis/scale_merge/FrenchWilson.cpp +++ b/image_analysis/scale_merge/FrenchWilson.cpp @@ -27,7 +27,16 @@ struct Posterior { Posterior integrate_posterior(double I, double sigma, double sigma_wilson, bool centric, int npts, std::vector &logw) { const double inv_2s2 = 1.0 / (2.0 * sigma * sigma); - const double j_max = std::max(I, 0.0) + 8.0 * sigma; + // The posterior is the Gaussian likelihood tilted by the exponential prior, so it peaks at + // I - sigma^2/S and decays over whichever of sigma and S is TIGHTER. Ranging to I + 8 sigma + // regardless is wrong once sigma greatly exceeds S: with npts fixed the whole prior then falls + // inside the first grid cell, the quadrature degenerates to that one point and returns + // F = sqrt(dj/2) with sigmaF -> 0 - i.e. a reflection we know nothing about comes back looking + // like the best measured one in the file. + const double prior_scale = centric ? 2.0 * sigma_wilson : sigma_wilson; + const double peak = std::max(I - sigma * sigma / prior_scale, 0.0); + const double width = peak > 0.0 ? sigma : std::min(sigma, prior_scale); + const double j_max = peak + 10.0 * width; const double dj = j_max / npts; double max_logw = -std::numeric_limits::infinity(); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 3693c421..9672cd8b 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -130,6 +130,7 @@ void print_usage() { std::cout << " --integration-high-resolution High resolution limit for prediction/integration. If omitted (or 0), integration extends as far as the detector reaches" << std::endl; std::cout << " --max-hkl 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 Monochromatic (rotation + still): high-side clip of the background ring at mean + n*sqrt(mean) (default 4; 0 = off). This is the default background estimator - it rejects neighbour cores and zingers without the symmetric trim's Poisson skew bias. Broadband data always clip, at 3 sigma; ignored by --integrator boxsum" << std::endl; + std::cout << " --background-radial[=on|off] Correct the background ring for the CURVATURE of the radial background (default off). The signal disk and the background ring are concentric, so a background linear in position cancels between them and only curvature survives - which on an ice ring reaches +26 counts on a single reflection. Costs one short dot product per reflection and no extra pixel reads" << std::endl; std::cout << " --background-trim Use the old symmetric trimmed mean for the background ring instead of the clip (0<=f<0.5; 0.10 was the former default). Switches --background-clip off. A symmetric trim is biased low on Poisson data and adds ~5 counts to every partial, so this is for back compatibility only; 0 = plain ring mean" << std::endl; std::cout << " --integrator Spot integrator boxsum|gaussian|empirical (default: gaussian profile-fit; boxsum is the classical fallback)" << std::endl; std::cout << " --simple-stills stills: treat every reflection as a full (p=1, single-pass scale/merge); disables the default physical partiality post-refinement" << std::endl; @@ -174,6 +175,7 @@ enum { OPT_FORCE_ROTATION_LATTICE, OPT_ROTATION_NO_POSTREFINE, OPT_BACKGROUND_CLIP, + OPT_BACKGROUND_RADIAL, OPT_REFINE_GEOMETRY, OPT_BANDWIDTH, OPT_INTEGRATION_RADIUS, @@ -283,6 +285,7 @@ static option long_options[] = { {"scaling-iterations", required_argument, nullptr, OPT_SCALING_ITERATIONS}, {"scaling-high-resolution", required_argument, nullptr, OPT_SCALING_HIGH_RESOLUTION}, {"background-clip", required_argument, nullptr, OPT_BACKGROUND_CLIP}, + {"background-radial", optional_argument, nullptr, OPT_BACKGROUND_RADIAL}, {"resolution-cutoff", required_argument, nullptr, OPT_RESOLUTION_CUTOFF}, {"resolution-cc-target", required_argument, nullptr, OPT_RESOLUTION_CC_TARGET}, {"resolution-shells", required_argument, nullptr, OPT_RESOLUTION_SHELLS}, @@ -561,6 +564,7 @@ static int RunRugnux(int argc, char **argv) { int64_t scaling_iter = 3; std::optional forced_rotation_lattice; std::optional background_clip_arg; // --background-clip: background-ring high-side sigma clip + std::optional background_radial_arg; // --background-radial: radial curvature correction std::optional refine_geometry; // --refine-geometry[=N]: stills global geometry-refinement pass bool refine_geometry_disabled = false; // --refine-geometry=off: opt out of the stills default-on @@ -897,6 +901,16 @@ static int RunRugnux(int argc, char **argv) { case OPT_INTEGRATION_HIGH_RES: integration_d_min_arg = parse_double_arg(optarg, "--integration-high-resolution", logger); break; + case OPT_BACKGROUND_RADIAL: + if (optarg == nullptr || strcmp(optarg, "on") == 0) + background_radial_arg = true; + else if (strcmp(optarg, "off") == 0) + background_radial_arg = false; + else { + logger.Error("Invalid --background-radial value: {} (expected on|off)", optarg); + return 1; + } + break; case OPT_BACKGROUND_CLIP: background_clip_arg = parse_double_arg(optarg, "--background-clip", logger); break; @@ -1634,6 +1648,14 @@ static int RunRugnux(int argc, char **argv) { logger.Info("Background ring: high-side clip at {:.1f} sigma", *background_clip_arg); } + if (background_radial_arg) { + BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings(); + bis.BackgroundRadialCorrection(*background_radial_arg); + experiment.ImportBraggIntegrationSettings(bis); + logger.Info("Background ring: radial curvature correction {}", + *background_radial_arg ? "on" : "off"); + } + SpotFindingSettings spot_settings; spot_settings.enable = true; spot_settings.indexing = true; diff --git a/tests/BraggIntegrationEngineGPUTest.cpp b/tests/BraggIntegrationEngineGPUTest.cpp index 78257b87..455620b2 100644 --- a/tests/BraggIntegrationEngineGPUTest.cpp +++ b/tests/BraggIntegrationEngineGPUTest.cpp @@ -82,6 +82,7 @@ Scene BuildScene(size_t width, size_t height, int spacing = 60) { // the CPU and GPU each implement separately are both covered. DiffractionExperiment MakeExperiment(IntegratorMode mode, std::optional bandwidth_fwhm, float clip_nsigma = 4.0f, + bool radial = false, const DetectorSetup &det = DetJF(2)) { DiffractionExperiment experiment(det); // DetJF(2) (small) keeps the correctness test fast experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV) @@ -93,13 +94,14 @@ DiffractionExperiment MakeExperiment(IntegratorMode mode, std::optional b settings.BackgroundClipNSigma(clip_nsigma); else settings.BackgroundTrimFraction(0.10f); + settings.BackgroundRadialCorrection(radial); experiment.ImportBraggIntegrationSettings(settings); return experiment; } void CompareCpuVsGpu(IntegratorMode mode, std::optional bandwidth_fwhm, - float clip_nsigma = 4.0f) { - const DiffractionExperiment experiment = MakeExperiment(mode, bandwidth_fwhm, clip_nsigma); + float clip_nsigma = 4.0f, bool radial = false) { + const DiffractionExperiment experiment = MakeExperiment(mode, bandwidth_fwhm, clip_nsigma, radial); const size_t width = experiment.GetXPixelsNum(); const size_t height = experiment.GetYPixelsNum(); const size_t npixel = experiment.GetPixelsNum(); @@ -154,6 +156,10 @@ TEST_CASE("BraggIntegrationEngineGPU_MatchesCPU") { SECTION("ProfileGaussian broadband") { CompareCpuVsGpu(IntegratorMode::ProfileGaussian, 0.03f); } 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); } } // Hidden ([.]) benchmark: the raison d'etre of the GPU port is < 2 ms/frame (vs ~142 ms on the CPU @@ -164,7 +170,7 @@ TEST_CASE("BraggIntegrationEngineGPU_Benchmark", "[.][bragg_bench]") { return; } const DiffractionExperiment experiment = MakeExperiment(IntegratorMode::ProfileGaussian, std::nullopt, - 4.0f, DetJF4M()); + 4.0f, false, DetJF4M()); const size_t width = experiment.GetXPixelsNum(); const size_t height = experiment.GetYPixelsNum(); const size_t npixel = experiment.GetPixelsNum();