diff --git a/image_analysis/bragg_integration/ProfileIntegrate2D.cpp b/image_analysis/bragg_integration/ProfileIntegrate2D.cpp index 39433081..29accbee 100644 --- a/image_analysis/bragg_integration/ProfileIntegrate2D.cpp +++ b/image_analysis/bragg_integration/ProfileIntegrate2D.cpp @@ -29,7 +29,8 @@ struct Rough { template Rough BoxSum(const T *image, size_t xpixel, size_t ypixel, int64_t special, int64_t saturation, - const Reflection &r, float r1_sq, float r2_sq, float r3_sq, float r3, float min_sigma_ratio) { + const Reflection &r, float r1_sq, float r2_sq, float r3_sq, float r3, float min_sigma_ratio, + bool apply_bkg_clip) { Rough out; const int64_t x0 = std::max(0, std::floor(r.predicted_x - r3 - 1.0)); const int64_t x1 = std::min(xpixel - 1, std::ceil(r.predicted_x + r3 + 1.0)); @@ -64,7 +65,9 @@ Rough BoxSum(const T *image, size_t xpixel, size_t ypixel, int64_t special, int6 // (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). - { + // Stills-only: on rotation (no bandwidth streak) it clips legitimate high background pixels and biases + // the mean low, slightly hurting weak (anomalous) intensities, so it is gated off there. + if (apply_bkg_clip) { const double thr = out.bkg + 3.0 * std::sqrt(std::max(out.bkg, 1.0)); double s = 0.0; int n = 0; @@ -95,6 +98,9 @@ std::vector ProfileIntegrateInternal(const DiffractionExperiment &ex const float r2 = settings.GetR2(), r2_sq = r2 * r2; const float r3 = settings.GetR3(), r3_sq = r3 * r3; const float min_sigma_ratio = settings.GetMinimumSigmaInRegardsToI(); + // The background-ring sigma-clip de-biases bandwidth-streaked stills but can clip legitimate + // background on rotation (no streaks) and bias it low; gate it to stills (a bandwidth is set). + const bool apply_bkg_clip = experiment.GetBandwidthFWHM().value_or(0.0f) > 0.0f; const int R = static_cast(std::ceil(r2)); // profile grid half-size const int G = 2 * R + 1, GG = G * G; auto idx = [G, R](int dx, int dy) { return (dy + R) * G + (dx + R); }; @@ -108,7 +114,7 @@ std::vector ProfileIntegrateInternal(const DiffractionExperiment &ex double inv_d2_min = std::numeric_limits::max(), inv_d2_max = 0.0; for (size_t i = 0; i < npredicted; ++i) { rough[i] = BoxSum(ptr, xpixel, ypixel, special, saturation, predicted[i], - r1_sq, r2_sq, r3_sq, r3, min_sigma_ratio); + r1_sq, r2_sq, r3_sq, r3, min_sigma_ratio, apply_bkg_clip); if (rough[i].ok && predicted[i].d > 0.0f) { const double inv_d2 = 1.0 / (predicted[i].d * predicted[i].d); inv_d2_min = std::min(inv_d2_min, inv_d2);