diff --git a/common/BraggIntegrationSettings.cpp b/common/BraggIntegrationSettings.cpp index 03b1a92bf..76b0752ae 100644 --- a/common/BraggIntegrationSettings.cpp +++ b/common/BraggIntegrationSettings.cpp @@ -166,6 +166,17 @@ float BraggIntegrationSettings::GetBackgroundClipNSigma() const { return bkg_clip_nsigma; } +BraggIntegrationSettings &BraggIntegrationSettings::BackgroundPixelGain(float input) { + check_finite("Background pixel gain", input); + check_min("Background pixel gain", input, 1.0); + bkg_pixel_gain = input; + return *this; +} + +float BraggIntegrationSettings::GetBackgroundPixelGain() const { + return bkg_pixel_gain; +} + BraggIntegrationSettings &BraggIntegrationSettings::BackgroundRadialCorrection(std::optional input) { bkg_radial_correction = input; return *this; diff --git a/common/BraggIntegrationSettings.h b/common/BraggIntegrationSettings.h index 3555ccbdc..9c8023add 100644 --- a/common/BraggIntegrationSettings.h +++ b/common/BraggIntegrationSettings.h @@ -97,6 +97,12 @@ class BraggIntegrationSettings { // engine applies (rugnux --background-clip); the front end picks the default, and rugnux lowers it // to 3 sigma for broadband (non-zero bandwidth) data, where longer spots leak further into the ring. float bkg_clip_nsigma = 4.0f; + // Variance of a background pixel per unit of its mean, the "n*sqrt(mean)" of the clip above being + // n*sqrt(this * mean). 1 for a detector that counts photons, which is what every HDF5 and CBF input + // is taken to be. A CCD reports ADU instead, and on a Rigaku Saturn a pixel scatters 4-5x its mean: + // there the photon clip cuts the background's upper tail at about 2 sigma and adds a positive offset + // to every reflection. Set by rugnux from the pre-scan for SMV and marCCD input (BackgroundGain.h). + float bkg_pixel_gain = 1.0f; // Radial background curvature correction. The signal disk and the background annulus are // concentric, so for ANY background linear in position their means are equal - a plane fit buys // nothing and the leading error is the CURVATURE of the radial background, which the flat annulus @@ -174,6 +180,7 @@ public: BraggIntegrationSettings& Integrator(IntegratorMode input); BraggIntegrationSettings& BackgroundTrimFraction(float input); BraggIntegrationSettings& BackgroundClipNSigma(float input); + BraggIntegrationSettings& BackgroundPixelGain(float input); BraggIntegrationSettings& BackgroundRadialCorrection(std::optional input); BraggIntegrationSettings& MaxHKL(std::optional input); BraggIntegrationSettings& Overlap(OverlapMode input); @@ -194,6 +201,7 @@ public: [[nodiscard]] float GetMinimumSigmaInRegardsToI() const; [[nodiscard]] float GetBackgroundTrimFraction() const; [[nodiscard]] float GetBackgroundClipNSigma() const; + [[nodiscard]] float GetBackgroundPixelGain() const; // Unset = auto (gate per image on the smooth-ice score); see bkg_radial_correction. [[nodiscard]] std::optional GetBackgroundRadialCorrection() const; [[nodiscard]] std::optional GetMaxHKL() const; diff --git a/image_analysis/bragg_integration/BackgroundGain.cpp b/image_analysis/bragg_integration/BackgroundGain.cpp new file mode 100644 index 000000000..b6732ec19 --- /dev/null +++ b/image_analysis/bragg_integration/BackgroundGain.cpp @@ -0,0 +1,95 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "BackgroundGain.h" + +#include +#include +#include +#include + +namespace background_gain { + +namespace { +constexpr int PAIR_DISTANCE = 2; // columns between the two pixels of a pair +constexpr size_t ROW_STEP = 4; // every fourth row: a quarter of a frame is ample for one number +constexpr double CLIP_NSIGMA = 4.0; +constexpr int CLIP_ROUNDS = 3; +const double LOG_BIN = std::log(1.25); +} + +Levels MeasureFrame(const std::vector &image, size_t width, size_t height, + const std::vector &mask, int64_t saturation) { + const auto usable = [&](size_t i) { + return mask[i] == 0 && image[i] >= 1 && image[i] < saturation; + }; + // The level of every usable pair, worked out once; -1 where the pair is not usable. + std::vector bin(width * height, -1); + for (size_t y = 0; y < height; y += ROW_STEP) + for (size_t x = 0; x + PAIR_DISTANCE < width; x++) { + const size_t i = y * width + x; + if (!usable(i) || !usable(i + PAIR_DISTANCE)) + continue; + const double m = 0.5 * (static_cast(image[i]) + image[i + PAIR_DISTANCE]); + bin[i] = static_cast(std::min(N_BINS - 1, static_cast(std::log(m) / LOG_BIN))); + } + + // Clip about zero at CLIP_NSIGMA of the level's own scatter, starting from no clip at all: a few + // rounds take the Bragg peaks and zingers out without needing any idea of the gain beforehand. + std::array limit2; + limit2.fill(std::numeric_limits::infinity()); + Levels out; + for (int round = 0; round <= CLIP_ROUNDS; round++) { + out = Levels{}; + for (size_t i = 0; i < bin.size(); i++) { + const int b = bin[i]; + if (b < 0) + continue; + const double d = static_cast(image[i]) - image[i + PAIR_DISTANCE]; + if (d * d > limit2[b]) + continue; + out.n[b]++; + out.sum_mean[b] += 0.5 * (static_cast(image[i]) + image[i + PAIR_DISTANCE]); + out.sum_diff2[b] += d * d; + } + for (int b = 0; b < N_BINS; b++) + if (out.n[b] > 0) + limit2[b] = CLIP_NSIGMA * CLIP_NSIGMA * out.sum_diff2[b] / static_cast(out.n[b]); + } + return out; +} + +std::optional Gain(const std::vector &frames) { + // (variance/mean, pairs) per level + std::vector> ratio; + uint64_t total = 0; + for (int b = 0; b < N_BINS; b++) { + uint64_t n = 0; + double sum_mean = 0.0, sum_diff2 = 0.0; + for (const auto &f : frames) { + n += f.n[b]; + sum_mean += f.sum_mean[b]; + sum_diff2 += f.sum_diff2[b]; + } + // var(a - b) = 2 var(pixel) + if (n >= MIN_PAIRS_PER_BIN && sum_mean > 0.0) { + ratio.emplace_back(0.5 * sum_diff2 / sum_mean, n); + total += n; + } + } + if (ratio.size() < MIN_BINS) + return std::nullopt; + // The median over PAIRS, not over levels: the background holds most of the pixels in a few levels, + // and the Bragg peaks spread a few pixels over many more - an unweighted median over the levels of + // a weak-background marCCD frame lands on the peaks and reads 5.7 where the background reads 0.3. + std::sort(ratio.begin(), ratio.end()); + uint64_t cum = 0; + for (const auto &[r, n] : ratio) { + cum += n; + if (2 * cum >= total) + return static_cast(r); + } + return static_cast(ratio.back().first); +} + +} // namespace background_gain diff --git a/image_analysis/bragg_integration/BackgroundGain.h b/image_analysis/bragg_integration/BackgroundGain.h new file mode 100644 index 000000000..de9c9ae75 --- /dev/null +++ b/image_analysis/bragg_integration/BackgroundGain.h @@ -0,0 +1,50 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include +#include +#include + +// The background noise of a detector that reports ADU rather than counted photons. For a photon +// counter a background pixel's variance equals its mean; a CCD's pixel value is a gain times the +// photons plus read noise, so its variance is g*mean (+ a constant). The gain is not in the headers +// (an SMV file states none, every marCCD file states the writer's default of 1 ADU/photon), so it is +// measured: the variance per unit mean of the difference between two pixels two columns apart, which +// share the background, but are far enough apart that the detector's point spread does not correlate +// them. Bragg peaks and zingers are clipped per level. Taken per level (bins of the pair's mean, a +// factor 1.25 wide) and summarised as the median ratio over all pairs, i.e. over the background that +// holds most of the pixels. +// +// The variance/mean of pixels in a neighbourhood is how XDS fills its GAIN table (INIT step, +// Kabsch (2010) Acta Cryst. D66, 125-132). +namespace background_gain { + +constexpr int N_BINS = 64; // levels 1.25^0 .. 1.25^63, i.e. up to ~1.3e6 ADU +constexpr size_t MIN_PAIRS_PER_BIN = 1000; +constexpr size_t MIN_BINS = 3; +// Frames of the pre-scan sample it is measured on, spread over the sample. The answer moves by under +// 1 % between frames of a sweep, so more only cost time. +constexpr size_t GAIN_MAX_FRAMES = 16; + +// Per level, over the pairs that survived the clip: the count, the sum of the pair means and the sum +// of the squared differences. +struct Levels { + std::array n{}; + std::array sum_mean{}; + std::array sum_diff2{}; +}; + +// One frame. A pixel takes part if the mask is zero there and 1 <= value < saturation (a CCD reads +// its bias in the dark, so a zero is a pixel outside the active area). +Levels MeasureFrame(const std::vector &image, size_t width, size_t height, + const std::vector &mask, int64_t saturation); + +// The gain over frames: the pair-weighted median of the per-level variance/mean, or nothing with +// fewer than MIN_BINS levels holding MIN_PAIRS_PER_BIN pairs. +std::optional Gain(const std::vector &frames); + +} // namespace background_gain diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp index eb453ad00..28a40d186 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp @@ -115,6 +115,7 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe // symmetric trimmed mean (rugnux --background-trim). The caller owns the choice - the engine no // longer overrides it for broadband data. bkg_clip_nsigma = settings.GetBackgroundClipNSigma(); + bkg_pixel_gain = settings.GetBackgroundPixelGain(); bkg_trim = bkg_clip_nsigma > 0.0f ? 0.0f : settings.GetBackgroundTrimFraction(); // Overlap treatment. Ownership is decided out to the fit grid's half size, which is where the diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index bc62af4b1..a5edef530 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -175,6 +175,7 @@ protected: double bw_sigma; // bandwidth sigma [dimensionless, * Rpx -> px] float bkg_clip_nsigma; // high-outlier background sigma-clip multiplier (0 = no clip) + float bkg_pixel_gain = 1.0f; // background pixel variance per unit mean the clip scales by bool use_ellipse; // radially elongate the per-reflection Gaussian double c_radial; // radial variance coefficient of tan^2(2theta): parallax + capture diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 3f57f5744..3bb37b74d 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -311,8 +311,10 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, } } else if (do_clip) { // One high-outlier sigma-clip pass on the background ring: reject pixels above - // mean + n*sqrt(mean) to strip a neighbour core or a zinger that biases the mean. - const double thr = out.bkg + bkg_clip_nsigma * std::sqrt(std::max(out.bkg, 1.0)); + // mean + n*sqrt(gain*mean) to strip a neighbour core or a zinger that biases the mean. + // The gain is 1 for counted photons, so this is then exactly the Poisson clip. + const double thr = out.bkg + bkg_clip_nsigma * std::sqrt(static_cast(bkg_pixel_gain) + * std::max(out.bkg, 1.0)); double s = 0.0; int n = 0; for (int y = y0; y <= y1; ++y) diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 0f458ce6e..62ecc50b5 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -19,6 +19,7 @@ struct BraggGpuParams { float r1_sq, r2, r2_sq, r3, r3_sq; int R, G, GG; float bkg_clip_nsigma; // high-side background sigma-clip multiplier (0 = no clip) + float bkg_pixel_gain; // background pixel variance per unit mean (1 = counted photons) int empirical; // ProfileEmpirical vs ProfileGaussian int use_ellipse; float bw_sigma; @@ -284,7 +285,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, if (s_nown < p.minpk * s_ndisk) atomicAdd(counts + COUNT_OVERLAPPED_MAJOR, 1ull); } s_bkg = s_accept ? ((double) (long long) s_bkgsum / (double) s_nbkg) : 0.0; - s_thr = s_bkg + (double) p.bkg_clip_nsigma * sqrt(fmax(s_bkg, 1.0)); + s_thr = s_bkg + (double) p.bkg_clip_nsigma * sqrt((double) p.bkg_pixel_gain * fmax(s_bkg, 1.0)); // The pixel is an integer, so comparing it against the floor of the threshold accepts // exactly the same set - and does it with an integer compare instead of a widening // conversion and a double comparison, per ring pixel. @@ -965,6 +966,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu .r1_sq = r1_sq, .r2 = r2, .r2_sq = r2_sq, .r3 = r3, .r3_sq = r3_sq, .R = R, .G = G, .GG = GG, .bkg_clip_nsigma = mode != IntegratorMode::BoxSum ? bkg_clip_nsigma : 0.0f, + .bkg_pixel_gain = bkg_pixel_gain, .empirical = empirical ? 1 : 0, .use_ellipse = use_ellipse ? 1 : 0, .bw_sigma = static_cast(bw_sigma), .c_radial = static_cast(c_radial), diff --git a/image_analysis/bragg_integration/CMakeLists.txt b/image_analysis/bragg_integration/CMakeLists.txt index 3aca9bfcc..ac34468cf 100644 --- a/image_analysis/bragg_integration/CMakeLists.txt +++ b/image_analysis/bragg_integration/CMakeLists.txt @@ -11,6 +11,8 @@ ADD_LIBRARY(JFJochBraggIntegration STATIC SystematicAbsence.h SpotWidth.cpp SpotWidth.h + BackgroundGain.cpp + BackgroundGain.h ) TARGET_LINK_LIBRARIES(JFJochBraggIntegration JFJochImagePreprocessing JFJochCommon) diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 1147a40f1..aa47acdde 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -137,6 +137,12 @@ struct ProcessConfig { // profile the model cannot represent. 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 + // measures the background noise per unit level and the background clip is scaled by it + // (BraggIntegrationSettings::bkg_pixel_gain, background_gain::Gain). Set by the front end from the + // file format; HDF5 and CBF input never has it. + bool counts_in_adu = false; + // Rotation two-pass geometry post-refinement (FullAnalysis, rotation only; on by default in the rugnux // CLI, --rotation-no-postrefine disables it). When set, a first pass integrates and post-refines the // detector distance + beam (from the observed spot positions) and the cell scale + rotation axis (from the @@ -754,6 +760,8 @@ class Rugnux { // The recorded spot width the pre-scan measured (config_.adaptive_integration_radius), kept so the // two-pass rotation run measures it once and both passes integrate at the same radius. bool spot_width_measured_ = false; + // The background gain (config_.counts_in_adu) is measured once, on the first pre-scan. + bool background_gain_measured_ = false; // The pixels the pre-scan's own frames show to be defective (HotPixelFinder) are measured once, // on the pre-scan that also measures the width, and stay in pixel_mask_ for every pass after it. bool hot_pixels_measured_ = false; diff --git a/rugnux/RugnuxPreScan.cpp b/rugnux/RugnuxPreScan.cpp index 30bc7ae46..c5020f201 100644 --- a/rugnux/RugnuxPreScan.cpp +++ b/rugnux/RugnuxPreScan.cpp @@ -7,6 +7,7 @@ #include "../writer/WriteModel.h" #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/SpotFootprint.h" #include "../image_analysis/bragg_integration/BraggStencil.h" // BRAGG_FOOTPRINT_NSIGMA #include "../image_analysis/hot_pixels/HotPixels.h" @@ -167,7 +168,9 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru const auto goniometer = experiment_.GetGoniometer(); const bool want_hot = !hot_pixels_measured_ && experiment_.IsRotationIndexing() && goniometer && goniometer->IsScanning(); - if (!want_shadow && !want_width && !want_hot) + // The background noise per unit level of a CCD, read off the same frames (BackgroundGain.h). + const bool want_gain = config_.counts_in_adu && !background_gain_measured_; + if (!want_shadow && !want_width && !want_hot && !want_gain) return; // The sample is taken from the frames given - on a rotation sweep its first part (PrepassEnd) - @@ -335,6 +338,9 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru const std::vector ordinals(shadow_set.begin(), shadow_set.end()); const size_t nworkers = std::min(std::max(config_.nthreads, 1), std::min(PRESCAN_MAX_WORKERS, ordinals.size())); + // Per frame, in sample order, so the gain does not depend on how the workers interleaved. + std::vector gain_levels(want_gain ? ordinals.size() : 0); + const size_t gain_frame_stride = (ordinals.size() + background_gain::GAIN_MAX_FRAMES - 1) / background_gain::GAIN_MAX_FRAMES; { std::atomic next{0}; std::vector> futures; @@ -342,6 +348,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru for (size_t t = 0; t < nworkers; t++) futures.emplace_back(std::async(std::launch::async, [&] { std::vector shadow_buffer; + std::vector gain_pixels; JFJochReaderRawImage raw_image; for (size_t i = next.fetch_add(1); i < ordinals.size(); i = next.fetch_add(1)) { const int ordinal = ordinals[i]; @@ -367,11 +374,56 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru msg.original_number = image_idx; finder.AddImage(msg, shadow_buffer); } + if (want_gain && i % gain_frame_stride == 0) { + const auto &img = raw_image.image; + const size_t n = img.GetWidth() * img.GetHeight(); + const uint8_t *ptr = img.GetUncompressedPtr(shadow_buffer); + gain_pixels.resize(n); + switch (img.GetMode()) { + case CompressedImageMode::Int32: + std::memcpy(gain_pixels.data(), ptr, n * sizeof(int32_t)); + break; + case CompressedImageMode::Uint16: + std::copy_n(reinterpret_cast(ptr), n, gain_pixels.begin()); + break; + case CompressedImageMode::Int16: + std::copy_n(reinterpret_cast(ptr), n, gain_pixels.begin()); + break; + default: + continue; // the CCD readers give 16- or 32-bit integers + } + gain_levels[i] = background_gain::MeasureFrame(gain_pixels, img.GetWidth(), + img.GetHeight(), pixel_mask_.GetMask(), + experiment_.GetSaturationLimit()); + } } })); for (auto &f : futures) f.get(); } + // Applied before anything below copies the integration settings (bragg_before_adaptive_), so every + // pass and every fallback integrates with it. At or below 1 - an ADSC or a marCCD, whose pixel + // scatters less than its mean - the photon clip is already looser than its nominal width and is left + // exactly as it is. + if (want_gain) { + background_gain_measured_ = true; + const auto g = background_gain::Gain(gain_levels); + if (g) { + BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings(); + bis.BackgroundPixelGain(std::max(1.0f, *g)); + experiment_.ImportBraggIntegrationSettings(bis); + } + if (!g) + logger.Info("Background gain: not measurable on {} frames; the background clip stays the photon " + "clip", (ordinals.size() + gain_frame_stride - 1) / gain_frame_stride); + else + logger.Info("Background gain: a background pixel scatters {:.2f}x its mean in ADU ({} frames){}", + *g, (ordinals.size() + gain_frame_stride - 1) / gain_frame_stride, + *g > 1.0f ? fmt::format(" => background clip at mean + {:.1f}*sqrt({:.2f}*mean)", + experiment_.GetBraggIntegrationSettings().GetBackgroundClipNSigma(), *g) + : std::string(" - at or below 1, the photon clip is kept")); + } + const auto spot_measurement = [&] { prescan_mapping = std::make_unique(prescan_x, prescan_mask, *azint_geometry_); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index f370f5be5..07cd8b7d0 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -3030,6 +3030,7 @@ static int RunRugnux(int argc, char **argv) { config.detect_beam_stop = detect_beam_stop; config.beam_center_check = beam_center_check; config.adaptive_integration_radius = adaptive_integration_radius; + config.counts_in_adu = input_is_smv || input_is_marccd; config.rotation_postrefine_geometry = rotation_postrefine_geometry; // Measure the spot budget from the data unless the user pinned it. config.measure_spot_budget = !max_spot_count_override.has_value(); diff --git a/tests/BackgroundGainTest.cpp b/tests/BackgroundGainTest.cpp new file mode 100644 index 000000000..ab9ac9db9 --- /dev/null +++ b/tests/BackgroundGainTest.cpp @@ -0,0 +1,75 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include + +#include +#include + +#include "../image_analysis/bragg_integration/BackgroundGain.h" +#include "../common/BraggIntegrationSettings.h" +#include "../common/JFJochException.h" + +namespace { + +constexpr size_t W = 512, H = 512; + +// A smooth background of 100..2000 photons per pixel, read out as bias + gain * photons + read noise, +// with a grid of bright 5x5 "Bragg" spots and a masked strip. +std::vector Frame(double gain, uint32_t seed, std::vector &mask) { + std::mt19937 rng(seed); + std::normal_distribution read_noise(0.0, 3.0); + std::vector img(W * H); + mask.assign(W * H, 0); + for (size_t y = 0; y < H; y++) + for (size_t x = 0; x < W; x++) { + const double photons = 100.0 + 1900.0 * (x + y) / static_cast(W + H); + std::poisson_distribution p(photons); + double v = (gain > 1.0 ? 20.0 : 0.0) + gain * p(rng) + (gain > 1.0 ? read_noise(rng) : 0.0); + if (x % 32 >= 14 && x % 32 < 19 && y % 32 >= 14 && y % 32 < 19) + v += 50.0 * gain * photons; + img[y * W + x] = static_cast(std::lround(v)); + if (y >= 200 && y < 210) + mask[y * W + x] = 1; + } + return img; +} + +} + +TEST_CASE("BackgroundGain_CCD_gain4") { + std::vector frames; + std::vector mask; + for (uint32_t seed = 1; seed <= 3; seed++) { + const auto img = Frame(4.0, seed, mask); + frames.push_back(background_gain::MeasureFrame(img, W, H, mask, 1 << 30)); + } + const auto g = background_gain::Gain(frames); + REQUIRE(g.has_value()); + CHECK_THAT(*g, Catch::Matchers::WithinAbs(4.0, 0.25)); +} + +TEST_CASE("BackgroundGain_photon_counting") { + std::vector mask; + const auto img = Frame(1.0, 7, mask); + const auto g = background_gain::Gain({background_gain::MeasureFrame(img, W, H, mask, 1 << 30)}); + REQUIRE(g.has_value()); + CHECK_THAT(*g, Catch::Matchers::WithinAbs(1.0, 0.1)); +} + +TEST_CASE("BackgroundGain_saturated_and_empty") { + std::vector mask; + const auto img = Frame(4.0, 11, mask); + // Everything at or above the saturation limit: nothing to measure. + CHECK_FALSE(background_gain::Gain({background_gain::MeasureFrame(img, W, H, mask, 1)}).has_value()); + CHECK_FALSE(background_gain::Gain({}).has_value()); +} + +TEST_CASE("BackgroundGain_settings") { + BraggIntegrationSettings s; + CHECK(s.GetBackgroundPixelGain() == 1.0f); + s.BackgroundPixelGain(4.5f); + CHECK(s.GetBackgroundPixelGain() == 4.5f); + CHECK_THROWS_AS(s.BackgroundPixelGain(0.5f), JFJochException); +} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index a8e0468d1..14af2bec2 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -79,6 +79,7 @@ ADD_EXECUTABLE(jfjoch_test BraggIntegrationEngineGPUTest.cpp BraggIntegrationEngineCompressedImageTest.cpp BraggIntegrationEngineCPUTest.cpp + BackgroundGainTest.cpp BraggStencilTest.cpp SpotFootprintTest.cpp LossyFilterTest.cpp