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,