From 61def00a1a0373eede2a4a4f129021929bd7f480 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 22 Oct 2026 09:20:52 +0200 Subject: [PATCH] Stills: background clip follows the radial curve; signal disk from the lattice spacing Two stills-only changes to the integration stencil. 1. Where the whole-frame radial curve is in use (the correction on), the ring's high-side clip follows it: a ring pixel is clipped against the threshold plus the curve's rise from the reflection's radius to its own (BraggIntegrationEngine::RadialClipOffset, rounded to whole counts so CPU, CUDA and Metal keep the same pixels). On a background falling steeply away from the beam the flat mean + 4 sqrt(mean) threshold cut the beam-side half of the ring as outliers, so the background read low and reflections near the beam stop read 3-150x too strong. The empty-aperture test's corrected column clips the same way, so the per-run correction decision sees it. Measured with r1 6 (dimple R_free, battery arguments, all frames): lcls222 0.196 -> 0.1856, pal86 0.216 -> 0.1997, sacla236 0.2559 -> 0.2559, alvra213_int16 0.1986 -> 0.1986. 2. The signal disk is taken from the lattice: it must not reach into the next node's spot. The closest two nodes can come on the detector is the shortest reciprocal-lattice vector the centring allows, s = F_px tan(2 asin(lambda / 2 d_max)), and a spot reaches 2 r80 (measured in the pre-scan); r1 = the largest of 4, 5, 6 px not above s - 2 r80, never below 4, never above the disk the run starts with (6). The empty apertures are measured for each candidate disk and the ring is chosen on the chosen disk's tally. Pixel detectors keep 6 (sacla236 s 22 px, alvra213 20 px); the binned CCDs with a long axis shrink to 4 (lcls222 5.4 px, pal86 6.7 px), which the radius scans there measured as the optimum; pal90 (s 21 px) and lcls242_s (15.5 px) keep 6. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu --- .../bragg_integration/BackgroundAperture.cpp | 47 ++++++++-- .../bragg_integration/BackgroundAperture.h | 11 +++ .../BraggIntegrationEngine.h | 14 +++ .../BraggIntegrationEngineCPU.cpp | 2 +- .../BraggIntegrationEngineGPU.cu | 20 ++-- .../BraggIntegrationEngineMetal.metal | 10 +- .../BraggIntegrationEngineMetal.mm | 9 ++ image_analysis/bragg_integration/BraggPassA.h | 1 + rugnux/Rugnux.h | 2 +- rugnux/RugnuxPreScan.cpp | 94 ++++++++++++------- tests/BackgroundApertureTest.cpp | 18 ++++ tests/BraggIntegrationEngineCPUTest.cpp | 43 +++++++++ tests/BraggIntegrationEngineGPUTest.cpp | 7 +- 13 files changed, 228 insertions(+), 50 deletions(-) diff --git a/image_analysis/bragg_integration/BackgroundAperture.cpp b/image_analysis/bragg_integration/BackgroundAperture.cpp index 4db9a6c94..781914c0d 100644 --- a/image_analysis/bragg_integration/BackgroundAperture.cpp +++ b/image_analysis/bragg_integration/BackgroundAperture.cpp @@ -7,6 +7,7 @@ #include #include "../../common/JFJochMath.h" +#include "SystematicAbsence.h" std::vector background_aperture::CandidateR3(float r2, float max_r3) { std::vector out; @@ -84,6 +85,12 @@ void background_aperture::MeasureFrame(ImagePreprocessorBuffer &image, int width } if (!readable || disk.empty()) continue; + // The ring background at the aperture's centre: what the clip with the correction follows the + // curve from (BraggIntegrationEngine::RadialClipOffset). + double centre_b = 0.0; + for (const auto &p : disk) + if (p.d2 == 0.0f) + centre_b = p.b; double disk_v = 0.0, disk_b = 0.0; for (const auto &p : disk) { disk_v += p.v; @@ -112,17 +119,26 @@ void background_aperture::MeasureFrame(ImagePreprocessorBuffer &image, int width } const double mean = sum / n; const double thr = mean + clip_nsigma * std::sqrt(gain * std::max(mean, 1.0)); - double kept = 0.0; - int n_kept = 0; - for (const auto &p : ring) - if (p.d2 < r3_sq && (clip_nsigma <= 0.0f || p.v <= thr)) { + // Plain, and - for the correction - with the clip following the curve, as the integrator + // clips where the correction is on. + double kept = 0.0, kept_radial = 0.0; + int n_kept = 0, n_kept_radial = 0; + for (const auto &p : ring) { + if (p.d2 >= r3_sq) + continue; + if (clip_nsigma <= 0.0f || p.v <= thr) { kept += p.v; ++n_kept; } - measured = n_kept > 0; + if (clip_nsigma <= 0.0f || p.v <= thr + std::round(p.b - centre_b)) { + kept_radial += p.v; + ++n_kept_radial; + } + } + measured = n_kept > 0 && n_kept_radial > 0; if (measured) { d_plain[k] = disk_v - kept / n_kept; - d_radial[k] = d_plain[k] - (disk_b - sum_b / n); + d_radial[k] = disk_v - kept_radial / n_kept_radial - (disk_b - sum_b / n); } } if (!measured) @@ -179,3 +195,22 @@ int background_aperture::Best(const std::vector &mean_square) { return -1; return static_cast(std::min_element(mean_square.begin(), mean_square.end()) - mean_square.begin()); } + +double background_aperture::NodeSpacing_px(const UnitCell &cell, char centring, double f_px, double wavelength_A) { + const gemmi::UnitCell uc = cell; + double d_max = 0.0; + for (int h = -2; h <= 2; ++h) + for (int k = -2; k <= 2; ++k) + for (int l = -2; l <= 2; ++l) + if (!systematic_absence(h, k, l, centring)) + d_max = std::max(d_max, uc.calculate_d({h, k, l})); + return f_px * std::tan(2.0 * std::asin(wavelength_A / (2.0 * d_max))); +} + +float background_aperture::DiskRadius(double node_spacing_px, double spot_extent_px, const std::vector &options) { + float r1 = options.front(); + for (const float o : options) + if (o <= node_spacing_px - spot_extent_px) + r1 = o; + return r1; +} diff --git a/image_analysis/bragg_integration/BackgroundAperture.h b/image_analysis/bragg_integration/BackgroundAperture.h index 44e0bb36e..1a69343cf 100644 --- a/image_analysis/bragg_integration/BackgroundAperture.h +++ b/image_analysis/bragg_integration/BackgroundAperture.h @@ -12,6 +12,7 @@ #include "../../common/DiffractionSpot.h" #include "../image_preprocessing/ImagePreprocessorBuffer.h" #include "RingBackground.h" +#include "../../common/UnitCell.h" // How far out the background ring should reach, measured on apertures with nothing in them. // @@ -72,6 +73,16 @@ void MeasureFrame(ImagePreprocessorBuffer &image, int width, int height, // The candidate with the smallest mean square difference, or -1 where nothing was measured. [[nodiscard]] int Best(const std::vector &mean_square); +// How close two lattice nodes can come on the detector [px]: the shortest reciprocal-lattice vector the +// centring allows, 1/d_max, seen from the sample at distance f_px (in pixels). In any one orientation +// neighbours may stand further apart; this is the closest they can be. +[[nodiscard]] double NodeSpacing_px(const UnitCell &cell, char centring, double f_px, double wavelength_A); + +// The signal disk from the lattice and the spot: it must not reach into the next node's spot, which +// extends spot_extent_px from its centre - r1 = s - extent, rounded down to one of `options` (ascending), +// and never below the first of them. +[[nodiscard]] float DiskRadius(double node_spacing_px, double spot_extent_px, const std::vector &options); + constexpr float SPOT_CLEARANCE_PX = 4.0f; // At most this many apertures a frame: a 16 Mpx detector would otherwise lay twenty thousand. constexpr size_t MAX_APERTURES_PER_FRAME = 2000; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index 3695579b8..2b884444f 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -31,6 +31,7 @@ // that the downstream scaling/merge consumes unchanged. // ============================================================================= +#include #include #include #include @@ -289,6 +290,19 @@ public: // a few dozen reflections, against every pixel of the frame here. An empty profile goes back to // accumulating the curve from the reflections. void RadialBackgroundFromRings(const std::vector &ring_bkg, int32_t q_bins, float low_q, float dq); + // Where that curve is in use (the correction on), the ring's high-side clip follows it as well: a ring + // pixel at radial offset `rad` from a reflection at radius r0 is clipped against the threshold plus + // the curve's rise from r0 to it, rounded to whole counts so that every engine keeps the same pixels. + // On a background falling steeply away from the beam the flat threshold otherwise cuts the beam-side + // half of the ring as outliers, and the background reads low. 0 where no curve is in use. + [[nodiscard]] int32_t RadialClipOffset(float r0, float rad) const { + if (!bkg_radial || radial_curve.empty()) + return 0; + const int n = static_cast(radial_curve.size()); + const int b = std::clamp(static_cast(std::lround(r0 + rad)), 0, n - 1); + const int b0 = std::clamp(static_cast(std::lround(r0)), 0, n - 1); + return static_cast(std::lround(radial_curve[b] - radial_curve[b0])); + } [[nodiscard]] const std::vector &GetRadialBackground() const { return radial_curve; } // What this engine has counted since it was built (BraggIntegrationCounts). The GPU engine keeps diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 726d355c4..fd607bec8 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -308,7 +308,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, if (ClipsRing(s)) { const int32_t thr = ClipThreshold(s); for (const auto &[px, rad] : clip_ring) - if (px <= thr) { + if (static_cast(px) <= static_cast(thr) + RadialClipOffset(st.r0, rad)) { clip.sum += px; ++clip.n; if (bkg_radial && radial_curve.empty()) { diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 3a592a138..c378275fd 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -169,7 +169,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p uint8_t *hasobs_o, uint8_t *overloaded_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, int accumulate_rad, - unsigned long long *counts) { + const float *rad_curve, unsigned long long *counts) { const int i = blockIdx.x; if (i >= n) return; @@ -358,7 +358,14 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p if (!BraggPanelHolds(pn, x, y) || mask[y * p.W + x] == MASK_FLUX) continue; const int32_t px = img[y * p.W + x]; if (!valid(px)) continue; - if ((long long) px <= s_thr_i) { + // Where the frame's radial curve is in use the clip follows it (RadialClipOffset). + long long lim = s_thr_i; + if (rad_curve) { + const int b = min(max((int) lroundf(s_r0 + d.rad), 0), n_rad - 1); + const int b0 = min(max((int) lroundf(s_r0), 0), n_rad - 1); + lim += lroundf(rad_curve[b] - rad_curve[b0]); + } + if ((long long) px <= lim) { c_l += px; ++cn_l; if (accumulate_rad) { // The radial curve is binned on the TRUE detector radius, so this is the @@ -1044,6 +1051,9 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu cuda_err(cudaMemsetAsync(d_rad_sum, 0, sizeof(unsigned long long) * n_rad, *stream)); cuda_err(cudaMemsetAsync(d_rad_cnt, 0, sizeof(int) * n_rad, *stream)); } + if (radial_given) + cuda_err(cudaMemcpyAsync(d_rad_curve, radial_curve.data(), sizeof(float) * n_rad, + cudaMemcpyHostToDevice, *stream)); reset<<<32, 256, 0, *stream>>>(d_shell_grid, d_global_grid, d_mom, d_shell_n, d_global_n, d_invd2, GG); cuda_err(cudaGetLastError()); dirty = true; @@ -1056,14 +1066,12 @@ 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_overloaded, d_invd2, d_isum, d_ninner, d_rbin, d_kbin, - d_rad_sum, d_rad_cnt, rad_n, rad_n > 0 && !radial_given, d_counts); + d_rad_sum, d_rad_cnt, rad_n, rad_n > 0 && !radial_given, + radial_given ? d_rad_curve.get() : nullptr, d_counts); cuda_err(cudaGetLastError()); // Correct the flat annulus background for the curvature of the radial background before anything // downstream (profile fit, variance) reads it. - if (radial_given) - cuda_err(cudaMemcpyAsync(d_rad_curve, radial_curve.data(), sizeof(float) * n_rad, - cudaMemcpyHostToDevice, *stream)); if (radial_on) { radial_correct<<<(n + threads - 1) / threads, threads, 0, *stream>>>( d_rad_sum, d_rad_cnt, radial_given ? d_rad_curve.get() : nullptr, n_rad, d_k_diff, k_len, k_off, diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineMetal.metal b/image_analysis/bragg_integration/BraggIntegrationEngineMetal.metal index 46d12341f..21750d674 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineMetal.metal +++ b/image_analysis/bragg_integration/BraggIntegrationEngineMetal.metal @@ -167,6 +167,7 @@ kernel void clip_sums(device const float2 *xy [[buffer(0)]], device atomic_uint *rad_sum [[buffer(8)]], device atomic_int *rad_cnt [[buffer(9)]], device const BraggPanel *panels [[buffer(10)]], + device const float *curve [[buffer(11)]], uint i [[thread_position_in_grid]]) { if (i >= f.n || !clips[i]) return; @@ -179,13 +180,20 @@ kernel void clip_sums(device const float2 *xy [[buffer(0)]], const int y0 = max(0, (int) floor(py - st.ey_out - 1.0f)); const int y1 = min(f.height - 1, (int) ceil(py + st.ey_out + 1.0f)); const int32_t t = thr[i]; + const int b0 = f.n_curve > 0 ? clamp((int) round(st.r0), 0, f.n_curve - 1) : 0; BraggClipSums c = {}; for (int y = y0; y <= y1; ++y) for (int x = x0; x <= x1; ++x) { const BraggStencilDist d = BraggStencilDistances(st, x - px, y - py); const int idx = y * f.width + x; const int32_t v = BraggPanelHolds(pn, x, y) ? img[idx] : PX_MASKED; - if (!(d.inner >= f.r2_sq && d.outer < f.r3_sq && valid_px(v) && mask[idx] != MASK_FLUX && v <= t)) + // Where the frame's radial curve is in use the clip follows it (RadialClipOffset). + int64_t lim = t; + if (f.n_curve > 0) { + const int b = clamp((int) round(st.r0 + d.rad), 0, f.n_curve - 1); + lim += (int64_t) round(curve[b] - curve[b0]); + } + if (!(d.inner >= f.r2_sq && d.outer < f.r3_sq && valid_px(v) && mask[idx] != MASK_FLUX && (int64_t) v <= lim)) continue; c.sum += v; ++c.n; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineMetal.mm b/image_analysis/bragg_integration/BraggIntegrationEngineMetal.mm index 2d5e7c4d2..e3c5d0ab4 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineMetal.mm +++ b/image_analysis/bragg_integration/BraggIntegrationEngineMetal.mm @@ -107,6 +107,7 @@ struct BraggIntegrationEngineMetal::Device { id unmark = nil; bool clear = false; id rad_sum = nil, rad_cnt = nil; // the radial background curve, as the clip kernel adds it up + id curve = nil; // the curve handed over, which the clip follows // Per reflection, grown with the frames. size_t capacity = 0; id xy = nil, panels = nil, level = nil, sums = nil, clips = nil, thr = nil, clip_out = nil; @@ -223,6 +224,13 @@ std::vector BraggIntegrationEngineMetal::Run(const ImagePreprocessor // The curve is added up here only where no ring profile was handed over (RadialBackgroundFromRings). const bool accumulate_radial = bkg_radial && radial_curve.empty(); f.n_rad = accumulate_radial ? a.n_rad : 0; + // The clip follows the handed-over curve where the correction uses it (RadialClipOffset). + f.n_curve = bkg_radial && !radial_curve.empty() ? static_cast(radial_curve.size()) : 0; + if (f.n_curve > 0) { + if (dev->curve == nil || dev->curve.length < radial_curve.size() * sizeof(float)) + dev->curve = buffer(radial_curve.size() * sizeof(float)); + std::copy(radial_curve.begin(), radial_curve.end(), static_cast(dev->curve.contents)); + } f.static_mask_words = static_cast(static_mask.size()); f.n = static_cast(npredicted); @@ -311,6 +319,7 @@ std::vector BraggIntegrationEngineMetal::Run(const ImagePreprocessor [enc setBuffer:(accumulate_radial ? dev->rad_sum : dev->clips) offset:0 atIndex:8]; // not written without the curve [enc setBuffer:(accumulate_radial ? dev->rad_cnt : dev->clips) offset:0 atIndex:9]; [enc setBuffer:dev->panels offset:0 atIndex:10]; + [enc setBuffer:(f.n_curve > 0 ? dev->curve : dev->thr) offset:0 atIndex:11]; // not read without a curve [enc dispatchThreads:grid threadsPerThreadgroup:group]; [enc endEncoding]; Wait(cb); diff --git a/image_analysis/bragg_integration/BraggPassA.h b/image_analysis/bragg_integration/BraggPassA.h index 8b1b10979..31df74e8e 100644 --- a/image_analysis/bragg_integration/BraggPassA.h +++ b/image_analysis/bragg_integration/BraggPassA.h @@ -40,6 +40,7 @@ struct BraggPassAFrame { int32_t use_owner; // an overlap treatment is on, so the owner map is built and read int32_t exclude; // ... and a pixel a neighbour owns leaves the signal disk int32_t n_rad; // radial background bins; 0 when the curve is not accumulated + int32_t n_curve; // bins of the handed-over curve the clip follows (RadialClipOffset); 0: none uint32_t static_mask_words; // length of the run's packed pixel mask uint32_t n; // reflections }; diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 2f02f1edf..1705dbec2 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -1334,7 +1334,7 @@ class Rugnux { std::vector footprint; std::vector spot_q; size_t width_images = 0; - background_aperture::Tally apertures; + std::vector apertures; // one per candidate signal disk }; struct Defects { size_t pixels = 0; diff --git a/rugnux/RugnuxPreScan.cpp b/rugnux/RugnuxPreScan.cpp index 9603581df..4796f81bd 100644 --- a/rugnux/RugnuxPreScan.cpp +++ b/rugnux/RugnuxPreScan.cpp @@ -369,20 +369,26 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru }; // Stills: how far the background ring should reach, measured on empty apertures of the same frames - // (BackgroundAperture.h) at the signal disk and gap the run integrates with. The candidates run up - // to the ring the run starts from, so where the noise wins that ring is what is kept. The clip gain - // is read once the gain below has been measured. + // (BackgroundAperture.h), for each signal disk the run may integrate with - from 4 px up to the one it + // starts from, the gap 2 px beyond it. Which disk is decided once the spots are measured; the ring + // candidates run up to the ring the run starts from, so where the noise wins that ring is what is + // kept. The clip gain is read once the gain below has been measured. const bool stills = !experiment_.IsRotationIndexing(); const BraggIntegrationSettings aperture_bis = prescan_x.GetBraggIntegrationSettings(); - const std::vector aperture_r3 = stills - ? background_aperture::CandidateR3(aperture_bis.GetR2(), aperture_bis.GetR3()) : std::vector(); + std::vector aperture_r1; + if (stills) + for (float r1 = std::min(4.0f, aperture_bis.GetR1()); r1 <= aperture_bis.GetR1(); r1 += 1.0f) + aperture_r1.push_back(r1); + const auto aperture_r3 = [&](size_t o) { + return background_aperture::CandidateR3(aperture_r1[o] + 2.0f, aperture_bis.GetR3()); + }; float aperture_gain = 1.0f; // The spots of one image, and what is taken off them. const auto find_spots = [&](PreScanWorker &w, const CompressedImage &image, int image_idx, bool for_width, std::vector &curves, std::vector &spot_q, std::vector &footprint, - background_aperture::Tally &apertures) { + std::vector &apertures) { try { // As in the image loops: a frame the device cannot decode goes to the host decoder. ImageStatistics stats; @@ -426,7 +432,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru MeasureSpotsOnTiles(*w.preprocessed, static_cast(prescan_x.GetXPixelsNumConv()), static_cast(prescan_x.GetYPixelsNumConv()), prescan_x.GetDiffractionGeometry(), prescan_x.GetModules(), spots, curves, footprint); - if (for_width && !aperture_r3.empty()) { + if (for_width && !aperture_r1.empty()) { // The frame's ring background, as the integration of this frame will read it - from the // adaptive finder, which has it whole-frame and peak-excluded. const std::vector no_rings; @@ -435,11 +441,14 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru prescan_mapping->GetQBinCount(), prescan_mapping->Settings().GetLowQ_recipA(), prescan_mapping->Settings().GetQSpacing_recipA()); - background_aperture::MeasureFrame(*w.preprocessed, static_cast(prescan_x.GetXPixelsNumConv()), - static_cast(prescan_x.GetYPixelsNumConv()), spots, - aperture_bis.GetR1(), aperture_bis.GetR2(), aperture_r3, - aperture_bis.GetBackgroundClipNSigma(), aperture_gain, rings, - prescan_x.GetDiffractionGeometry(), prescan_x.GetModules(), apertures); + apertures.resize(aperture_r1.size()); + for (size_t o = 0; o < aperture_r1.size(); o++) + background_aperture::MeasureFrame(*w.preprocessed, static_cast(prescan_x.GetXPixelsNumConv()), + static_cast(prescan_x.GetYPixelsNumConv()), spots, + aperture_r1[o], aperture_r1[o] + 2.0f, aperture_r3(o), + aperture_bis.GetBackgroundClipNSigma(), aperture_gain, rings, + prescan_x.GetDiffractionGeometry(), prescan_x.GetModules(), + apertures[o]); } }; @@ -450,8 +459,8 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru std::vector shape_curves; // Every spot's measured widths, for the spot footprint (SpotFootprint.h). std::vector footprint_spots; - // The empty apertures of the frames the width is measured on (stills). - background_aperture::Tally aperture_tally; + // The empty apertures of the frames the width is measured on, one tally per signal disk (stills). + std::vector aperture_tally(aperture_r1.size()); // Images the width ended up being measured on, for the log: the tiers below stop as soon as the // answer has settled, so this is a property of the crystal and worth reporting. size_t width_images = 0; @@ -625,7 +634,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru std::vector> curves_of(ordinals.size()); std::vector> spot_q_of(ordinals.size()); std::vector> footprint_of(ordinals.size()); - std::vector apertures_of(ordinals.size()); + std::vector> apertures_of(ordinals.size()); // Frames read for the bandwidth alone, past the tier the width settled on. Their spots stay // out of the powder measurement, which is taken on the frames it always was. std::vector shape_only(ordinals.size(), 0); @@ -743,7 +752,8 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // In sample order, over the frames the width was taken on. for (size_t i = 0; i < ordinals.size(); i++) if (position[i] < width_images) - aperture_tally.Add(apertures_of[i]); + for (size_t o = 0; o < apertures_of[i].size(); o++) + aperture_tally[o].Add(apertures_of[i][o]); // In sample order, so the table does not depend on how the workers interleaved. for (const auto &f : footprint_of) footprint_spots.insert(footprint_spots.end(), f.begin(), f.end()); @@ -783,21 +793,39 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // keeps the area it has at the shipped 4/6/13 - a ring that shrinks with the disk is what makes a // bare --integration-radius worse than the default it replaces. // - // Stills keep their signal disk and take the ring's outer radius from the empty apertures - // instead (BackgroundAperture.h): the candidate whose ring predicted the background under the - // disk best. Where the noise wins that is the widest, which is the ring the run started from. - // The same apertures decide the radial curvature correction, where it was left to auto: on - // where the corrected ring predicts the disk better. + // Stills take the signal disk from the lattice spacing and the spot width instead, and the ring's + // outer radius from the empty apertures (BackgroundAperture.h): the candidate whose ring + // predicted the background under the disk best. Where the noise wins that is the widest, which + // is the ring the run started from. The same apertures decide the radial curvature correction, + // where it was left to auto: on where the corrected ring predicts the disk better. if (want_width && stills) { spot_width_measured_ = true; BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings(); + // The signal disk: as large as the run starts with, but not reaching into the next lattice + // node's spot (DiskRadius). On a pixel detector the nodes stand 20 px apart and the disk is + // left alone; on a binned CCD with a long axis they are 5-7 px apart and it shrinks. Where the + // cell or the spot width is not known it stays as it was. + size_t o = aperture_r1.size() - 1; + const auto r80 = spot_width::R80AtReference(width_curves); + const auto cell = experiment_.GetUnitCell(); + double spacing = 0.0; + if (r80 && cell) { + const auto &g = experiment_.GetDiffractionGeometry(); + spacing = background_aperture::NodeSpacing_px(*cell, experiment_.GetCentering(), + g.GetDetectorDistance_mm() / g.GetPixelSize_mm(), + g.GetWavelength_A()); + const float r1 = background_aperture::DiskRadius(spacing, 2.0 * *r80, aperture_r1); + o = std::find(aperture_r1.begin(), aperture_r1.end(), r1) - aperture_r1.begin(); + } + const auto &tally = aperture_tally[o]; + const std::vector r3 = aperture_r3(o); const auto radial_setting = bis.GetBackgroundRadialCorrection(); // The ring is chosen on the plain differences, so its reach does not rest on the radial // model; the correction is then kept, where it was left to auto, if it makes that ring predict // the disk better. Read only where every frame gave a ring background. - const std::vector plain = aperture_tally.MeanSquare(false); - const std::vector radial = aperture_tally.apertures_radial == aperture_tally.apertures - ? aperture_tally.MeanSquare(true) : std::vector(); + const std::vector plain = tally.MeanSquare(false); + const std::vector radial = tally.apertures_radial == tally.apertures + ? tally.MeanSquare(true) : std::vector(); const int best = background_aperture::Best(plain); const bool use_radial = radial_setting.has_value() ? *radial_setting @@ -809,22 +837,24 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru const auto table = [&](const std::vector &m) { std::string t; for (size_t k = 0; k < m.size(); k++) - t += fmt::format(" {:.0f}:{:.4g}", aperture_r3[k], m[k]); + t += fmt::format(" {:.0f}:{:.4g}", r3[k], m[k]); return t; }; - bis.R3(aperture_r3[best]); + bis.R1(aperture_r1[o]).R2(aperture_r1[o] + 2.0f).R3(r3[best]); if (!radial_setting) bis.BackgroundRadialCorrection(use_radial ? std::nullopt : std::optional(false)); experiment_.ImportBraggIntegrationSettings(bis); + logger.Info("Signal disk: lattice nodes at least {:.1f} px apart on the detector, spots reaching " + "{:.1f} px => r1={:.1f}", spacing, r80 ? 2.0f * *r80 : 0.0f, bis.GetR1()); logger.Info("Background ring: on {} empty apertures from {} images the r{:.0f}..r3 ring predicts " "the r{:.0f} disk with a mean square error of (r3:value){}, and with the radial " "curvature correction{} => integration radii r1={:.1f} r2={:.1f} r3={:.1f}, radial " - "correction {} (z = {:.1f} at that ring)", aperture_tally.apertures, width_images, bis.GetR2(), bis.GetR1(), - table(plain), radial.empty() ? std::string(" (no ring background)") : table(radial), - bis.GetR1(), bis.GetR2(), bis.GetR3(), use_radial ? "on" : "off", - aperture_tally.RadialZ(best)); - Note(fmt::format("background ring {:.0f}..{:.0f} px, radial correction {}, measured on empty " - "apertures", bis.GetR2(), bis.GetR3(), use_radial ? "on" : "off")); + "correction {} (z = {:.1f} at that ring)", tally.apertures, width_images, bis.GetR2(), + bis.GetR1(), table(plain), radial.empty() ? std::string(" (no ring background)") : table(radial), + bis.GetR1(), bis.GetR2(), bis.GetR3(), use_radial ? "on" : "off", tally.RadialZ(best)); + Note(fmt::format("integration radii {:.0f}/{:.0f}/{:.0f} px, radial correction {}, measured on " + "the lattice and empty apertures", bis.GetR1(), bis.GetR2(), bis.GetR3(), + use_radial ? "on" : "off")); } } if (want_width && !stills) { diff --git a/tests/BackgroundApertureTest.cpp b/tests/BackgroundApertureTest.cpp index bedf72ee5..3487919d5 100644 --- a/tests/BackgroundApertureTest.cpp +++ b/tests/BackgroundApertureTest.cpp @@ -2,6 +2,7 @@ // SPDX-License-Identifier: GPL-3.0-only #include +#include #include #include @@ -116,3 +117,20 @@ TEST_CASE("BackgroundAperture_RadialCorrectionHelpsOnARadialBackground", "[Integ // With the curvature taken out, the noise wins again and a wide ring is best. CHECK(r3[background_aperture::Best(radial)] >= 12.0f); } + +// A 265 A axis at 1092 px and 1.3 A puts its nodes 5.4 px apart; with I-centring the (100) node of a +// 100 A cubic cell is absent and the closest allowed nodes are (110), 1/70.7 A apart. +TEST_CASE("BackgroundAperture_NodeSpacing", "[Integration]") { + const UnitCell hexagonal{78.6f, 78.6f, 265.3f, 90.0f, 90.0f, 120.0f}; + CHECK(background_aperture::NodeSpacing_px(hexagonal, 'P', 1092.0, 1.3) == Catch::Approx(5.35).margin(0.05)); + const UnitCell cubic{100.0f, 100.0f, 100.0f, 90.0f, 90.0f, 90.0f}; + CHECK(background_aperture::NodeSpacing_px(cubic, 'P', 1000.0, 1.0) == Catch::Approx(10.0).margin(0.01)); + CHECK(background_aperture::NodeSpacing_px(cubic, 'I', 1000.0, 1.0) == Catch::Approx(14.14).margin(0.02)); +} + +TEST_CASE("BackgroundAperture_DiskRadius", "[Integration]") { + const std::vector options{4.0f, 5.0f, 6.0f}; + CHECK(background_aperture::DiskRadius(5.4, 3.5, options) == 4.0f); // dense: never below the first + CHECK(background_aperture::DiskRadius(8.6, 3.5, options) == 5.0f); + CHECK(background_aperture::DiskRadius(20.0, 3.2, options) == 6.0f); // sparse: the disk the run starts with +} diff --git a/tests/BraggIntegrationEngineCPUTest.cpp b/tests/BraggIntegrationEngineCPUTest.cpp index 292ff9ee4..86e26f3d4 100644 --- a/tests/BraggIntegrationEngineCPUTest.cpp +++ b/tests/BraggIntegrationEngineCPUTest.cpp @@ -203,3 +203,46 @@ TEST_CASE("BraggIntegrationEngineCPU_RadialCurveFromRings", "[Integration][porta engine.RadialBackgroundFromRings({}, 500, low_q, dq); CHECK(engine.GetRadialBackground().empty()); } + +// A ring on a background falling steeply away from the beam: the flat high-side clip cuts the beam-side +// half of the ring and the background reads low. With the correction on, the clip follows the frame's +// radial curve (RadialClipOffset) and the empty aperture reads zero. +TEST_CASE("BraggIntegrationEngineCPU_ClipFollowsRadialCurve", "[Integration][portable]") { + DiffractionExperiment experiment(DetJF(2)); + experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(400.0f).BeamY_pxl(400.0f); + const float bx = 400.0f, by = 400.0f; + const auto background = [](double r) { return 40000.0 - 100.0 * r; }; + const size_t width = experiment.GetXPixelsNum(), npixel = experiment.GetPixelsNum(); + ImagePreprocessorBuffer image(npixel); + for (size_t i = 0; i < npixel; ++i) + image[i] = static_cast(std::lround(std::max(0.0, + background(std::hypot(static_cast(i % width) - bx, static_cast(i / width) - by))))); + const float low_q = 0.05f, dq = 0.01f; + const double F_px = 100.0 / 0.075, lambda = 1.0; + std::vector rings(500); + for (size_t i = 0; i < rings.size(); ++i) { + const double q = low_q + (i + 0.5) * dq; + rings[i] = static_cast(background(F_px * std::tan(2.0 * std::asin(q * lambda / (4.0 * PI))))); + } + Reflection r{}; + r.h = 1; r.k = 0; r.l = 0; + r.predicted_x = bx + 200.0f; + r.predicted_y = by; + r.d = 2.0f; + r.prescaling_corr = 1.0f; + r.partiality = 1.0f; + const std::vector predicted{r}; + auto background_read = [&](bool correct) { + experiment.ImportBraggIntegrationSettings(BraggIntegrationSettings().BackgroundRadialCorrection(correct)); + BraggIntegrationEngineCPU engine(experiment, PixelMask(experiment)); + engine.RadialBackgroundFromRings(rings, static_cast(rings.size()), low_q, dq); + const auto out = engine.Run(image, predicted, predicted.size(), 0); + REQUIRE(out.size() == 1); + return out[0].bkg; + }; + // The kernel models the radial offset only; the ring's tangential reach adds the circle's own + // curvature, slope * / (2 r) = 100 * 51 / 400 = 13 counts, which it leaves. + const double truth = background(200.0); + CHECK(background_read(false) < truth - 100.0); + CHECK(std::fabs(background_read(true) - truth) < 20.0); +} diff --git a/tests/BraggIntegrationEngineGPUTest.cpp b/tests/BraggIntegrationEngineGPUTest.cpp index 4f8788703..4ec667ec2 100644 --- a/tests/BraggIntegrationEngineGPUTest.cpp +++ b/tests/BraggIntegrationEngineGPUTest.cpp @@ -424,7 +424,9 @@ TEST_CASE("BraggIntegrationEngineGPU_RadialCurveFromRingsMatchesCPU", "[Integrat DiffractionExperiment experiment(DetJF(2)); experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(400.0f).BeamY_pxl(400.0f); const float bx = 400.0f, by = 400.0f; - const auto background = [](double r) { return 50.0 + 0.05 * (r - 300.0) * (r - 300.0); }; + // Curved, plus a slope steep enough that the ring's flat clip would cut its beam-side half: the + // GPU must follow the curve in the clip as the CPU does (RadialClipOffset). + const auto background = [](double r) { return 40000.0 - 100.0 * r + 0.05 * (r - 300.0) * (r - 300.0); }; const size_t width = experiment.GetXPixelsNum(), npixel = experiment.GetPixelsNum(); ImagePreprocessorBuffer cpu_image(npixel); ImagePreprocessorBufferGPU gpu_image(npixel); @@ -462,8 +464,7 @@ TEST_CASE("BraggIntegrationEngineGPU_RadialCurveFromRingsMatchesCPU", "[Integrat for (size_t i = 0; i < out_cpu.size(); ++i) { INFO("mode " << static_cast(mode) << " reflection " << i); CHECK(out_gpu[i].bkg == Catch::Approx(out_cpu[i].bkg).epsilon(0.001).margin(0.05)); - CHECK(std::fabs(out_gpu[i].I) < 10.0f); - CHECK(std::fabs(out_cpu[i].I) < 10.0f); + CHECK(out_gpu[i].I == Catch::Approx(out_cpu[i].I).epsilon(0.01).margin(2.0)); } } }