From 4e6eb18d937c313cfa06770fa5e6dde57772b327 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 26 Aug 2026 00:20:41 +0200 Subject: [PATCH] Integration: fall back to the fixed radius on a pattern too dense for it 48008e144 widens the signal radius on crystals with wide spots. On one battery crystal - simultaneously the widest-spot and among the highest in mosaicity - the wider radius left its neighbours' background rings with too few clean pixels and cost 28.5% of its observations. Pass 1 now measures how often that happens and, above a bound, hands pass 2 the settings from before the pre-scan widened them. The obvious quantity does not work. On the total rate of reflections dropped for a starved ring, the losing crystal reads 4.08% and the rule's four biggest winners read 1.89-2.04% - and a crystal at the shipped radius reads 2.37%, above all of them. Re-running the winners at r1 = 4 shows why: they read 2.20-2.32% there too, and widening moves them down. That floor is module gaps, the beam stop and the resolution mask, which are properties of the detector and do not move with the radius. So the counter separates the two. A ring is neighbour-starved when it would have kept more than five pixels but for the pixels a neighbouring reflection's signal region occupies. That is exact rather than estimated: the reflection mask marks the disk inside r2 and the ring is everything outside it, so a masked ring pixel always belongs to some other reflection's core. The separation goes from a factor of 2 to a factor of 13 - over the twelve crystals the radius moves, the rate is 0.000 five times, 0.001 three times, then 0.004, 0.235, 0.315 and 4.082 - and the bound is the log-space midpoint of that one gap, 0.0113, a factor 3.6 clear of the nearest measurement on either side. Predicted reflection spacing does not separate them at all: the losing crystal is 19th of 38, a winner sits at 21.9 px, and the loosest pattern in the battery starves 1.93% of its rings. Battery: the space group is identical on all 38 and the merged .hkl is byte-identical on 37, so it is inert wherever it does not fire. On the one crystal it fires on, is up 18.3%, R_meas down 29.4%, observations up 6.3%, CC1/2 0.944 to 0.974, and its two empty top shells come back as numbers. Its indexing rate, refined distance, beam centre and cell are bit-identical between the two arms, so this is the guard and not the two-pass gate. The counters are a shared channel through both engines, summed across workers and logged once per pass; on the GPU it is one atomic add per dropped reflection. The profile-fit runaway guard reports on the same channel, which is the first measurement of its trip rate. This does not recover that crystal fully. With the adaptive radius on, pass 1 reaches a different lattice and pass 2 indexes 21% fewer frames - which happens before the measurement this guard reads exists, and is unaffected by it. At matched indexing rate the guard recovers 96% of the baseline's observations against 90.5% without it. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01CHMmeM1d489zvNFT7ZMN2P --- docs/CHANGELOG.md | 1 + image_analysis/MXAnalysisWithoutFPGA.cpp | 4 ++ image_analysis/MXAnalysisWithoutFPGA.h | 4 ++ .../BraggIntegrationEngine.h | 27 +++++++++ .../BraggIntegrationEngineCPU.cpp | 22 ++++++- .../BraggIntegrationEngineGPU.cu | 50 +++++++++++++--- .../BraggIntegrationEngineGPU.h | 7 +++ rugnux/Rugnux.cpp | 58 +++++++++++++++++++ rugnux/Rugnux.h | 7 +++ rugnux/SpotWidth.h | 16 +++++ 10 files changed, 186 insertions(+), 10 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index cd5295cb..2e148d10 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -9,6 +9,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as * rugnux: a lattice centring the data could not test - the crystal was integrated on the primitive sub-cell, so the reflections the centring extinguishes were never measured - is marked `UNTESTED` in the space-group candidate table and, where it is adopted, is warned about as coming from the lattice metric rather than from the intensities. * rugnux: the space-group search prints the twin-law disagreement H for every operator it tested and the H ratio of the point group it adopted, on every run, instead of only when that ratio refuses a promotion. * rugnux: `--mode scale` works on a `_process.h5` whose space group came from re-seating the lattice; a file written before this stops with a message naming the two cells instead of failing inside the merge. +* rugnux: on a pattern too dense for the widened radius - where neighbouring reflections leave more than 1.1% of the reflections without a background ring - the second integration pass goes back to the fixed 4 px radius, and says so. * rugnux: on rotation data the integration signal radius is set from how wide the crystal's own spots are, measured in the pre-scan, instead of the fixed 4 px; `--adaptive-integration-radius=off` restores the fixed radius, and an explicit `--integration-radius` still overrides both. * rugnux: a directional diffraction limit that is the edge of the measured data rather than the crystal's own limit is marked as such - with a `<` in the report and in `ANISOTROPY_D_MIN_CENSORED`, and in the mmCIF - so `ANISOTROPY_D_MIN_SPREAD` is not read as a measurement when it is a lower bound. * rugnux: the anisotropy verdict line names which of `ANISOTROPY_DELTA_B` and `ANISOTROPY_DELTA_B_LINEAR` it is quoting, says which of the two to act on, and says why the second can be the larger. diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index e430407a..0c964db7 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -357,6 +357,10 @@ void MXAnalysisWithoutFPGA::RunROIOnly(DataMessage &output) { roi->Run(*preprocessor_buffer, output.roi); } +BraggIntegrationCounts MXAnalysisWithoutFPGA::BraggCounts() const { + return bragg_engine ? bragg_engine->Counts() : BraggIntegrationCounts{}; +} + void MXAnalysisWithoutFPGA::UpdateMaskResolution(const SpotFindingSettings &settings) { mask_low_res = settings.low_resolution_limit; mask_high_res = settings.high_resolution_limit; diff --git a/image_analysis/MXAnalysisWithoutFPGA.h b/image_analysis/MXAnalysisWithoutFPGA.h index e77955f1..54ec41a4 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.h +++ b/image_analysis/MXAnalysisWithoutFPGA.h @@ -93,6 +93,10 @@ public: void RebuildROI(); void AnalyzeROIOnly(DataMessage &output); void RunROIOnly(DataMessage &output); + + // What this worker's Bragg integrator counted (BraggIntegrationCounts). Each worker builds its own + // analysis, so a caller that wants the run's totals sums this over the workers it started. + [[nodiscard]] BraggIntegrationCounts BraggCounts() const; }; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index 7a982c82..3280ae62 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -104,6 +104,24 @@ 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 +// 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 + // 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 + // "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 + // costs data adds 4.1 percentage points on top of a floor of 0.02%. + uint64_t bkg_starved_by_neighbour = 0; + uint64_t profile_fallback = 0; // the profile fit disagreed with the box sum; the box sum was kept +}; + // One reflection's extracted intensity, produced by the derived engine and turned into a // Reflection by Finalize() (which owns the polarization correction and scale bookkeeping). struct BraggFitResult { @@ -199,6 +217,10 @@ protected: DiffractionGeometry geom; // kept for the per-reflection polarization correction std::optional polarization; + // Accumulated over every image this engine has run. An engine belongs to one worker thread, so + // these are plain counters and the caller sums over the engines it made. + BraggIntegrationCounts counts; + // Assemble output reflections from the per-reflection fit results (polarization + scale corr). std::vector Finalize(const std::vector &predicted, size_t npredicted, const std::vector &results, @@ -222,4 +244,9 @@ public: // radial buffers were never allocated, so the CPU would correct and the GPU would not. void BackgroundRadial(bool on) { bkg_radial = on && bkg_radial_built; } [[nodiscard]] bool IsBackgroundRadialAuto() const { return bkg_radial_auto; } + + // What this engine has counted since it was built (BraggIntegrationCounts). The GPU engine keeps + // its counts on the device and brings them back here, so this is a synchronising call - ask for it + // once a pass, not once an image. + [[nodiscard]] virtual BraggIntegrationCounts Counts() const { return counts; } }; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index a6d243d6..015053f6 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -53,6 +53,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, std::vector results(npredicted); if (npredicted == 0) return Finalize(predicted, npredicted, results, image_number); + counts.predicted += npredicted; const int W = static_cast(xpixel), H = static_cast(ypixel); const bool do_clip = bkg_clip_nsigma > 0.0f && mode != IntegratorMode::BoxSum; @@ -166,6 +167,11 @@ 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; bkg_vals.clear(); for (int y = y0; y <= y1; ++y) for (int x = x0; x <= x1; ++x) { @@ -188,7 +194,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, I_sum_y += static_cast(y) * px; ++n_inner_valid; } else if (d.inner >= r2_sq && d.outer < r3_sq) { - if (refl_mask[y * W + x]) continue; + if (refl_mask[y * W + x]) { ++n_bkg_neighbour; continue; } if (!valid(px)) continue; bkg_sum += static_cast(px); if (bkg_trim_frac > 0.0) bkg_vals.push_back(px); @@ -205,7 +211,18 @@ 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; - if ((full || mode != IntegratorMode::BoxSum) && n_bkg > 5) { + // 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 + // 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; + } + if (keep_partial && n_bkg > 5) { out.bkg = bkg_sum / n_bkg; if (bkg_trim_frac > 0.0 && bkg_vals.size() > 5 && bkg_vals.size() <= static_cast(bragg_engine::BKG_TRIM_MAX)) { @@ -555,6 +572,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, I = rh.I; sigma = rh.sigma; var_bkg = rh.var_bkg; + ++counts.profile_fallback; } // Carry the Pass-A box-sum intensity-weighted centroid (observed spot position) through the // profile path too - post-refinement uses it as the observed position (beam-centre / distance). diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 76d73877..2cc96b2a 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -42,6 +42,14 @@ __device__ inline bool valid(int32_t v) { return v != INT32_MIN && v != INT32_MA // Learned second moments, (sum v*rad^2, sum v*tan^2, sum v) per shell plus one global slot at N_SHELL. constexpr int MOM_STRIDE = 3; +// Slots of the device counter array behind BraggIntegrationCounts. They stay on the device for the +// engine's lifetime and come back only when Counts() is asked for, so an image costs nothing beyond +// the atomics themselves - and both are on paths taken by a small minority of reflections. +constexpr int COUNT_BKG_STARVED = 0; +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. // The same sweep fills the owner map when an overlap treatment is on: a pixel two signal regions @@ -149,12 +157,14 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, 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, float *isum_o, int *ninner_o, int *rbin_o, int *kbin_o, - unsigned long long *rad_sum, int *rad_cnt, int n_rad) { + unsigned long long *rad_sum, int *rad_cnt, int n_rad, + unsigned long long *counts) { const int i = blockIdx.x; if (i >= n) return; __shared__ unsigned long long s_Isum, s_Ix, s_Iy; __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__ unsigned long long s_bkgsum; __shared__ int s_accept, s_full; __shared__ double s_bkg, s_thr; @@ -179,7 +189,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; - s_ndisk = 0; s_nown = 0; + s_ndisk = 0; s_nown = 0; s_nbkg_nb = 0; } __syncthreads(); @@ -195,7 +205,7 @@ __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, l_nd = 0, l_no = 0; + int l_ni = 0, l_niv = 0, l_nb = 0, l_nd = 0, l_no = 0, l_nbnb = 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; @@ -214,7 +224,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, const int32_t px = img[y * p.W + x]; 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) { - if (mask[y * p.W + x]) continue; + if (mask[y * p.W + x]) { ++l_nbnb; continue; } const int32_t px = img[y * p.W + x]; if (!valid(px)) continue; l_bkg += px; ++l_nb; @@ -229,6 +239,7 @@ __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); __syncthreads(); for (int t = threadIdx.x; t < p.rad_w; t += blockDim.x) { s_radv[t] = 0; s_radn[t] = 0; } @@ -239,6 +250,13 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, // (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; + // 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); + } 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 @@ -567,7 +585,8 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x const float *bkgvar_a, const float *varbkg_seed, const uint8_t *ok_a, const float *shell_P, const float *global_P, const float *sigma2_r, const float *sigma2_t, const int *shell_n, - float *I_o, float *sigma_o, float *varbkg_o, uint8_t *ok_o, BraggGpuParams p, int n) { + float *I_o, float *sigma_o, float *varbkg_o, uint8_t *ok_o, BraggGpuParams p, int n, + unsigned long long *counts) { const int i = blockIdx.x; if (i >= n) return; extern __shared__ float Pbuf[]; @@ -737,6 +756,7 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x I = I_seed[i]; sigma = sigma_seed[i]; var_bkg = varbkg_seed[i]; + atomicAdd(counts + COUNT_PROFILE_FALLBACK, 1ull); } I_o[i] = I; sigma_o[i] = sigma; varbkg_o[i] = var_bkg; ok_o[i] = 1; } else ok_o[i] = 0; @@ -759,8 +779,10 @@ BraggIntegrationEngineGPU::BraggIntegrationEngineGPU(const DiffractionExperiment d_sigma2_t(bragg_engine::N_SHELL + 1), d_shell_n(bragg_engine::N_SHELL), d_global_n(1), - d_invd2(2) { + d_invd2(2), + d_counts(COUNT_SLOTS) { threads = 128; + cuda_err(cudaMemset(d_counts, 0, sizeof(unsigned long long) * COUNT_SLOTS)); // Fit profile grid: R for empirical / box, up to 3R (radially elongated) for the Gaussian. const int max_Rf = empirical ? R : 3 * R; @@ -856,6 +878,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu std::vector results(npredicted); if (image.size() != npixel || npredicted == 0) return Finalize(predicted, npredicted, results, image_number); + counts.predicted += npredicted; const int32_t *img = image.getGPUBuffer(); if (img == nullptr) @@ -925,7 +948,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu 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, - d_rad_sum, d_rad_cnt, rad_n); + d_rad_sum, d_rad_cnt, rad_n, d_counts); // Correct the flat annulus background for the curvature of the radial background before anything // downstream (profile fit, variance) reads it. @@ -945,7 +968,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu 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); + d_I, d_sigma, d_var_bkg, d_ok, p, n, d_counts); } // Nothing below reads the mask or the owner map, so put them back now, over the boxes that were @@ -983,3 +1006,14 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu } return Finalize(predicted, npredicted, results, image_number); } + +BraggIntegrationCounts BraggIntegrationEngineGPU::Counts() const { + BraggIntegrationCounts out = counts; // `predicted` is counted on the host, as on the CPU engine + unsigned long long host[COUNT_SLOTS] = {0, 0, 0}; + cuda_err(cudaMemcpyAsync(host, d_counts, sizeof(host), cudaMemcpyDeviceToHost, *stream)); + cuda_err(cudaStreamSynchronize(*stream)); + out.bkg_starved = host[COUNT_BKG_STARVED]; + out.bkg_starved_by_neighbour = host[COUNT_BKG_STARVED_NEIGHBOUR]; + out.profile_fallback = host[COUNT_PROFILE_FALLBACK]; + return out; +} diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h index 506d67ea..a95f78c4 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h @@ -64,6 +64,9 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { CudaDevicePtr d_sigma2_r, d_sigma2_t; // radial/tangential widths, N_SHELL + global CudaDevicePtr d_shell_n, d_global_n; CudaDevicePtr d_invd2; // [min,max] inv-d^2 as monotonic bit patterns + // BraggIntegrationCounts, accumulated on the device so an image costs no transfer; brought back + // only when Counts() is asked for. + CudaDevicePtr d_counts; // --- host staging (copied back once per frame) --- // Pinned, like every other engine's staging: a copy out of pageable memory does not return until the @@ -80,4 +83,8 @@ public: std::vector Run(const ImagePreprocessorBuffer &image, const std::vector &predicted, size_t npredicted, int64_t image_number) override; + + // Brings the two device counters back before answering. Synchronises the stream, so ask once a + // pass rather than once an image. + [[nodiscard]] BraggIntegrationCounts Counts() const override; }; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 3ab6108b..7c81a399 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -549,6 +549,9 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru const float r2 = r1 + 2.0f; const float r3 = std::sqrt(r2 * r2 + 133.0f); BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings(); + // Remembered only where the radius actually moved, so the guard in RunAllPasses has + // something to give back and cannot announce a fallback that changes nothing. + if (r1 != bis.GetR1()) bragg_before_adaptive_ = bis; bis.R1(r1).R2(r2).R3(r3); experiment_.ImportBraggIntegrationSettings(bis); logger.Info("Spot width: r80 = {:.2f} px at {:.0f} A ({} spots) => integration radii " @@ -929,6 +932,30 @@ ProcessResult Rugnux::RunAllPasses(RugnuxObserver *observer) { gonio_snapshot->GetIncrement_deg(), gonio_snapshot->GetIncrement_deg() * *prepass_rotation_scale_); } + // The integration radius the pre-scan chose from the recorded spot width, given back where pass 1 + // measured that the pattern could not take it. r2 = r1 + 2 is also the inner edge of the + // background ring, so widening the signal disk pushes that ring out into the neighbours; a + // reflection whose ring is left with five or fewer clean pixels has no background and is dropped + // whole. On a pattern dense enough in three dimensions at once that costs a quarter of the + // observations - and NOTHING available before integrating predicts it: the crystal it happens to + // has a mid-table predicted spot spacing, wider than crystals that survive r1 = 12. What does + // measure it is the integrator counting its own drops, so the decision is taken on pass 1's + // measurement and pass 2 - the canonical pass - integrates at the radius that keeps the data. + if (bragg_before_adaptive_ && bkg_starved_fraction_ + && *bkg_starved_fraction_ > spot_width::BKG_STARVED_MAX_FRACTION) { + logger.Info("Two-pass: at the widened integration radius r1={:.1f} the neighbouring " + "reflections left {:.2f}% of the predicted reflections without a background " + "ring (bound {:.2f}%), so the second pass integrates at r1={:.1f} r2={:.1f} " + "r3={:.2f}", + experiment_.GetBraggIntegrationSettings().GetR1(), + 100.0 * *bkg_starved_fraction_, + 100.0 * spot_width::BKG_STARVED_MAX_FRACTION, + bragg_before_adaptive_->GetR1(), bragg_before_adaptive_->GetR2(), + bragg_before_adaptive_->GetR3()); + experiment_.ImportBraggIntegrationSettings(*bragg_before_adaptive_); + bragg_before_adaptive_.reset(); + } + // Space group: the second pass RE-INDEXES DE NOVO (clear the group here) so the indexer's pseudo- // symmetry safeguards recover the true cell - reusing pass-1's group in the indexer forces build_sr // onto a wrong / doubled cell (a huge oblique cell collapses; a pseudo-centred cell doubles). Pass-1's @@ -1789,6 +1816,12 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b }; std::atomic total_uncompressed_bytes = 0; + // What the integrator did to the reflections it was handed, summed over the workers as each one + // finishes (BraggIntegrationCounts). Reading an engine's counts synchronises its stream, so it is + // done once per worker at the end rather than per image. + std::atomic bragg_predicted = 0, bragg_starved = 0, bragg_starved_nb = 0, + bragg_fallback = 0; + // Calibration by spots: the pooled spot list of the whole run. Ring clustering is O(n^2) in it and // the Hough circle centre O(n^3) in its first few hundred, so each image contributes a fair share of // a fixed budget instead of everything it found - a powder ring is over-determined either way. @@ -1908,6 +1941,12 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const int done = finished_count.fetch_add(1) + 1; if (observer) observer->OnProgress(done, images_to_process); } + + const auto bragg = analysis.BraggCounts(); + bragg_predicted += bragg.predicted; + bragg_starved += bragg.bkg_starved; + bragg_starved_nb += bragg.bkg_starved_by_neighbour; + bragg_fallback += bragg.profile_fallback; }; if (observer) @@ -1955,6 +1994,25 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.cancelled = cancelled_; result.images_processed = finished_count.load(); + // 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. + 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", + experiment_.GetBraggIntegrationSettings().GetR1(), + experiment_.GetBraggIntegrationSettings().GetR2(), + experiment_.GetBraggIntegrationSettings().GetR3(), + bragg_predicted.load(), frac(bragg_starved), frac(bragg_starved_nb), + frac(bragg_fallback)); + } + // Every image failing is a total failure, not a run that produced nothing: it used to be reported // only as per-image log lines while the process still exited 0 with no output file. if (!cancelled_ && images_to_process > 0 && result.images_processed == 0) diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 05094048..10fd7b11 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -309,6 +309,13 @@ class Rugnux { // The recorded spot width the pre-scan measured (config_.adaptive_integration_radius), kept so the // two-pass rotation run measures it once and both passes integrate at the same radius. bool spot_width_measured_ = false; + // The integration radii in force before the measured spot width widened them, so a widening that + // turns out to cost data can be given back exactly rather than to a written-down default. + std::optional bragg_before_adaptive_; + // Fraction of the predicted reflections the last completed pass dropped because NEIGHBOURING + // reflections left their background ring below six clean pixels. Measured, not predicted - it is + // the only thing that tells a radius the pattern can take from one it cannot (see RunAllPasses). + std::optional bkg_starved_fraction_; // Pre-scan: read a spread sample of frames and take three things off them - the shadow of the beam // stop and its holder, added to the pixel mask (config_.detect_beam_stop), the beam centre diff --git a/rugnux/SpotWidth.h b/rugnux/SpotWidth.h index 66946401..5b62b2ca 100644 --- a/rugnux/SpotWidth.h +++ b/rugnux/SpotWidth.h @@ -62,6 +62,22 @@ struct FluxCurve { // neighbour-ownership radius and the inner edge of the background ring, and past 6 px a dense // pattern starts losing reflections whose ring falls below six clean pixels. [[nodiscard]] float R1ForWidth(float r80); + +// Most of the predicted reflections a widened radius may leave without a background ring, measured +// (BraggIntegrationCounts::bkg_starved_by_neighbour) rather than predicted. Above this the widening +// has cost the pattern more than it can be worth and the shipped radius is restored. +// +// The bound is read off the rotation battery, not chosen. At the radius the rule picks, the eleven +// crystals whose observation count survives the widening sit between 0.000 % and 0.315 %; the one +// crystal that loses a quarter of its observations sits at 4.082 %. That is the only gap in the +// distribution, and this is its midpoint in log space - the maximum-margin split, a factor 3.6 +// clear of the nearest measurement on either side. +// +// It has to be the NEIGHBOUR count and not the total: the total also carries the detector's own +// floor (module gaps, the beam stop, the resolution mask), which reaches 2.3 % on some geometries, +// does not move with r1, and would put four of the rule's six wins within a factor of two of the +// loss. Splitting the two apart takes the separation from 2x to 13x. +constexpr double BKG_STARVED_MAX_FRACTION = 0.0113; } // namespace spot_width // Every isolated, strong, fully readable spot of one image, as a normalised encircled-flux curve,