diff --git a/image_analysis/bragg_integration/ProfileIntegrate2D.cpp b/image_analysis/bragg_integration/ProfileIntegrate2D.cpp index c43a5d2d..95cfeed1 100644 --- a/image_analysis/bragg_integration/ProfileIntegrate2D.cpp +++ b/image_analysis/bragg_integration/ProfileIntegrate2D.cpp @@ -39,6 +39,8 @@ Rough BoxSum(const T *image, size_t xpixel, size_t ypixel, int64_t special, int6 int64_t I_sum = 0, n_inner = 0, n_inner_valid = 0; double bkg_sum = 0.0; int n_bkg = 0; + std::vector bkg_vals; + bkg_vals.reserve(128); for (int64_t y = y0; y <= y1; ++y) for (int64_t x = x0; x <= x1; ++x) { const double d2 = (x - r.predicted_x) * (x - r.predicted_x) + (y - r.predicted_y) * (y - r.predicted_y); @@ -52,11 +54,23 @@ Rough BoxSum(const T *image, size_t xpixel, size_t ypixel, int64_t special, int6 if (px == special || px == special + 1 || px == saturation || px == saturation - 1) continue; bkg_sum += static_cast(px); ++n_bkg; + bkg_vals.push_back(static_cast(px)); } } if (n_inner_valid != n_inner || n_bkg <= 5) return out; out.bkg = bkg_sum / n_bkg; + // One high-outlier sigma-clip pass on the background ring. A bandwidth-streaked high-resolution spot + // (or a neighbour) leaks into the annulus and biases the mean high, which over-subtracts and drives + // the merged high-resolution intensities negative; rejecting pixels above mean + 3*sqrt(mean) removes + // that contamination (and barely touches a clean Poisson background, where ~0.1% exceed the cut). + { + const double thr = out.bkg + 3.0 * std::sqrt(std::max(out.bkg, 1.0)); + double s = 0.0; + int n = 0; + for (const double v : bkg_vals) if (v <= thr) { s += v; ++n; } + if (n > 5) out.bkg = s / n; + } out.I = static_cast(I_sum) - static_cast(n_inner) * out.bkg; out.sigma = std::max(1.0, out.I * min_sigma_ratio); if (I_sum > 0)