integration: neighbour-mask the profile-fit background ring, widen r3 8->10

ProfileIntegrate2D::BoxSum now excludes every predicted reflection's r2 disk
from the r2..r3 background annulus (mirroring BraggIntegrate2D), so a neighbour
Bragg peak can no longer bias a reflection's background high and over-subtract.
With neighbours excluded the annulus can safely widen, so the default r_3 goes
8 -> 10 (more background pixels, lower-variance estimate).

Measured (rotation lyso @1.0 A, external CCref/CCxds vs XDS): 1.05 A CCref +1.3
/ CCxds +1.3, R-meas -5 pts, low-res R-meas unchanged. Serial (Jet8 @0.0002
bandwidth): 1.68 A CC1/2 +3.9 / CCref +1.1, 1.58 A +4.0 / +1.6.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
2026-06-29 22:41:14 +02:00
co-authored by Claude Opus 4.8
parent bbd888dcc3
commit eccd10a0a7
2 changed files with 28 additions and 3 deletions
+1 -1
View File
@@ -16,7 +16,7 @@ class BraggIntegrationSettings {
IntegratorMode integrator_mode = IntegratorMode::ProfileGaussian;
float r_1 = 4;
float r_2 = 6;
float r_3 = 8;
float r_3 = 10;
float d_min_limit_A = 1.0;
std::optional<float> fixed_profile_radius;
float minimum_sigma_in_regards_to_i = 0.02;
@@ -17,6 +17,28 @@ constexpr int N_SHELL = 6;
constexpr double STRONG_I_OVER_SIGMA = 5.0;
constexpr int MIN_STRONG_PER_SHELL = 30;
// Mark the r2 signal disk of every predicted reflection. A neighbour whose disk falls inside this
// reflection's r2..r3 background ring would otherwise bias the background mean high and over-subtract;
// excluding the marked pixels from the ring removes that contamination (same as BraggIntegrate2D).
std::vector<uint8_t> BuildReflectionMask(const std::vector<Reflection> &predicted, size_t npredicted,
size_t xpixel, size_t ypixel, float r2) {
std::vector<uint8_t> mask(xpixel * ypixel, 0);
const float r2_sq = r2 * r2;
for (size_t i = 0; i < npredicted; ++i) {
const auto &r = predicted[i];
const int64_t x0 = std::max<int64_t>(0, std::floor(r.predicted_x - r2 - 1.0f));
const int64_t x1 = std::min<int64_t>(xpixel - 1, std::ceil(r.predicted_x + r2 + 1.0f));
const int64_t y0 = std::max<int64_t>(0, std::floor(r.predicted_y - r2 - 1.0f));
const int64_t y1 = std::min<int64_t>(ypixel - 1, std::ceil(r.predicted_y + r2 + 1.0f));
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);
if (d2 < r2_sq) mask[y * xpixel + x] = 1;
}
}
return mask;
}
// Rough box-sum of one reflection: total intensity over the inner disk (r1) minus the mean of the
// background ring (r2..r3), like BraggIntegrate2D. Used to seed the fit, pick strong spots for
// profile learning, and provide the per-reflection background. ok=false when the inner disk is
@@ -29,7 +51,8 @@ struct Rough {
template<class T>
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, const std::vector<uint8_t> &refl_mask,
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<int64_t>(0, std::floor(r.predicted_x - r3 - 1.0));
@@ -52,6 +75,7 @@ Rough BoxSum(const T *image, size_t xpixel, size_t ypixel, int64_t special, int6
I_sum += px;
++n_inner_valid;
} else if (d2 >= r2_sq && d2 < r3_sq) {
if (refl_mask[y * xpixel + x]) continue;
if (px == special || px == special + 1 || px == saturation || px == saturation - 1) continue;
bkg_sum += static_cast<double>(px);
++n_bkg;
@@ -110,10 +134,11 @@ std::vector<Reflection> ProfileIntegrateInternal(const DiffractionExperiment &ex
const size_t xpixel = image.GetWidth(), ypixel = image.GetHeight();
// --- Pass A: box-sum every reflection (rough I, background, centre, strong flag). ---
const auto refl_mask = BuildReflectionMask(predicted, npredicted, xpixel, ypixel, r2);
std::vector<Rough> rough(npredicted);
double inv_d2_min = std::numeric_limits<double>::max(), inv_d2_max = 0.0;
for (size_t i = 0; i < npredicted; ++i) {
rough[i] = BoxSum<T>(ptr, xpixel, ypixel, special, saturation, predicted[i],
rough[i] = BoxSum<T>(ptr, xpixel, ypixel, special, saturation, predicted[i], refl_mask,
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);