diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index ca860600b..6b4c101f7 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -235,22 +235,37 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, // soon as it has detected - so this must be called after detection and before integration. // Where no such background exists (no adaptive finder), auto leaves the correction off: the // plain profile carries the Bragg peaks and cannot support an absolute threshold. + // + // On stills the correction also takes its CURVE from that background, and auto applies it to + // every image that has one. A serial frame has a few dozen predicted reflections, so their own + // ring pixels - where the curve comes from otherwise - sample the radial background sparsely, + // while the finder's rings hold every pixel of the frame. And the correction is the measured + // curvature itself: where the background has none it is zero, so there is nothing for a gate + // to decide. The ice score is the wrong gate there: a smooth hump that is not ice - the + // solvent or a carrier medium - curves the background as much as an ice ring does, and the + // score does not see it. const auto decide_radial_background = [this, &spot_finding_settings, &output]() { - if (!bragg_engine->IsBackgroundRadialAuto()) - return; // Same condition as the score at the end of this function: the adaptive finder holds its // ring background from whenever it last ran, so requiring that it ran for THIS image is // what keeps a stale curve out. - if (!spot_finding_settings.adaptive_threshold) - return; - const std::vector &ring_bkg = adaptiveSpotFinder->GetRingBackground(); - if (ring_bkg.empty()) + const std::vector no_rings; + const std::vector &ring_bkg = spot_finding_settings.adaptive_threshold + ? adaptiveSpotFinder->GetRingBackground() : no_rings; + const bool stills = !experiment.IsRotationIndexing(); + if (stills) + bragg_engine->RadialBackgroundFromRings(ring_bkg, integration.GetQBinCount(), + integration.Settings().GetLowQ_recipA(), + integration.Settings().GetQSpacing_recipA()); + if (!bragg_engine->IsBackgroundRadialAuto() || ring_bkg.empty()) return; output.ice_ring_score = AzimuthalIntegrationProfile::IceRingScore( ring_bkg, integration.GetQBinCount(), integration.Settings(), spot_finding_settings.ice_ring_width_Q_recipA); - bragg_engine->BackgroundRadial(*output.ice_ring_score - >= experiment.GetScalingSettings().GetIceMinScore()); + if (stills) + bragg_engine->BackgroundRadial(!bragg_engine->GetRadialBackground().empty()); + else + bragg_engine->BackgroundRadial(*output.ice_ring_score + >= experiment.GetScalingSettings().GetIceMinScore()); }; // A missing min-pix (std::nullopt) means "choose it per image". This applies only to the stills diff --git a/image_analysis/bragg_integration/BackgroundAperture.cpp b/image_analysis/bragg_integration/BackgroundAperture.cpp new file mode 100644 index 000000000..a5b696756 --- /dev/null +++ b/image_analysis/bragg_integration/BackgroundAperture.cpp @@ -0,0 +1,122 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "BackgroundAperture.h" + +#include +#include + +std::vector background_aperture::CandidateR3(float r2, float max_r3) { + std::vector out; + for (float r3 = r2 + 2.0f; r3 <= max_r3 + 1e-3f; r3 += 1.0f) + out.push_back(r3); + return out; +} + +void background_aperture::MeasureFrame(ImagePreprocessorBuffer &image, int width, int height, + const std::vector &spots, float r1, float r2, + const std::vector &r3, float clip_nsigma, float gain, + Tally &tally) { + if (r3.empty()) + return; + tally.sum_d2.resize(r3.size(), 0.0); + const float r_max = r3.back(); + const int half = static_cast(std::ceil(r_max)); + const int side = 2 * half + 1; + + // A grid of centres, at least one aperture apart, kept clear of every found spot. + const double area = static_cast(width) * height; + const int step = std::max(side, static_cast(std::ceil(std::sqrt(area / MAX_APERTURES_PER_FRAME)))); + const float clear = r_max + SPOT_CLEARANCE_PX; + std::vector cx, cy; + for (int y = half; y + half < height; y += step) + for (int x = half; x + half < width; x += step) { + bool near_spot = false; + for (const auto &s : spots) { + const Coord c = s.RawCoord(); + if ((c.x - x) * (c.x - x) + (c.y - y) * (c.y - y) < clear * clear) { + near_spot = true; + break; + } + } + if (!near_spot) { + cx.push_back(x); + cy.push_back(y); + } + } + if (cx.empty()) + return; + const int32_t *tiles = image.GatherTiles(cx, cy, width, half); + + std::vector ring_sum(r3.size()); + std::vector ring_n(r3.size()); + std::vector ring_px; + std::vector ring_d2; + for (size_t a = 0; a < cx.size(); ++a) { + const int32_t *t = tiles + a * static_cast(side) * side; + double disk_sum = 0.0; + int disk_n = 0; + bool readable = true; + ring_px.clear(); + ring_d2.clear(); + for (int dy = -half; dy <= half && readable; ++dy) + for (int dx = -half; dx <= half; ++dx) { + const float d2 = static_cast(dx * dx + dy * dy); + if (d2 >= r_max * r_max) + continue; + const int32_t v = t[(dy + half) * side + dx + half]; + if (v == INT32_MIN || v == INT32_MAX) { + readable = false; + break; + } + if (d2 < r1 * r1) { + disk_sum += v; + ++disk_n; + } else if (d2 >= r2 * r2) { + ring_px.push_back(v); + ring_d2.push_back(d2); + } + } + if (!readable || disk_n == 0) + continue; + // Every candidate ring as the integrator reads it: its mean, then one high-side clip pass. + bool measured = true; + for (size_t k = 0; k < r3.size() && measured; ++k) { + double sum = 0.0; + int n = 0; + for (size_t p = 0; p < ring_px.size(); ++p) + if (ring_d2[p] < r3[k] * r3[k]) { + sum += ring_px[p]; + ++n; + } + if (n == 0) { + measured = false; + break; + } + const double mean = sum / n; + const double thr = mean + clip_nsigma * std::sqrt(gain * std::max(mean, 1.0)); + ring_sum[k] = 0.0; + ring_n[k] = 0; + for (size_t p = 0; p < ring_px.size(); ++p) + if (ring_d2[p] < r3[k] * r3[k] && (clip_nsigma <= 0.0f || ring_px[p] <= thr)) { + ring_sum[k] += ring_px[p]; + ++ring_n[k]; + } + measured = ring_n[k] > 0; + } + if (!measured) + continue; + const double disk = disk_sum / disk_n; + for (size_t k = 0; k < r3.size(); ++k) { + const double d = disk - ring_sum[k] / ring_n[k]; + tally.sum_d2[k] += d * d; + } + ++tally.apertures; + } +} + +int background_aperture::Best(const Tally &tally) { + if (tally.apertures == 0 || tally.sum_d2.empty()) + return -1; + return static_cast(std::min_element(tally.sum_d2.begin(), tally.sum_d2.end()) - tally.sum_d2.begin()); +} diff --git a/image_analysis/bragg_integration/BackgroundAperture.h b/image_analysis/bragg_integration/BackgroundAperture.h new file mode 100644 index 000000000..7227301d6 --- /dev/null +++ b/image_analysis/bragg_integration/BackgroundAperture.h @@ -0,0 +1,55 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include + +#include "../../common/DiffractionSpot.h" +#include "../image_preprocessing/ImagePreprocessorBuffer.h" + +// How far out the background ring should reach, measured on apertures with nothing in them. +// +// The ring mean stands in for the background under the signal disk. A wider ring averages more pixels, +// so its own noise falls; but it also reaches further from the reflection, so whatever the background +// does between the disk and the ring - a curved solvent ring, the tile structure of a fibre-optic CCD, +// the tails of reflections nobody predicted - enters the intensity as a bias. Which of the two wins is +// a property of the detector and the sample, not a constant: on a pixel detector counting a fraction of +// a photon per pixel the noise wins and the widest ring is best, on a binned CCD reading hundreds of +// ADU a pixel the structure wins and the narrowest is. +// +// So it is measured. Apertures are laid on a grid over pre-scan frames, away from every spot found; on +// each, the mean of the r1 disk and the mean of the r2..r3 ring are taken for every candidate r3, and +// the run keeps the r3 whose ring predicts the disk best - the smallest mean square of the difference. +// With nothing in the aperture the difference is exactly the error the ring makes in the background of +// a reflection: its noise plus its bias. The disk's own noise is the same for every candidate, so it +// does not move the choice. The ring is clipped as the integrator clips it (high side, at +// mean + nsigma * sqrt(gain * mean)). +namespace background_aperture { + +// The squared disk-minus-ring differences, summed per candidate r3, over every aperture measured. +struct Tally { + std::vector sum_d2; + size_t apertures = 0; +}; + +// The candidate outer radii: from r2 + 2 up to max_r3, one pixel apart. +[[nodiscard]] std::vector CandidateR3(float r2, float max_r3); + +// The apertures of one frame, appended to `tally`. Centres lie on a grid over the frame, and an aperture +// is used only where every pixel out to the largest candidate is readable and no found spot is within +// that radius plus SPOT_CLEARANCE_PX of its centre. +void MeasureFrame(ImagePreprocessorBuffer &image, int width, int height, + const std::vector &spots, float r1, float r2, + const std::vector &r3, float clip_nsigma, float gain, Tally &tally); + +// The candidate with the smallest mean square difference, or -1 where no aperture was measured. +[[nodiscard]] int Best(const Tally &tally); + +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; + +} // namespace background_aperture diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp index bfcb4401f..fec16cea1 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp @@ -218,6 +218,50 @@ void BraggIntegrationEngine::BuildRadialKernel(float grow) { k_diff.push_back(static_cast(hist_ann[k] / sa - hist_disk[k] / sd)); } +void BraggIntegrationEngine::RadialBackgroundFromRings(const std::vector &ring_bkg, int32_t q_bins, + float low_q, float dq) { + radial_curve.clear(); + if (q_bins <= 0 || !(dq > 0.0f)) + return; + // The rings averaged over their azimuthal bins, where the profile has any. + std::vector sum(q_bins, 0.0); + std::vector cnt(q_bins, 0); + for (size_t i = 0; i < ring_bkg.size(); ++i) + if (std::isfinite(ring_bkg[i])) { + sum[i % q_bins] += ring_bkg[i]; + ++cnt[i % q_bins]; + } + std::vector q, b; + for (int i = 0; i < q_bins; ++i) + if (cnt[i] > 0) { + q.push_back(low_q + (i + 0.5) * dq); + b.push_back(sum[i] / cnt[i]); + } + if (q.size() < 2) + return; + + // The same bins as the engines' own curve: one per pixel of distance from the beam, out to the far + // corner. Each is read off the rings by linear interpolation in q; beyond the first and the last + // ring the end value is held, which has no curvature to correct. + const int n = static_cast(std::ceil(r_max)) + 2; + const double two_k = 4.0 * PI / geom.GetWavelength_A(); + radial_curve.resize(n); + size_t j = 0; + for (int r = 0; r < n; ++r) { + const double qr = two_k * std::sin(0.5 * std::atan(r / F_px)); + while (j + 1 < q.size() && q[j + 1] < qr) + ++j; + double v; + if (qr <= q.front()) + v = b.front(); + else if (j + 1 >= q.size()) + v = b.back(); + else + v = b[j] + (b[j + 1] - b[j]) * (qr - q[j]) / (q[j + 1] - q[j]); + radial_curve[r] = static_cast(v); + } +} + std::vector BraggIntegrationEngine::Panels(const std::vector &predicted, size_t npredicted) const { std::vector out(npredicted); diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index bec850887..3695579b8 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -229,6 +229,10 @@ protected: int k_len = 0; // entries per kernel int n_kern = 1; // kernels in the table (1 = circular stencil) std::vector k_diff; // n_kern * k_len, annulus-minus-disk weight per radial offset + // The radial background B the kernels are applied to, one value per pixel of distance from the + // beam, where the caller handed one over (RadialBackgroundFromRings). Empty: the curve is + // accumulated from the clipped ring pixels of the predicted reflections themselves. + std::vector radial_curve; // One radial-offset kernel for a ring grown by `grow` px, appended to k_diff. Kept out of line // so the circular and elongated cases cannot drift apart. The signal disk does not change with @@ -277,6 +281,16 @@ public: void BackgroundRadial(bool on) { bkg_radial = on && bkg_radial_built; } [[nodiscard]] bool IsBackgroundRadialAuto() const { return bkg_radial_auto; } + // The radial background curve for the images that follow, from a peak-excluded background per + // resolution ring of the whole detector (the adaptive spot finder's GetRingBackground): q bin i + // holds the ring at low_q + (i + 0.5) * dq, entries past the first q_bins are further azimuthal + // bins and are averaged in, and a ring with no value (NAN) is bridged by its neighbours. On a + // serial frame the predicted reflections' own ring pixels are a sparse sample of the same curve - + // 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); + [[nodiscard]] const std::vector &GetRadialBackground() const { return radial_curve; } + // 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. diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index e99c37f81..726d355c4 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -311,7 +311,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, if (px <= thr) { clip.sum += px; ++clip.n; - if (bkg_radial) { + if (bkg_radial && radial_curve.empty()) { // The radial curve is binned on the TRUE detector radius, so the // offset here is the unshrunk radial projection. const int b = std::clamp(static_cast(std::lround(st.r0 + rad)), 0, a.n_rad - 1); @@ -335,7 +335,7 @@ BraggIntegrationEngineCPU::PassA BraggIntegrationEngineCPU::StartPassA(size_t np const float by0 = static_cast(plane_bounds.y0), by1 = static_cast(plane_bounds.y1); a.n_rad = static_cast(std::ceil(std::hypot(std::max(beam_x - bx0, bx1 - beam_x), std::max(beam_y - by0, by1 - beam_y)))) + 2; - if (bkg_radial) { + if (bkg_radial && radial_curve.empty()) { a.rad_sum.assign(a.n_rad, 0); a.rad_cnt.assign(a.n_rad, 0); } @@ -612,7 +612,11 @@ std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler & double corr = 0.0; for (int k = 0; k < k_len; ++k) { const int b = std::clamp(rh.r_bin + k - k_off, 0, a.n_rad - 1); - const double v = a.rad_cnt[b] > 0 ? static_cast(a.rad_sum[b]) / a.rad_cnt[b] : rh.bkg; + double v; + if (!radial_curve.empty()) + v = radial_curve[b]; + else + v = a.rad_cnt[b] > 0 ? static_cast(a.rad_sum[b]) / a.rad_cnt[b] : rh.bkg; corr += static_cast(kern[k]) * v; } rh.bkg -= corr; // annulus mean -> mean over the signal disk diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 2c1428552..3a592a138 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -168,7 +168,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p float *bkgvar_o, float *varbkg_o, float *obsx_o, float *obsy_o, uint8_t *ok_o, uint8_t *strong_o, 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, + unsigned long long *rad_sum, int *rad_cnt, int n_rad, int accumulate_rad, unsigned long long *counts) { const int i = blockIdx.x; if (i >= n) return; @@ -360,7 +360,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p if (!valid(px)) continue; if ((long long) px <= s_thr_i) { c_l += px; ++cn_l; - if (n_rad > 0) { + if (accumulate_rad) { // The radial curve is binned on the TRUE detector radius, so this is the // unshrunk radial projection - no per-pixel sqrt. const int idx = (int) lroundf(s_r0 + d.rad) - s_radbase; @@ -375,7 +375,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p } __syncthreads(); - if (n_rad > 0) { + if (accumulate_rad) { for (int t = threadIdx.x; t < p.rad_w; t += blockDim.x) { if (s_radn[t] == 0) continue; const int b = min(max(s_radbase + t, 0), n_rad - 1); @@ -602,7 +602,7 @@ __global__ void build_profiles(const unsigned long long *shell_grid, const unsig // correction because the kernel weights sum to zero. The kernel depends on how far this // reflection's ring was elongated, so k_diff is a table and kbin_a picks the row. // Mirrors BraggIntegrationEngineCPU. --- -__global__ void radial_correct(const unsigned long long *rad_sum, const int *rad_cnt, int n_rad, +__global__ void radial_correct(const unsigned long long *rad_sum, const int *rad_cnt, const float *curve, int n_rad, const float *k_diff, int k_len, int k_off, const float *isum_a, const int *ninner_a, const int *rbin_a, const int *kbin_a, @@ -615,7 +615,8 @@ __global__ void radial_correct(const unsigned long long *rad_sum, const int *rad for (int k = 0; k < k_len; ++k) { int b = rbin_a[i] + k - k_off; b = min(max(b, 0), n_rad - 1); - const float v = rad_cnt[b] > 0 ? (float) (long long) rad_sum[b] / (float) rad_cnt[b] : bkg; + const float v = curve ? curve[b] + : rad_cnt[b] > 0 ? (float) (long long) rad_sum[b] / (float) rad_cnt[b] : bkg; corr += kern[k] * v; } const float bkg_new = bkg - corr; @@ -906,6 +907,7 @@ BraggIntegrationEngineGPU::BraggIntegrationEngineGPU(const DiffractionExperiment n_rad = static_cast(std::ceil(r_max)) + 2; d_rad_sum = CudaDevicePtr(n_rad); d_rad_cnt = CudaDevicePtr(n_rad); + d_rad_curve = CudaDevicePtr(n_rad); d_k_diff = CudaDevicePtr(k_diff.size()); cuda_err(cudaMemcpyAsync(d_k_diff, k_diff.data(), sizeof(float) * k_diff.size(), cudaMemcpyHostToDevice, *this->stream)); @@ -1015,7 +1017,11 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu // Whether the radial correction runs for THIS image. n_rad only says the buffers exist - under // the auto mode they are allocated for every image and used for the ones the ice score selects. - const int rad_n = bkg_radial ? n_rad : 0; + // Where the caller handed over a curve (RadialBackgroundFromRings) it is read instead, and the + // box sums accumulate nothing. + const bool radial_on = bkg_radial && n_rad > 0; + const bool radial_given = radial_on && !radial_curve.empty(); + const int rad_n = radial_on ? n_rad : 0; // The mask and the owner map can be taken back either by clearing the whole frame or by revisiting // the boxes that were marked, and which is cheaper depends on the detector: on 18 Mpx the marks are @@ -1034,7 +1040,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu if (p.overlap) cuda_err(cudaMemsetAsync(d_owner, 0xff, sizeof(uint32_t) * npixel, *stream)); // BRAGG_OWNER_NONE } - if (rad_n > 0) { + if (rad_n > 0 && !radial_given) { 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)); } @@ -1050,14 +1056,17 @@ 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, d_counts); + d_rad_sum, d_rad_cnt, rad_n, rad_n > 0 && !radial_given, 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 (rad_n > 0) { + 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, rad_n, d_k_diff, k_len, k_off, + d_rad_sum, d_rad_cnt, radial_given ? d_rad_curve.get() : nullptr, n_rad, d_k_diff, k_len, k_off, d_isum, d_ninner, d_rbin, d_kbin, d_ok, d_bkg, d_I, n); cuda_err(cudaGetLastError()); } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h index 8a4de8253..c8d07c6ce 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h @@ -54,6 +54,7 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { CudaDevicePtr d_rad_sum; // integer pixel sums, see boxsum CudaDevicePtr d_k_diff; CudaDevicePtr d_rad_cnt; + CudaDevicePtr d_rad_curve; // the curve handed over (RadialBackgroundFromRings) // --- fixed-size device arrays --- // The learning/fit math is single precision: FP64 is heavily throttled on consumer GPUs and the diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineMetal.mm b/image_analysis/bragg_integration/BraggIntegrationEngineMetal.mm index d086742e2..2d5e7c4d2 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineMetal.mm +++ b/image_analysis/bragg_integration/BraggIntegrationEngineMetal.mm @@ -220,7 +220,9 @@ std::vector BraggIntegrationEngineMetal::Run(const ImagePreprocessor f.inv_claim = inv_claim; f.use_owner = use_owner; f.exclude = overlap == OverlapMode::Exclude && mode != IntegratorMode::BoxSum; - f.n_rad = bkg_radial ? a.n_rad : 0; + // 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; f.static_mask_words = static_cast(static_mask.size()); f.n = static_cast(npredicted); @@ -285,12 +287,12 @@ std::vector BraggIntegrationEngineMetal::Run(const ImagePreprocessor } } if (any_clip) { - if (bkg_radial && dev->rad_sum == nil) { + if (accumulate_radial && dev->rad_sum == nil) { dev->rad_sum = buffer(a.n_rad * 2 * sizeof(uint32_t)); dev->rad_cnt = buffer(a.n_rad * sizeof(int32_t)); } cb = [dev->queue commandBuffer]; - if (bkg_radial) { + if (accumulate_radial) { blit = [cb blitCommandEncoder]; [blit fillBuffer:dev->rad_sum range:NSMakeRange(0, a.n_rad * 2 * sizeof(uint32_t)) value:0]; [blit fillBuffer:dev->rad_cnt range:NSMakeRange(0, a.n_rad * sizeof(int32_t)) value:0]; @@ -306,13 +308,13 @@ std::vector BraggIntegrationEngineMetal::Run(const ImagePreprocessor [enc setBuffer:dev->clips offset:0 atIndex:5]; [enc setBuffer:dev->thr offset:0 atIndex:6]; [enc setBuffer:dev->clip_out offset:0 atIndex:7]; - [enc setBuffer:(bkg_radial ? dev->rad_sum : dev->clips) offset:0 atIndex:8]; // not written without the curve - [enc setBuffer:(bkg_radial ? dev->rad_cnt : dev->clips) offset:0 atIndex:9]; + [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 dispatchThreads:grid threadsPerThreadgroup:group]; [enc endEncoding]; Wait(cb); - if (bkg_radial) { + if (accumulate_radial) { const auto *w = static_cast(dev->rad_sum.contents); const auto *c = static_cast(dev->rad_cnt.contents); for (int b = 0; b < a.n_rad; ++b) { diff --git a/image_analysis/bragg_integration/CMakeLists.txt b/image_analysis/bragg_integration/CMakeLists.txt index a52ca23ea..3a632ffd7 100644 --- a/image_analysis/bragg_integration/CMakeLists.txt +++ b/image_analysis/bragg_integration/CMakeLists.txt @@ -18,6 +18,8 @@ ADD_LIBRARY(JFJochBraggIntegration STATIC SpotWidth.h BackgroundGain.cpp BackgroundGain.h + BackgroundAperture.cpp + BackgroundAperture.h ) TARGET_LINK_LIBRARIES(JFJochBraggIntegration JFJochImagePreprocessing JFJochCommon) diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index b0d485b71..2f02f1edf 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -45,6 +45,7 @@ #include "../image_analysis/WriteReflections.h" // OutputProvenance #include "../image_analysis/bragg_integration/SpotWidth.h" // FluxCurve #include "../image_analysis/bragg_integration/SpotFootprint.h" // FootprintSpot +#include "../image_analysis/bragg_integration/BackgroundAperture.h" // background_aperture::Tally #include "DiagnosticOutput.h" // PerImagePlot class JFJochReader; @@ -189,7 +190,9 @@ struct ProcessConfig { // aperture that owes nothing to the integrator (SpotWidth.h) - and sets r1 from it. r1 is the // aperture the integrator learns the profile WIDTH over, and a second moment over a disk of radius // a saturates at a^2/4, so at the shipped r1 = 4 a spot wider than sigma = 2 px is fitted with a - // profile the model cannot represent. Ignored where the radii were set by hand. + // profile the model cannot represent. On stills the signal disk is kept and the pre-scan measures + // the background ring's outer radius instead, on empty apertures (BackgroundAperture.h). Ignored + // where the radii were set by hand. bool adaptive_integration_radius = false; // The input's pixel values are ADU, not counted photons (SMV and marCCD: CCDs). The pre-scan then @@ -1331,6 +1334,7 @@ class Rugnux { std::vector footprint; std::vector spot_q; size_t width_images = 0; + background_aperture::Tally apertures; }; struct Defects { size_t pixels = 0; diff --git a/rugnux/RugnuxPreScan.cpp b/rugnux/RugnuxPreScan.cpp index 834c8546a..1eca770ad 100644 --- a/rugnux/RugnuxPreScan.cpp +++ b/rugnux/RugnuxPreScan.cpp @@ -7,6 +7,7 @@ #include "../image_analysis/indexing/SpindleCuspLoss.h" #include "../image_analysis/bragg_integration/SpotWidth.h" #include "../image_analysis/bragg_integration/BackgroundGain.h" +#include "../image_analysis/bragg_integration/BackgroundAperture.h" #include "../image_analysis/bragg_integration/SpotFootprint.h" #include "../image_analysis/bragg_integration/BraggStencil.h" // BRAGG_FOOTPRINT_NSIGMA #include "../image_analysis/hot_pixels/HotPixels.h" @@ -366,10 +367,21 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru return w; }; + // 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. + 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(); + 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) { + std::vector &spot_q, std::vector &footprint, + background_aperture::Tally &apertures) { try { // As in the image loops: a frame the device cannot decode goes to the host decoder. ImageStatistics stats; @@ -413,6 +425,11 @@ 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()) + 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, apertures); }; std::vector prescan_spot_q; @@ -422,6 +439,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; // 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; @@ -595,6 +614,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()); // 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); @@ -641,7 +661,8 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru const size_t i = visit[v]; const int image_idx = start_image + ordinals[i] * config_.stride; if (const CompressedImage *image = frame(i, w.raw_image)) - find_spots(w, *image, image_idx, true, curves_of[i], spot_q_of[i], footprint_of[i]); + find_spots(w, *image, image_idx, true, curves_of[i], spot_q_of[i], footprint_of[i], + apertures_of[i]); { std::lock_guard lock(m); done_in_tier[tier_of(v)]++; @@ -708,11 +729,21 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru for (size_t i = 0; i < ordinals.size(); i++) if (!shape_only[i]) prescan_spot_q.insert(prescan_spot_q.end(), spot_q_of[i].begin(), spot_q_of[i].end()); + // In sample order, over the frames the width was taken on. + aperture_tally.sum_d2.assign(aperture_r3.size(), 0.0); + for (size_t i = 0; i < ordinals.size(); i++) + if (position[i] < width_images && !apertures_of[i].sum_d2.empty()) { + for (size_t k = 0; k < aperture_r3.size(); k++) + aperture_tally.sum_d2[k] += apertures_of[i].sum_d2[k]; + aperture_tally.apertures += apertures_of[i].apertures; + } // 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()); }; std::future spots_measured; + // The gain is applied above, before the spots are measured, so the apertures clip as the run will. + aperture_gain = experiment_.GetBraggIntegrationSettings().GetBackgroundPixelGain(); if (want_spots) spots_measured = std::async(std::launch::async, spot_measurement); @@ -725,7 +756,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru if (!measured) measured.emplace(); measured->spots = PreScanMemo::Spots{width_curves, shape_curves, footprint_spots, prescan_spot_q, - width_images}; + width_images, aperture_tally}; } else if (want_width) { logger.Info("Spot width: the spots kept from an earlier pre-scan of this run on the same frames, " "geometry, mask and spot finding"); @@ -735,6 +766,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru footprint_spots = k.footprint; prescan_spot_q = k.spot_q; width_images = k.width_images; + aperture_tally = k.apertures; } // The integration radius from the recorded spot width. r1 is the aperture the integrator learns @@ -743,7 +775,33 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // the model cannot represent. r2 follows r1, and r3 is taken so that the r2..r3 background ring // 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. - if (want_width) { + // + // 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. + if (want_width && stills) { + spot_width_measured_ = true; + const int best = background_aperture::Best(aperture_tally); + if (best < 0) { + logger.Info("Background ring: no empty aperture measured on {} images, keeping r3={:.1f}", + width_images, experiment_.GetBraggIntegrationSettings().GetR3()); + } else { + std::string table; + for (size_t k = 0; k < aperture_r3.size(); k++) + table += fmt::format(" {:.0f}:{:.4g}", aperture_r3[k], + aperture_tally.sum_d2[k] / static_cast(aperture_tally.apertures)); + BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings(); + bis.R3(aperture_r3[best]); + experiment_.ImportBraggIntegrationSettings(bis); + 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){} => integration radii " + "r1={:.1f} r2={:.1f} r3={:.1f}", aperture_tally.apertures, width_images, + bis.GetR2(), bis.GetR1(), table, bis.GetR1(), bis.GetR2(), bis.GetR3()); + Note(fmt::format("background ring {:.0f}..{:.0f} px, measured on empty apertures", + bis.GetR2(), bis.GetR3())); + } + } + if (want_width && !stills) { spot_width_measured_ = true; const auto r80 = spot_width::R80AtReference(width_curves); if (!r80) { @@ -824,8 +882,9 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // The X-ray bandwidth, read off the shapes of the same spots (spot_width::EstimateBandwidth). A // value the file states, or --bandwidth, wins; otherwise a significant estimate is the bandwidth // of the run - prediction, partiality, the profile's radial width and the background stencil all - // take it - and anything short of significant leaves the beam monochromatic. - if (want_width) { + // take it - and anything short of significant leaves the beam monochromatic. Not on stills, where the + // estimate has not been shown to measure the beam (on the serial sets it reads not significant). + if (want_width && !stills) { const auto &det = experiment_.GetDetectorSetup(); const double lambda = experiment_.GetWavelength_A(); const auto estimate = spot_width::EstimateBandwidth( diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index b2e641876..3256a54d7 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -249,12 +249,12 @@ void print_usage() { std::cout << " Integration" << std::endl; std::cout << " --bandwidth Relative X-ray bandwidth FWHM (e.g. 0.01 for 1% DMM). Default: the file's value, else measured in the pre-scan from the radial elongation of strong spots and used where significant, else 0; --bandwidth 0 forces a monochromatic beam" << std::endl; std::cout << " --integration-radius Signal-box radius r1, or r1,r2,r3 (px). One value => r2=r1+2, r3=r1+4" << std::endl; - std::cout << " --adaptive-integration-radius[=on|off] Set the signal radius r1 from how wide this crystal's spots actually are, measured in the pre-scan (default: on for rotation, off for stills). r1 is the aperture the integrator learns the profile WIDTH over, and a second moment over a disk of radius a saturates at a^2/4, so at the shipped r1=4 the learned sigma can never exceed 2 px and a broader spot is fitted with a profile the model cannot represent. r80 - the radius holding 80% of a spot's flux, at a common resolution - is read off isolated strong spots over a fixed aperture that owes nothing to r1, and r1 = clamp(round(2*r80), 4, 6); the background ring keeps the area it has at the default 4,6,13. Ignored when --integration-radius is given" << std::endl; + std::cout << " --adaptive-integration-radius[=on|off] Set the signal radius r1 from how wide this crystal's spots actually are, measured in the pre-scan (default on). On stills it measures instead how far the background ring should reach: on empty apertures of the pre-scan frames, the r3 (up to 14) whose ring predicts the background under the signal disk best. r1 is the aperture the integrator learns the profile WIDTH over, and a second moment over a disk of radius a saturates at a^2/4, so at the shipped r1=4 the learned sigma can never exceed 2 px and a broader spot is fitted with a profile the model cannot represent. r80 - the radius holding 80% of a spot's flux, at a common resolution - is read off isolated strong spots over a fixed aperture that owes nothing to r1, and r1 = clamp(round(2*r80), 4, 6); the background ring keeps the area it has at the default 4,6,13. Ignored when --integration-radius is given" << std::endl; std::cout << " --integration-stencil Push the r2..r3 background ring out by k times the beam's radial streak (bandwidth*Rpx), per reflection (default 0 = a fixed circular ring). A fixed ring otherwise ends up on a streaked reflection's own tails at high resolution and measures them as background. Only the ring moves, and only radially - the r1 signal box stays a circle. On a monochromatic beam the streak is zero and this does nothing" << std::endl; std::cout << " --integration-high-resolution High resolution limit for prediction/integration. If omitted (or 0), integration extends as far as the detector reaches" << std::endl; std::cout << " --max-hkl Predict reflections with |h|,|k|,|l| <= n. Default: derived per crystal from the refined cell (ceil(longest axis / d_min) + 1), which is the exact bound - set it only to override that" << std::endl; std::cout << " --background-clip High-side clip of the background ring at mean + n*sqrt(mean) (default 4; 0 = off). This is the default background estimator - it rejects neighbour cores and zingers without the symmetric trim's Poisson skew bias. Ignored by --integrator boxsum" << std::endl; - std::cout << " --background-radial[=on|off|auto] Correct the background ring for the CURVATURE of the radial background (default auto). The signal disk and the background ring are concentric, so a background linear in position cancels between them and only curvature survives - which on a smooth ice ring reaches +26 counts on a single reflection. =auto applies it per image where that image's ice score shows a smooth powder ring, which is where a radius-only background model holds; on ice made of discrete crystallite spots there is no such ring and the correction makes the bias worse. Costs one short dot product per reflection and no extra pixel reads" << std::endl; + std::cout << " --background-radial[=on|off|auto] Correct the background ring for the CURVATURE of the radial background (default auto). The signal disk and the background ring are concentric, so a background linear in position cancels between them and only curvature survives - which on a smooth ice ring reaches +26 counts on a single reflection. =auto applies it per image where that image's ice score shows a smooth powder ring, which is where a radius-only background model holds; on ice made of discrete crystallite spots there is no such ring and the correction makes the bias worse. On stills =auto applies it to every image, with the curve taken from the spot finder's peak-excluded background of the whole frame: a smooth solvent or carrier hump curves the background as much as ice does. Costs one short dot product per reflection and no extra pixel reads" << std::endl; std::cout << " --background-trim Use the old symmetric trimmed mean for the background ring instead of the clip (0<=f<0.5; 0.10 was the former default). Switches --background-clip off. A symmetric trim is biased low on Poisson data and adds ~5 counts to every partial, so this is for back compatibility only; 0 = plain ring mean" << std::endl; std::cout << " --overlap What to do where two predicted reflections share signal pixels: off|reject|exclude (default exclude). A pixel inside two signal disks belongs to the nearer centre; before this, nothing kept a neighbour's flux out of a reflection's own disk, so on a dense pattern a crowded reflection read high. exclude drops the shared PIXELS from the profile fit, which renormalises itself, so the reflection is kept; reject instead drops the whole reflection when less than --overlap-minpk of its expected profile is cleanly its own (what XDS calls MINPK). --integrator boxsum has no profile to renormalise with, so exclude does nothing there and only reject acts" << std::endl; std::cout << " --overlap-minpk Least fraction of a reflection's expected profile that must be usable for the reflection to be kept (default 0.75, XDS MINPK). Governs both ways part of a profile is lost: the fraction that must be READABLE - not masked, untrusted, in a detector gap or overloaded - in every profile mode, and, under --overlap reject, the fraction that must be cleanly the reflection's own. --integrator boxsum has no profile to renormalise with, so there a signal disk with any unreadable pixel is still discarded outright and the reject fraction is by disk AREA, which cuts harder" << std::endl; @@ -3308,15 +3308,15 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { } else if (!rotation_indexing) { // An explicit --integration-radius always wins over the stills default. ApplyRugnuxIntegrationRadii(experiment, rotation_indexing); - logger.Info("Stills integration radii default to r1=6.0 r2=8.0 r3=14.0 (override with --integration-radius)"); + logger.Info("Stills integration radii start at r1=6.0 r2=8.0 r3=14.0, the ring's reach measured in the pre-scan (override with --integration-radius)"); } - // Adaptive signal radius: measured in the pre-scan and applied there (Rugnux::PreScan). Radii the - // user set by hand always win, and stills already default to a wide box, so the rule is on for - // rotation only. + // Adaptive radii: measured in the pre-scan and applied there (Rugnux::PreScan) - on rotation the + // signal radius from the spot width, on stills the background ring's reach from empty apertures. + // Radii the user set by hand always win. const bool adaptive_integration_radius = !integration_radius_arg - && adaptive_integration_radius_arg.value_or(rotation_indexing); + && adaptive_integration_radius_arg.value_or(true); if (forced_prediction_mosaicity_arg) { BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings(); @@ -3427,7 +3427,9 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { const auto radial = experiment.GetBraggIntegrationSettings().GetBackgroundRadialCorrection(); logger.Info("Background ring: radial curvature correction {}", radial.has_value() ? (*radial ? "on" : "off") - : "auto (per image, where the ice score sees a smooth ring)"); + : rotation_indexing + ? "auto (per image, where the ice score sees a smooth ring)" + : "auto (every image, from the radial background of the whole frame)"); } // Detection defaults differ by workflow (RugnuxDefaultSpotFinding); each flag always wins. diff --git a/tests/BackgroundApertureTest.cpp b/tests/BackgroundApertureTest.cpp new file mode 100644 index 000000000..3121546ed --- /dev/null +++ b/tests/BackgroundApertureTest.cpp @@ -0,0 +1,72 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include +#include + +#include "../image_analysis/bragg_integration/BackgroundAperture.h" +#include "../common/JFJochMath.h" + +namespace { + +constexpr int W = 600, H = 600; + +// A background of `level` counts a pixel with Poisson noise, modulated by `ripple` (relative) with a +// period of 20 px along x - the structure a ring further from the disk reads wrong. +ImagePreprocessorBuffer Frame(double level, double ripple, uint32_t seed) { + std::mt19937 rng(seed); + ImagePreprocessorBuffer img(static_cast(W) * H); + for (int y = 0; y < H; y++) + for (int x = 0; x < W; x++) { + std::poisson_distribution p(level * (1.0 + ripple * std::sin(2.0 * PI * x / 20.0))); + img[static_cast(y) * W + x] = p(rng); + } + return img; +} + +int BestR3(double level, double ripple) { + const auto r3 = background_aperture::CandidateR3(8.0f, 14.0f); + background_aperture::Tally tally; + for (uint32_t seed = 1; seed <= 4; seed++) { + auto img = Frame(level, ripple, seed); + background_aperture::MeasureFrame(img, W, H, {}, 6.0f, 8.0f, r3, 4.0f, 1.0f, tally); + } + REQUIRE(tally.apertures > 100); + return static_cast(r3[background_aperture::Best(tally)]); +} + +} // namespace + +TEST_CASE("BackgroundAperture_CandidateR3", "[Integration]") { + const auto r3 = background_aperture::CandidateR3(8.0f, 14.0f); + REQUIRE(r3.size() == 5); + CHECK(r3.front() == 10.0f); + CHECK(r3.back() == 14.0f); +} + +// A few counts a pixel on a flat background: the noise wins, and the widest ring is kept. +TEST_CASE("BackgroundAperture_FlatNoisyKeepsWidestRing", "[Integration]") { + CHECK(BestR3(0.5, 0.0) == 14); +} + +// Thousands of counts a pixel on a background that varies over the stencil's own scale: the structure +// wins, and the narrowest ring is kept. +TEST_CASE("BackgroundAperture_StructuredBrightTakesNarrowestRing", "[Integration]") { + CHECK(BestR3(2000.0, 0.05) == 10); +} + +// An aperture over an unreadable pixel is not used. +TEST_CASE("BackgroundAperture_SkipsUnreadablePixels", "[Integration]") { + const auto r3 = background_aperture::CandidateR3(8.0f, 14.0f); + auto img = Frame(1.0, 0.0, 7); + background_aperture::Tally all; + background_aperture::MeasureFrame(img, W, H, {}, 6.0f, 8.0f, r3, 4.0f, 1.0f, all); + for (size_t i = 0; i < static_cast(W) * H; i += 97) + img[i] = INT32_MIN; + background_aperture::Tally masked; + background_aperture::MeasureFrame(img, W, H, {}, 6.0f, 8.0f, r3, 4.0f, 1.0f, masked); + CHECK(masked.apertures == 0); + CHECK(all.apertures > 0); +} diff --git a/tests/BraggIntegrationEngineCPUTest.cpp b/tests/BraggIntegrationEngineCPUTest.cpp index 067723b69..292ff9ee4 100644 --- a/tests/BraggIntegrationEngineCPUTest.cpp +++ b/tests/BraggIntegrationEngineCPUTest.cpp @@ -10,6 +10,7 @@ #include "../common/BraggIntegrationSettings.h" #include "../common/DetectorSetup.h" #include "../common/DiffractionExperiment.h" +#include "../common/JFJochMath.h" #include "../common/Reflection.h" #include "../image_analysis/bragg_integration/BraggIntegrationEngineCPU.h" #include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h" @@ -148,3 +149,57 @@ TEST_CASE("BraggIntegrationEngineCPU_SaturatedPeakIsFlaggedNotDropped", "[Integr CHECK(saturated[0].overloaded); CHECK(saturated[0].clipped); } + +// An empty aperture on a background that is a curved function of the radius alone: the flat ring mean +// reads it wrong by the curvature, and the radial correction with the curve taken from a per-ring +// profile of the whole frame (as the spot finder hands it over on stills) takes that back out. +TEST_CASE("BraggIntegrationEngineCPU_RadialCurveFromRings", "[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 50.0 + 0.05 * (r - 300.0) * (r - 300.0); }; + 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( + background(std::hypot(static_cast(i % width) - bx, static_cast(i / width) - by)))); + + // The profile per q bin, as the spot finder's rings would hold it. + 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))))); + } + rings[100] = NAN; // a ring with no value is bridged + + Reflection r{}; + r.h = 1; r.k = 0; r.l = 0; + r.predicted_x = bx + 300.0f; + r.predicted_y = by; + r.d = 2.0f; + r.prescaling_corr = 1.0f; + r.partiality = 1.0f; + const std::vector predicted{r}; + + auto integrate = [&](bool correct) { + experiment.ImportBraggIntegrationSettings( + BraggIntegrationSettings().Integrator(IntegratorMode::BoxSum).BackgroundRadialCorrection(correct)); + BraggIntegrationEngineCPU engine(experiment, PixelMask(experiment)); + engine.RadialBackgroundFromRings(rings, static_cast(rings.size()), low_q, dq); + REQUIRE(!engine.GetRadialBackground().empty()); + const auto out = engine.Run(image, predicted, predicted.size(), 0); + REQUIRE(out.size() == 1); + return out[0].I; + }; + // The curvature is along the radius only, so the error is n_inner x 0.05 x (_ring - _disk) + // with u the radial offset: 49 px x 0.05 x ((6^2 + 13^2) / 4 - 4^2 / 4) = 116 counts read too low. + CHECK(integrate(false) < -80.0f); + CHECK(std::fabs(integrate(true)) < 5.0f); + + // An empty profile hands the curve back to the reflections' own rings. + BraggIntegrationEngineCPU engine(experiment, PixelMask(experiment)); + engine.RadialBackgroundFromRings({}, 500, low_q, dq); + CHECK(engine.GetRadialBackground().empty()); +} diff --git a/tests/BraggIntegrationEngineGPUTest.cpp b/tests/BraggIntegrationEngineGPUTest.cpp index c4030bc4f..4f8788703 100644 --- a/tests/BraggIntegrationEngineGPUTest.cpp +++ b/tests/BraggIntegrationEngineGPUTest.cpp @@ -15,6 +15,7 @@ #include "../common/BraggIntegrationSettings.h" #include "../common/DetectorSetup.h" #include "../common/DiffractionExperiment.h" +#include "../common/JFJochMath.h" #include "../common/Reflection.h" #include "../image_analysis/bragg_integration/BraggIntegrationEngineCPU.h" #include "../image_analysis/bragg_integration/BraggIntegrationEngineGPU.h" @@ -413,6 +414,60 @@ TEST_CASE("BraggIntegrationEngineGPU_ReusedEngineMatchesFresh") { } } +// The radial correction with its curve handed over from a per-ring profile (RadialBackgroundFromRings), as +// on stills: the GPU reads the same curve and corrects the same empty aperture as the CPU. +TEST_CASE("BraggIntegrationEngineGPU_RadialCurveFromRingsMatchesCPU", "[Integration]") { + if (get_gpu_count() == 0) { + WARN("No CUDA GPU present. Skipping"); + return; + } + 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); }; + const size_t width = experiment.GetXPixelsNum(), npixel = experiment.GetPixelsNum(); + ImagePreprocessorBuffer cpu_image(npixel); + ImagePreprocessorBufferGPU gpu_image(npixel); + for (size_t i = 0; i < npixel; ++i) { + cpu_image[i] = static_cast(std::lround( + background(std::hypot(static_cast(i % width) - bx, static_cast(i / width) - by)))); + gpu_image[i] = cpu_image[i]; + } + auto stream = std::make_shared(); + REQUIRE(cudaMemcpyAsync(gpu_image.getGPUBuffer(), gpu_image.getBuffer().data(), + npixel * sizeof(int32_t), cudaMemcpyHostToDevice, *stream) == cudaSuccess); + + 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))))); + } + + std::vector predicted; + for (int i = 0; i < 8; ++i) + predicted.push_back(MakeReflection(bx + 300.0f * std::cos(0.7f * i), by + 300.0f * std::sin(0.7f * i), 2.0f, i + 1)); + + for (const auto mode : {IntegratorMode::BoxSum, IntegratorMode::ProfileGaussian}) { + experiment.ImportBraggIntegrationSettings(BraggIntegrationSettings().Integrator(mode).BackgroundRadialCorrection(true)); + BraggIntegrationEngineCPU cpu(experiment, PixelMask(experiment)); + BraggIntegrationEngineGPU gpu(experiment, stream, PixelMask(experiment)); + cpu.RadialBackgroundFromRings(rings, static_cast(rings.size()), low_q, dq); + gpu.RadialBackgroundFromRings(rings, static_cast(rings.size()), low_q, dq); + const auto out_cpu = cpu.Run(cpu_image, predicted, predicted.size(), 0); + const auto out_gpu = gpu.Run(gpu_image, predicted, predicted.size(), 0); + REQUIRE(out_cpu.size() == predicted.size()); + REQUIRE(out_gpu.size() == out_cpu.size()); + 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); + } + } +} + // Hidden ([.]) benchmark: the raison d'etre of the GPU port is < 2 ms/frame (vs ~142 ms on the CPU // for ProfileIntegrate2D). Run explicitly with: ./jfjoch_test "[bragg_bench]" TEST_CASE("BraggIntegrationEngineGPU_Benchmark", "[.][bragg_bench]") { diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 3eac9e404..ff73e96f7 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -92,6 +92,7 @@ ADD_EXECUTABLE(jfjoch_test BraggIntegrationEngineCompressedImageTest.cpp BraggIntegrationEngineCPUTest.cpp BackgroundGainTest.cpp + BackgroundApertureTest.cpp BraggStencilTest.cpp SpotFootprintTest.cpp SpotTilesTest.cpp