diff --git a/common/BraggIntegrationSettings.cpp b/common/BraggIntegrationSettings.cpp index 2953582a..95b0cd39 100644 --- a/common/BraggIntegrationSettings.cpp +++ b/common/BraggIntegrationSettings.cpp @@ -160,3 +160,24 @@ BraggIntegrationSettings &BraggIntegrationSettings::BackgroundRadialCorrection(s std::optional BraggIntegrationSettings::GetBackgroundRadialCorrection() const { return bkg_radial_correction; } + +BraggIntegrationSettings &BraggIntegrationSettings::Overlap(OverlapMode input) { + overlap_mode = input; + return *this; +} + +OverlapMode BraggIntegrationSettings::GetOverlap() const { + return overlap_mode; +} + +BraggIntegrationSettings &BraggIntegrationSettings::OverlapMinPeak(float input) { + check_finite("Overlap minimum peak fraction", input); + check_min("Overlap minimum peak fraction", input, 0.0); + check_max("Overlap minimum peak fraction", input, 1.0); + overlap_min_peak = input; + return *this; +} + +float BraggIntegrationSettings::GetOverlapMinPeak() const { + return overlap_min_peak; +} diff --git a/common/BraggIntegrationSettings.h b/common/BraggIntegrationSettings.h index b3b64104..a2580ddd 100644 --- a/common/BraggIntegrationSettings.h +++ b/common/BraggIntegrationSettings.h @@ -12,6 +12,16 @@ // spots - see docs/CPU_DATA_ANALYSIS.md (Bragg integration). enum class IntegratorMode { BoxSum, ProfileGaussian, ProfileEmpirical }; +// What the integrator does about a signal region shared with a neighbouring reflection. Off is the +// historical behaviour: a neighbour's signal is kept out of this reflection's BACKGROUND ring, but +// nothing keeps it out of the reflection's own SIGNAL region, so on a dense pattern a crowded +// reflection reads high. Reject is XDS's MINPK - drop the reflection when too little of its expected +// profile is cleanly its own. Exclude - the default - drops only the shared PIXELS from the fit; a +// profile fit is the amplitude of a normalised profile, so leaving pixels out renormalises it by +// construction and the reflection is kept unbiased rather than discarded. A box sum has no profile +// to renormalise with, so Exclude does nothing there; only Reject acts on it. +enum class OverlapMode { Off, Reject, Exclude }; + // The hkl half-width the broker bootstraps when a config carries no bragg_integration block. Matches // the max_hkl default in broker/jfjoch_api.yaml, so an omitting client and an omitting config agree. constexpr int BRAGG_ONLINE_DEFAULT_MAX_HKL = 100; @@ -100,6 +110,13 @@ class BraggIntegrationSettings { // Offline that is what is wanted. ONLINE it is not: the broker bootstraps a concrete value // (BRAGG_ONLINE_DEFAULT_MAX_HKL) so per-image cost stays predictable across samples. std::optional max_hkl; + // Overlap treatment and, for OverlapMode::Reject, the least fraction of a reflection's expected + // profile that must be cleanly its own for the reflection to be kept (XDS calls it MINPK). + // Excluding the shared pixels is the default: over the rotation battery it costs 1.1% of the wall + // clock (23% on a genuinely crowded crystal, nothing where no two predictions touch) and buys ISa + // on 15 crystals against 5, cutting the summed shortfall against XDS by a third. + OverlapMode overlap_mode = OverlapMode::Exclude; + float overlap_min_peak = 0.75f; public: BraggIntegrationSettings& R1(float input); @@ -113,6 +130,8 @@ public: BraggIntegrationSettings& BackgroundClipNSigma(float input); BraggIntegrationSettings& BackgroundRadialCorrection(std::optional input); BraggIntegrationSettings& MaxHKL(std::optional input); + BraggIntegrationSettings& Overlap(OverlapMode input); + BraggIntegrationSettings& OverlapMinPeak(float input); [[nodiscard]] IntegratorMode GetIntegrator() const; @@ -129,4 +148,6 @@ public: // Unset = auto (gate per image on the smooth-ice score); see bkg_radial_correction. [[nodiscard]] std::optional GetBackgroundRadialCorrection() const; [[nodiscard]] std::optional GetMaxHKL() const; + [[nodiscard]] OverlapMode GetOverlap() const; + [[nodiscard]] float GetOverlapMinPeak() const; }; diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 46157b99..dfcdbd49 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -5,6 +5,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as * rugnux: The error-model **a** and **b** are reported in XDS's convention, and `_reflns.jfjoch_diffrn_ISa` now carries the whole-range `1/sqrt(a*b)` that XDS's ISa denotes; the strong-reflection asymptote moves to `_reflns.jfjoch_diffrn_ISa_asymptotic`. **A file written by an earlier version carries the asymptote under the old name.** * Bragg integration: the background ring's outer radius default changes from 10 px to **13 px**, which roughly doubles the pixels behind each background estimate; the signal disk is unchanged. +* Bragg integration: signal pixels shared with a neighbouring reflection are now dropped from the profile fit (`--overlap off|reject|exclude`, **default `exclude`**), so a crowded reflection no longer reads its neighbour's flux as its own; `reject` instead discards a reflection whose cleanly-observed profile fraction is below `--overlap-minpk` (default 0.75). * rugnux: The background ring can be **elongated radially per reflection** for broadband data (`--integration-stencil `, default 0 = the fixed circular ring), by `k` times the beam's radial streak; the `r1` signal box stays circular. * rugnux: New **beam-stop shadow detection**, **on by default** (`--detect-beam-stop[=N|off]`), finds the beam stop and its holder in a projection of N images (default 60) and adds them to the pixel mask as bit 9, which is cleared at the start of every run. * Viewer: the detected beam-stop shadow is drawn in coral, with a "Show beam stop" switch in the side panel. diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp index 07dbf6d3..3136cd17 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp @@ -88,6 +88,17 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe bkg_clip_nsigma = settings.GetBackgroundClipNSigma(); bkg_trim = bkg_clip_nsigma > 0.0f ? 0.0f : settings.GetBackgroundTrimFraction(); + // Overlap treatment. Ownership is decided out to the fit grid's half size, which is where the + // profile fit reads pixels; beyond it a pixel that nobody claims is this reflection's own. + // Excluding the shared pixels needs a profile to renormalise, so it cannot act on a box sum - + // drop it to Off there rather than build an owner map nothing will read. + overlap = settings.GetOverlap(); + if (overlap == OverlapMode::Exclude && mode == IntegratorMode::BoxSum) + overlap = OverlapMode::Off; + overlap_min_peak = settings.GetOverlapMinPeak(); + claim = static_cast(R); + inv_claim = 1.0f / claim; + // Radial-offset kernels for the background curvature correction. A stencil pixel at (dx, dy) // sits at radial offset dx*cos(phi) + dy*sin(phi) from the reflection, where phi is the // reflection's azimuth; averaging over phi makes the kernels position-independent, which is diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index 057e5143..9b05f6ad 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -134,6 +134,15 @@ protected: // the CPU and GPU engines. float bkg_trim = 0.0f; + // --- overlap treatment (BraggIntegrationSettings, rugnux --overlap) --- + // Off leaves the signal region exactly as it always was. Reject and Exclude both need the owner + // map (BraggStencil.h): which of two touching reflections a shared pixel belongs to. `claim` is + // how far a reflection claims pixels - the fit grid's half size, so ownership is decided + // everywhere the fit looks. + OverlapMode overlap = OverlapMode::Off; + float overlap_min_peak = 0.0f; // Reject: least clean profile fraction that is kept + float claim = 0.0f, inv_claim = 0.0f; + // --- radial background curvature correction (BraggIntegrationSettings) --- // The disk and the annulus are concentric, so any background LINEAR in position cancels between // them; what survives is the curvature of the radial background. mean_annulus(B) - mean_disk(B) diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 75be25ca..15201c36 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -85,6 +85,37 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, } } + // --- Signal-region ownership: a pixel inside two reflections' signal regions belongs to the + // NEARER predicted centre. The union mask above cannot answer that - it also marks a + // reflection's own core - so ownership gets its own map, one (distance, reflection) key per + // pixel (BraggStencil.h). Built only when an overlap treatment is asked for. --- + // A box sum has no profile to renormalise a disk it has taken pixels out of, so it never excludes. + const bool exclude = overlap == OverlapMode::Exclude && mode != IntegratorMode::BoxSum; + std::vector owner; + if (overlap != OverlapMode::Off) { + owner.assign(npixel, BRAGG_OWNER_NONE); + const float claim_sq = claim * claim; + for (size_t i = 0; i < npredicted; ++i) { + const auto &r = predicted[i]; + const int x0 = std::max(0, static_cast(std::floor(r.predicted_x - claim))); + const int x1 = std::min(W - 1, static_cast(std::ceil(r.predicted_x + claim))); + const int y0 = std::max(0, static_cast(std::floor(r.predicted_y - claim))); + const int y1 = std::min(H - 1, static_cast(std::ceil(r.predicted_y + claim))); + for (int y = y0; y <= y1; ++y) + for (int x = x0; x <= x1; ++x) { + const float dx = x - r.predicted_x, dy = y - r.predicted_y; + const float d2 = dx * dx + dy * dy; + if (d2 >= claim_sq) continue; + uint32_t &o = owner[y * W + x]; + o = std::min(o, BraggOwnerKey(std::sqrt(d2), inv_claim, static_cast(i))); + } + } + } + auto clean = [&](int x, int y, size_t i) { + return owner.empty() + || BraggOwnedBy(owner[static_cast(y) * W + x], static_cast(i)); + }; + // --- Pass A: box-sum every reflection (rough I, background, centroid, strong flag). --- struct Rough { double I = 0.0, sigma = NAN, bkg = 0.0, obs_x = 0.0, obs_y = 0.0; @@ -92,6 +123,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, double var_bkg = 0.0; // total non-signal variance carried to the merge int64_t I_sum = 0; // kept so I can be rebuilt after the radial background correction int n_inner = 0; + int n_disk = 0, n_own = 0; // signal-disk pixels, and how many of them are this reflection's int r_bin = 0; // rounded distance from the beam centre, indexes the radial curve int k_bin = 0; // which radial-background kernel this reflection's stencil needs int cx = 0, cy = 0, shell = -1; @@ -130,6 +162,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, out.k_bin = BraggStencilKernelIndex(st, n_kern); int64_t I_sum = 0, I_sum_x = 0, I_sum_y = 0, n_inner = 0, n_inner_valid = 0; + 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; bkg_vals.clear(); @@ -138,6 +171,15 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y); const int32_t px = img[y * W + x]; if (d.signal < r1_sq) { + // A pixel a nearer neighbour owns carries that neighbour's flux, so in Exclude + // mode it leaves the disk entirely - the sum, the pixel count the background is + // subtracted with, and the all-or-nothing validity gate alike. The box sum only + // counts how many were lost, which is all it can act on. + ++n_disk; + if (overlap != OverlapMode::Off) { + if (clean(x, y, i)) ++n_own; + else if (exclude) continue; + } ++n_inner; if (!valid(px)) continue; I_sum += px; @@ -207,6 +249,8 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, + static_cast(n_inner) * n_inner * out.bkg_var; out.I_sum = I_sum; out.n_inner = static_cast(n_inner); + out.n_disk = n_disk; + out.n_own = n_own; const double var_bkg_term = static_cast(n_inner) * n_inner * out.bkg_var; out.sigma = 1.0; if (I_sum > 0) { @@ -255,6 +299,14 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, for (size_t i = 0; i < npredicted; ++i) { const auto &rh = rough[i]; if (!rh.ok) continue; + // A box sum measures what is in the disk with no model of what should be there, so it can + // neither renormalise nor tell a neighbour's photon from its own: dropping the reflection + // is the only treatment it has, and it is applied only where that is what was ASKED for. + // Exclude means "leave the shared pixels out of the fit", and a box sum has no fit, so it + // is a no-op here rather than a rejection the caller never requested. The fraction is by + // AREA, not by profile mass, so the same threshold cuts harder here than in the profile + // modes. + if (overlap == OverlapMode::Reject && rh.n_own < overlap_min_peak * rh.n_disk) continue; results[i] = {static_cast(rh.I), static_cast(rh.sigma), static_cast(rh.bkg), static_cast(rh.obs_x), static_cast(rh.obs_y), static_cast(rh.var_bkg), true, rh.has_obs}; @@ -296,6 +348,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, if (x < 0 || y < 0 || x >= W || y >= H) continue; const int32_t px = img[y * W + x]; if (!valid(px)) continue; + if (exclude && !clean(x, y, i)) continue; const double v = (static_cast(px) - rh.bkg) / rh.I; shell_grid[rh.shell][grid_idx(dx, dy)] += v; global_grid[grid_idx(dx, dy)] += v; @@ -398,6 +451,32 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, } const int Gf = 2 * Rf + 1; + + // --- How much of the expected profile is cleanly this reflection's own. p_own is the whole + // grid's clean mass, i.e. XDS's MINPK quantity, and it is what Reject cuts on. m_own / + // m_all is the same fraction over the r1 disk alone, which is what the summation seed the + // runaway guard below compares against actually saw; with nothing excluded it is 1 and + // the guard is untouched. --- + double p_own = 1.0, m_all = 0.0, m_own = 0.0; + if (overlap != OverlapMode::Off) { + p_own = 0.0; + for (int dy = -Rf; dy <= Rf; ++dy) + for (int dx = -Rf; dx <= Rf; ++dx) { + const double Pp = (*Pvec)[(dy + Rf) * Gf + (dx + Rf)]; + if (Pp <= 0.0) continue; + const int x = rh.cx + dx, y = rh.cy + dy; + if (x < 0 || y < 0 || x >= W || y >= H) continue; + if (!valid(img[y * W + x])) continue; + const bool own = clean(x, y, i); + if (own) p_own += Pp; + if (dx * dx + dy * dy < r1_sq) { + m_all += Pp; + if (own) m_own += Pp; + } + } + } + if (overlap == OverlapMode::Reject && p_own < overlap_min_peak) continue; + const double B = std::max(rh.bkg, PIXEL_VARIANCE_FLOOR); double I = rh.I, den = 0.0, wsum = 0.0; for (int iter = 0; iter < 4; ++iter) { @@ -412,6 +491,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, if (x < 0 || y < 0 || x >= W || y >= H) continue; const int32_t px = img[y * W + x]; if (!valid(px)) continue; + if (exclude && !clean(x, y, i)) continue; const double v = std::max(B + I * Pp, WEIGHT_VARIANCE_MIN_FRACTION * B); num += Pp * (static_cast(px) - rh.bkg) / v; den += Pp * Pp / v; @@ -431,7 +511,11 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, double sigma = std::sqrt(1.0 / den + (wsum / den) * (wsum / den) * rh.bkg_var); double var_bkg = std::max(0.0, 1.0 / den - std::max(0.0, I) + (wsum / den) * (wsum / den) * rh.bkg_var); - if (std::abs(I - rh.I) > PROFILE_SUMMATION_MAX_NSIGMA * rh.sigma) { + // The seed is a sum over the disk the box sum actually read, so when Exclude has taken pixels + // out of both, the fit's full-profile intensity has to be scaled down to that same disk + // before the two are comparable. Nothing excluded gives exactly 1. + const double guard_scale = exclude && m_all > 0.0 ? m_own / m_all : 1.0; + if (std::abs(I * guard_scale - rh.I) > PROFILE_SUMMATION_MAX_NSIGMA * rh.sigma) { I = rh.I; sigma = rh.sigma; var_bkg = rh.var_bkg; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index b3560c12..38a35b16 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -28,6 +28,11 @@ struct BraggGpuParams { BraggStencilParams stencil; // per-reflection signal/background geometry (BraggStencil.h) int rad_w; // radial-background window held in shared memory, in bins of one pixel int n_kern; // radial-background kernels in the table + int overlap; // 0 = off, 1 = reject, 2 = exclude (OverlapMode) + int boxsum_reject; // BoxSum mode under Reject: drop on the disk-AREA fraction (see the CPU engine) + int exclude_px; // drop a neighbour's pixels from the disk: Exclude, and not a box sum + float claim_sq, inv_claim; // how far a reflection claims pixels in the owner map + float minpk; // Reject: least clean profile fraction that is kept }; __device__ inline bool valid(int32_t v) { return v != INT32_MIN && v != INT32_MAX; } @@ -36,8 +41,13 @@ __device__ inline bool valid(int32_t v) { return v != INT32_MIN && v != INT32_MA constexpr int MOM_STRIDE = 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. --- -__global__ void mark_mask(const float *px_x, const float *px_y, uint8_t *mask, BraggGpuParams p, int n) { +// 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 +// already walks, so ownership costs one atomic on the pixels already visited. --- +__global__ void mark_mask(const float *px_x, const float *px_y, uint8_t *mask, uint32_t *owner, + BraggGpuParams p, int n) { const int i = blockIdx.x; if (i >= n) return; const float cx = px_x[i], cy = px_y[i]; @@ -52,12 +62,15 @@ __global__ void mark_mask(const float *px_x, const float *px_y, uint8_t *mask, B 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] = 1; + if (p.overlap && d.signal < p.claim_sq) + atomicMin(&owner[y * p.W + x], BraggOwnerKey(sqrtf(d.signal), p.inv_claim, i)); } } // --- Pass A box-sum: rough I / background / centroid / strong flag, one block per reflection. --- __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, - const int32_t *img, const uint8_t *mask, BraggGpuParams p, int n, + const int32_t *img, const uint8_t *mask, const uint32_t *owner, + BraggGpuParams p, int n, int *cx_o, int *cy_o, float *I_o, float *sigma_o, float *bkg_o, float *bkgvar_o, float *varbkg_o, float *obsx_o, float *obsy_o, uint8_t *ok_o, uint8_t *strong_o, uint8_t *hasobs_o, unsigned long long *invd2mm, @@ -67,7 +80,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, if (i >= n) return; __shared__ unsigned long long s_Isum, s_Ix, s_Iy; - __shared__ int s_ninner, s_ninner_valid, s_nbkg; + __shared__ int s_ninner, s_ninner_valid, s_nbkg, s_ndisk, s_nown; __shared__ double s_bkgsum; __shared__ int s_accept; __shared__ double s_bkg, s_thr, s_clipsum; @@ -85,6 +98,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_ninner = 0; s_ninner_valid = 0; s_nbkg = 0; s_bkgsum = 0.0; + s_ndisk = 0; s_nown = 0; } __syncthreads(); @@ -100,13 +114,19 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, const int area = (bw > 0 && bh > 0) ? bw * bh : 0; long long l_Isum = 0, l_Ix = 0, l_Iy = 0; - int l_ni = 0, l_niv = 0, l_nb = 0; + int l_ni = 0, l_niv = 0, l_nb = 0, l_nd = 0, l_no = 0; double l_bkg = 0.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); const int32_t px = img[y * p.W + x]; if (d.signal < p.r1_sq) { + // A pixel a nearer neighbour owns carries that neighbour's flux; see the CPU engine. + ++l_nd; + if (p.overlap) { + if (BraggOwnedBy(owner[y * p.W + x], i)) ++l_no; + else if (p.exclude_px) continue; + } ++l_ni; if (valid(px)) { l_Isum += px; l_Ix += (long long) x * px; l_Iy += (long long) y * px; ++l_niv; } } else if (d.inner >= p.r2_sq && d.outer < p.r3_sq) { @@ -122,6 +142,8 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, atomicAdd(&s_ninner_valid, l_niv); atomicAdd(&s_nbkg, l_nb); atomicAdd(&s_bkgsum, l_bkg); + atomicAdd(&s_ndisk, l_nd); + atomicAdd(&s_nown, l_no); __syncthreads(); for (int t = threadIdx.x; t < p.rad_w; t += blockDim.x) { s_radv[t] = 0.0f; s_radn[t] = 0; } @@ -221,6 +243,11 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, if (threadIdx.x != 0) return; if (!s_accept) { ok_o[i] = 0; strong_o[i] = 0; hasobs_o[i] = 0; return; } + // A box sum has no profile to renormalise a disk it has taken pixels out of, so dropping the + // reflection is the only overlap treatment it has (see the CPU engine). + if (p.boxsum_reject && (float) s_nown < p.minpk * (float) s_ndisk) { + ok_o[i] = 0; strong_o[i] = 0; hasobs_o[i] = 0; return; + } double bkg = s_bkg; int n_bkg_used = s_nbkg; // pixels behind the FINAL background value (trim/clip shrink it) @@ -300,7 +327,8 @@ __global__ void reset(float *shell_grid, float *global_grid, float *mom, int *sh // the spot's own radial/tangential frame, because a detector-frame stack is azimuthally averaged // and cannot tell a radially smeared spot from a tangentially wide one (see the CPU engine). // One block per reflection. --- -__global__ void learn_profile(const int32_t *img, const float *px_x, const float *px_y, +__global__ void learn_profile(const int32_t *img, const uint32_t *owner, + const float *px_x, const float *px_y, const int *cx_a, const int *cy_a, const float *dd, const unsigned long long *invd2mm, const float *I_a, const float *bkg_a, const uint8_t *ok_a, const uint8_t *strong_a, @@ -326,6 +354,7 @@ __global__ void learn_profile(const int32_t *img, const float *px_x, const float if (x < 0 || y < 0 || x >= p.W || y >= p.H) continue; const int32_t px = img[y * p.W + x]; if (!valid(px)) continue; + if (p.overlap == 2 && !BraggOwnedBy(owner[y * p.W + x], i)) continue; const float v = ((float) px - bkg) / I; atomicAdd(&sg[k], v); atomicAdd(&global_grid[k], v); @@ -414,7 +443,7 @@ __global__ void radial_correct(const float *rad_sum, const int *rad_cnt, int n_r // --- Pass B Kabsch profile fit: I = sum P(c-B)/v over sum P^2/v, v = B + max(I,0)P (iterate). // One block per reflection; the (possibly elongated) profile is built in shared memory. --- -__global__ void fit(const int32_t *img, const float *px_x, const float *px_y, +__global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x, const float *px_y, const int *cx_a, const int *cy_a, const float *dd, const unsigned long long *invd2mm, const float *I_seed, const float *sigma_seed, const float *bkg_a, const float *bkgvar_a, const float *varbkg_seed, const uint8_t *ok_a, @@ -425,6 +454,7 @@ __global__ void fit(const int32_t *img, const float *px_x, const float *px_y, if (i >= n) return; extern __shared__ float Pbuf[]; __shared__ float s_gs, s_num, s_den, s_I, s_wsum; + __shared__ float s_pown, s_mall, s_mown; __shared__ int s_Rf, s_Gf; if (!ok_a[i]) { if (threadIdx.x == 0) ok_o[i] = 0; return; } @@ -476,6 +506,40 @@ __global__ void fit(const int32_t *img, const float *px_x, const float *px_y, const int Rf = s_Rf, Gf = s_Gf, GfGf = Gf * Gf; const float B = fmaxf(bkg, (float) PIXEL_VARIANCE_FLOOR); + + // How much of the expected profile is cleanly this reflection's own: s_pown over the whole grid + // (XDS's MINPK quantity, what Reject cuts on) and s_mown / s_mall over the r1 disk alone, which + // is what the summation seed the runaway guard compares against actually saw. See the CPU engine. + if (p.overlap) { + if (threadIdx.x == 0) { s_pown = 0.0f; s_mall = 0.0f; s_mown = 0.0f; } + __syncthreads(); + float l_pown = 0.0f, l_mall = 0.0f, l_mown = 0.0f; + for (int k = threadIdx.x; k < GfGf; k += blockDim.x) { + const float Pp = Pbuf[k]; + if (Pp <= 0.0f) continue; + const int dx = k % Gf - Rf, dy = k / Gf - Rf; + const int x = cx + dx, y = cy + dy; + if (x < 0 || y < 0 || x >= p.W || y >= p.H) continue; + if (!valid(img[y * p.W + x])) continue; + const bool own = BraggOwnedBy(owner[y * p.W + x], i); + // Zeroing the profile here is how Exclude drops the pixel: the fit skips any cell with + // P <= 0 already, and P is not renormalised, so the fitted amplitude comes out on the + // scale of the WHOLE profile - the renormalisation is the estimator's own doing. + if (!own && p.overlap == 2) Pbuf[k] = 0.0f; + if (own) l_pown += Pp; + if ((float) (dx * dx + dy * dy) < p.r1_sq) { + l_mall += Pp; + if (own) l_mown += Pp; + } + } + atomicAdd(&s_pown, l_pown); atomicAdd(&s_mall, l_mall); atomicAdd(&s_mown, l_mown); + __syncthreads(); + if (p.overlap == 1 && s_pown < p.minpk) { + if (threadIdx.x == 0) ok_o[i] = 0; + return; + } + } + if (threadIdx.x == 0) s_I = I_seed[i]; __syncthreads(); @@ -510,8 +574,11 @@ __global__ void fit(const int32_t *img, const float *px_x, const float *px_y, float I = s_I, sigma = sqrtf(1.0f / s_den + wr * wr * bkgvar_a[i]); float var_bkg = fmaxf(0.0f, 1.0f / s_den - fmaxf(0.0f, I) + wr * wr * bkgvar_a[i]); // Guard against profile-fit runaways (see the CPU engine): fall back to the summation seed - // when the profile result diverges from it. - if (fabsf(I - I_seed[i]) > (float) PROFILE_SUMMATION_MAX_NSIGMA * sigma_seed[i]) { + // when the profile result diverges from it. Exclude has taken pixels out of both, so the + // full-profile intensity is scaled back to the disk the seed read; nothing excluded gives + // exactly 1. + const float gs = (p.overlap == 2 && s_mall > 0.0f) ? s_mown / s_mall : 1.0f; + if (fabsf(I * gs - I_seed[i]) > (float) PROFILE_SUMMATION_MAX_NSIGMA * sigma_seed[i]) { I = I_seed[i]; sigma = sigma_seed[i]; var_bkg = varbkg_seed[i]; @@ -569,6 +636,9 @@ BraggIntegrationEngineGPU::BraggIntegrationEngineGPU(const DiffractionExperiment // Radial background curve: one bin per pixel of distance from the beam, out to the far corner. // Allocated whenever the correction COULD run, so the auto mode can turn it on for an individual // image; whether it runs for a given image is decided in Run() from the per-image bkg_radial. + if (overlap != OverlapMode::Off) + d_owner = CudaDevicePtr(npixel); + if (bkg_radial || bkg_radial_auto) { n_rad = static_cast(std::ceil(r_max)) + 2; d_rad_sum = CudaDevicePtr(n_rad); @@ -649,6 +719,10 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu .beam_x = beam_x, .beam_y = beam_y, .bkg_trim = bkg_trim, // effective trim fraction (0 for stills), set by the base ctor from settings .stencil = stencil, .rad_w = rad_w, .n_kern = n_kern, + .overlap = overlap == OverlapMode::Reject ? 1 : (overlap == OverlapMode::Exclude ? 2 : 0), + .boxsum_reject = (overlap == OverlapMode::Reject && mode == IntegratorMode::BoxSum) ? 1 : 0, + .exclude_px = (overlap == OverlapMode::Exclude && mode != IntegratorMode::BoxSum) ? 1 : 0, + .claim_sq = claim * claim, .inv_claim = inv_claim, .minpk = overlap_min_peak, }; // Whether the radial correction runs for THIS image. n_rad only says the buffers exist - under @@ -657,13 +731,15 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu // Pass A: reset accumulators, mask, then box-sum. cuda_err(cudaMemsetAsync(d_mask, 0, npixel, *stream)); + if (p.overlap) + cuda_err(cudaMemsetAsync(d_owner, 0xff, sizeof(uint32_t) * npixel, *stream)); // BRAGG_OWNER_NONE if (rad_n > 0) { cuda_err(cudaMemsetAsync(d_rad_sum, 0, sizeof(float) * n_rad, *stream)); cuda_err(cudaMemsetAsync(d_rad_cnt, 0, sizeof(int) * n_rad, *stream)); } reset<<<32, 256, 0, *stream>>>(d_shell_grid, d_global_grid, d_mom, d_shell_n, d_global_n, d_invd2, GG); - mark_mask<<>>(d_px_x, d_px_y, d_mask, p, n); - boxsum<<>>(d_px_x, d_px_y, d_d, img, d_mask, p, n, + mark_mask<<>>(d_px_x, d_px_y, d_mask, d_owner, p, n); + 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, d_ok, d_strong, d_has_obs, d_invd2, d_isum, d_ninner, d_rbin, d_kbin, @@ -678,13 +754,13 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu if (mode != IntegratorMode::BoxSum) { // Pass B: learn (shell computed inline) -> build -> fit. - learn_profile<<>>(img, d_px_x, d_px_y, d_cx, d_cy, d_d, d_invd2, + learn_profile<<>>(img, d_owner, d_px_x, d_px_y, d_cx, d_cy, d_d, d_invd2, d_I, d_bkg, d_ok, d_strong, d_shell_grid, d_global_grid, d_mom, d_shell_n, d_global_n, p, n); build_profiles<<>>( d_shell_grid, d_global_grid, d_mom, d_shell_n, d_global_n, d_shell_P, d_global_P, d_sigma2_r, d_sigma2_t, p); - fit<<>>(img, d_px_x, d_px_y, d_cx, d_cy, d_d, d_invd2, + fit<<>>(img, d_owner, d_px_x, d_px_y, d_cx, d_cy, d_d, d_invd2, d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_ok, d_shell_P, d_global_P, d_sigma2_r, d_sigma2_t, d_shell_n, d_I, d_sigma, d_var_bkg, d_ok, p, n); diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h index f0a4f4df..f8e71f60 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h @@ -45,6 +45,9 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { // The learning/fit math is single precision: FP64 is heavily throttled on consumer GPUs and the // extraction is Poisson-noise limited, so float reproduces the double CPU path to ~1e-4. CudaDevicePtr d_mask; // per-pixel inner-stencil reflection mask + // Per-pixel (distance, reflection) key naming the nearest predicted centre; allocated only when + // an overlap treatment is on, so the default path costs no extra device memory. + CudaDevicePtr d_owner; CudaDevicePtr d_shell_grid, d_global_grid; // learned profile accumulators (N_SHELL*GG, GG) CudaDevicePtr d_shell_P, d_global_P; // normalised profiles (empirical mode) CudaDevicePtr d_mom; // learned 2nd moments, 3 per shell + global diff --git a/image_analysis/bragg_integration/BraggStencil.h b/image_analysis/bragg_integration/BraggStencil.h index d4191988..f87f859d 100644 --- a/image_analysis/bragg_integration/BraggStencil.h +++ b/image_analysis/bragg_integration/BraggStencil.h @@ -46,6 +46,7 @@ // ============================================================================= #include +#include #ifdef __CUDACC__ #define BRAGG_STENCIL_HD __host__ __device__ inline @@ -132,6 +133,38 @@ struct BraggStencilDist { float signal, inner, outer, rad; }; +// --- Ownership of a shared signal region --------------------------------------------------------- +// Two reflections whose centres are closer than 2*r1 have signal disks that intersect, and a pixel in +// the intersection carries both. It belongs to the NEARER centre. The union mask above cannot say +// that (it also marks a reflection's own core), so the decision is made once per image into an owner +// map: every reflection writes (quantised distance << 24) | index over its claim disk with an atomic +// minimum, so the nearest centre wins whatever order the writes arrive in and the lowest index breaks +// a tie. Eight bits of distance over the claim radius is 1/40 px at the default r2 = 6 px, far finer +// than the prediction itself. A pixel nobody claimed is nobody's neighbour, hence the owner's. +constexpr uint32_t BRAGG_OWNER_NONE = 0xffffffffu; +constexpr uint32_t BRAGG_OWNER_IDX_MASK = 0x00ffffffu; + +BRAGG_STENCIL_HD int BraggOwnerQuant(float d, float inv_claim) { + int q = (int) (d * inv_claim * 255.0f); + return q > 255 ? 255 : q; +} + +BRAGG_STENCIL_HD uint32_t BraggOwnerKey(float d, float inv_claim, int index) { + return ((uint32_t) BraggOwnerQuant(d, inv_claim) << 24) + | ((uint32_t) index & BRAGG_OWNER_IDX_MASK); +} + +BRAGG_STENCIL_HD bool BraggOwnedBy(uint32_t owner, int index) { + return owner == BRAGG_OWNER_NONE + || (owner & BRAGG_OWNER_IDX_MASK) == ((uint32_t) index & BRAGG_OWNER_IDX_MASK); +} + +// Widening that split - keeping a pixel only where no other centre is within its distance PLUS a +// margin - was built and measured, and it is worse, monotonically: on the crowded crystal the residual +// bias of the pixels that were kept grew from +0.072 at no margin to +0.087, +0.144 and +0.209 in ln +// intensity at 1, 2 and 3 px, and the merge fell apart with it. What the margin removes is the +// reflection's own profile, not the neighbour's tail, so the plain nearest-centre split is the rule. + BRAGG_STENCIL_HD BraggStencilDist BraggStencilDistances(const BraggStencil &s, float ddx, float ddy) { BraggStencilDist d; d.rad = ddx * s.ux + ddy * s.uy; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 8d801f88..d7161066 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -164,6 +164,8 @@ void print_usage() { std::cout << " --background-clip High-side clip of the background ring at mean + n*sqrt(mean) (default 4, or 3 when --bandwidth is set; 0 = off). This is the default background estimator - it rejects neighbour cores and zingers without the symmetric trim's Poisson skew bias. Ignored by --integrator boxsum" << std::endl; std::cout << " --background-radial[=on|off|auto] Correct the background ring for the CURVATURE of the radial background (default off). The signal disk and the background ring are concentric, so a background linear in position cancels between them and only curvature survives - which on a smooth ice ring reaches +26 counts on a single reflection. =auto applies it per image where that image's ice score shows a smooth powder ring, which is where a radius-only background model holds; on ice made of discrete crystallite spots there is no such ring and the correction makes the bias worse. Costs one short dot product per reflection and no extra pixel reads" << std::endl; std::cout << " --background-trim Use the old symmetric trimmed mean for the background ring instead of the clip (0<=f<0.5; 0.10 was the former default). Switches --background-clip off. A symmetric trim is biased low on Poisson data and adds ~5 counts to every partial, so this is for back compatibility only; 0 = plain ring mean. Applies whatever --bandwidth is set to" << std::endl; + std::cout << " --overlap What to do where two predicted reflections share signal pixels: off|reject|exclude (default exclude). A pixel inside two signal disks belongs to the nearer centre; before this, nothing kept a neighbour's flux out of a reflection's own disk, so on a dense pattern a crowded reflection read high. exclude drops the shared PIXELS from the profile fit, which renormalises itself, so the reflection is kept; reject instead drops the whole reflection when less than --overlap-minpk of its expected profile is cleanly its own (what XDS calls MINPK). --integrator boxsum has no profile to renormalise with, so exclude does nothing there and only reject acts" << std::endl; + std::cout << " --overlap-minpk Least fraction of a reflection's expected profile that must be cleanly its own for --overlap reject to keep it (default 0.75, XDS MINPK). With --integrator boxsum the fraction is by disk AREA instead, which cuts harder" << std::endl; std::cout << " --integrator Spot integrator boxsum|gaussian|empirical (default: gaussian profile-fit; boxsum is the classical fallback)" << std::endl; std::cout << " --simple-stills stills: treat every reflection as a full (p=1, single-pass scale/merge); disables the default physical partiality post-refinement" << std::endl; std::cout << " -q, --azim-q-spacing Azimuthal-integration Q bin spacing (1/A) (default: 0.01)" << std::endl; @@ -214,6 +216,8 @@ enum { OPT_INTEGRATION_RADIUS, OPT_INTEGRATION_STENCIL, OPT_BACKGROUND_TRIM, + OPT_OVERLAP, + OPT_OVERLAP_MINPK, OPT_MAX_HKL, OPT_INTEGRATION_HIGH_RES, OPT_REJECT_OUTLIERS, @@ -332,6 +336,8 @@ static option long_options[] = { {"integration-radius", required_argument, nullptr, OPT_INTEGRATION_RADIUS}, {"integration-stencil", required_argument, nullptr, OPT_INTEGRATION_STENCIL}, {"background-trim", required_argument, nullptr, OPT_BACKGROUND_TRIM}, + {"overlap", required_argument, nullptr, OPT_OVERLAP}, + {"overlap-minpk", required_argument, nullptr, OPT_OVERLAP_MINPK}, {"max-hkl", required_argument, nullptr, OPT_MAX_HKL}, {"integration-high-resolution", required_argument, nullptr, OPT_INTEGRATION_HIGH_RES}, {"integrator", required_argument, nullptr, OPT_INTEGRATOR}, @@ -630,6 +636,8 @@ static int RunRugnux(int argc, char **argv) { std::optional integration_radius_arg; std::optional integration_stencil_arg; // --integration-stencil: ring elongation, in sigma std::optional background_trim_arg; // --background-trim: background-ring trimmed-mean fraction + std::optional overlap_arg; // --overlap: treatment of shared signal pixels + std::optional overlap_minpk_arg; // --overlap-minpk: XDS-like MINPK threshold std::optional max_hkl_arg; // --max-hkl: half-width of the predicted hkl box std::optional integration_d_min_arg; // --integration-high-resolution; unset = detector reach std::optional integrator_mode; // --integrator boxsum|gaussian|empirical @@ -981,6 +989,15 @@ static int RunRugnux(int argc, char **argv) { case OPT_BACKGROUND_TRIM: background_trim_arg = parse_double_arg(optarg, "--background-trim", logger); break; + case OPT_OVERLAP: + if (strcmp(optarg, "off") == 0) overlap_arg = OverlapMode::Off; + else if (strcmp(optarg, "reject") == 0) overlap_arg = OverlapMode::Reject; + else if (strcmp(optarg, "exclude") == 0) overlap_arg = OverlapMode::Exclude; + else { logger.Error("--overlap expects off|reject|exclude"); return 1; } + break; + case OPT_OVERLAP_MINPK: + overlap_minpk_arg = parse_double_arg(optarg, "--overlap-minpk", logger); + break; case OPT_MAX_HKL: max_hkl_arg = parse_number_arg(optarg, "--max-hkl", logger, 1, 511); break; @@ -1913,6 +1930,20 @@ static int RunRugnux(int argc, char **argv) { BROADBAND_BACKGROUND_CLIP_NSIGMA); } + if (overlap_arg || overlap_minpk_arg) { + BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings(); + if (overlap_arg) + bis.Overlap(*overlap_arg); + if (overlap_minpk_arg) + bis.OverlapMinPeak(static_cast(*overlap_minpk_arg)); + experiment.ImportBraggIntegrationSettings(bis); + const auto ovl = bis.GetOverlap(); + logger.Info("Overlapping signal regions: {}", ovl == OverlapMode::Reject + ? fmt::format("reject below {:.2f} of the profile", bis.GetOverlapMinPeak()) + : (ovl == OverlapMode::Exclude ? std::string("exclude the shared pixels") + : std::string("off"))); + } + if (background_radial_given) { BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings(); bis.BackgroundRadialCorrection(background_radial_arg); diff --git a/tests/BraggIntegrationEngineGPUTest.cpp b/tests/BraggIntegrationEngineGPUTest.cpp index 020cec72..e869db18 100644 --- a/tests/BraggIntegrationEngineGPUTest.cpp +++ b/tests/BraggIntegrationEngineGPUTest.cpp @@ -38,7 +38,10 @@ Reflection MakeReflection(float x, float y, float d, int hkl) { return r; } -Scene BuildScene(size_t width, size_t height, int spacing = 60) { +// companion_dx > 0 puts a second spot that many pixels beside every grid spot, so their r1 signal +// disks share pixels while the background rings still see clean sky - which is what a dense pattern +// actually looks like (crowded along one reciprocal axis, sparse across it). +Scene BuildScene(size_t width, size_t height, int spacing = 60, float companion_dx = 0.0f) { Scene s; s.width = width; s.height = height; @@ -66,6 +69,19 @@ Scene BuildScene(size_t width, size_t height, int spacing = 60) { } const float d = 1.4f + 0.12f * static_cast((gx + gy) % 12); // 1.4..2.72 A s.predicted.push_back(MakeReflection(cx, cy, d, hkl++)); + if (companion_dx > 0.0f) { + const float ccx = cx + companion_dx; + for (int dy = -6; dy <= 6; ++dy) + for (int dx = -6; dx <= 6; ++dx) { + const int x = static_cast(std::lround(ccx)) + dx; + const int y = static_cast(std::lround(cy)) + dy; + if (x < 0 || y < 0 || x >= static_cast(width) || y >= static_cast(height)) continue; + const double ex = x - ccx, ey = y - cy; + const double g = 0.6 * amp * std::exp(-(ex * ex + ey * ey) / (2.0 * sigma * sigma)); + s.image[y * width + x] += static_cast(std::lround(g)); + } + s.predicted.push_back(MakeReflection(ccx, cy, d, hkl++)); + } } } @@ -85,7 +101,8 @@ DiffractionExperiment MakeExperiment(IntegratorMode mode, std::optional b bool radial = false, const DetectorSetup &det = DetJF(2), float stencil_k = 0.0f, - float r1 = 0.0f, float r2 = 0.0f, float r3 = 0.0f) { + float r1 = 0.0f, float r2 = 0.0f, float r3 = 0.0f, + OverlapMode overlap = OverlapMode::Off) { DiffractionExperiment experiment(det); // DetJF(2) (small) keeps the correctness test fast experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV) .BeamX_pxl(400.0f).BeamY_pxl(400.0f); @@ -100,6 +117,7 @@ DiffractionExperiment MakeExperiment(IntegratorMode mode, std::optional b settings.BackgroundTrimFraction(0.10f); settings.BackgroundRadialCorrection(radial); settings.StencilKSigma(stencil_k); + settings.Overlap(overlap); experiment.ImportBraggIntegrationSettings(settings); return experiment; } @@ -107,15 +125,17 @@ DiffractionExperiment MakeExperiment(IntegratorMode mode, std::optional b void 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) { + float r1 = 0.0f, float r2 = 0.0f, float r3 = 0.0f, + OverlapMode overlap = OverlapMode::Off, float companion_dx = 0.0f) { const DiffractionExperiment experiment = - MakeExperiment(mode, bandwidth_fwhm, clip_nsigma, radial, DetJF(2), stencil_k, r1, r2, r3); + MakeExperiment(mode, bandwidth_fwhm, clip_nsigma, radial, DetJF(2), stencil_k, r1, r2, r3, + overlap); const size_t width = experiment.GetXPixelsNum(); const size_t height = experiment.GetYPixelsNum(); const size_t npixel = experiment.GetPixelsNum(); REQUIRE(npixel == width * height); - const Scene scene = BuildScene(width, height, spacing); + const Scene scene = BuildScene(width, height, spacing, companion_dx); REQUIRE(scene.image.size() == npixel); REQUIRE(scene.predicted.size() > 60); @@ -197,6 +217,27 @@ TEST_CASE("BraggIntegrationEngineGPU_MatchesCPU") { CompareCpuVsGpu(IntegratorMode::ProfileGaussian, 0.04f, 0.0f, false, 120, 4.0f, 6.0f, 8.0f, 12.0f); } SECTION("ProfileEmpirical") { CompareCpuVsGpu(IntegratorMode::ProfileEmpirical, std::nullopt); } + // Overlap treatment: companions 4 px apart put each reflection's centre inside its neighbour's + // signal disk, so the owner map, the excluded pixels and the profile fraction the two modes act on + // all have to come out the same in both engines - the ownership atomic in particular is settled by + // an atomicMin on the GPU and a serial minimum on the CPU. + SECTION("ProfileGaussian overlap exclude") { + CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 60, 0.0f, + 0.0f, 0.0f, 0.0f, OverlapMode::Exclude, 4.0f); + } + SECTION("ProfileGaussian overlap reject") { + CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 60, 0.0f, + 0.0f, 0.0f, 0.0f, OverlapMode::Reject, 4.0f); + } + SECTION("BoxSum overlap reject") { + CompareCpuVsGpu(IntegratorMode::BoxSum, std::nullopt, 4.0f, false, 60, 0.0f, + 0.0f, 0.0f, 0.0f, OverlapMode::Reject, 4.0f); + } + // Nothing shares a pixel at this spacing, so an overlap treatment has to leave the result alone. + SECTION("ProfileGaussian overlap inert") { + CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 60, 0.0f, + 0.0f, 0.0f, 0.0f, OverlapMode::Exclude); + } SECTION("ProfileGaussian mono trim") { CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 0.0f); } // The radial background curvature correction is computed independently in the two engines // (host loop vs radial_correct kernel), so it needs its own parity coverage. @@ -217,8 +258,12 @@ TEST_CASE("BraggIntegrationEngineGPU_Benchmark", "[.][bragg_bench]") { WARN("No CUDA GPU present. Skipping benchmark"); return; } + // The overlap treatment is priced here too: it adds an owner map over the whole frame plus one + // atomic per claimed pixel, so what it costs is a property of the frame more than of the crowding. + for (OverlapMode ovl : {OverlapMode::Off, OverlapMode::Reject, OverlapMode::Exclude}) { const DiffractionExperiment experiment = MakeExperiment(IntegratorMode::ProfileGaussian, std::nullopt, - 4.0f, false, DetJF4M()); + 4.0f, false, DetJF4M(), 0.0f, 0.0f, 0.0f, 0.0f, + ovl); const size_t width = experiment.GetXPixelsNum(); const size_t height = experiment.GetYPixelsNum(); const size_t npixel = experiment.GetPixelsNum(); @@ -226,7 +271,7 @@ TEST_CASE("BraggIntegrationEngineGPU_Benchmark", "[.][bragg_bench]") { auto stream = std::make_shared(); BraggIntegrationEngineGPU gpu(experiment, stream); - for (int spacing : {28, 40, 60, 90}) { + for (int spacing : {28, 60}) { const Scene scene = BuildScene(width, height, spacing); const size_t nrefl = scene.predicted.size(); @@ -254,9 +299,10 @@ TEST_CASE("BraggIntegrationEngineGPU_Benchmark", "[.][bragg_bench]") { const auto c1 = std::chrono::steady_clock::now(); const double cpu_ms = std::chrono::duration(c1 - c0).count(); - WARN(width << "x" << height << " | " << nrefl << " refl (" << observed / iters - << " obs) | GPU " << ms << " ms | CPU " << cpu_ms << " ms (" << cpu_observed - << " obs) | speedup " << cpu_ms / ms << "x"); + WARN((int) ovl << " | " << width << "x" << height << " | " << nrefl << " refl (" + << observed / iters << " obs) | GPU " << ms << " ms | CPU " << cpu_ms << " ms (" + << cpu_observed << " obs) | speedup " << cpu_ms / ms << "x"); + } } }