diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 9fc33481c..3fdfd34bc 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -15,6 +15,7 @@ * Rugnux's `--model` rigid-body refinement computes the model's density and solvent mask on all threads, with results unchanged. * Rugnux's `--model` rigid-body refinement computes the model's structure factors from one copy of the model and derives the translation part of its Jacobian exactly, so it runs faster; placements can differ slightly. * `rugnux --model` also draws a null replicate again when its rigid-body placement ends near an orientation equivalent to the model's, not only when it starts near one. +* Rugnux keeps a reflection's background ring clear only of the neighbouring reflections that put at least 5% of their flux on the frame, so finely sliced dense rotation data no longer lose most of their observations for want of a background. ### 1.0.0-rc.173 diff --git a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md index e11387285..cf9aed468 100644 --- a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md +++ b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md @@ -87,7 +87,7 @@ Pixels are classified by their squared distance $r^2=(x-x_p)^2+(y-y_p)^2$: - **signal region:** $r^2 < r_1^2$, - **background annulus:** $r_2^2 \le r^2 < r_3^2$. -Invalid pixels (masked/bad/saturated) are excluded from both sums. In addition, pixels lying inside the signal disk ($r`, default 0).** The three radii above are one triple for the whole run, identical for every reflection at every resolution. A reflection is not round, though: a finite bandwidth streaks it radially by $\sigma_\mathrm{bw}=\text{bandwidth}\cdot R_\mathrm{px}$, with $R_\mathrm{px}$ the distance from the beam centre — the same physical smear as §8.2's and §11.1's $\sigma_\mathrm{bw}$, expressed here in detector pixels where those sections use reciprocal units; the two forms are never mixed in one formula. Throughout, $\text{bandwidth}$ is the **rms** relative energy spread: the user-facing `--bandwidth` takes a FWHM (a DMM's usual specification) and it is divided by 2.355 on input. On a radially smeared spot the fixed $6\ldots13$ px ring therefore sits only $\approx1.3$–$2.2$ radial $\sigma$ from the centre — on the reflection's own tails, which it then measures as background. @@ -189,9 +189,9 @@ The factor 2 is not fitted: for a Gaussian $r_{80} = 1.794\,\sigma$, so $r_1 = 2 **It applies to the final pass only.** A rotation run integrates twice (§7.5), and the widened radius is handed to the canonical second pass, not to the geometry pre-pass. The reason is that post-refinement takes its observed positions from the integrator, and an observed position is a first moment over the signal disk with the background still in it: a flat background adds nothing to the numerator but adds $n\,b$ to the denominator, so every measured offset is pulled toward its prediction by $I/(I + n b)$, and $n$ nearly doubles between $r_1 = 4$ and $r_1 = 6$. A wider disk therefore *under*-corrects the geometry — enough, on a crystal whose metric is half a degree off orthorhombic, to flip the second pass's de-novo Bravais choice. -**The density guard.** Widening $r_1$ pushes $r_2$, the ring's inner edge, into the neighbours; a reflection whose ring is left with five or fewer clean pixels has no background and is dropped whole. The integrator counts these, and separates the ones lost to **neighbouring reflections** from the ones lost to the detector itself (module gaps, the beam stop, the resolution mask) — a floor that reaches a couple of percent on some geometries and does not move with $r_1$. Where the neighbour-driven loss exceeds 1.13 % of the predicted reflections, the pattern is too dense for the widened radius and the final pass is integrated again at the fixed $4,6,13$, reported as pass 3 of 3 with the reason in `PASS_DECISION`. +**The density guard.** Widening $r_1$ pushes $r_2$, the ring's inner edge, into the neighbours; a reflection whose ring is left with five or fewer clean pixels has no background and is dropped whole. The integrator counts the rings the **neighbouring reflections** would starve — every predicted neighbour, the rocking-curve tails the background itself does not exclude included, so the count measures the pattern's density — apart from the ones the detector itself starves (module gaps, the beam stop, the resolution mask), a floor that reaches a couple of percent on some geometries and does not move with $r_1$. Where the neighbour-driven count exceeds 1.13 % of the predicted reflections, the pattern is too dense for the widened radius and the final pass is integrated again at the fixed $4,6,13$, reported as pass 3 of 3 with the reason in `PASS_DECISION`. -Every integration pass, adaptive or not, now logs the radii it used together with the fraction of predicted reflections that lost their background ring, split into the neighbour and detector parts, and the profile-fit fallback rate. +Every integration pass, adaptive or not, now logs the radii it used together with the fraction of predicted reflections that lost their background ring, the fraction of rings the neighbouring predictions crowd, and the profile-fit fallback rate. --- diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index 613f21625..83e89b583 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -102,21 +102,34 @@ constexpr double PROFILE_SUMMATION_MAX_NSIGMA = 10.0; // needs nothing the fit does not already have, so it costs one max-reduction in the loop that // measures the readable fraction, and it applies to the learned empirical profile unchanged. constexpr double MINPK_MAX_MISSING_PEAK = 0.9; + +// A predicted reflection sets its r2 region aside from its neighbours' background rings only where it +// puts at least this fraction of its flux on the frame. Prediction reaches +-4 sigma of the rocking +// curve, so on a finely sliced frame most predictions are the tails of reflections recorded on the +// frames either side: on a dense pattern those fill every ring while the frame shows nothing there, +// and a reflection whose ring they take is dropped. A still, whose partiality is 1, masks as before. +constexpr float NEIGHBOUR_MASK_MIN_PARTIALITY = 0.05f; +// The two marks a neighbour leaves in the reflection mask: its signal region holds a tail of a rocking +// curve only, or flux the ring has to be kept clear of. The larger wins where two regions meet. +constexpr uint8_t MASK_TAIL = 1; +constexpr uint8_t MASK_FLUX = 2; } // namespace bragg_engine -// How often the integrator lost a reflection's clean background ring, or silently declined its own fit, over +// How often the integrator silently dropped a reflection, or silently declined its own fit, over // every image an engine has run. Both engines keep the same two counts, so a caller sees the same // numbers whichever one it got. Nothing inside the engine reads them: they exist so a caller that // CHOSE the stencil can find out what that choice cost, which no quantity available before // integration measures. rugnux reads them out of its first pass (see Rugnux::RunAllPasses). struct BraggIntegrationCounts { uint64_t predicted = 0; // reflections offered to the engine - uint64_t bkg_starved = 0; // the r2..r3 ring kept 5 or fewer clean pixels - // Of those, the ones the NEIGHBOURS starved: their ring would have kept more than five pixels but - // for the ones a neighbouring reflection's signal region occupies. These are integrated against - // the whole ring rather than dropped; the rest of bkg_starved is dropped whole. This is the count that answers - // "is the aperture too wide for this pattern", because the rest of bkg_starved is module gaps, the - // beam stop and the resolution mask - a property of the detector that a wider r1 does not change. + uint64_t bkg_starved = 0; // dropped whole: the r2..r3 ring kept 5 or fewer clean pixels + // The rings the NEIGHBOURS would starve: more than five readable pixels, but five or fewer that no + // predicted neighbour's signal region covers - counting every prediction, the tails of rocking + // curves the background does not exclude (NEIGHBOUR_MASK_MIN_PARTIALITY) included, so it measures + // the pattern's density whatever the ring then keeps. This is the count that answers + // "is the aperture too wide for this pattern", because a ring the detector alone starves is module + // gaps, the beam stop and the resolution mask - a property of the detector that a wider r1 does not + // change. // Measured over the battery, that floor reaches 2.3% of all reflections while the widening that // costs data adds 4.1 percentage points on top of a floor of 0.02%. uint64_t bkg_starved_by_neighbour = 0; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index d4fc3dc93..9c5d509d9 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -88,14 +88,18 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, auto grid_idx = [this](int dx, int dy) { return (dy + R) * G + (dx + R); }; - // --- Reflection mask: mark the r2 signal region of every predicted reflection so a neighbour's - // signal is excluded from this reflection's r2..r3 background ring. The region is the INNER - // stencil ellipse, taken in each neighbour's OWN frame - an elongated neighbour whose streak - // is still marked as a disk would leak its tails into this reflection's ring. --- + // --- Reflection mask: mark the r2 signal region of every predicted reflection, so a neighbour's + // signal is excluded from this reflection's r2..r3 background ring. MASK_FLUX where the + // neighbour puts flux on this frame (NEIGHBOUR_MASK_MIN_PARTIALITY) - only those are excluded + // - and MASK_TAIL where only the tail of a rocking curve does, which is counted for the density + // guard but read as background. The region is the INNER stencil ellipse, taken in each + // neighbour's OWN frame - an elongated neighbour whose streak is still marked as a disk would + // leak its tails into this reflection's ring. --- refl_mask.Clear(); owner.Clear(); for (size_t i = 0; i < npredicted; ++i) { const auto &r = predicted[i]; + const uint8_t level = r.partiality >= NEIGHBOUR_MASK_MIN_PARTIALITY ? MASK_FLUX : MASK_TAIL; const BraggStencil st = MakeBraggStencil(r.predicted_x, r.predicted_y, stencil); const int x0 = std::max(0, static_cast(std::floor(r.predicted_x - st.ex_in - 1.0f))); const int x1 = std::min(W - 1, static_cast(std::ceil(r.predicted_x + st.ex_in + 1.0f))); @@ -104,7 +108,10 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, 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); - if (d.inner < r2_sq) refl_mask.At(x, y) = 1; + if (d.inner < r2_sq) { + uint8_t &m = refl_mask.At(x, y); + m = std::max(m, level); + } } } @@ -155,7 +162,6 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, std::vector rough(npredicted); double inv_d2_min = std::numeric_limits::max(), inv_d2_max = 0.0; std::vector bkg_vals; // reused per reflection for the trimmed-mean background (idea 1) - std::vector 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, @@ -207,17 +213,13 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, int n_disk = 0, n_own = 0; // pixels in the signal disk, and how many are this reflection's double bkg_sum = 0.0; int n_bkg = 0; - // Ring pixels a NEIGHBOUR's signal region occupies. refl_mask marks d.inner < r2_sq and the - // ring is d.inner >= r2_sq, so the two are complementary and a masked ring pixel is always - // 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; + // Readable ring pixels, and those of them no predicted neighbour's signal region covers. + // refl_mask marks d.inner < r2_sq and the ring is d.inner >= r2_sq, so the two are + // complementary and a masked ring pixel is always someone else's - never this reflection's + // own core. That makes the difference an exact count of what the pattern's density took, + // separable from what the detector took. + int n_ring = 0, n_clean = 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); @@ -243,16 +245,11 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, y_sum += y; ++n_inner_valid; } else if (d.inner >= r2_sq && d.outer < r3_sq) { - if (refl_mask.Get(x, y)) { - ++n_bkg_neighbour; - if (valid(px)) { - bkg_sum_neighbour += static_cast(px); - ++n_bkg_neighbour_valid; - if (bkg_trim_frac > 0.0) bkg_vals_neighbour.push_back(px); - } - continue; - } if (!valid(px)) continue; + ++n_ring; + const uint8_t m = refl_mask.Get(x, y); + if (m == 0) ++n_clean; + if (m == MASK_FLUX) continue; bkg_sum += static_cast(px); if (bkg_trim_frac > 0.0) bkg_vals.push_back(px); ++n_bkg; @@ -273,24 +270,11 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, // 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) { - ++counts.bkg_starved; - // Would the ring have been enough without the neighbours? Then it is the pattern, not the - // 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) ++counts.bkg_starved; + // Would every predicted neighbour, the tails included, leave the ring starved that the + // detector alone would not? That is the pattern's density, what the guard on a widened + // radius reads (BraggIntegrationCounts). + if (keep_partial && n_clean <= 5 && n_ring > 5) ++counts.bkg_starved_by_neighbour; if (keep_partial && n_bkg > 5) { out.bkg = bkg_sum / n_bkg; if (bkg_trim_frac > 0.0 && bkg_vals.size() > 5 @@ -317,7 +301,7 @@ std::vector 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.Get(x, y) && !ring_unmasked) continue; + if (refl_mask.Get(x, y) == MASK_FLUX) continue; const int32_t px = img[y * W + x]; if (!valid(px)) continue; if (static_cast(px) <= thr) { diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index dea04ac83..4ff0d148b 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -50,8 +50,8 @@ constexpr int COUNT_BKG_STARVED_NEIGHBOUR = 1; constexpr int COUNT_PROFILE_FALLBACK = 2; constexpr int COUNT_SLOTS = 3; -// --- Mark the r2 signal region of every predicted reflection (race-free: all writes are 1). The -// region is the INNER stencil ellipse in the neighbour's own frame; see the CPU engine. +// --- Mark the r2 signal region of every predicted reflection (race-free: one launch writes one +// value). The region is the INNER stencil ellipse in the neighbour's own frame; see the CPU engine. // The same sweep fills the owner map when an overlap treatment is on: a pixel two signal regions // share belongs to the nearer centre, which an atomicMin over (distance, index) settles without // any ordering or sorting (BraggStencil.h). The claim disk sits inside the ellipse this loop @@ -60,10 +60,15 @@ constexpr int COUNT_SLOTS = 3; // pixels back the way they were. The two directions share one kernel because they have to visit // exactly the same set: what is written here is what has to be unwritten, and a second kernel that // recomputed the boxes for itself would be free to drift from this one. -__global__ void mark_mask(const float *px_x, const float *px_y, uint8_t *mask, uint32_t *owner, - BraggGpuParams p, int n, int clear) { +// `value` is what the launch writes: MASK_TAIL from every reflection, which also claims the owner map; +// then MASK_FLUX from the ones that put flux on this frame (`marks`), a later launch so that it wins +// over the first where two regions meet, as the CPU engine's maximum does; and 0 to clear. +__global__ void mark_mask(const float *px_x, const float *px_y, const uint8_t *marks, uint8_t *mask, + uint32_t *owner, BraggGpuParams p, int n, int value) { const int i = blockIdx.x; if (i >= n) return; + if (value == MASK_FLUX && !marks[i]) return; + const bool clear = value == 0; const float cx = px_x[i], cy = px_y[i]; const BraggStencil st = MakeBraggStencil(cx, cy, p.stencil); const int x0 = max(0, (int) floorf(cx - st.ex_in - 1.0f)); @@ -75,8 +80,8 @@ __global__ void mark_mask(const float *px_x, const float *px_y, uint8_t *mask, u for (int t = threadIdx.x; t < bw * bh; t += blockDim.x) { const int x = x0 + t % bw, y = y0 + t / bw; const BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy); - if (d.inner < p.r2_sq) mask[y * p.W + x] = clear ? (uint8_t) 0 : (uint8_t) 1; - if (p.overlap && d.signal < p.claim_sq) { + if (d.inner < p.r2_sq) mask[y * p.W + x] = (uint8_t) value; + if (p.overlap && value != MASK_FLUX && d.signal < p.claim_sq) { // Clearing writes one fixed value, so the blocks that share a pixel cannot disagree and // no atomic is needed; marking keeps the nearest centre and does. if (clear) @@ -165,10 +170,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, __shared__ unsigned long long s_Isum, s_Ix, s_Iy; __shared__ unsigned long long s_x, s_y; // positions behind s_Isum, to take the background out of the centroid __shared__ int s_ninner, s_ninner_valid, s_nbkg, s_ndisk, s_nown; - __shared__ int s_nbkg_nb; // ring pixels a NEIGHBOUR's signal region holds; see the CPU engine - __shared__ int s_nbkg_nbv; // the readable ones among them, and their sum - __shared__ unsigned long long s_bkgsum_nb; - __shared__ int s_unmasked; // the ring was starved by neighbours and is taken whole; see the CPU engine + __shared__ int s_nring, s_nclean; // readable ring pixels, and those no neighbour covers; see the CPU engine __shared__ unsigned long long s_bkgsum; __shared__ int s_accept, s_full; __shared__ double s_bkg, s_thr; @@ -193,7 +195,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, s_r0 = sqrtf(rx * rx + ry * ry); s_Isum = 0; s_Ix = 0; s_Iy = 0; s_x = 0; s_y = 0; s_ninner = 0; s_ninner_valid = 0; s_nbkg = 0; s_bkgsum = 0; - s_ndisk = 0; s_nown = 0; s_nbkg_nb = 0; s_nbkg_nbv = 0; s_bkgsum_nb = 0; + s_ndisk = 0; s_nown = 0; s_nring = 0; s_nclean = 0; } __syncthreads(); @@ -210,8 +212,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, long long l_Isum = 0, l_Ix = 0, l_Iy = 0; long long l_x = 0, l_y = 0; // positions behind l_Isum, for the background-free centroid - int l_ni = 0, l_niv = 0, l_nb = 0, l_nd = 0, l_no = 0, l_nbnb = 0, l_nbnbv = 0; - long long l_bkg = 0, l_bkg_nb = 0; + int l_ni = 0, l_niv = 0, l_nb = 0, l_nd = 0, l_no = 0, l_nring = 0, l_nclean = 0; + long long l_bkg = 0; for (int t = threadIdx.x; t < area; t += blockDim.x) { const int x = x0 + t % bw, y = y0 + t / bw; const BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy); @@ -231,12 +233,11 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, l_x += x; l_y += y; ++l_niv; } } else if (d.inner >= p.r2_sq && d.outer < p.r3_sq) { const int32_t px = img[y * p.W + x]; - if (mask[y * p.W + x]) { - ++l_nbnb; - if (valid(px)) { l_bkg_nb += px; ++l_nbnbv; } - continue; - } if (!valid(px)) continue; + ++l_nring; + const uint8_t m = mask[y * p.W + x]; + if (m == 0) ++l_nclean; + if (m == MASK_FLUX) continue; l_bkg += px; ++l_nb; } } @@ -251,9 +252,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, WARP_ATOMIC_ADD(s_bkgsum, (unsigned long long) l_bkg); WARP_ATOMIC_ADD(s_ndisk, l_nd); WARP_ATOMIC_ADD(s_nown, l_no); - WARP_ATOMIC_ADD(s_nbkg_nb, l_nbnb); - WARP_ATOMIC_ADD(s_nbkg_nbv, l_nbnbv); - WARP_ATOMIC_ADD(s_bkgsum_nb, (unsigned long long) l_bkg_nb); + WARP_ATOMIC_ADD(s_nring, l_nring); + WARP_ATOMIC_ADD(s_nclean, l_nclean); __syncthreads(); for (int t = threadIdx.x; t < p.rad_w; t += blockDim.x) { s_radv[t] = 0; s_radn[t] = 0; } @@ -263,21 +263,15 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, // renormalises to the pixels it can read and Pass B cuts on how much of the profile survived // (XDS's MINPK). A box sum has no profile to renormalise with. See the CPU engine. s_full = (s_ninner_valid == s_ninner) ? 1 : 0; - s_unmasked = 0; - // A reflection the disk would have kept but the ring cannot support: dropped whole for want of - // a background. One atomic per DROPPED reflection, so nothing is paid where none is dropped. - // See BraggIntegrationCounts; the CPU engine counts the same condition. - if ((s_full || p.partial_ok) && s_nbkg <= 5) { - atomicAdd(counts + COUNT_BKG_STARVED, 1ull); - if (s_nbkg + s_nbkg_nb > 5) atomicAdd(counts + COUNT_BKG_STARVED_NEIGHBOUR, 1ull); - // A ring the neighbours starved is taken whole instead of dropping the reflection. - if (s_nbkg + s_nbkg_nbv > 5) { - s_unmasked = 1; - s_nbkg += s_nbkg_nbv; - s_bkgsum += s_bkgsum_nb; - } - } s_accept = ((s_full || p.partial_ok) && s_nbkg > 5) ? 1 : 0; + // A reflection the disk would have kept but the ring cannot support: dropped whole for want of + // a background. And a ring every predicted neighbour would starve: the density the guard on a + // widened radius reads. One atomic per such reflection, so nothing is paid where there is none. + // See BraggIntegrationCounts; the CPU engine counts the same conditions. + if ((s_full || p.partial_ok) && s_nbkg <= 5) + atomicAdd(counts + COUNT_BKG_STARVED, 1ull); + if ((s_full || p.partial_ok) && s_nclean <= 5 && s_nring > 5) + atomicAdd(counts + COUNT_BKG_STARVED_NEIGHBOUR, 1ull); s_bkg = s_accept ? ((double) (long long) s_bkgsum / (double) s_nbkg) : 0.0; s_thr = s_bkg + (double) p.bkg_clip_nsigma * sqrt(fmax(s_bkg, 1.0)); // The pixel is an integer, so comparing it against the floor of the threshold accepts @@ -301,7 +295,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, const int x = x0 + t % bw, y = y0 + t / bw; const BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy); if (!(d.inner >= p.r2_sq && d.outer < p.r3_sq)) continue; - if (mask[y * p.W + x] && !s_unmasked) continue; + if (mask[y * p.W + x] == MASK_FLUX) continue; const int32_t px = img[y * p.W + x]; if (!valid(px)) continue; const int slot = atomicAdd(&s_bn, 1); @@ -343,7 +337,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, const int x = x0 + t % bw, y = y0 + t / bw; const BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy); if (!(d.inner >= p.r2_sq && d.outer < p.r3_sq)) continue; - if (mask[y * p.W + x] && !s_unmasked) continue; + if (mask[y * p.W + x] == MASK_FLUX) continue; const int32_t px = img[y * p.W + x]; if (!valid(px)) continue; if ((long long) px <= s_thr_i) { @@ -876,6 +870,7 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) { d_px_x = CudaDevicePtr(new_capacity); d_px_y = CudaDevicePtr(new_capacity); d_d = CudaDevicePtr(new_capacity); + d_mark = CudaDevicePtr(new_capacity); d_cx = CudaDevicePtr(new_capacity); d_cy = CudaDevicePtr(new_capacity); d_I = CudaDevicePtr(new_capacity); @@ -896,6 +891,7 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) { h_px_x = CudaHostPtr(new_capacity); h_px_y = CudaHostPtr(new_capacity); h_d = CudaHostPtr(new_capacity); + h_mark = CudaHostPtr(new_capacity); h_I = CudaHostPtr(new_capacity); h_sigma = CudaHostPtr(new_capacity); h_bkg = CudaHostPtr(new_capacity); @@ -926,10 +922,12 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu h_px_x[i] = predicted[i].predicted_x; h_px_y[i] = predicted[i].predicted_y; h_d[i] = predicted[i].d; + h_mark[i] = predicted[i].partiality >= NEIGHBOUR_MASK_MIN_PARTIALITY ? 1 : 0; } cuda_err(cudaMemcpyAsync(d_px_x, h_px_x.get(), sizeof(float) * npredicted, cudaMemcpyHostToDevice, *stream)); cuda_err(cudaMemcpyAsync(d_px_y, h_px_y.get(), sizeof(float) * npredicted, cudaMemcpyHostToDevice, *stream)); cuda_err(cudaMemcpyAsync(d_d, h_d.get(), sizeof(float) * npredicted, cudaMemcpyHostToDevice, *stream)); + cuda_err(cudaMemcpyAsync(d_mark, h_mark.get(), sizeof(uint8_t) * npredicted, cudaMemcpyHostToDevice, *stream)); BraggGpuParams p{ .W = static_cast(xpixel), .H = static_cast(ypixel), @@ -979,7 +977,9 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu reset<<<32, 256, 0, *stream>>>(d_shell_grid, d_global_grid, d_mom, d_shell_n, d_global_n, d_invd2, GG); cuda_err(cudaGetLastError()); dirty = true; - mark_mask<<>>(d_px_x, d_px_y, d_mask, d_owner, p, n, 0); + mark_mask<<>>(d_px_x, d_px_y, d_mark, d_mask, d_owner, p, n, MASK_TAIL); + cuda_err(cudaGetLastError()); + mark_mask<<>>(d_px_x, d_px_y, d_mark, d_mask, d_owner, p, n, MASK_FLUX); cuda_err(cudaGetLastError()); boxsum<<>>(d_px_x, d_px_y, d_d, img, d_mask, d_owner, p, n, d_cx, d_cy, d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_obs_x, d_obs_y, @@ -1017,7 +1017,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu // Nothing below reads the mask or the owner map, so put them back now, over the boxes that were // marked rather than over the frame - where that is the cheaper of the two. if (clear_by_box) { - mark_mask<<>>(d_px_x, d_px_y, d_mask, d_owner, p, n, 1); + mark_mask<<>>(d_px_x, d_px_y, d_mark, d_mask, d_owner, p, n, 0); cuda_err(cudaGetLastError()); } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h index a95f78c47..e6aa5b497 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h @@ -37,6 +37,7 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { // --- per-reflection device arrays (grown by EnsureCapacity) --- CudaDevicePtr d_px_x, d_px_y, d_d; + CudaDevicePtr d_mark; // the reflection marks its signal region in d_mask CudaDevicePtr d_cx, d_cy; CudaDevicePtr d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_obs_x, d_obs_y; CudaDevicePtr d_isum; // box-sum raw sum, for the radial correction @@ -73,6 +74,7 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { // driver has staged it through a bounce buffer, so eight of them in a row are eight serialised // round-trips rather than eight queued transfers. CudaHostPtr h_px_x, h_px_y, h_d; + CudaHostPtr h_mark; CudaHostPtr h_I, h_sigma, h_bkg, h_var_bkg, h_obs_x, h_obs_y; CudaHostPtr h_ok, h_has_obs; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 0fbb27295..d5be30ec0 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -6126,17 +6126,18 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.harmonic.null_fraction, result.harmonic.mean_d_A); // What the stencil cost, as the integrator measured it: the reflections dropped whole for a - // background ring left below six clean pixels, and - of those - the ones the NEIGHBOURS took - // rather than the detector. The second is what a wider r1 can change, and it is what the two-pass - // guard reads; the first is reported beside it because the difference is the detector's own floor. + // background ring left below six clean pixels, and the rings the NEIGHBOURS' predicted regions + // would starve where the detector alone would not - counting the tails of rocking curves the ring + // is not kept clear of, so it measures the pattern's density. The second is what a wider r1 can + // change, and it is what the two-pass guard reads. if (bragg_predicted > 0) { const auto frac = [&](uint64_t n) { return 100.0 * static_cast(n) / static_cast(bragg_predicted); }; bkg_starved_fraction_ = static_cast(bragg_starved_nb) / static_cast(bragg_predicted); logger.Info("Integration at r1={:.1f} r2={:.1f} r3={:.2f}: of {} predicted reflections, " - "{:.3f}% lost their background ring - {:.3f}% to neighbouring reflections and the " - "rest to the detector; {:.3f}% profile-fit fallbacks", + "{:.3f}% lost their background ring; the neighbouring predictions crowd {:.3f}% of the " + "rings; {:.3f}% profile-fit fallbacks", experiment_.GetBraggIntegrationSettings().GetR1(), experiment_.GetBraggIntegrationSettings().GetR2(), experiment_.GetBraggIntegrationSettings().GetR3(), diff --git a/tests/BraggIntegrationEngineCPUTest.cpp b/tests/BraggIntegrationEngineCPUTest.cpp new file mode 100644 index 000000000..f00690f42 --- /dev/null +++ b/tests/BraggIntegrationEngineCPUTest.cpp @@ -0,0 +1,78 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include +#include +#include + +#include "../common/BraggIntegrationSettings.h" +#include "../common/DetectorSetup.h" +#include "../common/DiffractionExperiment.h" +#include "../common/Reflection.h" +#include "../image_analysis/bragg_integration/BraggIntegrationEngineCPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h" + +// One spot on a flat background, inside a grid of predictions 8 px apart that put no flux on the +// frame - the tails of reflections recorded on the frames either side of a finely sliced one. Their +// r2 regions cover every background ring, so whether the spot keeps a clean ring depends only on +// whether those predictions are allowed to mask it. +TEST_CASE("BraggIntegrationEngineCPU_NeighbourMaskFollowsPartiality", "[Integration][portable]") { + DiffractionExperiment experiment(DetJF(2)); + experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(400.0f).BeamY_pxl(400.0f); + experiment.ImportBraggIntegrationSettings(BraggIntegrationSettings()); + const size_t width = experiment.GetXPixelsNum(), npixel = experiment.GetPixelsNum(); + + const float cx = 600.3f, cy = 300.2f, amp = 800.0f, sigma = 1.3f; + ImagePreprocessorBuffer image(npixel); + for (size_t i = 0; i < npixel; ++i) + image[i] = 12; + for (int dy = -6; dy <= 6; ++dy) + for (int dx = -6; dx <= 6; ++dx) { + const int x = static_cast(std::lround(cx)) + dx, y = static_cast(std::lround(cy)) + dy; + const float ex = x - cx, ey = y - cy; + image[y * width + x] += static_cast(std::lround(amp * std::exp(-(ex * ex + ey * ey) / (2 * sigma * sigma)))); + } + + auto run = [&](float neighbour_partiality, BraggIntegrationCounts &counts) { + std::vector predicted; + for (int gy = -4; gy <= 4; ++gy) + for (int gx = -4; gx <= 4; ++gx) { + Reflection r{}; + r.h = gx; r.k = gy; r.l = 1; + r.predicted_x = cx + 8.0f * gx; + r.predicted_y = cy + 8.0f * gy; + r.d = 2.0f; + r.prescaling_corr = 1.0f; + r.partiality = (gx == 0 && gy == 0) ? 0.5f : neighbour_partiality; + predicted.push_back(r); + } + BraggIntegrationEngineCPU engine(experiment); + const auto out = engine.Run(image, predicted, predicted.size(), 0); + counts = engine.Counts(); + std::vector spot; + for (const auto &r : out) + if (r.h == 0 && r.k == 0) + spot.push_back(r); + return spot; + }; + + // Tails: no reflection is dropped, and the spot is measured against the clean flat background - + // while the density the widened-radius guard reads still sees every prediction. + BraggIntegrationCounts tails; + const auto spot = run(0.01f, tails); + CHECK(tails.bkg_starved == 0); + CHECK(tails.bkg_starved_by_neighbour > 0); + REQUIRE(spot.size() == 1); + CHECK(spot[0].bkg == Catch::Approx(12.0f)); + CHECK(spot[0].I == Catch::Approx(6.2831853 * sigma * sigma * amp).epsilon(0.05)); + + // The same predictions with their flux on this frame do take the rings, the spot's among them, and + // it is the neighbours, not the detector, that starve them - the same rings the density counted. + BraggIntegrationCounts on_frame; + CHECK(run(1.0f, on_frame).empty()); + CHECK(on_frame.bkg_starved > 0); + CHECK(on_frame.bkg_starved_by_neighbour == on_frame.bkg_starved); + CHECK(on_frame.bkg_starved_by_neighbour == tails.bkg_starved_by_neighbour); +} diff --git a/tests/BraggIntegrationEngineGPUTest.cpp b/tests/BraggIntegrationEngineGPUTest.cpp index cc74813cb..2f5c45425 100644 --- a/tests/BraggIntegrationEngineGPUTest.cpp +++ b/tests/BraggIntegrationEngineGPUTest.cpp @@ -146,7 +146,7 @@ double CompareCpuVsGpu(IntegratorMode mode, std::optional 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 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 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 diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 7ce420525..85458f40a 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -18,6 +18,7 @@ ADD_EXECUTABLE(jfjoch_portable_test EXCLUDE_FROM_ALL CalcBraggPredictionTest.cpp BraggStencilTest.cpp BraggIntegrationEngineCompressedImageTest.cpp + BraggIntegrationEngineCPUTest.cpp WriteReflectionsTest.cpp LoadReferenceMtzFreeFlagsTest.cpp LoadReferenceSfCifTest.cpp @@ -76,6 +77,7 @@ ADD_EXECUTABLE(jfjoch_test ROIIntegrationGPUTest.cpp BraggIntegrationEngineGPUTest.cpp BraggIntegrationEngineCompressedImageTest.cpp + BraggIntegrationEngineCPUTest.cpp BraggStencilTest.cpp LossyFilterTest.cpp ImageBufferTest.cpp