mx/Jungfraujoch
3
2
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 81 Wiki Activity

Bragg integration: a ring the neighbours starved is taken whole, not dropped

Browse Source 
The r2..r3 background ring excludes every pixel inside any predicted
reflection's r2 region, and a reflection whose ring kept five or fewer
clean pixels was dropped. The mask marks every prediction in the +-4 sigma
rocking window, most of which put a few percent of their flux on the
frame. On a finely sliced, dense pattern (0.05 deg frames, sigma_M 0.067
deg, 0.35 A, 7200 frames) they fill every ring: ~9400 predictions per frame
but only ~900 with partiality above 0.2, and 88.9% of all partials were
dropped at the fixed radius - including the peak frames. The merge was 17.9%
complete.

Such a ring now falls back on all its readable pixels; the existing
high-side clip removes the neighbour cores that are really there. The
starvation counts are unchanged, so the adaptive-radius guard reads the
same numbers. Detector-starved rings are still dropped.

On that set (fixed radius, no --model): 7.4M -> 67M partials ingested,
completeness 17.9 -> 79.7% (88% to 0.65 A), d_min 0.69 -> 0.61 A,
ISa 12.8 -> 14.3, reflections in common with the deposition 60k -> 252k,
rank CC vs |Fc|^2 minus the depositor's, outer shells -0.10 -> +0.07.
A dense 360-frame set: 14.9M -> 17.6M partials, dCC outer -0.116 -> -0.096.
myob_x10sa MTZ byte-identical.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C
This commit is contained in:
leonarski_f
2026-09-28 00:11:01 +02:00
co-authored by Claude Opus 5.5
parent f4f9157cc0
commit 9ac2ca677f
4 changed files with 67 additions and 16 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp
+31 -5
View File
@@ -148,6 +148,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
std::vector<Rough> rough(npredicted);
double inv_d2_min = std::numeric_limits<double>::max(), inv_d2_max = 0.0;
std::vector<int32_t> bkg_vals; // reused per reflection for the trimmed-mean background (idea 1)
std::vector<int32_t> bkg_vals_neighbour; // the ring pixels a neighbour's region holds, likewise
// 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,
@@ -187,7 +188,12 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
// someone else's - never this reflection's own core. That makes this an exact count of what
// the pattern's density took, separable from what the detector took.
int n_bkg_neighbour = 0;
// The readable pixels among those, kept apart so a ring the neighbours starved can fall back
// on them (below).
double bkg_sum_neighbour = 0.0;
int n_bkg_neighbour_valid = 0;
bkg_vals.clear();
bkg_vals_neighbour.clear();
for (int y = y0; y <= y1; ++y)
for (int x = x0; x <= x1; ++x) {
const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y);
@@ -213,7 +219,15 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
y_sum += y;
++n_inner_valid;
} else if (d.inner >= r2_sq && d.outer < r3_sq) {
if (refl_mask[y * W + x]) { ++n_bkg_neighbour; continue; }
if (refl_mask[y * W + x]) {
++n_bkg_neighbour;
if (valid(px)) {
bkg_sum_neighbour += static_cast<double>(px);
++n_bkg_neighbour_valid;
if (bkg_trim_frac > 0.0) bkg_vals_neighbour.push_back(px);
}
continue;
}
if (!valid(px)) continue;
bkg_sum += static_cast<double>(px);
if (bkg_trim_frac > 0.0) bkg_vals.push_back(px);
@@ -230,9 +244,9 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
// survived to constrain the amplitude (XDS's MINPK, dials' valid_foreground_threshold). A box
// sum has no profile to renormalise with, so there it stays all or nothing.
const bool full = n_inner_valid == n_inner;
// A ring left with five or fewer clean pixels cannot estimate a background, so the reflection
// is dropped whole - the one thing the stencil geometry does to the DATA rather than to a
// measurement. Counted here because it is the only direct evidence of a radius that has
// A ring left with five or fewer readable pixels cannot estimate a background, so the
// reflection is dropped whole - the one thing the stencil geometry does to the DATA rather
// than to a measurement. Counted here because it is the only direct evidence of a radius that has
// outgrown the pattern it is integrating (BraggIntegrationCounts).
const bool keep_partial = full || mode != IntegratorMode::BoxSum;
if (keep_partial && n_bkg <= 5) {
@@ -241,6 +255,18 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
// detector, that took it.
if (n_bkg + n_bkg_neighbour > 5) ++counts.bkg_starved_by_neighbour;
}
// A ring the NEIGHBOURS starved still holds readable pixels - only the mask set them aside.
// Most of those neighbours are predictions in the tail of their rocking curve, which put a
// few percent of their flux on this frame: on a finely sliced, dense pattern they fill the
// whole ring while the frame shows no spot there. Dropping the reflection loses a measurement
// to protect a background from flux that is mostly not there, so the ring is taken whole
// instead, and the high-side clip below removes the neighbour cores that are.
const bool ring_unmasked = keep_partial && n_bkg <= 5 && n_bkg + n_bkg_neighbour_valid > 5;
if (ring_unmasked) {
bkg_sum += bkg_sum_neighbour;
n_bkg += n_bkg_neighbour_valid;
bkg_vals.insert(bkg_vals.end(), bkg_vals_neighbour.begin(), bkg_vals_neighbour.end());
}
if (keep_partial && n_bkg > 5) {
out.bkg = bkg_sum / n_bkg;
if (bkg_trim_frac > 0.0 && bkg_vals.size() > 5
@@ -267,7 +293,7 @@ std::vector<Reflection> BraggIntegrationEngineCPU::RunImpl(const Sampler &img,
for (int x = x0; x <= x1; ++x) {
const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y);
if (!(d.inner >= r2_sq && d.outer < r3_sq)) continue;
if (refl_mask[y * W + x]) continue;
if (refl_mask[y * W + x] && !ring_unmasked) continue;
const int32_t px = img[y * W + x];
if (!valid(px)) continue;
if (static_cast<double>(px) <= thr) {