Bragg integration: keep the background ring clear only of neighbours with flux on the frame

The r2 regions set aside from a reflection's r2..r3 background ring were
those of every prediction in the +-4 sigma rocking window. On a finely
sliced dense pattern most of them are the tails of reflections recorded on
the frames either side, which fill every ring while the frame shows nothing
there. 9ac2ca677 kept the reflections those rings starved by taking the ring
whole, neighbour pixels included, and relying on the high-side clip.

A prediction now masks the ring only where it puts at least 5% of its flux
on the frame (partiality >= 0.05). A ring still starved by those neighbours
has real flux in it, and the reflection is dropped, as before 9ac2ca677:
taking it whole let the neighbours' wings into the background. The CPU mask
holds two levels (tail, flux); the GPU writes them in two launches so the
flux mark wins, with no atomics.

The count the widened-radius guard reads keeps measuring the pattern's
density over every prediction, tails included, so the guard decides on
the quantity its 1.13% bound was read off. Read on the gated mask it let a
cubic set keep r1=6, and R_meas went 36.9 -> 43.2%. Masked but
unreadable ring pixels are no longer counted as the neighbours' doing.

Fixed radius, same code base, --model, the 0.05 deg / 7200-frame set:
                      pre-9ac2ca677  9ac2ca677  this
  rings starved            88.9%       88.9%    21.5%
  partials ingested         7.4M        67M      53M
  completeness             17.9%       79.7%    85.1%
  ISa                      12.8        14.3     14.7
  R_meas                    6.8%        8.4%     7.6%
  R_free                   0.139       0.176    0.171
  radial misfit            0.11        0.185    0.066
  R vs model, common hkl to 0.69 A (61k):
                           0.133       0.130    0.129
  (the R_free rise over pre-9ac2ca677 is composition: 63k -> 319k
  reflections to 0.69 A, the added ones weaker.) Gate at 0.2 instead:
  0.1% starved but radial misfit 0.35, R_free 0.180.

A crystal whose header-geometry pass indexes a 7x supercell:
9ac2ca677's whole rings raised that pass's I/sigma >= 2 count past the
refined pass's by more than 10%. RefinedPassIsWorse then sent the run
back to the supercell. Now: the true cell and space group as before
9ac2ca677, same ISa, 104M partials in the header pass against 219M.

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_fandClaude Opus 5.5 committed 2026-09-28 17:13:16 +02:00
1 parent de352fdaf3
commit 5803dee2ce
10 files changed
+194 -103

No files matched your search

+15 -5
View File
@@ -146,7 +146,7 @@ double CompareCpuVsGpu(IntegratorMode mode, std::optional<float> bandwidth_fwhm,
float stencil_k = 0.0f,
float r1 = 0.0f, float r2 = 0.0f, float r3 = 0.0f,
OverlapMode overlap = OverlapMode::Off, float companion_dx = 0.0f,
bool clip_spots = false) {
bool clip_spots = false, float odd_partiality = 1.0f) {
const DiffractionExperiment experiment =
MakeExperiment(mode, bandwidth_fwhm, clip_nsigma, radial, DetJF(2), stencil_k, r1, r2, r3,
overlap);
@@ -155,7 +155,9 @@ double CompareCpuVsGpu(IntegratorMode mode, std::optional<float> bandwidth_fwhm,
const size_t npixel = experiment.GetPixelsNum();
REQUIRE(npixel == width * height);
const Scene scene = BuildScene(width, height, spacing, companion_dx, clip_spots);
Scene scene = BuildScene(width, height, spacing, companion_dx, clip_spots);
for (size_t i = 1; i < scene.predicted.size(); i += 2)
scene.predicted[i].partiality = odd_partiality;
REQUIRE(scene.image.size() == npixel);
REQUIRE(scene.predicted.size() > 60);
@@ -181,6 +183,8 @@ double CompareCpuVsGpu(IntegratorMode mode, std::optional<float> bandwidth_fwhm,
// atomic summation of the learned profile, so compare up to a small tolerance.
REQUIRE(out_gpu.size() == out_cpu.size());
REQUIRE(out_cpu.size() > 40);
CHECK(gpu.Counts().bkg_starved == cpu.Counts().bkg_starved);
CHECK(gpu.Counts().bkg_starved_by_neighbour == cpu.Counts().bkg_starved_by_neighbour);
if (clip_spots) {
// Guard against the coverage going vacuous: the punched pixels have to actually cost some
// reflections, or the two engines are being compared on a case neither of them meets.
@@ -272,10 +276,16 @@ TEST_CASE("BraggIntegrationEngineGPU_MatchesCPU") {
0.0f, 0.0f, 0.0f, OverlapMode::Exclude);
}
// Spots 8 px apart: the neighbours' r2 regions cover every background ring, so only the edge of
// the grid keeps a clean ring pixel. Both engines have to fall back to the whole ring and keep the
// reflections rather than drop them for want of a background.
// the grid keeps a clean ring pixel, and both engines have to drop the same reflections for want
// of a background - and count them as the neighbours' doing.
SECTION("ProfileGaussian neighbour-starved rings") {
CHECK(CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 8) > 0.95);
CHECK(CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 8) < 0.5);
}
// The same crowding with every other prediction in the tail of its rocking curve: those do not
// mask the rings, and both engines have to agree on which ones do.
SECTION("ProfileGaussian tail neighbours") {
CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 8, 0.0f,
0.0f, 0.0f, 0.0f, OverlapMode::Off, 0.0f, false, 0.01f);
}
SECTION("ProfileGaussian mono trim") { CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 0.0f); }
// Unreadable pixels inside the signal disks themselves: the MINPK rescue keeps the reflection and