From 711f902aa72a8abd9fc6fdfb9379103da0f4f9fb Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 23 Oct 2026 16:44:09 +0200 Subject: [PATCH] One gain mechanism: frames analysed in their own units, the gain applied where a photon count is assumed The gain of an input in detector units (SMV, marCCD, XTC, the XFEL formats) is measured once, before anything reads the frames (RugnuxMeasureGain: photon-peak spacing and background variance per unit mean), or given with --gain, and handed to the integration and the spot finding by one function (ApplyRugnuxPixelGain). The rotation pre-scan no longer measures it, and the stills reader no longer converts frames of a detector whose photons are not resolved: SACLA MPCCD, Rayonix and CCD frames are integrated in ADU at the gain, rotation and stills alike. Where a photon count was assumed, it is now taken at the gain (bit-identical at gain 1): - integration: the background clip's one-photon floor, the sigma floor (one photon), the profile fit's background floor (0.01 photon), and a still's overload kept as a lower bound, which was left in ADU while every other reflection left the integrator in photons - the cause of the earlier loss of the engine route on clipped CCD stills (merged low resolution inflated by the gain, the rocking-width fit's bandwidth collapsing); CPU and CUDA, Metal through the shared CPU code; - the adaptive spot threshold (Poisson arm and its one-photon floor), CPU and GPU; - the pre-scan spot-width measurement (strong-spot SNR, the ring clip) behind the stills bandwidth and the rotation signal radius; - the hot/dead pixel tests, CPU and GPU; the empty-aperture ring clip. BackgroundPixelGain is gone: one pixel_gain. A detector whose background resolves single photons (photon peaks: JUNGFRAU in keV, AGIPD) is a photon counter and its frames are still counted in whole photons as read: analysed in keV at the gain instead, a photon-counting still set lost 0.008 R-free (fewer spots, a noisier beam-centre trial) - the rounding is the detector's own discrimination. --gain sets the value, the background still says which kind. Validation, against rc178: photon-counting and g<=1 rotation p.hkl md5-identical (lyso, thau, cytc, myob, insu, a small-molecule set, 5reo, 6cee, 6i3j, 2xfw, 6s1u, 8v4j); CCD g>1 at or within noise of rc178 (5uth R_meas 0.113 -> 0.108). Stills R-free (0-cycle REFMAC): MPCCD sets 0.184 -> 0.185 and 0.199 -> 0.196, Rayonix identical, photon-resolving sets back on their rc178 numbers. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu (cherry picked from commit da8dbfc8147354b63e375f2d9521dffb80938b6d) --- common/BraggIntegrationSettings.cpp | 11 -- common/BraggIntegrationSettings.h | 26 ++-- docs/CPU_DATA_ANALYSIS_INTEGRATION.md | 18 ++- docs/RUGNUX_ADVANCED.md | 2 +- docs/RUGNUX_FORMATS.md | 13 +- .../bragg_integration/BackgroundAperture.cpp | 2 +- .../BraggIntegrationEngine.cpp | 7 +- .../BraggIntegrationEngine.h | 5 +- .../BraggIntegrationEngineCPU.cpp | 12 +- .../BraggIntegrationEngineGPU.cu | 10 +- .../bragg_integration/SpotWidth.cpp | 25 +-- image_analysis/bragg_integration/SpotWidth.h | 10 +- image_analysis/hot_pixels/HotPixels.cpp | 11 +- image_analysis/hot_pixels/HotPixels.h | 4 +- image_analysis/hot_pixels/HotPixelsGPU.cu | 2 +- image_analysis/hot_pixels/HotPixelsGPU.h | 3 +- .../spot_finding/AdaptiveSpotFinderCPU.cpp | 6 +- .../spot_finding/AdaptiveSpotFinderGPU.cu | 6 +- .../spot_finding/AdaptiveThreshold.h | 12 +- .../spot_finding/SpotFindingSettings.h | 4 + rugnux/Rugnux.h | 9 -- rugnux/RugnuxDefaults.cpp | 42 ++++++ rugnux/RugnuxDefaults.h | 18 +++ rugnux/RugnuxPreScan.cpp | 94 +----------- rugnux/rugnux_cli.cpp | 142 ++++++++---------- tests/AdaptiveThresholdTest.cpp | 24 ++- tests/BackgroundGainTest.cpp | 8 +- tests/BraggIntegrationEngineCPUTest.cpp | 71 +++++++++ tests/SpotTilesTest.cpp | 2 +- viewer/JFJochProcessController.cpp | 9 +- 30 files changed, 341 insertions(+), 267 deletions(-) diff --git a/common/BraggIntegrationSettings.cpp b/common/BraggIntegrationSettings.cpp index 71494dc0b..4843482a8 100644 --- a/common/BraggIntegrationSettings.cpp +++ b/common/BraggIntegrationSettings.cpp @@ -166,17 +166,6 @@ 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::PixelGain(float input) { check_finite("Pixel gain", input); check_min("Pixel gain", input, 1.0); diff --git a/common/BraggIntegrationSettings.h b/common/BraggIntegrationSettings.h index 18425d85b..fb42d0a3c 100644 --- a/common/BraggIntegrationSettings.h +++ b/common/BraggIntegrationSettings.h @@ -97,18 +97,18 @@ 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; - // Detector units per photon of the integrated intensities. Above 1 the variances the integrator - // forms - the box sum's, the profile fit's and the background's - are taken at this many ADU per - // photon, so the strong-spot selection and the profile/summation guard compare in photons, and - // each reflection leaves the integrator converted to photons (I, sigma, bkg, var_bkg), which is - // what the merge's own variance model (var_bkg + I) assumes. Set by rugnux for rotation sweeps of a - // CCD from the same measurement as bkg_pixel_gain; 1 for a detector that counts photons. + // Detector units per photon of the frames being integrated (ADU/photon): 1 for a detector that counts + // photons, which is what every HDF5 and CBF input is taken to be. The frames are integrated in the + // units they are stored in, and the gain enters only where a photon count is assumed: every variance + // the integrator forms - the box sum's, the profile fit's and the background's, a pixel's variance being + // g * value in ADU^2 - so that the strong-spot selection and the profile/summation guard compare in + // photons; the background clip at mean + n*sqrt(g*mean), with its one-photon floor; and each reflection + // leaves the integrator converted to photons (I, sigma, bkg, var_bkg), which is what the merge's own + // variance model (var_bkg + I) assumes. On a Rigaku Saturn a pixel scatters 4-5x its mean: the photon + // clip would cut the background's upper tail at about 2 sigma and add a positive offset to every + // reflection. Set by rugnux from the gain it measures on the frames (BackgroundGain.h) or is given + // (--gain); a gain below 1 - a CCD read out binned, whose pixels scatter less than their mean through + // the point spread, not through a gain - is taken as 1. float 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 @@ -187,7 +187,6 @@ public: BraggIntegrationSettings& Integrator(IntegratorMode input); BraggIntegrationSettings& BackgroundTrimFraction(float input); BraggIntegrationSettings& BackgroundClipNSigma(float input); - BraggIntegrationSettings& BackgroundPixelGain(float input); BraggIntegrationSettings& PixelGain(float input); BraggIntegrationSettings& BackgroundRadialCorrection(std::optional input); BraggIntegrationSettings& MaxHKL(std::optional input); @@ -209,7 +208,6 @@ public: [[nodiscard]] float GetMinimumSigmaInRegardsToI() const; [[nodiscard]] float GetBackgroundTrimFraction() const; [[nodiscard]] float GetBackgroundClipNSigma() const; - [[nodiscard]] float GetBackgroundPixelGain() const; [[nodiscard]] float GetPixelGain() const; // Unset = auto (gate per image on the smooth-ice score); see bkg_radial_correction. [[nodiscard]] std::optional GetBackgroundRadialCorrection() const; diff --git a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md index f2eb05d06..8d195964f 100644 --- a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md +++ b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md @@ -112,6 +112,20 @@ What a circular $r_1$ loses is flux, and that loss is **not** a function of reso ### 9.2 Box summation (seed and fallback) +**Detector units and the gain.** Frames are integrated in the units the file stores them in — photons +for a photon-counting detector (and for a detector whose background resolves single photons, whose +frames are read as whole photons), ADU for a CCD or an analogue XFEL detector whose file states no gain — with +$g$ detector units to a photon ($g = 1$ for photon counting; measured on a few frames, or given with +`--gain`; a $g$ below 1, a CCD read out binned whose pixels scatter less than their mean through the +point spread, is taken as 1). The gain enters only where a photon count is assumed: a pixel's variance +is $g$ times its value, so every variance below is taken at $g$; the background clip is +$\hat b + n\sqrt{g\,\max(\hat b, g)}$; the variance floors are one photon ($g$ ADU for $\sigma$, $0.01g$ +for the fit's background term); and every reflection leaves the integrator divided by $g$ — $I$, $\sigma$ +and the background in photons — which is what the merge's variance model assumes. Integrated this way a +frame in ADU gives the reflections its photons would, the same stencil, clip and profile, rotation and +stills alike. The spot finder's threshold, the spot-width measurement and the hot- and dead-pixel tests +take their photon terms at the same $g$. + Let: - $S = \sum I(x,y)$ over signal pixels, - $n_S$ = number of valid signal pixels, @@ -123,9 +137,9 @@ $ \hat{b} = \frac{B}{n_B},\qquad \hat{I} = S - n_S \hat{b}, $ -with a Poisson-like uncertainty $\sigma(\hat{I})=\max\!\big(1,\ r_\sigma\hat{I},\ \sqrt{S + n_S^2\,\mathrm{var}(\hat{b})}\big)$, i.e. $\sqrt{S}$ floored both at 1 count (pixel values are photon counts) and at a small fraction $r_\sigma$ of the intensity. The second term under the root is the **uncertainty of the background estimate itself**: $\hat b$ is measured from a finite number of ring pixels, $\mathrm{var}(\hat b)=\hat b/n_B$, and it is subtracted $n_S$ times over, so it enters squared. Omitting it understates the **variance** by $1+n_S/n_B$ — 1.11 with the shipped circular stencil ($n_S = 45$, $n_B = 408$) — and so understates $\sigma$ by up to $\sqrt{1+n_S/n_B} \approx 1.05$, a bound attained on background-limited (weak) reflections and falling towards 1 on strong ones, where $S$ dominates; with an elongated ring $n_B$ grows with resolution, so the factor is no longer one number for a run. The same term is carried into the profile fit (§9.3), where it adds $\big(\sum P/v \,\big/ \sum P^2/v\big)^2\,\mathrm{var}(\hat b)$ — the square of $\partial I/\partial\hat b$ for that fit; $n_B$ is the count of pixels behind the *final* background value, so a clip or trim that discards ring pixels raises it. A box sum is accepted as “observed” only if all signal pixels were valid and $n_B$ exceeds a minimum — it measures what is in the disk with no model of what should be there, so it cannot renormalise a disk it has lost pixels out of. The profile modes can, and do (§9.3). This box sum is the classical estimator; it is used directly with `--integrator boxsum`, and otherwise seeds the profile fit below, where $S$ and $n_S$ then count only the pixels that were actually read. +with a Poisson-like uncertainty $\sigma(\hat{I})=\max\!\big(1,\ r_\sigma\hat{I},\ \sqrt{S + n_S^2\,\mathrm{var}(\hat{b})}\big)$, i.e. $\sqrt{S}$ floored both at 1 photon (the formula is in photons; on frames in ADU it is taken at the gain, above) and at a small fraction $r_\sigma$ of the intensity. The second term under the root is the **uncertainty of the background estimate itself**: $\hat b$ is measured from a finite number of ring pixels, $\mathrm{var}(\hat b)=\hat b/n_B$, and it is subtracted $n_S$ times over, so it enters squared. Omitting it understates the **variance** by $1+n_S/n_B$ — 1.11 with the shipped circular stencil ($n_S = 45$, $n_B = 408$) — and so understates $\sigma$ by up to $\sqrt{1+n_S/n_B} \approx 1.05$, a bound attained on background-limited (weak) reflections and falling towards 1 on strong ones, where $S$ dominates; with an elongated ring $n_B$ grows with resolution, so the factor is no longer one number for a run. The same term is carried into the profile fit (§9.3), where it adds $\big(\sum P/v \,\big/ \sum P^2/v\big)^2\,\mathrm{var}(\hat b)$ — the square of $\partial I/\partial\hat b$ for that fit; $n_B$ is the count of pixels behind the *final* background value, so a clip or trim that discards ring pixels raises it. A box sum is accepted as “observed” only if all signal pixels were valid and $n_B$ exceeds a minimum — it measures what is in the disk with no model of what should be there, so it cannot renormalise a disk it has lost pixels out of. The profile modes can, and do (§9.3). This box sum is the classical estimator; it is used directly with `--integrator boxsum`, and otherwise seeds the profile fit below, where $S$ and $n_S$ then count only the pixels that were actually read. -**High-side clipped background (default on).** Because $\hat{I}=S-n_S\hat{b}$ is a small difference of large numbers for weak reflections, a per-pixel background bias $\delta\hat{b}$ becomes a *fractional* intensity bias $\approx n_S\,\delta\hat{b}/\hat{I}$ that grows as $\hat{I}$ shrinks — worst at the resolution edge. A plain ring mean reads high there, because neighbour-spot wings that survive the signal-disk mask, tails and zingers are one-sided (positive) contaminants. The ring mean is therefore made robust: pixels above $\hat{b}+n\sqrt{\hat{b}}$ are rejected and the mean recomputed, with $n=4$ (`--background-clip`; $n=0$ disables), whatever the bandwidth. A clean Poisson ring is essentially unchanged by the cut (measured false-rejection rate 0.04–0.39 % at $4\sigma$), while a 40-pixel neighbour core at $+100$ counts shifts the estimate by $+0.009$ ct/px. +**High-side clipped background (default on).** Because $\hat{I}=S-n_S\hat{b}$ is a small difference of large numbers for weak reflections, a per-pixel background bias $\delta\hat{b}$ becomes a *fractional* intensity bias $\approx n_S\,\delta\hat{b}/\hat{I}$ that grows as $\hat{I}$ shrinks — worst at the resolution edge. A plain ring mean reads high there, because neighbour-spot wings that survive the signal-disk mask, tails and zingers are one-sided (positive) contaminants. The ring mean is therefore made robust: pixels above $\hat{b}+n\sqrt{\hat{b}}$ (in photons; above) are rejected and the mean recomputed, with $n=4$ (`--background-clip`; $n=0$ disables), whatever the bandwidth. A clean Poisson ring is essentially unchanged by the cut (measured false-rejection rate 0.04–0.39 % at $4\sigma$), while a 40-pixel neighbour core at $+100$ counts shifts the estimate by $+0.009$ ct/px. The clip cuts only the high tail, which matters: the **symmetric** trimmed mean it replaced (drop the lowest and highest fraction $f$ of ring pixels, $f=0.10$; still reachable with `--background-trim`, which switches the clip off) is *not* a consistent estimator of the mean of a right-skewed Poisson sample. It sits $\approx0.1$ ct/px **below** the true mean at every level, and with $n_S = 45$ signal pixels in the $r_1$ disk that under-estimate adds $\approx4.5$ counts to **every** partial — negligible at low resolution, but a large fraction of a partial in the outermost shell. The trim also collapses once contamination exceeds $\approx10\,\%$ of the ring, where the clip does not. Note that removing a positive background bias *lowers* $\langle I/\sigma\rangle$ and *raises* edge $R_\text{meas}$, because both are inflated by information-free counts — so neither may be read as evidence against the change. diff --git a/docs/RUGNUX_ADVANCED.md b/docs/RUGNUX_ADVANCED.md index 308f6a945..70c0c85a3 100644 --- a/docs/RUGNUX_ADVANCED.md +++ b/docs/RUGNUX_ADVANCED.md @@ -493,7 +493,7 @@ Geometry overrides (defaults are taken from the input file; override them to rep | `--detector-distance ` | Detector distance (mm) | | `--wavelength ` | Wavelength (Å) | | `--pixel-size ` | Pixel size (µm), for a file that does not state it (an XTC stream of an unknown camera) | -| `--gain ` | Detector units per photon, for SACLA, CXI, XTC and `--crystfel-geometry` input: the frames are converted to photons with it. Default: for SACLA and CXI input, measured on the background of a few frames; XTC frames stay in detector units ([Serial XFEL inputs](RUGNUX_FORMATS.md#serial-xfel-inputs)) | +| `--gain ` | Detector units per photon, for input in detector units (SACLA, CXI, XTC, SMV, marCCD and `--crystfel-geometry` input). The frames are analysed in their own units and the gain enters where a photon count is assumed; the integrated intensities come out in photons. Default: measured on the background of a few frames ([Serial XFEL inputs](RUGNUX_FORMATS.md#serial-xfel-inputs)) | | `--saturation ` | Pixel value from which a pixel is overloaded, in the units the frames are read in. Replaces the file's limit and the one stills read off the frames | | `--crystfel-geometry ` | Place the panels of an HDF5/CXI input by a CrystFEL geometry file instead of the file's own metadata: module geometry, camera length, mask arrays and bad regions (optional; [Serial XFEL inputs](RUGNUX_FORMATS.md#serial-xfel-inputs)) | | `--module-geometry ` | Where the detector's modules are: an NXmx file with one NXdetector_module per module (what a run writes as `_modules.h5`) or a CrystFEL geometry file whose panels are laid on the image. Overrides the input's own | diff --git a/docs/RUGNUX_FORMATS.md b/docs/RUGNUX_FORMATS.md index 8e3140066..a93f06f8a 100644 --- a/docs/RUGNUX_FORMATS.md +++ b/docs/RUGNUX_FORMATS.md @@ -73,9 +73,14 @@ the file does not state is either measured by rugnux or given on the command lin at **one wavelength** (the mean of the per-shot photon energies where the file states them), the **beam centre is measured** (the file's centre, the background-isotropy measurement and the point-symmetry capture are indexed on a sample of frames and the best is kept), and a SACLA or CXI -file that states **no gain** is converted to photons with the gain its background shows - where the -pixel values resolve single photons, the spacing of their photon peaks, so a baseline the file left in -does not hide it (`--gain` overrides; XTC frames stay in detector units unless `--gain` is given). +file that states **no gain** is read at the gain its background shows. Where the pixel values resolve +single photons (a JUNGFRAU file in keV, an AGIPD file) the gain is the spacing of their photon peaks, so a +baseline the file left in does not hide it, and the detector is a photon counter: its frames are counted in +whole photons as they are read. Otherwise (a CCD, an MPCCD, a Rayonix) the gain is the background's +variance per unit mean and the frames are analysed in the units they are stored in: the gain enters only +where a photon count is assumed (the spot finder's threshold, the background clip, the integrated +intensities' variances, which come out in photons). `--gain` gives the value; the background still says +which of the two the detector is. Float frames analysed in their own units are rounded to whole ones. The saturation level is read off the frames where the file states none. | Format | Read from the file | Measured by rugnux | Given on the command line | @@ -123,7 +128,7 @@ has any of ``, e.g. `@0xfeff`), in the frames' own layout, as stored. With read from each data file - a per-frame mask stored beside the frames. Both may be given more than once, and add to the masks the data files or the geometry file state (the CXI `mask_shared` and per-shot mask, the `/metadata/pixelmask` of a panel-table file, a `--crystfel-geometry` file's masks). A gain the background shows at or below one detector unit -per photon is not applied: such frames stay in detector units. +per photon is taken as one. ### Module geometry diff --git a/image_analysis/bragg_integration/BackgroundAperture.cpp b/image_analysis/bragg_integration/BackgroundAperture.cpp index 4db9a6c94..4861edc8f 100644 --- a/image_analysis/bragg_integration/BackgroundAperture.cpp +++ b/image_analysis/bragg_integration/BackgroundAperture.cpp @@ -111,7 +111,7 @@ void background_aperture::MeasureFrame(ImagePreprocessorBuffer &image, int width break; } const double mean = sum / n; - const double thr = mean + clip_nsigma * std::sqrt(gain * std::max(mean, 1.0)); + const double thr = mean + clip_nsigma * std::sqrt(gain * std::max(mean, static_cast(gain))); double kept = 0.0; int n_kept = 0; for (const auto &p : ring) diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp index 417febd8c..9f63e44a8 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp @@ -117,7 +117,6 @@ 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(); pixel_gain = settings.GetPixelGain(); bkg_trim = bkg_clip_nsigma > 0.0f ? 0.0f : settings.GetBackgroundTrimFraction(); @@ -290,9 +289,11 @@ std::vector BraggIntegrationEngine::Finalize(const std::vector 0) { - refl.I = fr.I + static_cast(fr.saturated * (saturation_limit - fr.bkg)); - refl.sigma = std::sqrt(fr.var_bkg + std::max(refl.I, 1.0f)); + const float I = fr.I + static_cast(fr.saturated * (saturation_limit - fr.bkg)); + refl.I = I / pixel_gain; + refl.sigma = std::sqrt(fr.var_bkg + pixel_gain * std::max(I, pixel_gain)) / pixel_gain; refl.overloaded = false; ++counts.overloads_kept; } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index 8d968748b..3a7d84f9b 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -60,6 +60,7 @@ constexpr double C_CAPTURE = 2.5; // weak-spot radial capture term ( // variance by floor/bkg for every pixel below it, so the previous 1/12 inflated sigma by 1.3x at // 0.05 ct/px and 1.7x at 0.03 - exactly where the weakest high-resolution data live. Digitisation // noise, where a detector has it, is additive on top of the background and does not belong here. +// In photons: on frames in ADU it is taken at the gain, as the background it bounds is. constexpr double PIXEL_VARIANCE_FLOOR = 0.01; // The plug-in signal term of the fit weights may lower the per-pixel variance as well as raise it, // but not below this fraction of the background. Half-wave rectifying it (max(0, I)) instead makes @@ -180,9 +181,9 @@ 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 // ADU per photon (BraggIntegrationSettings::pixel_gain): every variance the engine forms is taken - // at it, and Finalize converts the reflections to photons with it. 1 changes nothing, bit for bit. + // at it, the background clip and the variance floors are a photon wide in it, and Finalize converts + // the reflections to photons with it. 1 changes nothing, bit for bit. float pixel_gain = 1.0f; bool use_ellipse; // radially elongate the per-reflection Gaussian diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 11a1459e0..22bc416c7 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -354,10 +354,12 @@ bool BraggIntegrationEngineCPU::ClipsRing(const BraggBoxSums &s) const { int32_t BraggIntegrationEngineCPU::ClipThreshold(const BraggBoxSums &s) const { // Reject pixels above 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. A pixel - // value is an integer, so it is at or below the threshold exactly when it is at or below its floor. + // the mean, the mean taken at no less than one photon. The gain is 1 for counted photons, so this is + // then exactly the Poisson clip. A pixel value is an integer, so it is at or below the threshold + // exactly when it is at or below its floor. const double bkg = static_cast(s.bkg_sum) / s.n_bkg; - const double thr = bkg + bkg_clip_nsigma * std::sqrt(static_cast(bkg_pixel_gain) * std::max(bkg, 1.0)); + const double g = pixel_gain; + const double thr = bkg + bkg_clip_nsigma * std::sqrt(g * std::max(bkg, g)); return static_cast(std::clamp(std::floor(thr), static_cast(INT32_MIN), static_cast(INT32_MAX))); } @@ -440,7 +442,7 @@ void BraggIntegrationEngineCPU::FinishPassA(PassA &a, size_t i, const Reflection out.n_disk = s.n_disk; out.n_own = s.n_own; const double var_bkg_term = static_cast(s.n_inner_valid) * s.n_inner_valid * out.bkg_var; - out.sigma = 1.0; + out.sigma = g; // a photon if (s.I_sum > 0) { out.sigma = std::max(out.sigma, std::sqrt(g * (static_cast(s.I_sum) + var_bkg_term))); // The centroid of the SIGNAL, not of the disk. Weighting by the raw counts weights by @@ -825,7 +827,7 @@ std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler & fs.inv_2s2t = 0.5f / static_cast(s2t); } fs.bkg = static_cast(rh.bkg); - fs.B = static_cast(std::max(rh.bkg, PIXEL_VARIANCE_FLOOR)); + fs.B = static_cast(std::max(rh.bkg, pixel_gain * PIXEL_VARIANCE_FLOOR)); fs.v_min = static_cast(WEIGHT_VARIANCE_MIN_FRACTION) * fs.B; fs.I0 = static_cast(rh.I); } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index c802cf73b..412562901 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -19,8 +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) - float pixel_gain; // ADU per photon every variance is taken at (see the CPU engine) + float pixel_gain; // ADU per photon every variance and floor is taken at (see the CPU engine) int empirical; // ProfileEmpirical vs ProfileGaussian int use_ellipse; float bw_sigma; @@ -294,7 +293,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p 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((double) p.bkg_pixel_gain * fmax(s_bkg, 1.0)); + s_thr = s_bkg + (double) p.bkg_clip_nsigma * sqrt((double) p.pixel_gain * fmax(s_bkg, (double) p.pixel_gain)); // 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. @@ -419,7 +418,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p // var(I) = Isum + n_inner^2 * bkg/n_bkg_used. Both engines must agree. const double bkg_var = bkg / (double) n_bkg_used; const double var_bkg_term = (double) s_ninner_valid * (double) s_ninner_valid * bkg_var; - double sigma = 1.0; + double sigma = p.pixel_gain; // a photon uint8_t hasobs = 0; double ox = 0.0, oy = 0.0; if (Isum > 0) { sigma = fmax(sigma, sqrt((double) p.pixel_gain * ((double) Isum + var_bkg_term))); @@ -736,7 +735,7 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x } const int Rf = s_Rf, Gf = s_Gf, GfGf = Gf * Gf; - const float B = fmaxf(bkg, (float) PIXEL_VARIANCE_FLOOR); + const float B = fmaxf(bkg, (float) (p.pixel_gain * PIXEL_VARIANCE_FLOOR)); // How much of the expected profile the fit can actually see: s_pvalid over the whole grid against // s_pgrid, the mass that falls on the detector at all (XDS's MINPK), s_pown the same over @@ -1013,7 +1012,6 @@ 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, .pixel_gain = pixel_gain, .empirical = empirical ? 1 : 0, .use_ellipse = use_ellipse ? 1 : 0, diff --git a/image_analysis/bragg_integration/SpotWidth.cpp b/image_analysis/bragg_integration/SpotWidth.cpp index 7c83fc243..24a48a4f1 100644 --- a/image_analysis/bragg_integration/SpotWidth.cpp +++ b/image_analysis/bragg_integration/SpotWidth.cpp @@ -108,15 +108,16 @@ constexpr int R_MOMENT = 9; // the background ring's counts. The background is the ring's mean with anything 5 sigma above it // clipped: a median is biased low on sparse Poisson counts, and a bias here adds the same flux at // every radius of the disk - an extra second moment that does not cancel between the two directions -// once they are weighted by the obliquity. +// once they are weighted by the obliquity. The counts are `gain` to a photon, the clip taken in photons. SpotShape shape_of(const int32_t *centre_px, int stride, int cx, int cy, double mx, double my, const std::vector &ring, const DiffractionGeometry &geometry, - const DetectorModules &modules, int64_t panel, const BraggPanel &pn) { + const DetectorModules &modules, int64_t panel, const BraggPanel &pn, double gain) { SpotShape shape; double sum = 0.0; for (const int32_t v : ring) sum += v; const double clip = sum / static_cast(ring.size()) - + 5.0 * std::sqrt(std::abs(sum / static_cast(ring.size())) + 0.5) + 2.0; + + 5.0 * std::sqrt(gain * (std::abs(sum / static_cast(ring.size())) + 0.5 * gain)) + + 2.0 * gain; sum = 0.0; size_t n = 0; for (const int32_t v : ring) @@ -228,20 +229,20 @@ double median_of(std::vector &v) { void MeasureSpotFluxCurves(const ImagePreprocessorBuffer &image, int width, int height, const DiffractionGeometry &geometry, const DetectorModules &modules, const std::vector &spots, - std::vector &out) { + std::vector &out, float gain) { const int32_t *pixels = image.data(); if (pixels == nullptr) return; const auto centre_px = [&](size_t i) { const Coord c = spots[i].RawCoord(); return pixels + static_cast(std::lround(c.y)) * width + std::lround(c.x); }; - MeasureSpotFluxCurves(centre_px, width, width, height, geometry, modules, spots, out); + MeasureSpotFluxCurves(centre_px, width, width, height, geometry, modules, spots, out, gain); } void MeasureSpotFluxCurves(const std::function ¢re_px, int stride, int width, int height, const DiffractionGeometry &geometry, const DetectorModules &modules, const std::vector &spots, - std::vector &out) { + std::vector &out, float gain) { if (spots.empty()) return; const float beam_x = geometry.GetBeamX_pxl(), beam_y = geometry.GetBeamY_pxl(); @@ -361,8 +362,9 @@ void MeasureSpotFluxCurves(const std::function ¢re_ } } if (core <= 0.0) continue; - const double noise = std::sqrt(core + n_core * std::max(bkg, 0.05) - * (1.0 + static_cast(n_core) / n_ring)); + const double g = gain; + const double noise = std::sqrt(g * (core + n_core * std::max(bkg, 0.05 * g) + * (1.0 + static_cast(n_core) / n_ring))); if (core / noise < SNR_MIN) continue; mx /= core; my /= core; @@ -400,7 +402,8 @@ void MeasureSpotFluxCurves(const std::function ¢re_ if (!(curve.c[R_NORM - 1] > 0.0f) || !(curve.c[R_MAX - 1] > 0.0f)) continue; const float norm = curve.c[R_NORM - 1]; for (float &v : curve.c) v /= norm; - curve.shape = shape_of(centre, stride, cx, cy, mx, my, ring, geometry, modules, cand[k].panel, cand[k].pn); + curve.shape = shape_of(centre, stride, cx, cy, mx, my, ring, geometry, modules, cand[k].panel, cand[k].pn, + gain); out.push_back(curve); } } @@ -410,7 +413,7 @@ SpotTileMisses MeasureSpotsOnTiles(ImagePreprocessorBuffer &image, int width, in const DiffractionGeometry &geometry, const DetectorModules &modules, const std::vector &spots, std::vector &curves, std::vector &footprint, - int half, int footprint_half) { + float gain, int half, int footprint_half) { SpotTileMisses misses; // The pixel each spot rounds to, which both measurements centre on. std::vector tx(spots.size()), ty(spots.size()); @@ -428,7 +431,7 @@ SpotTileMisses MeasureSpotsOnTiles(ImagePreprocessorBuffer &image, int width, in MeasureSpotFluxCurves([&](size_t s) { return tiled ? tiles + s * square(half) + static_cast(half) * (side + 1) : whole + static_cast(ty[s]) * width + tx[s]; - }, tiled ? side : width, width, height, geometry, modules, spots, curves); + }, tiled ? side : width, width, height, geometry, modules, spots, curves, gain); const float beam_x = geometry.GetBeamX_pxl(), beam_y = geometry.GetBeamY_pxl(); const auto tile_of = [&](const int32_t *t, size_t k, size_t s, int h) { diff --git a/image_analysis/bragg_integration/SpotWidth.h b/image_analysis/bragg_integration/SpotWidth.h index 1635968a3..2268d6c4b 100644 --- a/image_analysis/bragg_integration/SpotWidth.h +++ b/image_analysis/bragg_integration/SpotWidth.h @@ -188,11 +188,13 @@ struct BandwidthEstimate { // Every isolated, strong, fully readable spot of one image, as a normalised encircled-flux curve // with its shape, appended to `out`. Called on one worker's own vector, so the pool is independent of the order the // workers finish in. The spots are in the image; `modules` places its panels on the plane the geometry -// describes (empty: the image is the plane), and a spot's aperture has to lie on its own panel. +// describes (empty: the image is the plane), and a spot's aperture has to lie on its own panel. The pixel +// values are `gain` to a photon (BraggIntegrationSettings::pixel_gain): "strong" and the background's clip +// are judged in photons. void MeasureSpotFluxCurves(const ImagePreprocessorBuffer &image, int width, int height, const DiffractionGeometry &geometry, const DetectorModules &modules, const std::vector &spots, - std::vector &out); + std::vector &out, float gain = 1.0f); // The same, with spot i's pixels read about centre_px(i) - the pixel its position rounds to - with // rows `stride` apart: R_BKG_OUT pixels of them every way that the image itself has. The frame is @@ -200,7 +202,7 @@ void MeasureSpotFluxCurves(const ImagePreprocessorBuffer &image, int width, int void MeasureSpotFluxCurves(const std::function ¢re_px, int stride, int width, int height, const DiffractionGeometry &geometry, const DetectorModules &modules, const std::vector &spots, - std::vector &out); + std::vector &out, float gain = 1.0f); // The squares MeasureSpotsOnTiles reads about each spot. The flux curves stay within // spot_width::R_BKG_OUT of their spot; the footprint's window re-centres without a bound, so a spot @@ -224,4 +226,4 @@ SpotTileMisses MeasureSpotsOnTiles(ImagePreprocessorBuffer &image, int width, in const DiffractionGeometry &geometry, const DetectorModules &modules, const std::vector &spots, std::vector &curves, std::vector &footprint, - int half = SPOT_TILE_HALF, int footprint_half = FOOTPRINT_TILE_HALF); + float gain = 1.0f, int half = SPOT_TILE_HALF, int footprint_half = FOOTPRINT_TILE_HALF); diff --git a/image_analysis/hot_pixels/HotPixels.cpp b/image_analysis/hot_pixels/HotPixels.cpp index 79736c80b..2a516db88 100644 --- a/image_analysis/hot_pixels/HotPixels.cpp +++ b/image_analysis/hot_pixels/HotPixels.cpp @@ -78,6 +78,7 @@ HotPixelFinder::HotPixelFinder(const DiffractionExperiment &experiment, const Pi : geometry(experiment.GetDiffractionGeometry()), modules(experiment.GetModules()), axis(experiment.GetGoniometer() ? experiment.GetGoniometer()->GetAxis().Normalize() : Coord(0, 0, 0)), + gain(experiment.GetBraggIntegrationSettings().GetPixelGain()), width(experiment.GetXPixelsNumConv()), height(experiment.GetYPixelsNumConv()), key(std::make_unique_for_overwrite(width * height)), n_lit(std::make_unique_for_overwrite(width * height)), @@ -230,10 +231,10 @@ void HotPixelFinder::AddLevels(const std::vector §or_level, const s for (size_t k = 0; k < nkeys; k++) { const int r = static_cast(k / SECTORS); level[k] = std::max(ring_level[r], sector_level[k]); - const float noise = std::max(std::sqrt(static_cast(std::max(level[k], 0))), ring_spread[r]); + const float noise = std::max(std::sqrt(gain * static_cast(std::max(level[k], 0))), ring_spread[r]); // One rounding for level + nsigma * noise and one for the offset, written out so that every // compiler takes the same two - the device takes them too (HotPixelsGPU.cu). - threshold[k] = std::fma(LIT_NSIGMA, noise, static_cast(level[k])) + LIT_OFFSET; + threshold[k] = std::fma(LIT_NSIGMA, noise, static_cast(level[k])) + LIT_OFFSET * gain; } // Every sum is an integer, so the result does not depend on the order the frames arrive in. @@ -269,7 +270,7 @@ void HotPixelFinder::PrepareDevice() { set_gpu(c); return std::make_unique(key.get(), width * height, key_begin, nrings, SECTORS, HotPixelLevelRules{MIN_SECTOR_PIXELS, MIN_RING_PIXELS, - LIT_NSIGMA, LIT_OFFSET}); + LIT_NSIGMA, LIT_OFFSET * gain, gain}); })); for (auto &b : built) gpus.push_back(b.get()); @@ -422,7 +423,7 @@ HotPixelFinder::Result HotPixelFinder::GetMask(double oscillation_deg, double sp && 2 * sum_value[i] <= key_local_level_sum[key[i]]) { const double mean_local = static_cast(key_local_level_sum[key[i]]) / n_valid(i); const double deficit = mean_local - static_cast(sum_value[i]) / n_valid(i); - if (deficit >= LIT_NSIGMA * std::sqrt(mean_local) + LIT_OFFSET) { + if (deficit >= LIT_NSIGMA * std::sqrt(gain * mean_local) + LIT_OFFSET * gain) { ret.mask[i] = 3; continue; } @@ -431,7 +432,7 @@ HotPixelFinder::Result HotPixelFinder::GetMask(double oscillation_deg, double sp const double mean_level = static_cast(sum_level(i)) / n_valid(i); const double excess = static_cast(sum_value[i] - sum_level(i)) / n_valid(i); if (static_cast(sum_value[i]) >= STRONG_RATIO * static_cast(sum_level(i)) - && excess >= LIT_NSIGMA * std::sqrt(std::max(mean_level, 0.0)) + LIT_OFFSET) + && excess >= LIT_NSIGMA * std::sqrt(gain * std::max(mean_level, 0.0)) + LIT_OFFSET * gain) ret.mask[i] = 1; } }); diff --git a/image_analysis/hot_pixels/HotPixels.h b/image_analysis/hot_pixels/HotPixels.h index 7320979bc..c2357c553 100644 --- a/image_analysis/hot_pixels/HotPixels.h +++ b/image_analysis/hot_pixels/HotPixels.h @@ -69,7 +69,8 @@ public: static constexpr int MIN_SECTOR_PIXELS = 8; // Lit: above level + LIT_NSIGMA * noise + LIT_OFFSET (noise = the larger of sqrt(level) and the // ring's robust spread; the offset keeps a pixel in a ring with no background from being lit by - // one or two photons). + // one or two photons). In photons: on frames in detector units the Poisson noise is sqrt(gain * + // level) and the offset LIT_OFFSET photons, gain of them each (BraggIntegrationSettings::pixel_gain). static constexpr float LIT_NSIGMA = 3.3f; static constexpr float LIT_OFFSET = 2.0f; // How long one reflection can stay on a pixel, degrees of rotation beyond the oscillation: +-3 @@ -119,6 +120,7 @@ private: const DiffractionGeometry geometry; const DetectorModules modules; // where the pixels are on the geometry's plane; empty = as imaged const Coord axis; // rotation axis; zero where there is none + const float gain; // detector units per photon of the frames (BraggIntegrationSettings::pixel_gain) const size_t width, height; int nrings = 0; diff --git a/image_analysis/hot_pixels/HotPixelsGPU.cu b/image_analysis/hot_pixels/HotPixelsGPU.cu index 61dbc3474..8b6e3ea4b 100644 --- a/image_analysis/hot_pixels/HotPixelsGPU.cu +++ b/image_analysis/hot_pixels/HotPixelsGPU.cu @@ -169,7 +169,7 @@ __global__ void levels_kernel(size_t nkeys, int sectors, HotPixelLevelRules rule const int32_t sector_level = sector_ok ? sector_median[k] : 0; const int32_t local_level = sector_ok ? sector_median[k] : ring_level; const int32_t lv = max(ring_level, sector_level); - const float root = sqrtf(static_cast(max(lv, 0))); + const float root = sqrtf(rules.gain * static_cast(max(lv, 0))); const float noise = root < ring_spread ? ring_spread : root; level[k] = lv; threshold[k] = __fadd_rn(__fmaf_rn(rules.lit_nsigma, noise, static_cast(lv)), rules.lit_offset); diff --git a/image_analysis/hot_pixels/HotPixelsGPU.h b/image_analysis/hot_pixels/HotPixelsGPU.h index 7d9cad18b..c490c9438 100644 --- a/image_analysis/hot_pixels/HotPixelsGPU.h +++ b/image_analysis/hot_pixels/HotPixelsGPU.h @@ -16,7 +16,8 @@ struct HotPixelLevelRules { int min_sector_pixels; int min_ring_pixels; float lit_nsigma; - float lit_offset; + float lit_offset; // in the frames' units + float gain; // the frames' units per photon }; // The device half of HotPixelFinder, for frames already preprocessed on the GPU. Everything a frame adds diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp index 3c53ba9af..4e799720f 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp @@ -203,12 +203,14 @@ void AdaptiveSpotFinderCPU::Detect(const ImagePreprocessorBuffer &image, const double g_mean = g_sum / n_total; const double g_sigma = std::sqrt(std::max(0.0, g_sum2 / n_total - g_mean * g_mean)); const float g_thr = adaptive_threshold::RingThreshold(static_cast(g_mean), - static_cast(g_sigma), p, z); + static_cast(g_sigma), p, z, + settings.adu_per_photon); for (size_t b = 0; b < nbins; ++b) ring_thr[b] = (ring_cnt[b] < adaptive_threshold::MIN_RING_PIXELS) ? g_thr - : adaptive_threshold::RingThreshold(ring_mean[b], ring_sigma[b], p, z); + : adaptive_threshold::RingThreshold(ring_mean[b], ring_sigma[b], p, z, + settings.adu_per_photon); // --- Stage C: the ring threshold, intersected with the classic local-box SNR test --- std::fill(ring_bits.begin(), ring_bits.end(), 0); diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu index 7dc5e36b0..5d3f2c576 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu @@ -566,7 +566,8 @@ void AdaptiveSpotFinderGPU::ComputeThresholds(const SpotFindingSettings &setting const double g_mean = g_sum / n_total; const double g_sigma = std::sqrt(std::max(0.0, g_sum2 / n_total - g_mean * g_mean)); const float g_thr = adaptive_threshold::RingThreshold(static_cast(g_mean), - static_cast(g_sigma), p, z); + static_cast(g_sigma), p, z, + settings.adu_per_photon); host_thr.assign(nbins, 0.0f); for (int b = 0; b < nbins; ++b) { @@ -577,7 +578,8 @@ void AdaptiveSpotFinderGPU::ComputeThresholds(const SpotFindingSettings &setting const double m = static_cast(static_cast(host_sum[b])) / host_count[b]; const double var = std::max(0.0, static_cast(host_sum2[b]) / host_count[b] - m * m); host_thr[b] = adaptive_threshold::RingThreshold(static_cast(m), - static_cast(std::sqrt(var)), p, z); + static_cast(std::sqrt(var)), p, z, + settings.adu_per_photon); host_bkg[b] = static_cast(m); } } diff --git a/image_analysis/spot_finding/AdaptiveThreshold.h b/image_analysis/spot_finding/AdaptiveThreshold.h index 745c64d27..6184bd175 100644 --- a/image_analysis/spot_finding/AdaptiveThreshold.h +++ b/image_analysis/spot_finding/AdaptiveThreshold.h @@ -98,10 +98,14 @@ inline float PoissonThreshold(double mu, double p, double z) { // A ring's threshold is background mean + z sigmas, computed two ways and max'd: Poisson significance // (correct where the background is countable) floored by a read-noise-aware Gaussian arm (which alone // survives mean -> 0, where Poisson degenerates to "one photon is significant" and would flood the -// empty high-resolution rings). p, z are the frame-wide operating point (p = E / N_pixels). -inline float RingThreshold(float mean, float sigma, double p, float z) { - const float gauss = mean + z * std::sqrt(sigma * sigma + READ * READ); - const float poisson = PoissonThreshold(static_cast(mean), p, static_cast(z)); +// empty high-resolution rings). p, z are the frame-wide operating point (p = E / N_pixels). mean and +// sigma are in the frame's own units, gain of them to a photon (SpotFindingSettings::adu_per_photon): +// the Poisson arm counts photons and the floor is a photon, so both are taken in photons and handed +// back in the frame's units. A gain of 1 is exactly the photon-counting threshold. +inline float RingThreshold(float mean, float sigma, double p, float z, float gain) { + const float read = READ * gain; + const float gauss = mean + z * std::sqrt(sigma * sigma + read * read); + const float poisson = gain * PoissonThreshold(static_cast(mean / gain), p, static_cast(z)); return std::max(gauss, poisson); } diff --git a/image_analysis/spot_finding/SpotFindingSettings.h b/image_analysis/spot_finding/SpotFindingSettings.h index 281840394..ae4f7c47c 100644 --- a/image_analysis/spot_finding/SpotFindingSettings.h +++ b/image_analysis/spot_finding/SpotFindingSettings.h @@ -105,4 +105,8 @@ struct SpotFindingSettings { // per frame (the threshold's operating point), ~100 for a multi-megapixel detector. bool adaptive_threshold = false; float false_pixels_per_frame = 100.0f; + // Detector units per photon of the frames (BraggIntegrationSettings::pixel_gain, set from the same + // number): the adaptive threshold's Poisson arm and its one-photon floor are taken in photons. + // 1 for a detector that counts photons. + float adu_per_photon = 1.0f; }; diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 962558a1a..651db9faf 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -201,12 +201,6 @@ struct ProcessConfig { // 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 @@ -1077,8 +1071,6 @@ 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; @@ -1356,7 +1348,6 @@ class Rugnux { size_t frames_read = 0; std::optional spots; std::optional defects; - std::optional> gain; }; std::mutex m; std::vector entries; diff --git a/rugnux/RugnuxDefaults.cpp b/rugnux/RugnuxDefaults.cpp index f827c8675..352e11050 100644 --- a/rugnux/RugnuxDefaults.cpp +++ b/rugnux/RugnuxDefaults.cpp @@ -3,6 +3,9 @@ #include "RugnuxDefaults.h" +#include +#include + namespace { // Smoothing range for the per-frame scale G on rotation data (XDS DELPHI-like), in degrees. constexpr double SMOOTH_G_DEG = 5.0; @@ -109,3 +112,42 @@ void ApplyRugnuxIntegrationRadii(DiffractionExperiment &experiment, bool rotatio bool RugnuxDetectIceRings(std::optional file_detect_ice_rings, bool rotation) { return file_detect_ice_rings.value_or(rotation); } + +RugnuxMeasuredGain RugnuxMeasureGain(JFJochReader &reader) { + const auto dataset = reader.GetDataset(); + const size_t n = reader.GetNumberOfImages(); + const size_t frames = std::min(n, background_gain::GAIN_MAX_FRAMES); + std::vector levels; + std::vector histogram(background_gain::PEAK_HISTOGRAM_BINS, 0); + JFJochReaderRawImage raw; + std::vector buffer; + std::vector pixels; + for (size_t i = 0; i < frames; i++) { + if (!reader.ReadRawImage(static_cast(i * n / frames), raw)) + continue; + const size_t w = raw.image.GetWidth(), h = raw.image.GetHeight(); + pixels.resize(w * h); + const uint8_t *ptr = raw.image.GetUncompressedPtr(buffer); + if (raw.image.GetMode() == CompressedImageMode::Int32) + std::memcpy(pixels.data(), ptr, w * h * sizeof(int32_t)); + else if (raw.image.GetMode() == CompressedImageMode::Int16) + std::copy_n(reinterpret_cast(ptr), w * h, pixels.begin()); + else if (raw.image.GetMode() == CompressedImageMode::Uint16) + std::copy_n(reinterpret_cast(ptr), w * h, pixels.begin()); + else + continue; + levels.push_back(background_gain::MeasureFrame(pixels, w, h, dataset->pixel_mask->GetMask(), + dataset->experiment.GetSaturationLimit())); + background_gain::AddToHistogram(pixels, dataset->pixel_mask->GetMask(), + dataset->experiment.GetSaturationLimit(), histogram); + } + return {background_gain::PhotonPeak(histogram), background_gain::Gain(levels)}; +} + +void ApplyRugnuxPixelGain(DiffractionExperiment &experiment, SpotFindingSettings &spot_finding, float adu_per_photon) { + const float g = std::max(1.0f, adu_per_photon); + BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings(); + bis.PixelGain(g); + experiment.ImportBraggIntegrationSettings(bis); + spot_finding.adu_per_photon = g; +} diff --git a/rugnux/RugnuxDefaults.h b/rugnux/RugnuxDefaults.h index a94b687eb..8d53f694b 100644 --- a/rugnux/RugnuxDefaults.h +++ b/rugnux/RugnuxDefaults.h @@ -6,6 +6,8 @@ #include "../common/DiffractionExperiment.h" #include "../common/ScalingSettings.h" #include "../image_analysis/spot_finding/SpotFindingSettings.h" +#include "../image_analysis/bragg_integration/BackgroundGain.h" +#include "../reader/JFJochReader.h" // What "process this dataset with nothing set" means, in one place. // @@ -73,3 +75,19 @@ void ApplyRugnuxIntegrationRadii(DiffractionExperiment &experiment, bool rotatio // and the ice-presence gate keeps it inert on a clean crystal; a serial stills run has too few spots per // image to spend any of them on flagging. A file that states it is followed. bool RugnuxDetectIceRings(std::optional file_detect_ice_rings, bool rotation); + +// The gain of an input in detector units whose file states none, ADU per photon, measured on a few frames +// spread over the run both ways (BackgroundGain.h): the spacing of the photon peaks of the pixel values, +// with the baseline they sit on, where the background resolves single photons, and the background's +// variance per unit mean. +struct RugnuxMeasuredGain { + std::optional peaks; + std::optional variance; +}; +RugnuxMeasuredGain RugnuxMeasureGain(JFJochReader &reader); + +// The gain the frames are analysed at, into everything that takes a photon count from them: the +// integration (BraggIntegrationSettings::pixel_gain, which the defective-pixel tests read too) and the +// spot finding. A gain below 1 - a CCD read out binned, whose pixels scatter less than their mean through +// the point spread, not through a gain - is taken as 1. +void ApplyRugnuxPixelGain(DiffractionExperiment &experiment, SpotFindingSettings &spot_finding, float adu_per_photon); diff --git a/rugnux/RugnuxPreScan.cpp b/rugnux/RugnuxPreScan.cpp index e84061b90..a112a9d6c 100644 --- a/rugnux/RugnuxPreScan.cpp +++ b/rugnux/RugnuxPreScan.cpp @@ -6,7 +6,6 @@ #include "../image_analysis/structure_refinement/ModelValidation.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/BackgroundAperture.h" #include "../image_analysis/bragg_integration/RingBackground.h" #include "../image_analysis/bragg_integration/SpotFootprint.h" @@ -203,9 +202,7 @@ 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())); - // 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) + if (!want_shadow && !want_width && !want_hot) return; // The sample is taken from the frames given - on a rotation sweep its first part (PrepassEnd) - @@ -263,7 +260,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru } const bool find_shadow = want_shadow && !shadow_memo; - // The spot measurement, the defective pixels and the background gain are kept the same way, under + // The spot measurement and the defective pixels are kept the same way, under // the same key and everything else they read: the settings the spots are found with, the // azimuthal binning the adaptive finder thresholds on, and whether the shadow is detected first // (the defective pixels are compared with it masked). A copy of the run that pre-scans again on @@ -297,7 +294,6 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru } const bool measure_width = want_width && !(measured && measured->spots); const bool measure_hot = want_hot && !(measured && measured->defects); - const bool measure_gain = want_gain && !(measured && measured->gain); std::optional finder; if (find_shadow) @@ -370,13 +366,13 @@ 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. + // to the ring the run starts from, so where the noise wins that ring is what is kept. The ring is + // clipped at the gain the run integrates at. 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; + const float aperture_gain = aperture_bis.GetPixelGain(); // The spots of one image, and what is taken off them. const auto find_spots = [&](PreScanWorker &w, const CompressedImage &image, int image_idx, @@ -425,7 +421,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru if (for_width) MeasureSpotsOnTiles(*w.preprocessed, static_cast(prescan_x.GetXPixelsNumConv()), static_cast(prescan_x.GetYPixelsNumConv()), prescan_x.GetDiffractionGeometry(), - prescan_x.GetModules(), spots, curves, footprint); + prescan_x.GetModules(), spots, curves, footprint, aperture_gain); if (for_width && !aperture_r3.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. @@ -474,10 +470,6 @@ 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(PrescanWorkers(), ordinals.size())); - // Per frame, in sample order, so the gain does not depend on how the workers interleaved. - std::vector gain_levels(measure_gain ? ordinals.size() : 0); - std::vector> gain_histograms(measure_gain ? ordinals.size() : 0); - const size_t gain_frame_stride = (ordinals.size() + background_gain::GAIN_MAX_FRAMES - 1) / background_gain::GAIN_MAX_FRAMES; // The sample as read here, by position, for the spot measurement and the defective pixels below // (PRESCAN_KEPT_BYTES). @@ -517,12 +509,11 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru }); } - if (find_shadow || measure_gain || keep_frames) { + if (find_shadow || keep_frames) { TimingMark mark(logger, "pre-scan read and projection"); std::atomic next{0}; OnWorkers(nworkers, [&](size_t) { std::vector shadow_buffer; - std::vector gain_pixels; JFJochReaderRawImage scratch; for (size_t i = next.fetch_add(1); i < ordinals.size(); i = next.fetch_add(1)) { const int ordinal = ordinals[i]; @@ -557,34 +548,6 @@ 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 (measure_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); - bool measurable = true; - 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: - measurable = false; // the CCD readers give 16- or 32-bit integers - } - if (measurable) { - gain_levels[i] = background_gain::MeasureFrame(gain_pixels, img.GetWidth(), - img.GetHeight(), pixel_mask_.GetMask(), - experiment_.GetSaturationLimit()); - gain_histograms[i].assign(background_gain::PEAK_HISTOGRAM_BINS, 0); - background_gain::AddToHistogram(gain_pixels, pixel_mask_.GetMask(), - experiment_.GetSaturationLimit(), gain_histograms[i]); - } - } if (keep) { kept_bytes.fetch_add(raw_image.image_buffer.capacity()); is_kept[i] = 1; @@ -596,47 +559,6 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru } else frames_read = measured ? measured->frames_read : shadow_memo->frames_read; - // 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; - if (measure_gain) { - if (!measured) - measured.emplace(); - // The estimate a stills file is converted with (BackgroundGain.h); a CCD's background shows no - // photon peaks, so this is its variance per unit mean. - std::vector histogram(background_gain::PEAK_HISTOGRAM_BINS, 0); - for (const auto &h: gain_histograms) - for (size_t b = 0; b < h.size(); b++) - histogram[b] += h[b]; - const auto estimate = background_gain::Estimate(gain_levels, histogram); - measured->gain = estimate ? std::optional(estimate->gain) : std::nullopt; - } - const std::optional g = *measured->gain; - if (g) { - BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings(); - bis.BackgroundPixelGain(std::max(1.0f, *g)); - // On a rotation sweep the same gain sets the integrated intensities' variances (and their - // unit, photons). Stills keep their sigmas in ADU for now. - if (experiment_.IsRotationIndexing()) - bis.PixelGain(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, - experiment_.IsRotationIndexing() - ? "; intensities and sigmas in photons at this gain" : "") - : std::string(" - at or below 1, the photon clip is kept")); - } - const auto spot_measurement = [&] { TimingMark mark(logger, "pre-scan spot measurement"); engines_built.get(); @@ -767,8 +689,6 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru 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); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 3ef431b71..c64f6c68c 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -34,7 +34,6 @@ #include "../reader/CrystFELGeometry.h" #include "../reader/NXmxModules.h" #include "../reader/JFJochStartMessageReader.h" -#include "../image_analysis/bragg_integration/BackgroundGain.h" #include "BrokerFeed.h" #include "../writer/HDF5NXmx.h" #include "../writer/NXmxModules.h" @@ -161,7 +160,7 @@ void print_usage() { std::cout << std::endl; std::cout << " Detector geometry" << std::endl; - std::cout << " (Serial XFEL stills are read without a geometry file: SACLA Cheetah files are placed by their own /metadata; Cheetah CXI files of a single-panel detector (square frame) are one flat panel with the pixel size, photon energy and camera length they state; LCLS-I XTC streams of a Rayonix take the photon energy and pixel size from the XTC and need --detector-distance, or -C to have it searched. The beam centre is measured, the run is read at the mean photon energy, and SACLA or CXI files that state no gain are converted to photons with the gain their background shows)" << std::endl; + std::cout << " (Serial XFEL stills are read without a geometry file: SACLA Cheetah files are placed by their own /metadata; Cheetah CXI files of a single-panel detector (square frame) are one flat panel with the pixel size, photon energy and camera length they state; LCLS-I XTC streams of a Rayonix take the photon energy and pixel size from the XTC and need --detector-distance, or -C to have it searched. The beam centre is measured, the run is read at the mean photon energy, and frames that state no gain are counted in photons where their background resolves single photons, and otherwise analysed in the units they are stored in at the gain their background shows - see --gain)" << std::endl; std::cout << " --crystfel-geometry Read the HDF5/CXI files with their panels placed by this CrystFEL geometry file: the module geometry (panel ranges, fast/slow directions, corners, pixel size, coffset, the frames' data path), if stated the camera length, and the masks (mask arrays, bad regions, mask_edge_pixels, no_index) are taken from it; photon energy, camera length and - where the geometry names no mask array - masks otherwise come from the data files (Cheetah CXI, SACLA) or the command line, gain and saturation are measured. Tilted panels, dimN layouts, per-panel data and detector shifts are refused" << std::endl; std::cout << std::endl; @@ -277,7 +276,7 @@ void print_usage() { std::cout << " --detector-distance Detector distance (mm). On stills it is held: neither searched nor refined (the cell absorbs any error). Without it, stills with a reference cell (-C) search the distance and refine it together with the cell, the file's value being one hypothesis" << std::endl; std::cout << " --wavelength Wavelength (A)" << std::endl; std::cout << " --pixel-size Pixel size (um)" << std::endl; - std::cout << " --gain Detector units per photon, for a file that states none (SACLA, CXI, XTC and --crystfel-geometry input; the frames are converted to photons with it). Default: for SACLA, CXI and --crystfel-geometry input, measured on a few frames - the single-photon peak of the pixel values where they resolve photons (for a JUNGFRAU file in keV, the photon energy), else the background's variance per unit mean; XTC frames stay in detector units" << std::endl; + std::cout << " --gain Detector units per photon, for input in detector units (SACLA, CXI, XTC, SMV, marCCD and --crystfel-geometry input). Default: measured on a few frames - the spacing of the photon peaks of the pixel values where they resolve single photons (for a JUNGFRAU file in keV, the photon energy), else the background's variance per unit mean. Where the background resolves single photons (SACLA, CXI and --crystfel-geometry input) the frames are counted in whole photons with it as they are read; otherwise they are analysed in their own units, and the gain enters where a photon count is assumed - the spot finder's threshold, the background clip, the spot-width and defective-pixel tests, and the integrated intensities' variances, which come out in photons. At or below 1 it is taken as 1. --gain replaces the value, not the kind: the background still says whether the photons are resolved" << std::endl; std::cout << " --saturation Pixel value from which a pixel is overloaded, in the units the frames are read in (photons where they are converted). Replaces the file's limit and the one stills read off the frames" << std::endl; std::cout << " --rot1 PONI rotation 1 (rad). On stills, --rot1 or --rot2 holds the tilt: it is not refined. Without them stills refine the tilt where frames the fit has not seen agree" << std::endl; std::cout << " --rot2 PONI rotation 2 (rad)" << std::endl; @@ -716,39 +715,6 @@ void configure_offline_output(DiffractionExperiment &experiment, const std::stri namespace { -// The gain of a detector whose file states none - ADU per photon, on a few frames spread over the run - -// so a foreign stills file can be converted to photons as it is read: the spacing of the photon peaks of -// the pixel values where the background resolves photons, with the baseline they sit on, else the -// background's variance per unit mean - the estimate the pre-scan takes for a CCD (BackgroundGain.h). -std::optional MeasureGain(JFJochReader &reader) { - const auto dataset = reader.GetDataset(); - const size_t n = reader.GetNumberOfImages(); - const size_t frames = std::min(n, background_gain::GAIN_MAX_FRAMES); - std::vector levels; - std::vector histogram(background_gain::PEAK_HISTOGRAM_BINS, 0); - JFJochReaderRawImage raw; - std::vector buffer; - std::vector pixels; - for (size_t i = 0; i < frames; i++) { - if (!reader.ReadRawImage(static_cast(i * n / frames), raw)) - continue; - const size_t w = raw.image.GetWidth(), h = raw.image.GetHeight(); - pixels.resize(w * h); - const uint8_t *ptr = raw.image.GetUncompressedPtr(buffer); - if (raw.image.GetMode() == CompressedImageMode::Int32) - std::memcpy(pixels.data(), ptr, w * h * sizeof(int32_t)); - else if (raw.image.GetMode() == CompressedImageMode::Int16) - std::copy_n(reinterpret_cast(ptr), w * h, pixels.begin()); - else - continue; - levels.push_back(background_gain::MeasureFrame(pixels, w, h, dataset->pixel_mask->GetMask(), - dataset->experiment.GetSaturationLimit())); - background_gain::AddToHistogram(pixels, dataset->pixel_mask->GetMask(), - dataset->experiment.GetSaturationLimit(), histogram); - } - return background_gain::Estimate(levels, histogram); -} - std::atomic g_active_process{nullptr}; // --simulate-live: a stored sweep read as the collection it was while being written - image i is @@ -1947,7 +1913,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { } xfel_reader.ExtraMasks(extra_masks); if (input_is_crystfel) { - xfel_reader.ReadFiles(input_files, crystfel_geometry, gain_per_photon); + xfel_reader.ReadFiles(input_files, crystfel_geometry); reader_ptr = &xfel_reader; if (!xfel_reader.HasPhotonEnergy() && !wavelength_A) { logger.Error("Nothing read states the photon energy of {}; give it with --wavelength", input_file); @@ -1956,11 +1922,11 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { logger.Info("Read {} frames from {} files with the CrystFEL geometry {}", xfel_reader.GetNumberOfImages(), input_files.size(), crystfel_geometry); } else if (input_is_panel_table) { - xfel_reader.ReadPanelTable(input_files, gain_per_photon); + xfel_reader.ReadPanelTable(input_files); reader_ptr = &xfel_reader; logger.Info("Read {} frames from {} HDF5 files placed by their /metadata panel table", xfel_reader.GetNumberOfImages(), input_files.size()); } else if (input_is_cxi) { - xfel_reader.ReadCXI(input_files, gain_per_photon); + xfel_reader.ReadCXI(input_files); reader_ptr = &xfel_reader; logger.Info("Read {} frames from {} CXI files", xfel_reader.GetNumberOfImages(), input_files.size()); } else if (const auto refusal = JFJochXFELReader::CXIRefusal(input_file)) { @@ -1974,9 +1940,9 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { && JFJochXTCReader::CanRead(input_file); if (input_is_xtc && !all_inputs(JFJochXTCReader::CanRead)) refuse_mixed(); - if (gain_per_photon && !reader_ptr && !input_is_xtc) { - logger.Error("--gain is for SACLA, CXI, XTC or --crystfel-geometry input; {} is read in the units it " - "states", input_file); + if (gain_per_photon && !reader_ptr && !input_is_xtc && !input_is_smv && !input_is_marccd) { + logger.Error("--gain is for input in detector units (SACLA, CXI, XTC, SMV, marCCD or " + "--crystfel-geometry input); {} is read in the units it states", input_file); exit(EXIT_FAILURE); } if (reader_ptr) { @@ -2016,8 +1982,6 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { xtc_reader.HasWavelength() ? "" : ", the wavelength", options); exit(EXIT_FAILURE); } - if (gain_per_photon) - xtc_reader.ApplyGain(*gain_per_photon); reader_ptr = &xtc_reader; logger.Info("Read {} XTC still images from {} streams", xtc_reader.GetNumberOfImages(), input_files.size()); @@ -2066,40 +2030,63 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { } JFJochReader &reader = *reader_ptr; - // A SACLA or Cheetah CXI file that states no gain is converted to photons with the one its background - // shows, unless --gain gave it. Not an XTC Rayonix: there the conversion cost R-free (one LCLS set, 0.28 -> - // 0.30) where it gained it on the PAL Rayonix, so XTC frames stay in detector units, with the - // background clip widened by the gain the pre-scan measures, as for any CCD. - // A gain at or below one detector unit per photon (a CCD read out with a low gain: one PAL Rayonix - // set lost most of its indexing when its frames were scaled up by it) is left as it is, as for any CCD. - // Float frames are measured at a sixteenth of their unit, as they would otherwise be rounded to whole - // units first: a photon of a few keV in a file in keV is a peak only a few units wide. - if (!gain_per_photon && (input_is_panel_table || input_is_cxi || input_is_crystfel) && !xfel_reader.CountsInPhotons()) { - const float unit = xfel_reader.StoredFloat() ? 1.0f / 16.0f : 1.0f; - if (unit != 1.0f) + // The gain of an input in detector units, which states none (SMV, marCCD, XTC and the XFEL formats): + // --gain, else measured on a few frames spread over the run. Two kinds of detector, told apart by + // their background: + // - one whose pixel values resolve single photons (a JUNGFRAU file in keV, an AGIPD file in ADU) + // counts photons: its frames are read as whole photons, the spacing of the photon peaks to one, as + // a photon-counting detector hands them out. The rounding is the detector's own discrimination - it + // takes the sub-photon noise of the empty pixels out - and analysing such frames in their stored + // units instead, at the gain, cost a photon-counting still set 0.008 R-free (fewer spots found, a + // noisier beam-centre trial). + // - one whose photons are not resolved (a CCD, an MPCCD, a Rayonix) is analysed in the units it + // stores, and the gain - its background's variance per unit mean - enters only where a photon count + // is assumed: the spot finder's threshold, the spot-width and defective-pixel tests and the + // integration's variances, clip and floors, whose results come out in photons + // (BraggIntegrationSettings::pixel_gain). A gain at or below one detector unit per photon (a CCD + // read out binned, whose pixels scatter less than their mean through the point spread) is taken as one. + // A baseline the photon peaks sit on is left in: taken off, the spot finder's threshold falls with the + // background and picks up noise. Float frames are measured at a sixteenth of their unit, as they + // would otherwise be rounded to whole units first: a photon of a few keV in a file in keV is a peak + // only a few units wide; analysed in their own units, they are rounded to whole ones. + // --gain replaces the value measured, not the kind of detector, which its background still says. + std::optional detector_gain; + const bool counts_in_adu = input_is_smv || input_is_marccd || input_is_xtc || input_is_panel_table + || input_is_cxi || input_is_crystfel; + const bool xfel_input = input_is_panel_table || input_is_cxi || input_is_crystfel; + if (counts_in_adu) { + const bool float_frames = xfel_input && xfel_reader.StoredFloat(); + const float unit = float_frames ? 1.0f / 16.0f : 1.0f; + if (float_frames) xfel_reader.ApplyGain(unit); - auto g = MeasureGain(reader); - if (g) { - g->gain *= unit; - g->pedestal *= unit; - } - if (g && g->gain > 1.0f) { - logger.Info("Gain: {:.3f} detector units per photon, measured on the background{}; the frames are " - "converted to photons with it", g->gain, - g->pedestal != 0.0f ? fmt::format(", its zero-photon peak on a baseline of {:.1f}, left in", g->pedestal) - : std::string()); - xfel_reader.ApplyGain(g->gain); - } else { - if (unit != 1.0f) - xfel_reader.ClearGain(); - if (g) - logger.Info("Gain: {:.3f} detector units per photon, measured on the background - at or below 1, " - "the frames stay in detector units", g->gain); - else - logger.Info("Gain: the file states none and the background does not show it; the frames stay in " - "detector units"); - } + const RugnuxMeasuredGain g = RugnuxMeasureGain(reader); + if (float_frames) + xfel_reader.ClearGain(); + if (g.peaks) { + const float peak_gain = gain_per_photon.value_or(g.peaks->gain * unit); + logger.Info("Gain: {:.3f} detector units per photon, {}{}", peak_gain, + gain_per_photon ? "as given; the background shows photon peaks" + : "the spacing of the photon peaks of the background", + g.peaks->pedestal != 0.0f + ? fmt::format(", its zero-photon peak on a baseline of {:.1f}, left in", g.peaks->pedestal * unit) + : std::string()); + if (xfel_input) { + xfel_reader.ApplyGain(peak_gain); + logger.Info("Gain: the detector resolves single photons, so its frames are counted in whole photons as read"); + } else + detector_gain = peak_gain; + } else if (gain_per_photon || g.variance) { + detector_gain = gain_per_photon ? *gain_per_photon : *g.variance; + logger.Info("Gain: {:.3f} detector units per photon, {}", *detector_gain, + gain_per_photon ? "as given" : "the variance per unit mean of the background"); + } else + logger.Info("Gain: the input states none and its background does not show it; its counts are taken as photons"); } + const float adu_per_photon = std::max(1.0f, detector_gain.value_or(1.0f)); // as ApplyRugnuxPixelGain takes it + if (counts_in_adu && !(xfel_input && xfel_reader.CountsInPhotons())) + logger.Info("Gain: the frames are analysed in their own units{}", adu_per_photon > 1.0f + ? fmt::format(", {:.3f} of them to a photon; intensities and sigmas in photons", adu_per_photon) + : std::string(", one to a photon")); const auto dataset = reader.GetDataset(); if (!dataset) { @@ -3467,6 +3454,7 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { spot_settings.high_resolution_limit = d_min_spot_finding; spot_settings.false_pixels_per_frame = false_pixels_per_frame; spot_settings.low_resolution_limit = d_max_spot_finding; + ApplyRugnuxPixelGain(experiment, spot_settings, adu_per_photon); // Validate the assembled spot-finding settings the same way the online receivers do (broker and // receiver call this same function). It enforces the cross-field constraints that per-argument @@ -3503,8 +3491,6 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { && (input_is_cbf || input_is_marccd || input_is_smv || input_is_crystfel || input_is_panel_table || input_is_cxi || input_is_xtc); config.adaptive_integration_radius = adaptive_integration_radius; - config.counts_in_adu = input_is_smv || input_is_marccd || (input_is_xtc && !xtc_reader.CountsInPhotons()) - || ((input_is_crystfel || input_is_panel_table || input_is_cxi) && !xfel_reader.CountsInPhotons()); 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/AdaptiveThresholdTest.cpp b/tests/AdaptiveThresholdTest.cpp index dc5c88257..fa944a9b2 100644 --- a/tests/AdaptiveThresholdTest.cpp +++ b/tests/AdaptiveThresholdTest.cpp @@ -93,19 +93,31 @@ TEST_CASE("AdaptiveThreshold_RingThreshold", "[SpotFinding]") { for (const float mean: {0.0f, 0.5f, 5.0f, 50.0f}) { for (const float sigma: {0.0f, 1.0f, 7.0f}) { const float gauss = mean + z * std::sqrt(sigma * sigma + READ * READ); - CHECK(RingThreshold(mean, sigma, p, z) >= Catch::Approx(gauss).epsilon(1e-6)); + CHECK(RingThreshold(mean, sigma, p, z, 1.0f) >= Catch::Approx(gauss).epsilon(1e-6)); } } - CHECK(RingThreshold(0.0f, 0.0f, p, z) >= z * READ); + CHECK(RingThreshold(0.0f, 0.0f, p, z, 1.0f) >= z * READ); // Non-decreasing in the background mean and in the background scatter. - CHECK(RingThreshold(20.0f, 4.0f, p, z) > RingThreshold(2.0f, 4.0f, p, z)); - CHECK(RingThreshold(5.0f, 9.0f, p, z) > RingThreshold(5.0f, 1.0f, p, z)); + CHECK(RingThreshold(20.0f, 4.0f, p, z, 1.0f) > RingThreshold(2.0f, 4.0f, p, z, 1.0f)); + CHECK(RingThreshold(5.0f, 9.0f, p, z, 1.0f) > RingThreshold(5.0f, 1.0f, p, z, 1.0f)); // Where the background is countable and quiet, Poisson significance is the binding arm: a ring // with mean 1 and no measured scatter must still demand several photons. - CHECK(RingThreshold(1.0f, 0.0f, p, z) > 1.0f + z * READ); + CHECK(RingThreshold(1.0f, 0.0f, p, z, 1.0f) > 1.0f + z * READ); // A ring whose scatter is far above Poisson (flat-field / read excess) is set by the Gaussian arm. - CHECK(RingThreshold(10.0f, 30.0f, p, z) == Catch::Approx(10.0f + z * std::sqrt(900.0f + READ * READ)).epsilon(1e-6)); + CHECK(RingThreshold(10.0f, 30.0f, p, z, 1.0f) == Catch::Approx(10.0f + z * std::sqrt(900.0f + READ * READ)).epsilon(1e-6)); +} + +TEST_CASE("AdaptiveThreshold_RingThreshold_Gain", "[SpotFinding]") { + const double p = 1e-5; + const float z = static_cast(NormalQuantile(1.0 - p)); + // A frame in detector units, g of them to a photon, is thresholded as its photons would be: the same + // threshold, g times over, for the background and scatter g times over. + for (const float g: {3.0f, 12.0f}) + for (const float mean: {0.0f, 0.05f, 1.0f, 5.0f, 300.0f}) + for (const float sigma: {0.0f, 1.0f, 7.0f}) + CHECK(RingThreshold(g * mean, g * sigma, p, z, g) + == Catch::Approx(g * RingThreshold(mean, sigma, p, z, 1.0f)).epsilon(1e-5)); } diff --git a/tests/BackgroundGainTest.cpp b/tests/BackgroundGainTest.cpp index 6234176c3..5322973f0 100644 --- a/tests/BackgroundGainTest.cpp +++ b/tests/BackgroundGainTest.cpp @@ -68,10 +68,10 @@ TEST_CASE("BackgroundGain_saturated_and_empty") { 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); + CHECK(s.GetPixelGain() == 1.0f); + s.PixelGain(4.5f); + CHECK(s.GetPixelGain() == 4.5f); + CHECK_THROWS_AS(s.PixelGain(0.5f), JFJochException); } namespace { diff --git a/tests/BraggIntegrationEngineCPUTest.cpp b/tests/BraggIntegrationEngineCPUTest.cpp index 153ddb7bb..7b68a0014 100644 --- a/tests/BraggIntegrationEngineCPUTest.cpp +++ b/tests/BraggIntegrationEngineCPUTest.cpp @@ -394,3 +394,74 @@ TEST_CASE("BraggIntegrationEngineCPU_PixelGainGivesPhotons", "[Integration][port CHECK(out_ph[0].var_bkg == Catch::Approx(out_adu[0].var_bkg / 4.0f)); } } + +// The same frame, once in photons and once in a detector's units at 7 ADU per photon, integrates to the +// same reflections in photons when the engine is given the gain: the background clip, its one-photon +// floor, the variance floors, the profile fit's weights and a still's overload kept as a lower bound +// are all taken at the gain. A sparse background of 0-2 photons per pixel, where the floors bind, weak +// and strong spots, and one spot clipped at the saturation limit. +TEST_CASE("BraggIntegrationEngineCPU_AduIntegratesAsPhotons", "[Integration][portable]") { + constexpr int32_t g = 7, limit = 400; + DiffractionExperiment experiment(DetJF(2)); + experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(400.0f).BeamY_pxl(400.0f); + const size_t width = experiment.GetXPixelsNum(), npixel = experiment.GetPixelsNum(); + + ImagePreprocessorBuffer photons(npixel); + uint32_t state = 12345; + for (size_t i = 0; i < npixel; ++i) { + state = state * 1664525u + 1013904223u; + const uint32_t u = state >> 24; + photons[i] = u < 200 ? 0 : (u < 245 ? 1 : 2); + } + std::vector predicted; + int n = 0; + for (int cy = 60; cy < 480; cy += 40) + for (int cx = 60; cx < 1000; cx += 40, ++n) { + const double amplitude = n == 7 ? 3000.0 : 2.0 + 15.0 * (n % 23); + for (int dy = -6; dy <= 6; ++dy) + for (int dx = -6; dx <= 6; ++dx) { + int32_t &v = photons[(cy + dy) * width + cx + dx]; + v += static_cast(std::lround(amplitude * std::exp(-(dx * dx + dy * dy) / (2 * 1.3 * 1.3)))); + if (v >= limit) + v = INT32_MAX; + } + Reflection r{}; + r.h = n; r.k = 1; r.l = 1; + r.predicted_x = static_cast(cx); + r.predicted_y = static_cast(cy); + r.d = 2.0f; + r.prescaling_corr = 1.0f; + r.partiality = 1.0f; + predicted.push_back(r); + } + ImagePreprocessorBuffer adu(npixel); + for (size_t i = 0; i < npixel; ++i) + adu[i] = photons[i] == INT32_MAX ? INT32_MAX : g * photons[i]; + + for (const auto mode : {IntegratorMode::BoxSum, IntegratorMode::ProfileGaussian}) { + BraggIntegrationSettings settings; + settings.Integrator(mode); + experiment.ImportBraggIntegrationSettings(settings); + experiment.Detector().SaturationLimit(limit); + BraggIntegrationEngineCPU in_photons(experiment, PixelMask(experiment)); + const auto ref = in_photons.Run(photons, predicted, predicted.size(), 0); + settings.PixelGain(static_cast(g)); + experiment.ImportBraggIntegrationSettings(settings); + experiment.Detector().SaturationLimit(g * limit); + BraggIntegrationEngineCPU in_adu(experiment, PixelMask(experiment)); + const auto out = in_adu.Run(adu, predicted, predicted.size(), 0); + REQUIRE(out.size() == ref.size()); + REQUIRE(ref.size() > 100); + CHECK(in_adu.Counts().overloads_kept == 1); + CHECK(in_photons.Counts().overloads_kept == 1); + for (size_t i = 0; i < ref.size(); ++i) { + CHECK(out[i].h == ref[i].h); + CHECK(out[i].I == Catch::Approx(ref[i].I).epsilon(1e-4).margin(1e-3)); + CHECK(out[i].sigma == Catch::Approx(ref[i].sigma).epsilon(1e-4)); + CHECK(out[i].bkg == Catch::Approx(ref[i].bkg).epsilon(1e-4).margin(1e-5)); + // 1/den - I cancels: float rounding of den shows at 1e-4 of the difference. + CHECK(out[i].var_bkg == Catch::Approx(ref[i].var_bkg).epsilon(1e-3).margin(1e-3)); + CHECK(out[i].overloaded == ref[i].overloaded); + } + } +} diff --git a/tests/SpotTilesTest.cpp b/tests/SpotTilesTest.cpp index 0cb1d0eb2..5696fa911 100644 --- a/tests/SpotTilesTest.cpp +++ b/tests/SpotTilesTest.cpp @@ -92,7 +92,7 @@ SpotTileMisses CheckFrames(std::vector &frames, Image &&image, int half, } MeasureFootprintSpots(f.image.data(), width, height, g.GetBeamX_pxl(), g.GetBeamY_pxl(), sx, sy, f_frame); - const auto misses = MeasureSpotsOnTiles(image(f), width, height, g, DetectorModules(), f.spots, c_tiles, f_tiles, + const auto misses = MeasureSpotsOnTiles(image(f), width, height, g, DetectorModules(), f.spots, c_tiles, f_tiles, 1.0f, half, footprint_half); CHECK(Same(c_frame, c_tiles)); CHECK(Same(f_frame, f_tiles)); diff --git a/viewer/JFJochProcessController.cpp b/viewer/JFJochProcessController.cpp index 5eda0868f..d9fd334fd 100644 --- a/viewer/JFJochProcessController.cpp +++ b/viewer/JFJochProcessController.cpp @@ -8,6 +8,7 @@ #include "../reader/JFJochHDF5Reader.h" #include "../rugnux/ResultReport.h" #include "../rugnux/RugnuxOutputFiles.h" +#include "../rugnux/RugnuxDefaults.h" #include #include @@ -67,16 +68,20 @@ void JFJochProcessController::run_(QString file_path, DiffractionExperiment expe } else if (JFJochMarCCDReader::CanRead(path)) { marccd_reader.ReadFiles(path); reader_ptr = &marccd_reader; - config.counts_in_adu = true; // a CCD: pixel values are ADU, as for rugnux } else if (JFJochSMVReader::CanRead(path)) { smv_reader.ReadFiles(path); reader_ptr = &smv_reader; - config.counts_in_adu = true; } else { hdf5_reader.ReadFile(path); reader_ptr = &hdf5_reader; } JFJochReader &reader = *reader_ptr; + // A CCD's pixel values are ADU: analysed at the gain its background shows, as rugnux does. + if (reader_ptr == &marccd_reader || reader_ptr == &smv_reader) { + const RugnuxMeasuredGain g = RugnuxMeasureGain(reader); + const std::optional gain = g.peaks ? std::optional(g.peaks->gain) : g.variance; + ApplyRugnuxPixelGain(experiment, config.spot_finding, gain.value_or(1.0f)); + } // Seed the live dataset with the experiment so the chart has geometry context; per-image // results are filled in by OnImageProcessed as the run progresses.