diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index 92b21bf8..38bf8339 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -49,12 +49,20 @@ constexpr int N_SHELL = 6; // resolution shells for per-shell constexpr double STRONG_I_OVER_SIGMA = 5.0; // strong-spot threshold that seeds the profile constexpr int MIN_STRONG_PER_SHELL = 30; // below this a shell falls back to the global profile constexpr double C_CAPTURE = 2.5; // weak-spot radial capture term (monochromatic only) -// Per-pixel variance floor for the Kabsch fit weights (v = floor + signal). The detector noise floor is -// the quantization noise from rounding the charge-spread deposited energy to an integer: a uniform -// rounding error has variance 1/12. Electronic noise is far below this for both EIGER and JUNGFRAU. A -// larger floor (the previous 1.0) silently over-regularizes — it inflates weak-reflection sigma and -// pins the scaling error model's `a` term at its floor. -constexpr double PIXEL_VARIANCE_FLOOR = 1.0 / 12.0; +// Lower bound on the background term of the Kabsch fit weights (v = max(bkg, floor) + signal). It +// guards the background ESTIMATE, not the detector: the r2..r3 ring mean of a high-angle reflection +// can come out exactly zero, and v = 0 makes the weights P^2/v diverge. A ring of n pixels cannot +// resolve a background below ~1/n (0.005..0.02 for the default r2=6/r3=10 stencil), so that is the +// scale the floor has to work at. Anything larger over-regularizes: the floor multiplies the reported +// variance by floor/bkg for every pixel below it, so the previous 1/12 inflated sigma by 1.3x at +// 0.05 ct/px and 1.7x at 0.03 - exactly where the weakest high-resolution data live. Digitisation +// noise, where a detector has it, is additive on top of the background and does not belong here. +constexpr double PIXEL_VARIANCE_FLOOR = 0.01; +// The plug-in signal term of the fit weights may lower the per-pixel variance as well as raise it, +// but not below this fraction of the background. Half-wave rectifying it (max(0, I)) instead makes +// the weights - and so the reported 1/den - respond only to upward fluctuations of a noisy intensity +// 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; // 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 diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 2bde3fca..d7ce3f0a 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -284,13 +284,15 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, } // Isotropic width (2nd moment) of a learned grid: over the r1 disk (monochromatic) or the full - // grid (broadband); = 2 sigma^2 in 2D. + // grid (broadband); = 2 sigma^2 in 2D. The cells are signed: away from the peak a learned + // cell is pure background noise centred on zero, and clamping it at zero turns that noise into a + // positive pedestal spread over the whole domain, which the r^2 weight then reads as extra width. auto measure_sigma2 = [&](const std::vector &grid) { double m2 = 0.0, m2w = 0.0; for (int dy = -R; dy <= R; ++dy) for (int dx = -R; dx <= R; ++dx) { if (!broadband && dx * dx + dy * dy >= r1_sq) continue; - const double g = std::max(0.0, grid[grid_idx(dx, dy)]); + const double g = grid[grid_idx(dx, dy)]; m2 += g * (dx * dx + dy * dy); m2w += g; } @@ -396,7 +398,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, if (x < 0 || y < 0 || x >= W || y >= H) continue; const int32_t px = img[y * W + x]; if (!valid(px)) continue; - const double v = B + std::max(0.0, I) * Pp; + const double v = std::max(B + I * Pp, WEIGHT_VARIANCE_MIN_FRACTION * B); num += Pp * (static_cast(px) - rh.bkg) / v; den += Pp * Pp / v; wsum += Pp / v; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index c6046ac2..c132c6b3 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -334,9 +334,12 @@ __global__ void build_profiles(const float *shell_grid, const float *global_grid float l_m2 = 0.0f, l_m2w = 0.0f, l_sum = 0.0f; for (int k = threadIdx.x; k < p.GG; k += blockDim.x) { const int dx = k % p.G - p.R, dy = k / p.G - p.R; - const float g = fmaxf(0.0f, grid[k]); + // The moment reads the SIGNED cell (see the CPU engine): clamping pure-background noise at + // zero leaves a positive pedestal that the r^2 weight reads as extra width. The empirical + // profile still takes the clamped grid - a profile has to be non-negative. + const float raw = grid[k], g = fmaxf(0.0f, raw); const int r2i = dx * dx + dy * dy; - if (p.broadband || (float) r2i < p.r1_sq) { l_m2 += g * (float) r2i; l_m2w += g; } + if (p.broadband || (float) r2i < p.r1_sq) { l_m2 += raw * (float) r2i; l_m2w += raw; } l_sum += g; if (p.empirical) P[k] = g; // pre-store clamped grid for in-place normalisation below } @@ -454,7 +457,7 @@ __global__ void fit(const int32_t *img, const float *px_x, const float *px_y, if (x < 0 || y < 0 || x >= p.W || y >= p.H) continue; const int32_t px = img[y * p.W + x]; if (!valid(px)) continue; - const float v = B + fmaxf(0.0f, Ihere) * Pp; + const float v = fmaxf(B + Ihere * Pp, (float) WEIGHT_VARIANCE_MIN_FRACTION * B); l_num += Pp * ((float) px - bkg) / v; l_den += Pp * Pp / v; l_wsum += Pp / v;