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)); } } }