diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index 3280ae622..613f21625 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -104,16 +104,17 @@ constexpr double PROFILE_SUMMATION_MAX_NSIGMA = 10.0; constexpr double MINPK_MAX_MISSING_PEAK = 0.9; } // namespace bragg_engine -// How often the integrator silently dropped a reflection, or silently declined its own fit, over +// How often the integrator lost a reflection's clean background ring, 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; // dropped whole: the r2..r3 ring kept 5 or fewer clean pixels + 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. This is the count that answers + // 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. // Measured over the battery, that floor reaches 2.3% of all reflections while the widening that diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 1adc41900..d4fc3dc93 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -155,6 +155,7 @@ 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, @@ -211,7 +212,12 @@ std::vector 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); @@ -237,7 +243,15 @@ 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; continue; } + 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; bkg_sum += static_cast(px); if (bkg_trim_frac > 0.0) bkg_vals.push_back(px); @@ -254,9 +268,9 @@ std::vector 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) { @@ -265,6 +279,18 @@ std::vector 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 @@ -291,7 +317,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)) continue; + if (refl_mask.Get(x, y) && !ring_unmasked) 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 84cd2305d..dea04ac83 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -166,6 +166,9 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, __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__ unsigned long long s_bkgsum; __shared__ int s_accept, s_full; __shared__ double s_bkg, s_thr; @@ -190,7 +193,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_ndisk = 0; s_nown = 0; s_nbkg_nb = 0; s_nbkg_nbv = 0; s_bkgsum_nb = 0; } __syncthreads(); @@ -207,8 +210,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; - long long l_bkg = 0; + 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; 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); @@ -227,8 +230,12 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, if (valid(px)) { l_Isum += px; l_Ix += (long long) x * px; l_Iy += (long long) y * px; l_x += x; l_y += y; ++l_niv; } } else if (d.inner >= p.r2_sq && d.outer < p.r3_sq) { - if (mask[y * p.W + x]) { ++l_nbnb; continue; } 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_bkg += px; ++l_nb; } @@ -245,6 +252,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, 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); __syncthreads(); for (int t = threadIdx.x; t < p.rad_w; t += blockDim.x) { s_radv[t] = 0; s_radn[t] = 0; } @@ -254,14 +263,21 @@ __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_accept = ((s_full || p.partial_ok) && s_nbkg > 5) ? 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; 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 @@ -285,7 +301,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]) continue; + if (mask[y * p.W + x] && !s_unmasked) continue; const int32_t px = img[y * p.W + x]; if (!valid(px)) continue; const int slot = atomicAdd(&s_bn, 1); @@ -327,7 +343,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]) continue; + if (mask[y * p.W + x] && !s_unmasked) continue; const int32_t px = img[y * p.W + x]; if (!valid(px)) continue; if ((long long) px <= s_thr_i) { diff --git a/tests/BraggIntegrationEngineGPUTest.cpp b/tests/BraggIntegrationEngineGPUTest.cpp index d2df89cb4..cc74813cb 100644 --- a/tests/BraggIntegrationEngineGPUTest.cpp +++ b/tests/BraggIntegrationEngineGPUTest.cpp @@ -140,7 +140,8 @@ DiffractionExperiment MakeExperiment(IntegratorMode mode, std::optional b return experiment; } -void CompareCpuVsGpu(IntegratorMode mode, std::optional bandwidth_fwhm, +// Returns the fraction of the predicted reflections the engines kept. +double CompareCpuVsGpu(IntegratorMode mode, std::optional bandwidth_fwhm, float clip_nsigma = 4.0f, bool radial = false, int spacing = 60, float stencil_k = 0.0f, float r1 = 0.0f, float r2 = 0.0f, float r3 = 0.0f, @@ -200,6 +201,7 @@ void CompareCpuVsGpu(IntegratorMode mode, std::optional bandwidth_fwhm, CHECK(out_gpu[i].I == Catch::Approx(out_cpu[i].I).epsilon(0.03).margin(2.0)); CHECK(out_gpu[i].sigma == Catch::Approx(out_cpu[i].sigma).epsilon(0.03).margin(0.5)); } + return static_cast(out_cpu.size()) / static_cast(scene.predicted.size()); } } // namespace @@ -269,6 +271,12 @@ TEST_CASE("BraggIntegrationEngineGPU_MatchesCPU") { CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 60, 0.0f, 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. + SECTION("ProfileGaussian neighbour-starved rings") { + CHECK(CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 8) > 0.95); + } 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 // fits it over what is left, and the peak-loss rule throws back the ones that lost the profile's