From a29c36600f8a59665da8b454aba2b1708be3cad4 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 9 Aug 2026 01:05:31 +0200 Subject: [PATCH] Beam-stop shadow detection, and a low-resolution limit for scaling rugnux finds the beam stop and its holder in a projection of 60 images and marks them in the pixel mask as bit 9 (--detect-beam-stop[=N|off], on by default). Reflections behind the stop are attenuated but not flagged, so they integrate low with a plausible sigma and nothing downstream catches them: the signal-box gate requires 100% valid pixels and shadow pixels are valid, the background clip is high-side only, and the |zeta| cut applies only to the space-group search merge. The detection compares each pixel's background against the typical background at the same radius on two channels. An azimuthal one (the ring median) finds the holder arm, which is a minority of its ring; a radial one (the background just outside) finds the disk, which the ring median cannot see because inside a fully blocked ring the median is the shadow itself. Pixels are pooled over a 5x5 box and tested only where the background has actually been counted, so low-background data no longer masks the whole detector. Recorded reflections are carved back out - a beam stop cannot block a reflection that was measured. Bit 9 belongs to the run that found it, not to the dataset: it is cleared when a run starts, so a mask read back from a file that carries one starts clear. The user mask (bit 8) is left alone. Scaling and merging gain a low-resolution limit, default 50 A (--scaling-low-resolution , 0 removes it), applied per observation before scaling so it also protects the per-frame scale fit and the space-group search. 50 A is the value XDS configurations use; rugnux_vs_xds.py now matches both of XDS's resolution limits instead of only the high one, so the lowest shell is the same shell in the two programs. The viewer draws the detected shadow in coral with a "Show beam stop" switch in the side panel, exposes the low-resolution limit in the settings dock, and offers detection in its processing jobs. Adding an image marker meant giving the reader a MIN_REAL_PXL_VALUE, because several places classify a pixel by range rather than by equality and would otherwise read the new marker as a very negative intensity. Co-Authored-By: Claude Opus 5 (1M context) --- common/ColorScale.cpp | 2 + common/ColorScale.h | 5 +- common/PixelMask.cpp | 14 + common/PixelMask.h | 5 + common/ScalingSettings.cpp | 11 + common/ScalingSettings.h | 7 + docs/CHANGELOG.md | 3 + docs/PIXEL_MASK.md | 5 + docs/RUGNUX.md | 7 + image_analysis/CMakeLists.txt | 2 + image_analysis/beam_stop/ShadowFinder.cpp | 288 ++++++++++-------- image_analysis/beam_stop/ShadowFinder.h | 74 ++--- image_analysis/scale_merge/HKLKey.cpp | 8 +- image_analysis/scale_merge/HKLKey.h | 4 +- image_analysis/scale_merge/Merge.cpp | 13 +- image_analysis/scale_merge/Merge.h | 2 + .../scale_merge/ReindexAmbiguity.cpp | 2 +- .../scale_merge/RotationScaleMerge.cpp | 2 + .../scale_merge/RotationScaleMerge.h | 1 + image_analysis/scale_merge/ScaleOnTheFly.cpp | 4 +- .../scale_merge/StillsPartialityRefine.cpp | 5 +- .../scale_merge/StillsPartialityRefine.h | 1 + preview/PreviewImage.cpp | 7 +- reader/JFJochReaderImage.cpp | 6 +- reader/JFJochReaderImage.h | 8 +- rugnux/Rugnux.cpp | 59 +++- rugnux/Rugnux.h | 9 + rugnux/RugnuxCommandLine.cpp | 12 + rugnux/rugnux_cli.cpp | 35 +++ rugnux_vs_xds.py | 17 +- tests/JFJochReaderTest.cpp | 3 + tests/MergeScaleTest.cpp | 26 +- tests/PixelMaskTest.cpp | 26 +- viewer/JFJochImageReadingWorker.cpp | 2 +- viewer/JFJochViewerSidePanel.cpp | 7 + viewer/JFJochViewerSidePanel.h | 1 + viewer/JFJochViewerWindow.cpp | 2 + .../image_viewer/JFJochDiffractionImage.cpp | 18 +- viewer/image_viewer/JFJochDiffractionImage.h | 1 + viewer/image_viewer/JFJochFollowerImage.cpp | 2 + viewer/image_viewer/JFJochImage.cpp | 1 + viewer/image_viewer/JFJochImage.h | 4 +- viewer/widgets/JFJochViewerSettingsDock.cpp | 16 + viewer/windows/JFJochProcessingJobsWindow.cpp | 29 ++ viewer/windows/JFJochProcessingJobsWindow.h | 2 + 45 files changed, 562 insertions(+), 196 deletions(-) diff --git a/common/ColorScale.cpp b/common/ColorScale.cpp index 65a88fd0..4d6e4085 100644 --- a/common/ColorScale.cpp +++ b/common/ColorScale.cpp @@ -97,6 +97,8 @@ rgb ColorScale::Apply(ColorScaleSpecial input) const { switch (input) { case ColorScaleSpecial::Gap: return gap; + case ColorScaleSpecial::BeamStop: + return beam_stop; default: case ColorScaleSpecial::BadPixel: return bad; diff --git a/common/ColorScale.h b/common/ColorScale.h index cb8c33e8..0daac58b 100644 --- a/common/ColorScale.h +++ b/common/ColorScale.h @@ -29,7 +29,8 @@ enum class ColorScaleEnum : int { enum class ColorScaleSpecial { Gap, - BadPixel + BadPixel, + BeamStop }; class ColorScale { @@ -72,6 +73,8 @@ class ColorScale { rgb gap = {.r = 190, .g = 190, .b = 190}; // Gray rgb bad = {.r = 255, .g = 0, .b = 255}; // Magenta + // Coral: distinct from the gray gap and the magenta bad pixel, and from every colormap + rgb beam_stop = {.r = 255, .g = 127, .b = 80}; static rgb Apply(float input, const std::vector &map); diff --git a/common/PixelMask.cpp b/common/PixelMask.cpp index 2b2dfd81..8b34b2ec 100644 --- a/common/PixelMask.cpp +++ b/common/PixelMask.cpp @@ -218,6 +218,20 @@ void PixelMask::LoadUserMask(const DiffractionExperiment& experiment, const std: "Size of input user mask invalid"); } +void PixelMask::LoadBeamStopMask(const DiffractionExperiment& experiment, const std::vector &in_mask) { + if (in_mask.size() != mask.size()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Size of input beam stop mask invalid"); + LoadMask(in_mask, BeamStopPixelBit); + UpdateRawMask(experiment); +} + +void PixelMask::ClearBeamStopMask(const DiffractionExperiment& experiment) { + for (auto &i: mask) + i &= ~(1u << BeamStopPixelBit); + UpdateRawMask(experiment); +} + void PixelMask::LoadUserMask(const DiffractionExperiment& experiment, const CompressedImage& image) { const size_t width = image.GetWidth(); const size_t height = image.GetHeight(); diff --git a/common/PixelMask.h b/common/PixelMask.h index f151c4d1..90bf8bcd 100644 --- a/common/PixelMask.h +++ b/common/PixelMask.h @@ -32,6 +32,7 @@ public: constexpr static const uint8_t ErrorPixelBit = 1; constexpr static const uint8_t NoisyPixelBit = 4; constexpr static const uint8_t UserMaskedPixelBit = 8; + constexpr static const uint8_t BeamStopPixelBit = 9; constexpr static const uint8_t ChipGapPixelBit = 31; constexpr static const uint8_t ModuleEdgePixelBit = 30; @@ -43,6 +44,10 @@ public: void CalcEdgePixels(const DiffractionExperiment& experiment); void LoadUserMask(const DiffractionExperiment& experiment, const std::vector& mask); void LoadUserMask(const DiffractionExperiment& experiment, const CompressedImage& image); + void LoadBeamStopMask(const DiffractionExperiment& experiment, const std::vector& mask); + // The beam-stop shadow belongs to the run that found it, not to the dataset, so a mask read back + // from a file that carries one starts clear. The user mask (bit 8) is deliberately left alone. + void ClearBeamStopMask(const DiffractionExperiment& experiment); void LoadDECTRISBadPixelMask(const std::vector& mask); void LoadDarkBadPixelMask(const DiffractionExperiment& experiment, const std::vector& mask); void LoadDetectorBadPixelMask(const DiffractionExperiment& experiment, const JFCalibration *calib); diff --git a/common/ScalingSettings.cpp b/common/ScalingSettings.cpp index fb7d27f8..773ed9fa 100644 --- a/common/ScalingSettings.cpp +++ b/common/ScalingSettings.cpp @@ -24,6 +24,13 @@ ScalingSettings& ScalingSettings::HighResolutionLimit_A(std::optional li return *this; } +ScalingSettings& ScalingSettings::LowResolutionLimit_A(std::optional limit) { + if (limit.has_value() && limit.value() <= 0.0) + throw JFJochException(JFJochExceptionCategory::InputParameterBelowMin, "Low resolution limit must be positive"); + low_resolution_limit_A = limit; + return *this; +} + bool ScalingSettings::GetMergeFriedel() const { return merge_friedel; } @@ -41,6 +48,10 @@ std::optional ScalingSettings::GetHighResolutionLimit_A() const { return high_resolution_limit_A; } +std::optional ScalingSettings::GetLowResolutionLimit_A() const { + return low_resolution_limit_A; +} + double ScalingSettings::GetMinMosaicity() const { return 0.001; } diff --git a/common/ScalingSettings.h b/common/ScalingSettings.h index c5beadf0..56f41159 100644 --- a/common/ScalingSettings.h +++ b/common/ScalingSettings.h @@ -16,6 +16,11 @@ class ScalingSettings { bool merge_friedel = true; std::optional high_resolution_limit_A; + // Low-resolution limit for scaling and merging. Reflections coarser than this are behind or + // beside the beam stop and are measured on a background the stop has eaten into; past 50 A a + // large share of them come out negative. 50 A is what XDS's own configurations use, so keeping + // it here is also what makes the two comparable at the coarse end. + std::optional low_resolution_limit_A = 50.0; std::optional wedge_for_scaling; std::optional forced_mosaicity; // diagnostic: fix the scaling mosaicity (deg) instead of the per-image seed double min_partiality = 0.02; @@ -107,6 +112,7 @@ public: ScalingSettings& MergeFriedel(bool input); ScalingSettings& HighResolutionLimit_A(double limit); ScalingSettings& HighResolutionLimit_A(std::optional limit); // nullopt clears the limit + ScalingSettings& LowResolutionLimit_A(std::optional limit); // nullopt clears the limit ScalingSettings& MinPartiality(double min_partiality); ScalingSettings& ForcedMosaicity(std::optional input); ScalingSettings& CaptureUncertaintyCoeff(double input); @@ -143,6 +149,7 @@ public: [[nodiscard]] bool GetMergeFriedel() const; [[nodiscard]] std::optional GetHighResolutionLimit_A() const; + [[nodiscard]] std::optional GetLowResolutionLimit_A() const; [[nodiscard]] double GetMinPartiality() const; [[nodiscard]] std::optional GetForcedMosaicity() const; diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index d391fdf0..d5d4ee96 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,6 +3,9 @@ ### 1.0.0-rc.161 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. +* rugnux: New **beam-stop shadow detection**, **on by default** (`--detect-beam-stop[=N|off]`), finds the beam stop and its holder in a projection of N images (default 60) and adds them to the pixel mask as bit 9, which is cleared at the start of every run. +* Viewer: the detected beam-stop shadow is drawn in coral, with a "Show beam stop" switch in the side panel. +* rugnux: Scaling and merging now apply a **low-resolution limit** of 50 Å (`--scaling-low-resolution `, 0 removes it). * Spot finding: Self-calibrating **adaptive detection** added, now the **default** for rotation as well as stills; a fused GPU engine runs it together with azimuthal integration in one image pass. * Spot finding: Serial stills pick `--min-pix-per-spot` per image, connected components run on the GPU, and detection is configurable over the API. * Spot finding: `--spot-sigma` now defaults to 4.0 (was 3.0) and `--max-spots` to 1000 spots per image (was 250). diff --git a/docs/PIXEL_MASK.md b/docs/PIXEL_MASK.md index 143f0fd5..033809ca 100644 --- a/docs/PIXEL_MASK.md +++ b/docs/PIXEL_MASK.md @@ -16,6 +16,11 @@ Bit 4 - noisy pixel (for PSI JUNGFRAU: pixel pedestal G0 RMS is over threshold, Bit 8 - user defined mask +Bit 9 - beam stop shadow (found by `rugnux --detect-beam-stop`, on by default; see [rugnux](RUGNUX.md)). +Unlike the other bits this one belongs to the run that found it, not to the dataset: rugnux clears it +at the start of every run, so a mask read back from a file that carries one starts clear. The user +mask (bit 8) is left alone. + Bit 30 - module edge (only for PSI systems) Bit 31 - chip edge interpolated pixel (multipixel) diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index ddf8f9f0..c49fa9c7 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -250,6 +250,12 @@ Calibration (`--mode calibration`): | `--calibrant ` | Powder standard: `lab6` \| `agbh` \| `ceo2` \| `si` \| `ice` (default `lab6`, case-insensitive) | | `--calibration ` | How the rings are measured: `rings` \| `spots` (default `rings`; see above). `rings` defaults `--azim-phi-bins` to 32 | +Detector mask: + +| Option | Description | +| --- | --- | +| `--detect-beam-stop[=N\|off]` | Find the beam stop and its holder in a projection of N images and add them to the pixel mask as bit 9, so nothing shadowed by them is integrated. **On by default** (60 images); `=off` disables. Reflections behind the stop are attenuated but not flagged, so they integrate low with a plausible sigma and no existing rejection catches them | + Spot finding: | Option | Description | @@ -316,6 +322,7 @@ Scaling and merging: | `--capture-uncertainty ` | rot3d: systematic sigma on under-captured fulls, ~num·(1−captured_fraction)·I (default: 1.0 for rotation, 0 otherwise) | | `--min-captured-fraction ` | rot3d: drop a combined full whose rocking curve was captured below this fraction — edge-of-sweep truncated fulls (default: 0.7 for rotation, 0 otherwise; 0 = off) | | `--scaling-high-resolution ` | High-resolution limit for scaling, Å — manual override (default: no limit; disables the automatic cutoff below) | +| `--scaling-low-resolution ` | Low-resolution limit for scaling and merging, Å (default: 50, the value XDS configurations use; 0 removes the limit). Reflections coarser than this sit behind or beside the beam stop and are measured on a background it has eaten into | | `--resolution-cutoff ` | Automatic high-resolution cutoff for the written reflections and reported shells: `cc-logistic` \| `off` (default: `cc-logistic`; ignored when `--scaling-high-resolution` is set) | | `--resolution-cc-target ` | CC1/2 target defining the `cc-logistic` fall-off (default: 0.30) | | `--resolution-shells ` | Number of resolution shells in the reported statistics table (default: 10) | diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index 7056c0d1..0ce15825 100644 --- a/image_analysis/CMakeLists.txt +++ b/image_analysis/CMakeLists.txt @@ -26,6 +26,8 @@ ADD_LIBRARY(JFJochImageAnalysis STATIC IndexAndRefine.h dark_mask_analysis/DarkMaskAnalysis.cpp dark_mask_analysis/DarkMaskAnalysis.h + beam_stop/ShadowFinder.cpp + beam_stop/ShadowFinder.h rotation_indexer/RotationIndexer.cpp rotation_indexer/RotationIndexer.h WriteReflections.cpp diff --git a/image_analysis/beam_stop/ShadowFinder.cpp b/image_analysis/beam_stop/ShadowFinder.cpp index c2044337..4687ce0c 100644 --- a/image_analysis/beam_stop/ShadowFinder.cpp +++ b/image_analysis/beam_stop/ShadowFinder.cpp @@ -11,15 +11,38 @@ #include "../../common/JFJochException.h" -// --------------------------------------------------------------------------------- -// Small binary-image helpers on a width*height frame stored row-major as char (0/1). -// All run once, at GetMask() time. The BFS forms keep them O(pixels) rather than -// O(pixels * radius), so a radius-14 dilation is still a single sweep. -// --------------------------------------------------------------------------------- +// A pixel is shadow when its background is below this fraction of the background it is +// compared against. +constexpr float SHADOW_RATIO = 0.35f; + +// The boundary grows outward into partially shadowed pixels down to this fraction, but no +// further than PENUMBRA_MAX_PX from the core. +constexpr float PENUMBRA_RATIO = 0.72f; +constexpr int PENUMBRA_MAX_PX = 14; + +// Bridge module gaps and small breaks that the holder arm crosses. +constexpr int BRIDGE_PX = 6; + +// A pixel whose maximum reaches this recorded a real reflection and is never masked - a +// beam stop cannot block a reflection that was measured. +constexpr int64_t MIN_REFLECTION = 25; + +// Counts the background must have accumulated over the frames and the pooled pixels before +// a dip in it is believable. Below this a Poisson hole is indistinguishable from a shadow, +// and testing anyway masks whole detectors on low-background data. +constexpr double MIN_EXPECTED_COUNTS = 60; + +// Side of the box the background is pooled over before testing, and how far out the radial +// comparison looks for unshadowed background. +constexpr int POOL_PX = 5; +constexpr float ENVELOPE_MM = 6.0f; + +// Binary-image helpers on a width*height frame stored row-major as char (0/1). All run once, +// at GetMask() time; the BFS forms keep them O(pixels) rather than O(pixels * radius). namespace { // 8-connected dilation by `r` pixels (Chebyshev), via a multi-source BFS. -std::vector Dilate(const std::vector &in, int W, int H, int r) { +std::vector dilate(const std::vector &in, int W, int H, int r) { if (r <= 0) return in; std::vector dist(in.size(), -1); @@ -46,12 +69,12 @@ std::vector Dilate(const std::vector &in, int W, int H, int r) { return out; } -// Erosion by `r` = dilation of the complement (image border counts as outside). -std::vector Erode(const std::vector &in, int W, int H, int r) { +// Erosion by `r` = dilation of the complement; outside the frame counts as complement. +std::vector erode(const std::vector &in, int W, int H, int r) { std::vector comp(in.size()); for (size_t i = 0; i < in.size(); i++) comp[i] = !in[i]; - const auto grown = Dilate(comp, W, H, r); + const auto grown = dilate(comp, W, H, r); std::vector out(in.size()); for (size_t i = 0; i < out.size(); i++) out[i] = !grown[i]; @@ -59,13 +82,11 @@ std::vector Erode(const std::vector &in, int W, int H, int r) { } // Pixels of `passable` reachable from any of `seeds` (8-connected flood). -std::vector Flood(const std::vector &passable, int W, int H, const std::vector &seeds) { +std::vector flood(const std::vector &passable, int W, int H, const std::vector &seeds) { std::vector visited(passable.size(), 0); std::queue q; for (const int s : seeds) - if (s >= 0 && s < static_cast(passable.size()) && passable[s] && !visited[s]) { - visited[s] = 1; q.push(s); - } + if (passable[s] && !visited[s]) { visited[s] = 1; q.push(s); } while (!q.empty()) { const int i = q.front(); q.pop(); const int y = i / W, x = i % W; @@ -82,7 +103,7 @@ std::vector Flood(const std::vector &passable, int W, int H, const s } // Fill holes: background not reachable from the image border becomes region. -std::vector FillHoles(const std::vector ®ion, int W, int H) { +std::vector fill_holes(const std::vector ®ion, int W, int H) { std::vector bg_visited(region.size(), 0); std::queue q; auto push = [&](int i) { if (!region[i] && !bg_visited[i]) { bg_visited[i] = 1; q.push(i); } }; @@ -107,9 +128,36 @@ std::vector FillHoles(const std::vector ®ion, int W, int H) { return out; } +// Sum of `in` over the k x k box centred on each pixel, zero outside the frame. +std::vector box_sum(const std::vector &in, int W, int H, int k) { + const int half = k / 2; + std::vector row(in.size(), 0.0), out(in.size(), 0.0); + for (int y = 0; y < H; y++) { + double s = 0; + for (int x = 0; x <= std::min(half, W - 1); x++) + s += in[y * W + x]; + for (int x = 0; x < W; x++) { + row[y * W + x] = s; + if (x + half + 1 < W) s += in[y * W + x + half + 1]; + if (x - half >= 0) s -= in[y * W + x - half]; + } + } + for (int x = 0; x < W; x++) { + double s = 0; + for (int y = 0; y <= std::min(half, H - 1); y++) + s += row[y * W + x]; + for (int y = 0; y < H; y++) { + out[y * W + x] = s; + if (y + half + 1 < H) s += row[(y + half + 1) * W + x]; + if (y - half >= 0) s -= row[(y - half) * W + x]; + } + } + return out; +} + // Median of `values` per integer radius, over the pixels flagged in `use`. -std::vector RingMedian(const std::vector &values, const std::vector &use, - const std::vector &radius, int max_radius) { +std::vector ring_median(const std::vector &values, const std::vector &use, + const std::vector &radius, int max_radius) { std::vector> bins(max_radius + 1); for (size_t i = 0; i < values.size(); i++) if (use[i]) @@ -126,39 +174,42 @@ std::vector RingMedian(const std::vector &values, const std::vecto return median; } -// Fraction of each integer-radius ring that is flagged in `blocked`. -std::vector RingFraction(const std::vector &blocked, const std::vector &radius, int max_radius) { - std::vector num(max_radius + 1, 0), den(max_radius + 1, 0); - for (size_t i = 0; i < blocked.size(); i++) { - den[radius[i]]++; - if (blocked[i]) num[radius[i]]++; +// Largest baseline over [r, r + win] - the background just outside radius r. +std::vector outer_envelope(const std::vector &baseline, int win) { + const int n = static_cast(baseline.size()); + std::vector out(n, 0.0f); + for (int r = 0; r < n; r++) { + float v = baseline[r]; + for (int k = 1; k <= win; k++) + v = std::max(v, baseline[std::min(r + k, n - 1)]); + out[r] = v; } - std::vector frac(max_radius + 1, 0.0f); - for (int r = 0; r <= max_radius; r++) - frac[r] = den[r] ? static_cast(num[r]) / static_cast(den[r]) : 0.0f; - return frac; + return out; } } // namespace -// --------------------------------------------------------------------------------- - -ShadowFinder::ShadowFinder(const DiffractionExperiment &experiment, ShadowFinderSettings in_settings) +ShadowFinder::ShadowFinder(const DiffractionExperiment &experiment, const PixelMask &mask) : width(static_cast(experiment.GetXPixelsNumConv())), height(static_cast(experiment.GetYPixelsNumConv())), beam_x(experiment.GetBeamX_pxl()), beam_y(experiment.GetBeamY_pxl()), - settings(in_settings), + envelope_px(std::max(4, static_cast(std::lround(ENVELOPE_MM / experiment.GetPixelSize_mm())))), + pixel_mask(mask.GetMask(experiment)), max_value(static_cast(width) * height, 0), sum_value(static_cast(width) * height, 0), - valid_count(static_cast(width) * height, 0) {} + valid_count(static_cast(width) * height, 0) { + if (pixel_mask.size() != max_value.size()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "ShadowFinder: pixel mask does not match the detector"); +} template void ShadowFinder::Add(const T *ptr) { // The pixel type's sentinel extreme marks "no data" (module gap / masked): the - // preprocessor/writer stores INT*_MIN for signed and UINT*_MAX for unsigned. For - // signed types the opposite extreme (INT*_MAX) is a genuine saturated value and is - // kept, so a saturated reflection still registers as bright. + // preprocessor/writer stores INT*_MIN for signed and UINT*_MAX for unsigned. For signed + // types the opposite extreme is a genuine saturated value and is kept, so a saturated + // reflection still registers as bright. T masked; if constexpr (std::is_signed_v) masked = std::numeric_limits::min(); @@ -170,7 +221,7 @@ void ShadowFinder::Add(const T *ptr) { const T v = ptr[i]; if (v == masked) continue; - const int32_t vi = static_cast(v); + const int64_t vi = static_cast(v); if (valid_count[i] == 0 || vi > max_value[i]) max_value[i] = vi; sum_value[i] += vi; @@ -207,70 +258,96 @@ std::vector ShadowFinder::GetMask() const { std::unique_lock ul(m); const int W = width, H = height; - const int N = W * H; - const ShadowFinderSettings &S = settings; + const int n_pixels = W * H; - std::vector mask(N, 0); + std::vector mask(n_pixels, 0); if (frames == 0) return mask; - // --- mean projection, per-pixel validity and radius from the beam centre --- - std::vector mean(N, 0.0f); - std::vector valid(N, 0); - std::vector radius(N, 0); + // mean projection, usable pixels and radius from the beam centre + std::vector mean(n_pixels, 0.0f); + std::vector valid(n_pixels, 0); + std::vector radius(n_pixels, 0); int max_radius = 0; for (int y = 0; y < H; y++) for (int x = 0; x < W; x++) { const int i = y * W + x; - if (valid_count[i] > 0) { + if (valid_count[i] > 0 && pixel_mask[i] == 0) { mean[i] = static_cast(static_cast(sum_value[i]) / valid_count[i]); valid[i] = 1; } - const double dx = x - beam_x, dy = y - beam_y; - const int r = static_cast(std::lround(std::sqrt(dx * dx + dy * dy))); - radius[i] = r; - if (r > max_radius) max_radius = r; + const float dx = x - beam_x, dy = y - beam_y; + radius[i] = static_cast(std::lround(std::sqrt(dx * dx + dy * dy))); + max_radius = std::max(max_radius, radius[i]); } - // --- robust radial baseline; iterate to keep the shadow out of its own baseline --- - std::vector ratio(N, 1.0f); - std::vector excluded(N, 0); + // Pool the background over a small box before testing it. A background of a fraction of + // a count per pixel per frame gives no single pixel enough counts to tell a shadow from + // a Poisson hole; the stop and its arm are wider than the box, so pooling costs no + // resolution that matters and multiplies the statistics by the pixels in the box. + std::vector num(n_pixels), den(n_pixels); + for (int i = 0; i < n_pixels; i++) { + num[i] = valid[i] ? mean[i] : 0.0; + den[i] = valid[i] ? 1.0 : 0.0; + } + const auto pooled_sum = box_sum(num, W, H, POOL_PX); + const auto pooled_count = box_sum(den, W, H, POOL_PX); + std::vector pooled(n_pixels, 0.0f); + for (int i = 0; i < n_pixels; i++) + if (pooled_count[i] > 0) + pooled[i] = static_cast(pooled_sum[i] / pooled_count[i]); + + // Azimuthal comparison: the median of the ring, iterated so the shadow stays out of the + // baseline it is measured against. + std::vector ratio(n_pixels, 1.0f); + std::vector excluded(n_pixels, 0); + std::vector baseline; for (int iter = 0; iter < 3; iter++) { - std::vector use(N); - for (int i = 0; i < N; i++) + std::vector use(n_pixels); + for (int i = 0; i < n_pixels; i++) use[i] = valid[i] && !excluded[i]; - const auto baseline = RingMedian(mean, use, radius, max_radius); - for (int i = 0; i < N; i++) + baseline = ring_median(pooled, use, radius, max_radius); + for (int i = 0; i < n_pixels; i++) if (valid[i]) - ratio[i] = mean[i] / std::max(baseline[radius[i]], 1e-6f); - for (int i = 0; i < N; i++) - excluded[i] = valid[i] && ratio[i] < S.shadow_ratio; + ratio[i] = pooled[i] / std::max(baseline[radius[i]], 1e-6f); + for (int i = 0; i < n_pixels; i++) + excluded[i] = valid[i] && ratio[i] < SHADOW_RATIO; } - // --- shadow core: low-ratio pixels connected to the beam centre (bridging gaps) --- - std::vector low(N); - for (int i = 0; i < N; i++) - low[i] = valid[i] && ratio[i] < S.shadow_ratio; + // Radial comparison: the background just outside this radius. The disk blocks its rings + // completely, so their median is the shadow itself and only this comparison sees it. + const auto envelope = outer_envelope(baseline, envelope_px); + std::vector ratio_radial(n_pixels, 1.0f); + for (int i = 0; i < n_pixels; i++) + if (valid[i]) + ratio_radial[i] = pooled[i] / std::max(envelope[radius[i]], 1e-6f); - const std::vector grown = Dilate(low, W, H, S.bridge_px); - std::vector seeds; // a small disk at the beam centre - for (int y = 0; y < H; y++) - for (int x = 0; x < W; x++) { - const double dx = x - beam_x, dy = y - beam_y; - if (dx * dx + dy * dy < 4.0 * 4.0) - seeds.push_back(y * W + x); - } - const std::vector connected = Flood(grown, W, H, seeds); - std::vector core(N); - for (int i = 0; i < N; i++) - core[i] = low[i] && connected[i]; + // A dip counts only where the background it is compared against was actually counted. + std::vector low(n_pixels, 0); + for (int i = 0; i < n_pixels; i++) { + if (!valid[i]) + continue; + const double counted = frames * pooled_count[i]; + low[i] = (ratio[i] < SHADOW_RATIO && baseline[radius[i]] * counted >= MIN_EXPECTED_COUNTS) + || (ratio_radial[i] < SHADOW_RATIO && envelope[radius[i]] * counted >= MIN_EXPECTED_COUNTS); + } - // --- real reflections: any pixel that recorded signal is never masked. Require a - // small cluster so a single-frame zinger does not count as a reflection. --- - std::vector lit(N, 0); - for (int i = 0; i < N; i++) - lit[i] = (valid_count[i] > 0) && (max_value[i] >= static_cast(S.min_reflection)); - std::vector reflection(N, 0); + // The shadow is the low region connected to the beam centre, bridging the gaps it crosses. + const std::vector bridged = dilate(low, W, H, BRIDGE_PX); + std::vector seeds; + for (int i = 0; i < n_pixels; i++) + if (radius[i] < 4) + seeds.push_back(i); + const std::vector connected = flood(bridged, W, H, seeds); + std::vector region(n_pixels); + for (int i = 0; i < n_pixels; i++) + region[i] = low[i] && connected[i]; + + // Recorded reflections. A small cluster is required so a single-frame zinger does not count. + std::vector lit(n_pixels, 0); + for (int i = 0; i < n_pixels; i++) + lit[i] = (valid_count[i] > 0) && (max_value[i] >= MIN_REFLECTION); + std::vector reflection(n_pixels, 0); for (int y = 0; y < H; y++) for (int x = 0; x < W; x++) { const int i = y * W + x; @@ -285,54 +362,23 @@ std::vector ShadowFinder::GetMask() const { reflection[i] = (neighbours >= 2); } - // --- central low-res disk: the fully-blocked region about the beam centre. Sized by - // the azimuthal blocked fraction (a disk blocks ~every azimuth; a thin arm or - // gap does not), and capped just inside the innermost reflection. --- - std::vector blocked(N); - for (int i = 0; i < N; i++) - blocked[i] = (valid_count[i] == 0) || low[i]; - const auto blocked_frac = RingFraction(blocked, radius, max_radius); - - int disk_radius = 0; - { - float head = 0.0f; int head_n = 0; - for (int r = 0; r <= std::min(5, max_radius); r++) { head += blocked_frac[r]; head_n++; } - if (head_n > 0 && head / head_n >= 0.65f) { // the beam centre is behind a disk - disk_radius = max_radius; - for (int r = 1; r <= max_radius; r++) - if (blocked_frac[r] < 0.55f) { disk_radius = r; break; } - } - } - int reflection_radius = max_radius + 1; // innermost reflection (ignore the very centre) - for (int i = 0; i < N; i++) - if (reflection[i] && radius[i] > 12 && radius[i] < reflection_radius) - reflection_radius = radius[i]; - if (disk_radius > reflection_radius - 4) - disk_radius = reflection_radius - 4; - if (disk_radius < 0) - disk_radius = 0; - - // --- assemble: core + disk, grow the soft penumbra, round, fill the disk interior --- - std::vector region(N); - for (int i = 0; i < N; i++) - region[i] = core[i] || (disk_radius > 0 && radius[i] < disk_radius); - - const std::vector near = Dilate(region, W, H, S.penumbra_max_px); - for (int i = 0; i < N; i++) - if (near[i] && valid[i] && ratio[i] < S.penumbra_ratio) + // Grow the soft boundary, round it and fill the disk interior. + const std::vector penumbra = dilate(region, W, H, PENUMBRA_MAX_PX); + for (int i = 0; i < n_pixels; i++) + if (penumbra[i] && valid[i] && std::min(ratio[i], ratio_radial[i]) < PENUMBRA_RATIO) region[i] = 1; - region = Erode(Dilate(region, W, H, 2), W, H, 2); // close: round the boundary - region = FillHoles(region, W, H); + region = erode(dilate(region, W, H, 2), W, H, 2); + region = fill_holes(region, W, H); - // Expose recorded reflections - done last, with no fill afterwards, so a spot the - // geometry still covered is given back rather than re-enclosed. - const std::vector reflection_grown = Dilate(reflection, W, H, 1); - for (int i = 0; i < N; i++) + // Expose recorded reflections - done last, with no fill afterwards, so a spot the shadow + // still covered is given back rather than re-enclosed. + const std::vector reflection_grown = dilate(reflection, W, H, 1); + for (int i = 0; i < n_pixels; i++) if (reflection_grown[i]) region[i] = 0; - for (int i = 0; i < N; i++) + for (int i = 0; i < n_pixels; i++) mask[i] = region[i] ? 1 : 0; return mask; } diff --git a/image_analysis/beam_stop/ShadowFinder.h b/image_analysis/beam_stop/ShadowFinder.h index 14ee4adc..ca7aea34 100644 --- a/image_analysis/beam_stop/ShadowFinder.h +++ b/image_analysis/beam_stop/ShadowFinder.h @@ -9,67 +9,51 @@ #include "../../common/CompressedImage.h" #include "../../common/DiffractionExperiment.h" -#include "../../common/JFJochMessages.h" // DataMessage +#include "../../common/JFJochMessages.h" +#include "../../common/PixelMask.h" -// Tunable parameters for ShadowFinder. Plain struct with sensible defaults; when the -// finder is wired into the workflow these can move onto DiffractionExperiment the way -// DarkMaskSettings does. See SHADOW_FINDER.md for what each one does. -struct ShadowFinderSettings { - // A pixel is "shadow core" when its mean is below this fraction of the typical - // (azimuthal-median) background at the same radius. - float shadow_ratio = 0.35f; - - // The soft boundary grows outward into partially-shadowed pixels down to this - // fraction of the background, but no further than penumbra_max_px from the core. - float penumbra_ratio = 0.72f; - int penumbra_max_px = 14; - - // Bridge module gaps / small breaks that the holder arm crosses (pixels). - int bridge_px = 6; - - // A pixel whose max-projection reaches this value recorded a real reflection and is - // never masked - a beam stop cannot block a reflection that was measured. This also - // caps the central disk just inside the innermost such reflection. - float min_reflection = 25.0f; -}; - -// Detects the beam-stop shadow (central disk + holder arm) from a small number of -// images, mirroring the accumulate-then-finalize shape of DarkMaskAnalysis: feed frames -// with AddImage(), then read the mask once with GetMask(). The returned mask is in -// converted geometry and is 1 where the beam stop shadows the detector. +// Finds the beam-stop shadow - the central disk and the holder arm - from a set of images, +// mirroring the accumulate-then-finalize shape of DarkMaskAnalysis: feed frames with +// AddImage(), then read the mask once with GetMask(). The mask is in converted geometry +// and is 1 where the beam stop shadows the detector. // -// The shadow is treated as an azimuthal anomaly: a per-radius background baseline is -// robust to the shadow, so a localized dip connected to the beam centre is the beam -// stop. See SHADOW_FINDER.md for the full algorithm and the (deferred) wiring plan. +// The shadow is a place where the background is missing, so it is found by comparing each +// pixel's mean against the typical background at the same radius. Two comparisons are +// needed: an azimuthal one (median over the ring) finds the arm, which is a minority of +// its ring, and a radial one (the background just outside) finds the disk, which is not - +// inside a fully blocked ring the ring median is itself the shadow. +// +// Frames are chosen by the caller; the detection needs enough of them that the background +// is counted rather than guessed (see MIN_EXPECTED_COUNTS in the .cpp). // Thread-safe: AddImage may be called from several worker threads. class ShadowFinder { mutable std::mutex m; - const int width; - const int height; - const double beam_x; - const double beam_y; - const ShadowFinderSettings settings; + const int width; + const int height; + const float beam_x; + const float beam_y; + const int envelope_px; + + std::vector pixel_mask; // pixels already masked carry no background to test // Per-pixel projection over the frames added so far (converted geometry). - std::vector max_value; // maximum over frames - std::vector sum_value; // sum of valid values - std::vector valid_count; // number of frames the pixel carried data + std::vector max_value; + std::vector sum_value; + std::vector valid_count; uint32_t frames = 0; template void Add(const T *ptr); public: - ShadowFinder(const DiffractionExperiment &experiment, ShadowFinderSettings settings = {}); + ShadowFinder(const DiffractionExperiment &experiment, const PixelMask &mask); - // Accumulate one full converted-geometry image into the projection. Gap / masked - // pixels (the pixel type's sentinel extreme) are skipped. `buffer` is scratch space - // for decompression (mirrors DarkMaskAnalysis::AnalyzeImage). + // Accumulate one full converted-geometry image. Gap / masked pixels (the pixel type's + // sentinel extreme) are skipped. `buffer` is scratch space for decompression. void AddImage(const DataMessage &data, std::vector buffer); - // Compute the beam-stop shadow mask (1 = shadow, 0 = keep). Size is the converted - // pixel count. Recomputed from the accumulators on each call - meant to be called - // once at the end; not cheap (see SHADOW_FINDER.md). + // Compute the shadow mask (1 = shadow, 0 = keep), of the converted pixel count. + // Recomputed from the accumulators on each call - meant to be called once at the end. [[nodiscard]] std::vector GetMask() const; [[nodiscard]] uint32_t GetFrameCount() const; diff --git a/image_analysis/scale_merge/HKLKey.cpp b/image_analysis/scale_merge/HKLKey.cpp index 6929b804..ba60ab77 100644 --- a/image_analysis/scale_merge/HKLKey.cpp +++ b/image_analysis/scale_merge/HKLKey.cpp @@ -76,13 +76,15 @@ bool HKLKeyGenerator::IsSystematicallyAbsent(const Reflection &r) const { return IsSystematicallyAbsent(r.h, r.k, r.l); } -bool AcceptReflection(const Reflection &r, std::optional d_min_limit) { +bool AcceptReflection(const Reflection &r, std::optional d_min_limit, std::optional d_max_limit) { if (!std::isfinite(r.I)) return false; if (!std::isfinite(r.d) || r.d <= 0.0f) return false; if (d_min_limit && r.d < d_min_limit) return false; + if (d_max_limit && r.d > d_max_limit) + return false; if (!std::isfinite(r.rlp) || r.rlp == 0.0f) return false; if (!std::isfinite(r.sigma) || r.sigma <= 0.0) @@ -90,13 +92,15 @@ bool AcceptReflection(const Reflection &r, std::optional d_min_limit) { return true; } -bool AcceptReflection(const Reflection &r, double d_min_limit) { +bool AcceptReflection(const Reflection &r, double d_min_limit, double d_max_limit) { if (!std::isfinite(r.I)) return false; if (!std::isfinite(r.d) || r.d <= 0.0f) return false; if (d_min_limit > 0.0 && r.d < d_min_limit) return false; + if (d_max_limit > 0.0 && r.d > d_max_limit) + return false; if (!std::isfinite(r.rlp) || r.rlp == 0.0f) return false; if (!std::isfinite(r.sigma) || r.sigma <= 0.0) diff --git a/image_analysis/scale_merge/HKLKey.h b/image_analysis/scale_merge/HKLKey.h index 70044c74..26fc8bad 100644 --- a/image_analysis/scale_merge/HKLKey.h +++ b/image_analysis/scale_merge/HKLKey.h @@ -48,5 +48,5 @@ HKLKey CanonicalHKL(const Reflection &r, bool merge_friedel, const std::optional HKLKey CanonicalHKL(const MergedReflection &r, bool merge_friedel, const std::optional &sg); HKLKey CanonicalHKL(int32_t h, int32_t k, int32_t l, bool merge_friedel, const std::optional &sg); -bool AcceptReflection(const Reflection &r, std::optional d_min_limit); -bool AcceptReflection(const Reflection &r, double d_min_limit); \ No newline at end of file +bool AcceptReflection(const Reflection &r, std::optional d_min_limit, std::optional d_max_limit); +bool AcceptReflection(const Reflection &r, double d_min_limit, double d_max_limit); \ No newline at end of file diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index aa65734f..4a99ea95 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -37,6 +37,7 @@ MergeOnTheFly::MergeOnTheFly(const DiffractionExperiment &x) scaling_settings(x.GetScalingSettings()), indexing_settings(x.GetIndexingSettings()), high_resolution_limit(scaling_settings.GetHighResolutionLimit_A()), + low_resolution_limit(scaling_settings.GetLowResolutionLimit_A()), // A min-image-CC of 0 (the default) means "no limit": leave the optional // empty so the per-image CC cut is inactive. Otherwise a 0.0 threshold // would silently drop every image with a non-positive per-image CC. @@ -68,7 +69,7 @@ void MergeOnTheFly::AddImage(const IntegrationOutcome &outcome, int64_t image_id if (r.image_scale_corr <= 0.0 || !std::isfinite(r.image_scale_corr)) continue; - if (!AcceptReflection(r, high_resolution_limit)) + if (!AcceptReflection(r, high_resolution_limit, low_resolution_limit)) continue; if (exclude_ice_rings && r.on_ice_ring) continue; @@ -186,7 +187,7 @@ void MergeOnTheFly::RefineErrorModel(const std::vector &outc continue; if (r.image_scale_corr <= 0.0 || !std::isfinite(r.image_scale_corr)) continue; - if (!AcceptReflection(r, high_resolution_limit)) + if (!AcceptReflection(r, high_resolution_limit, low_resolution_limit)) continue; if (exclude_ice_rings && r.on_ice_ring) continue; @@ -407,6 +408,7 @@ std::pair ImageReferenceCC(const std::vector &reflec const std::map &reference, const HKLKeyGenerator &generator, std::optional d_min_limit, + std::optional d_max_limit, double min_partiality) { constexpr size_t MIN_REFLECTIONS = 20; @@ -420,7 +422,7 @@ std::pair ImageReferenceCC(const std::vector &reflec for (const auto &r: reflections) { if (r.on_ice_ring) continue; - if (!AcceptReflection(r, d_min_limit)) + if (!AcceptReflection(r, d_min_limit, d_max_limit)) continue; if (r.partiality < min_partiality) continue; @@ -510,6 +512,7 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m auto d_min_limit_A = d_min_override.has_value() ? d_min_override : scaling_settings.GetHighResolutionLimit_A(); + const auto d_max_limit_A = scaling_settings.GetLowResolutionLimit_A(); std::unordered_map reference_intensities; if (!reference.empty()) { @@ -531,6 +534,8 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m continue; if (d_min_limit_A && m.d < d_min_limit_A) continue; + if (d_max_limit_A && m.d > d_max_limit_A) + continue; d_min = std::min(d_min, m.d); d_max = std::max(d_max, m.d); @@ -607,7 +612,7 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m continue; if (r.image_scale_corr <= 0.0 || !std::isfinite(r.image_scale_corr)) continue; - if (!AcceptReflection(r, d_min_limit_A)) + if (!AcceptReflection(r, d_min_limit_A, d_max_limit_A)) continue; if (r.partiality < min_partiality) continue; diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 32e168fd..7adf10eb 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -98,6 +98,7 @@ class MergeOnTheFly { std::optional reference_cell; std::optional high_resolution_limit; + std::optional low_resolution_limit; std::optional image_cc_limit; // Apply image_cc_limit in Mask(). One flag for the whole engine, not a per-call argument, so the // merge, the error model and MergeStats can never disagree about which images are in. @@ -186,4 +187,5 @@ std::pair ImageReferenceCC(const std::vector &reflec const std::map &reference, const HKLKeyGenerator &generator, std::optional d_min_limit, + std::optional d_max_limit, double min_partiality); diff --git a/image_analysis/scale_merge/ReindexAmbiguity.cpp b/image_analysis/scale_merge/ReindexAmbiguity.cpp index 8e13266d..46753047 100644 --- a/image_analysis/scale_merge/ReindexAmbiguity.cpp +++ b/image_analysis/scale_merge/ReindexAmbiguity.cpp @@ -128,7 +128,7 @@ ReindexAmbiguityResolver::ReindexAmbiguityResolver(const DiffractionExperiment & bool ReindexAmbiguityResolver::Accept(const Reflection &r) const { if (r.on_ice_ring) // ice-contaminated intensity would bias the correlation; keep it out return false; - return AcceptReflection(r, s.GetHighResolutionLimit_A()); + return AcceptReflection(r, s.GetHighResolutionLimit_A(), s.GetLowResolutionLimit_A()); } double ReindexAmbiguityResolver::ReferenceCC(const std::vector &reflections, diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 927f0027..225d8cbe 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -174,6 +174,7 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment, const auto s = x.GetScalingSettings(); min_partiality = s.GetMinPartiality(); d_min_limit = s.GetHighResolutionLimit_A(); + d_max_limit = s.GetLowResolutionLimit_A(); merge_friedel = s.GetMergeFriedel(); capture_uncertainty_coeff = s.GetCaptureUncertaintyCoeff(); min_captured_fraction = s.GetMinCapturedFraction(); @@ -574,6 +575,7 @@ int RotationScaleMerge::ComputeAsuGroups(const HKLKeyGenerator &keygen) { const float d = rawrun_d[r]; // resolution is a per-raw-hkl property (all its partials share d) if (!std::isfinite(d) || d <= 0.0f) continue; if (d_min_limit && d < *d_min_limit) continue; + if (d_max_limit && d > *d_max_limit) continue; key[r] = keygen(rawrun_h[r], rawrun_k[r], rawrun_l[r]).pack(); eligible[r] = 1; } diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 1f9e3a78..15677ee5 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -100,6 +100,7 @@ private: int n_frames = 0; double min_partiality = 0.02; std::optional d_min_limit; + std::optional d_max_limit; bool merge_friedel = true; double capture_uncertainty_coeff = 0.0; double min_captured_fraction = 0.0; diff --git a/image_analysis/scale_merge/ScaleOnTheFly.cpp b/image_analysis/scale_merge/ScaleOnTheFly.cpp index de9db0be..b1dc1887 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.cpp +++ b/image_analysis/scale_merge/ScaleOnTheFly.cpp @@ -88,7 +88,7 @@ bool ScaleOnTheFly::Accept(const Reflection &r) const { if (r.on_ice_ring) // ice-contaminated intensity would drag the per-image scale; keep it out of the fit return false; - return AcceptReflection(r, s.GetHighResolutionLimit_A()); + return AcceptReflection(r, s.GetHighResolutionLimit_A(), s.GetLowResolutionLimit_A()); } void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const { @@ -135,7 +135,7 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const { const auto [cc, cc_n] = ImageReferenceCC(integration_outcome.reflections, reference_data, hkl_key_generator, s.GetHighResolutionLimit_A(), - s.GetMinPartiality()); + s.GetLowResolutionLimit_A(), s.GetMinPartiality()); result.cc = cc; result.cc_n = cc_n; diff --git a/image_analysis/scale_merge/StillsPartialityRefine.cpp b/image_analysis/scale_merge/StillsPartialityRefine.cpp index f01bae07..6396d35e 100644 --- a/image_analysis/scale_merge/StillsPartialityRefine.cpp +++ b/image_analysis/scale_merge/StillsPartialityRefine.cpp @@ -143,6 +143,7 @@ StillsPartialityRefine::StillsPartialityRefine(const DiffractionExperiment &x) : experiment_(x), hkl_key_generator_(x.GetScalingSettings().GetMergeFriedel(), x.GetSpaceGroupNumber().value_or(1)), d_min_limit_(x.GetScalingSettings().GetHighResolutionLimit_A()), + d_max_limit_(x.GetScalingSettings().GetLowResolutionLimit_A()), min_partiality_(x.GetScalingSettings().GetMinPartiality()), bandwidth_sigma_(x.GetBandwidthFWHM().value_or(0.0f) / 2.3548f) {} @@ -171,7 +172,7 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, double m_n = 0.0, m_x = 0.0, m_xx = 0.0, m_y = 0.0, m_xy = 0.0; size_t n_de = 0; for (const Reflection &r: outcome.reflections) { - if (r.on_ice_ring || !AcceptReflection(r, d_min_limit_)) + if (r.on_ice_ring || !AcceptReflection(r, d_min_limit_, d_max_limit_)) continue; if (!std::isfinite(r.I) || !std::isfinite(r.sigma) || r.sigma <= 0.0f) continue; @@ -322,7 +323,7 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, // Refresh it here: it is reported per image and --min-image-cc drops images by it, so it has to be // the CC of the data that is actually merged. const auto [cc, cc_n] = ImageReferenceCC(outcome.reflections, reference, hkl_key_generator_, - d_min_limit_, min_partiality_); + d_min_limit_, d_max_limit_, min_partiality_); // Adopt the refined model only if it correlates with the reference at least as well as the model it // replaces. Rejecting puts the crystal back exactly as it arrived, which is the same state a diff --git a/image_analysis/scale_merge/StillsPartialityRefine.h b/image_analysis/scale_merge/StillsPartialityRefine.h index 4db99990..fbce8c1d 100644 --- a/image_analysis/scale_merge/StillsPartialityRefine.h +++ b/image_analysis/scale_merge/StillsPartialityRefine.h @@ -50,6 +50,7 @@ private: const Settings settings_{}; const HKLKeyGenerator hkl_key_generator_; const std::optional d_min_limit_; + const std::optional d_max_limit_; const double min_partiality_; const float bandwidth_sigma_; diff --git a/preview/PreviewImage.cpp b/preview/PreviewImage.cpp index be38a7bb..b7f94172 100644 --- a/preview/PreviewImage.cpp +++ b/preview/PreviewImage.cpp @@ -331,8 +331,11 @@ void PreviewImage::ConfigurePixel(const std::vector &mask_tmp, size_t | (1u << PixelMask::ChipGapPixelBit) | (1u << PixelMask::ModuleEdgePixelBit); - constexpr uint32_t det_bits = 0xFEFEu; // bits 1-7 and 9-15 - constexpr uint32_t usr_bits = (1u << PixelMask::UserMaskedPixelBit); + // Bits 1-7 and 10-15. The beam-stop shadow (bit 9) is a deliberate exclusion rather than + // a detector defect, so it is shown the same way as the user mask. + constexpr uint32_t det_bits = 0xFCFEu; + constexpr uint32_t usr_bits = (1u << PixelMask::UserMaskedPixelBit) + | (1u << PixelMask::BeamStopPixelBit); for (size_t i = pixel_begin; i < pixel_end; i++) { const auto pixel_val = mask_tmp[i]; diff --git a/reader/JFJochReaderImage.cpp b/reader/JFJochReaderImage.cpp index 109d7413..57428ca6 100644 --- a/reader/JFJochReaderImage.cpp +++ b/reader/JFJochReaderImage.cpp @@ -119,6 +119,8 @@ void JFJochReaderImage::ProcessInputImage(const void *data, size_t npixel, int64 | (1<(i)); @@ -211,6 +213,8 @@ void JFJochReaderImage::AddImage(const JFJochReaderImage &other) { for (size_t i = 0; i < image.size(); i++) { if (image[i] == GAP_PXL_VALUE || other.image[i] == GAP_PXL_VALUE) { image[i] = GAP_PXL_VALUE; + } else if (image[i] == BEAM_STOP_PXL_VALUE || other.image[i] == BEAM_STOP_PXL_VALUE) { + image[i] = BEAM_STOP_PXL_VALUE; } else if (image[i] == ERROR_PXL_VALUE || other.image[i] == ERROR_PXL_VALUE) { image[i] = ERROR_PXL_VALUE; error_pixel.emplace(static_cast(i)); @@ -219,7 +223,7 @@ void JFJochReaderImage::AddImage(const JFJochReaderImage &other) { saturated_pixel.emplace(static_cast(i)); } else { int64_t sum = static_cast(image[i]) + static_cast(other.image[i]); - if (sum <= INT32_MIN + 5) [[unlikely]] { + if (sum < MIN_REAL_PXL_VALUE) [[unlikely]] { image[i] = ERROR_PXL_VALUE; error_pixel.emplace(static_cast(i)); } else if (sum > dataset->experiment.GetSaturationLimit()) [[unlikely]] { diff --git a/reader/JFJochReaderImage.h b/reader/JFJochReaderImage.h index bb7ebbe6..113270b8 100644 --- a/reader/JFJochReaderImage.h +++ b/reader/JFJochReaderImage.h @@ -14,8 +14,14 @@ #include "../common/CrystalLattice.h" #include "../common/Histogram.h" -constexpr static int32_t GAP_PXL_VALUE = INT32_MIN + 1; +// Markers stored in place of an intensity. They occupy the bottom of the int32 range and +// INT32_MAX at the top, so anything in between is a real count. Add a marker at the BOTTOM and move +// MIN_REAL_PXL_VALUE with it - several places classify a pixel by range rather than by equality, +// and they all test against MIN_REAL_PXL_VALUE. constexpr static int32_t ERROR_PXL_VALUE = INT32_MIN; +constexpr static int32_t GAP_PXL_VALUE = INT32_MIN + 1; +constexpr static int32_t BEAM_STOP_PXL_VALUE = INT32_MIN + 2; +constexpr static int32_t MIN_REAL_PXL_VALUE = INT32_MIN + 3; constexpr static int32_t SATURATED_PXL_VALUE = INT32_MAX; struct JFJochReaderRawImage { diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index d2b79e3f..52d1d535 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -24,6 +24,7 @@ #include "../common/time_utc.h" #include "../writer/FileWriter.h" #include "../image_analysis/MXAnalysisWithoutFPGA.h" +#include "../image_analysis/beam_stop/ShadowFinder.h" #include "../image_analysis/IndexAndRefine.h" #include "../image_analysis/geom_refinement/GeometryRefiner.h" #include "../image_analysis/indexing/IndexerThreadPool.h" @@ -79,7 +80,58 @@ namespace { Rugnux::Rugnux(JFJochHDF5Reader &reader, DiffractionExperiment experiment, PixelMask pixel_mask, ProcessConfig config) : reader_(reader), experiment_(std::move(experiment)), - pixel_mask_(std::move(pixel_mask)), config_(std::move(config)) {} + pixel_mask_(std::move(pixel_mask)), config_(std::move(config)) { + // Bit 9 describes where THIS run found the beam stop, so a mask read back from a file that + // already carries one starts clear; the user mask (bit 8) is left as it was loaded. + pixel_mask_.ClearBeamStopMask(experiment_); +} + +void Rugnux::FindBeamStop(int start_image, int images_to_process, int frame_count) { + Logger logger("Rugnux"); + + // Two-pass rotation runs this twice. The shadow does not move, and re-detecting would find + // nothing (its pixels are masked by now) and so clear the mask the first pass established. + const auto ¤t = pixel_mask_.GetMask(); + if (std::any_of(current.begin(), current.end(), + [](uint32_t v) { return (v & (1u << PixelMask::BeamStopPixelBit)) != 0; })) + return; + + const auto sample = select_equally_spaced_image_ordinals(images_to_process, frame_count); + if (sample.empty()) + return; + + ShadowFinder finder(experiment_, pixel_mask_); + std::vector buffer; + for (const int ordinal : sample) { + const int image_idx = start_image + ordinal * config_.stride; + std::shared_ptr img; + try { + img = reader_.GetRawImage(image_idx); + } catch (const std::exception &e) { + logger.Warning("Beam stop detection: failed to load image {}: {}", image_idx, e.what()); + continue; + } + if (!img) continue; + + DataMessage msg{}; + msg.image = img->image; + msg.number = ordinal; + msg.original_number = image_idx; + finder.AddImage(msg, buffer); + } + + if (finder.GetFrameCount() == 0) { + logger.Warning("Beam stop detection: no image could be read. Skipping."); + return; + } + + const auto shadow = finder.GetMask(); + const auto shadowed = std::count(shadow.begin(), shadow.end(), 1u); + pixel_mask_.LoadBeamStopMask(experiment_, shadow); + logger.Info("Beam stop shadow: {} pixels ({:.2f}% of the detector) found in {} images", + shadowed, 100.0 * static_cast(shadowed) / static_cast(shadow.size()), + finder.GetFrameCount()); +} void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_process, RugnuxObserver *observer) { @@ -502,6 +554,11 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b if (full && config_.refine_geometry.has_value()) RefineStillsGeometry(start_image, end_image, images_to_process, observer); + // After any geometry refinement, so the shadow is found about the beam centre actually used, + // and before the azimuthal mapping and the output mask below, which both read pixel_mask_. + if (config_.detect_beam_stop.has_value()) + FindBeamStop(start_image, images_to_process, config_.detect_beam_stop.value()); + AzimuthalIntegrationMapping mapping(experiment_, pixel_mask_); JFJochReceiverPlots plots; diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 91aed5dc..dcf5ff64 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -75,6 +75,11 @@ struct ProcessConfig { // re-indexes with it. The value is the number of strong frames fed to the bundle (--refine-geometry). std::optional refine_geometry; + // Beam-stop shadow detection (--detect-beam-stop). When set, a pre-pass projects this many + // frames, finds where the beam stop and its holder shadow the detector and marks them in the + // pixel mask, so nothing behind them is integrated. The value is the number of frames projected. + std::optional detect_beam_stop; + // 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 @@ -214,6 +219,10 @@ class Rugnux { // frames, bundle-adjust the shared beam/distance/cell from the strongest ones, and apply the result // to experiment_ so the main pass re-indexes with it. No-op (leaves experiment_ unchanged) if it // cannot run. Stills only; rotation has its own two-pass. + // Beam-stop shadow pre-pass (config_.detect_beam_stop): project a spread sample of frames and + // add the shadow of the stop and its holder to the pixel mask. + void FindBeamStop(int start_image, int images_to_process, int frame_count); + void RefineStillsGeometry(int start_image, int end_image, int images_to_process, RugnuxObserver *observer); diff --git a/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index 05a10cab..368f248f 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/rugnux/RugnuxCommandLine.cpp @@ -79,6 +79,12 @@ std::string RugnuxCommandLine(const ProcessConfig &config, add("-e", std::to_string(config.end_image)); if (config.stride != 1) add("-t", std::to_string(config.stride)); + // Beam-stop detection is ON by default in the CLI, so emit only a deviation from that. getopt takes + // an optional argument only when attached (=N), never as a separate token. + if (!config.detect_beam_stop.has_value()) + args.emplace_back("--detect-beam-stop=off"); + else if (*config.detect_beam_stop != 60) + args.push_back("--detect-beam-stop=" + std::to_string(*config.detect_beam_stop)); if (calibration) { if (!calibrant_name.empty()) @@ -191,6 +197,12 @@ std::string RugnuxCommandLine(const ProcessConfig &config, // no CLI equivalent - the CLI always writes the .mtz/.cif when merging.) if (config.write_process_h5) args.emplace_back("--write-process-h5"); + // The low-resolution limit is on by default at the same value in both, so emit it only when + // it differs - including 0, which is how the CLI spells "no limit". + const auto d_max = sc.GetLowResolutionLimit_A(); + const auto d_max_default = ScalingSettings().GetLowResolutionLimit_A(); + if (d_max != d_max_default) + args.push_back("--scaling-low-resolution=" + num(d_max.value_or(0.0))); } else { args.emplace_back("--no-merge"); } diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index ee0e63d5..e7a39f82 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -89,6 +89,10 @@ void print_usage() { std::cout << " --calibration How the rings are measured: rings|spots (default: rings). rings sums the (q x azimuth) azimuthal profile over every processed image and fits the ring arcs in it; spots pools the found spots and fits those. -s/-e/--stride select the images; rings defaults --azim-phi-bins to 32" << std::endl; std::cout << std::endl; + std::cout << " Detector mask" << std::endl; + std::cout << " --detect-beam-stop[=N|off] Find the beam stop and its holder in a projection of N images and add them to the pixel mask (bit 9), so nothing shadowed by them is integrated. ON by default (60 images); =off disables. Reflections behind the stop are attenuated but not flagged, so they are integrated low with a plausible sigma and no existing rejection catches them" << std::endl; + std::cout << std::endl; + std::cout << " Spot finding" << std::endl; std::cout << " --spot-sigma Noise sigma level for spot finding (default: 4.0)" << std::endl; std::cout << " --spot-threshold Photon count threshold for spot finding (default: 10)" << std::endl; @@ -128,6 +132,7 @@ void print_usage() { std::cout << " --no-expected-variance-merge stills: disable the default expected-variance merge weighting (which rebuilds each weak observation's signal variance at the reflection mean to de-bias the inverse-variance merge); restores observed-sigma weighting" << std::endl; std::cout << " -A, --anomalous Anomalous mode (don't merge Friedel pairs)" << std::endl; std::cout << " --scaling-high-resolution High resolution limit for scaling/merging (manual override; default: no limit)" << std::endl; + std::cout << " --scaling-low-resolution Low resolution limit for scaling/merging, in A (default: 50, the value XDS configurations use; 0 = no limit). Reflections coarser than this sit behind or beside the beam stop and are measured on a background it has eaten into" << std::endl; std::cout << " --resolution-cutoff Automatic high-resolution cutoff for the written reflections + reported shells: cc-logistic|off (default: cc-logistic; ignored when --scaling-high-resolution is set)" << std::endl; std::cout << " --resolution-cc-target CC1/2 target defining the cc-logistic fall-off (default: 0.30)" << std::endl; std::cout << " --resolution-shells Number of resolution shells in the reported statistics table (default: 10)" << std::endl; @@ -189,6 +194,7 @@ enum { OPT_SEARCH_MIN_ZETA, OPT_SCALING_ITERATIONS, OPT_SCALING_HIGH_RESOLUTION, + OPT_SCALING_LOW_RESOLUTION, OPT_RESOLUTION_CUTOFF, OPT_RESOLUTION_CC_TARGET, OPT_RESOLUTION_SHELLS, @@ -198,6 +204,7 @@ enum { OPT_BACKGROUND_CLIP, OPT_BACKGROUND_RADIAL, OPT_REFINE_GEOMETRY, + OPT_DETECT_BEAM_STOP, OPT_BANDWIDTH, OPT_INTEGRATION_RADIUS, OPT_BACKGROUND_TRIM, @@ -289,6 +296,7 @@ static option long_options[] = { {"force-rotation-lattice", required_argument, nullptr, OPT_FORCE_ROTATION_LATTICE}, {"rotation-no-postrefine", no_argument, nullptr, OPT_ROTATION_NO_POSTREFINE}, {"refine-geometry", optional_argument, nullptr, OPT_REFINE_GEOMETRY}, + {"detect-beam-stop", optional_argument, nullptr, OPT_DETECT_BEAM_STOP}, {"spot-sigma", required_argument, nullptr, OPT_SPOT_SIGMA}, @@ -308,6 +316,7 @@ static option long_options[] = { {"search-min-zeta", required_argument, nullptr, OPT_SEARCH_MIN_ZETA}, {"scaling-iterations", required_argument, nullptr, OPT_SCALING_ITERATIONS}, {"scaling-high-resolution", required_argument, nullptr, OPT_SCALING_HIGH_RESOLUTION}, + {"scaling-low-resolution", required_argument, nullptr, OPT_SCALING_LOW_RESOLUTION}, {"background-clip", required_argument, nullptr, OPT_BACKGROUND_CLIP}, {"background-radial", optional_argument, nullptr, OPT_BACKGROUND_RADIAL}, {"resolution-cutoff", required_argument, nullptr, OPT_RESOLUTION_CUTOFF}, @@ -593,6 +602,7 @@ static int RunRugnux(int argc, char **argv) { std::optional background_clip_arg; // --background-clip: background-ring high-side sigma clip bool background_radial_given = false; // --background-radial seen at all (unset => auto) std::optional background_radial_arg; // when given: set = force on/off, unset = auto + std::optional detect_beam_stop = 60; // --detect-beam-stop[=N|off]; on by default std::optional refine_geometry; // --refine-geometry[=N]: stills global geometry-refinement pass bool refine_geometry_disabled = false; // --refine-geometry=off: opt out of the stills default-on @@ -606,6 +616,7 @@ static int RunRugnux(int argc, char **argv) { // it, and --spot-low-resolution 0 resets it to unset, i.e. no limit at that end. std::optional d_max_spot_finding = SpotFindingSettings{}.low_resolution_limit; std::optional d_min_scale_merge; + std::optional d_max_scale_merge; // --scaling-low-resolution; 0 removes the default limit std::optional resolution_cutoff_method; // --resolution-cutoff cc-logistic|off std::optional resolution_cc_target; // --resolution-cc-target std::optional report_shell_count; // --resolution-shells @@ -671,6 +682,17 @@ static int RunRugnux(int argc, char **argv) { case OPT_ROTATION_NO_POSTREFINE: rotation_postrefine_geometry = false; break; + case OPT_DETECT_BEAM_STOP: + // Frames projected to find the shadow. The default is what the detection was validated + // on; fewer leaves the background too sparsely counted to tell a shadow from noise. + if (optarg && std::string(optarg) == "off") { + detect_beam_stop = std::nullopt; + break; + } + detect_beam_stop = optarg + ? parse_number_arg(optarg, "--detect-beam-stop", logger, 1, 1000000) + : 60; + break; case OPT_REFINE_GEOMETRY: { if (optarg && std::string(optarg) == "off") { refine_geometry = std::nullopt; @@ -992,6 +1014,10 @@ static int RunRugnux(int argc, char **argv) { d_min_scale_merge = parse_number_arg(optarg, "--scaling-high-resolution", logger, 0.1f, 1000.0f); break; + case OPT_SCALING_LOW_RESOLUTION: + d_max_scale_merge = parse_number_arg(optarg, "--scaling-low-resolution", logger, + 0.0f, 100000.0f); + break; case OPT_RESOLUTION_CUTOFF: if (strcmp(optarg, "cc-logistic") == 0) resolution_cutoff_method = ResolutionCutoffMethod::CCHalfLogistic; @@ -1228,6 +1254,9 @@ static int RunRugnux(int argc, char **argv) { ScalingSettings scaling_settings = RugnuxDefaultScalingSettings(rot); if (d_min_scale_merge) scaling_settings.HighResolutionLimit_A(d_min_scale_merge.value()); + if (d_max_scale_merge) + scaling_settings.LowResolutionLimit_A(*d_max_scale_merge > 0.0f + ? std::optional(*d_max_scale_merge) : std::nullopt); if (resolution_cutoff_method) scaling_settings.ResolutionCutoff(*resolution_cutoff_method); if (resolution_cc_target) scaling_settings.ResolutionCCTarget(*resolution_cc_target); if (report_shell_count) scaling_settings.ReportShellCount(*report_shell_count); @@ -1508,6 +1537,7 @@ static int RunRugnux(int argc, char **argv) { config.stride = image_stride; config.nthreads = nthreads; config.output_prefix = output_prefix; + config.detect_beam_stop = detect_beam_stop; Rugnux process(reader, experiment, *dataset->pixel_mask, config); g_active_process = &process; @@ -1549,6 +1579,7 @@ static int RunRugnux(int argc, char **argv) { config.stride = image_stride; config.nthreads = nthreads; config.output_prefix = output_prefix; + config.detect_beam_stop = detect_beam_stop; config.write_process_h5 = false; // the .poni below is the output of this mode // Spot finding for --calibration spots. Indexing is off: a calibration wants the spot positions @@ -1725,6 +1756,9 @@ static int RunRugnux(int argc, char **argv) { scaling_settings.IceMinSpotRatio(static_cast(*ice_min_spot_ratio_arg)); if (d_min_scale_merge) scaling_settings.HighResolutionLimit_A(d_min_scale_merge.value()); + if (d_max_scale_merge) + scaling_settings.LowResolutionLimit_A(*d_max_scale_merge > 0.0f + ? std::optional(*d_max_scale_merge) : std::nullopt); if (resolution_cutoff_method) scaling_settings.ResolutionCutoff(*resolution_cutoff_method); if (resolution_cc_target) scaling_settings.ResolutionCCTarget(*resolution_cc_target); if (report_shell_count) scaling_settings.ReportShellCount(*report_shell_count); @@ -1886,6 +1920,7 @@ static int RunRugnux(int argc, char **argv) { config.spot_finding = spot_settings; config.rotation_indexing = rotation_indexing; config.two_pass_rotation = two_pass_rotation; + config.detect_beam_stop = detect_beam_stop; config.rotation_postrefine_geometry = rotation_postrefine_geometry; config.rotation_indexing_image_count = rotation_indexing_image_count; config.forced_rotation_lattice = forced_rotation_lattice; diff --git a/rugnux_vs_xds.py b/rugnux_vs_xds.py index 3423bf16..42aae0a6 100755 --- a/rugnux_vs_xds.py +++ b/rugnux_vs_xds.py @@ -80,8 +80,9 @@ def parse_xds(correct_lp): m = re.search(r"FRIEDEL'S_LAW=\s*(TRUE|FALSE)", txt) r["anomalous"] = (m.group(1) == "FALSE") if m else False - m = re.search(r"INCLUDE_RESOLUTION_RANGE=\s*[\d.]+\s+([\d.]+)", txt) - include_high = float(m.group(1)) if m else 0.0 + m = re.search(r"INCLUDE_RESOLUTION_RANGE=\s*([\d.]+)\s+([\d.]+)", txt) + include_low = float(m.group(1)) if m else 0.0 + include_high = float(m.group(2)) if m else 0.0 m = re.search(r"^\s*a\s+b\s+ISa\s*\n\s*[\d.Ee+-]+\s+[\d.Ee+-]+\s+([\d.]+)", txt, re.M) r["isa"] = float(m.group(1)) if m else None @@ -111,6 +112,10 @@ def parse_xds(correct_lp): # (XDS.INP often leaves the high limit at 0.0 = "use the full detector range"). table_high = float(shells[-1][0]) if shells else 0.0 r["dmin"] = (include_high if include_high > 0 else table_high) or None + # Same for the LOW-resolution limit. XDS configurations normally cut at 50 A while rugnux + # cuts at its own default, so without matching this the "lowest shell" columns below are not + # the same shell in the two programs and R_meas_lo is not comparable. + r["dmax"] = include_low or None # Mosaicity for comparison only. XDS refines this per run and the last value is the one it # settled on; it is a useful sanity check on rugnux's own estimate but NOT ground truth -- @@ -243,6 +248,8 @@ def run_rugnux(master, workdir, name, xds, rugnux_bin, threads, timeout, reuse, cmd.append("-A") if xds.get("dmin"): cmd += ["--scaling-high-resolution", f"{xds['dmin']:.3f}"] + if xds.get("dmax"): + cmd += ["--scaling-low-resolution", f"{xds['dmax']:.3f}"] if threads: cmd += ["-N", str(threads)] if extra_args: @@ -405,7 +412,8 @@ def main(): if not args.xds_only: if args.progress: print(f"[{i}/{len(crystals)}] rugnux {name} " - f"(dmin={xds.get('dmin')}, anom={xds.get('anomalous')}) ...", + f"(dmin={xds.get('dmin')}, dmax={xds.get('dmax')}, " + f"anom={xds.get('anomalous')}) ...", file=sys.stderr, flush=True) cif, err, elapsed = run_rugnux(master, workdir / name, name, xds, rugnux_bin, args.threads, args.timeout, args.reuse, @@ -426,7 +434,8 @@ def main(): print(f" rugnux vs XDS · {root} · {today}") if not args.xds_only: print(f" rugnux: {rugnux_bin}") - print(f" rugnux run de-novo, resolution + Friedel matched to XDS. XDS SG = symmetry") + print(f" rugnux run de-novo, resolution range (both limits) + Friedel matched to XDS.") + print(f" XDS SG = symmetry") print(f" CORRECT merged in; match is point-group level (screw/enantiomorph ignored). [ - = missing ]") print(f" Mosaic (deg): XDS cell = post-refined (CORRECT.LP) | integration-stage MLE (INTEGRATE.LP);") print(f" rugnux cell = median per-image mosaicity. Compare rugnux against the INTEGRATION-STAGE") diff --git a/tests/JFJochReaderTest.cpp b/tests/JFJochReaderTest.cpp index 243ce349..c3b5db93 100644 --- a/tests/JFJochReaderTest.cpp +++ b/tests/JFJochReaderTest.cpp @@ -210,6 +210,7 @@ TEST_CASE("JFJochReader_PixelMask", "[HDF5][Full]") { pixel_mask[5767] = 1; pixel_mask[x.GetPixelsNum() - 1] = 4; pixel_mask[0] = 256; + pixel_mask[3] = 1u << PixelMask::BeamStopPixelBit; ScanResultGenerator generator(x); @@ -250,6 +251,8 @@ TEST_CASE("JFJochReader_PixelMask", "[HDF5][Full]") { CHECK(reader_image->Image().at(0) == ERROR_PXL_VALUE); CHECK(reader_image->Image().at(1) == 0); CHECK(reader_image->Image().at(2) == 0); + // The beam-stop shadow reads back as its own marker, not as a bad pixel + CHECK(reader_image->Image().at(3) == BEAM_STOP_PXL_VALUE); CHECK(reader_image->Image().at(x.GetPixelsNum() - 1) == ERROR_PXL_VALUE); } remove("test16_master.h5"); diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index 0fd5d971..45e33624 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -87,4 +87,28 @@ TEST_CASE("HKLKey_sys_absence_P212121") { CHECK(hkl_key_gen.IsSystematicallyAbsent(0,0,5)); CHECK(!hkl_key_gen.IsSystematicallyAbsent(0,4,0)); CHECK(!hkl_key_gen.IsSystematicallyAbsent(5,5,5)); -} \ No newline at end of file +} +TEST_CASE("AcceptReflection_ResolutionLimits") { + Reflection r{}; + r.I = 100.0f; + r.sigma = 5.0f; + r.rlp = 1.0f; + r.d = 20.0f; + + // No limits: only the finiteness checks apply. + CHECK(AcceptReflection(r, std::nullopt, std::nullopt)); + + // Low-resolution limit rejects anything coarser than the limit, and is exclusive at it. + CHECK_FALSE(AcceptReflection(r, std::nullopt, std::optional(15.0))); + CHECK(AcceptReflection(r, std::nullopt, std::optional(20.0))); + CHECK(AcceptReflection(r, std::nullopt, std::optional(50.0))); + + // High-resolution limit still rejects anything finer, in the same direction as before. + CHECK_FALSE(AcceptReflection(r, std::optional(25.0), std::nullopt)); + CHECK(AcceptReflection(r, std::optional(2.0), std::optional(50.0))); + + // The plain-double overload treats 0 as "no limit" at both ends. + CHECK(AcceptReflection(r, 0.0, 0.0)); + CHECK_FALSE(AcceptReflection(r, 0.0, 15.0)); + CHECK(AcceptReflection(r, 2.0, 50.0)); +} diff --git a/tests/PixelMaskTest.cpp b/tests/PixelMaskTest.cpp index d8342eca..51be8586 100644 --- a/tests/PixelMaskTest.cpp +++ b/tests/PixelMaskTest.cpp @@ -350,4 +350,28 @@ TEST_CASE("PixelMask_GetMaskRaw_ThrowsForDECTRIS", "[PixelMask]") { REQUIRE_THROWS(mask.GetMaskRaw()); REQUIRE(mask.GetMask(experiment).size() == experiment.GetPixelsNum()); -} \ No newline at end of file +} +TEST_CASE("PixelMask_LoadBeamStopMask","[PixelMask]") { + DiffractionExperiment experiment(DetJF(1, 1)); + experiment.MaskModuleEdges(false).MaskChipEdges(false); + + PixelMask mask(experiment); + + std::vector shadow(experiment.GetPixelsNumConv(), 0); + shadow[345] = 1; + REQUIRE_NOTHROW(mask.LoadBeamStopMask(experiment, shadow)); + + CHECK((mask.GetMask()[345] & (1u << PixelMask::BeamStopPixelBit)) != 0); + CHECK(mask.GetMask()[346] == 0); + + // A user mask on the same pixel is carried alongside, not overwritten. + std::vector user(experiment.GetPixelsNumConv(), 0); + user[345] = 1; + REQUIRE_NOTHROW(mask.LoadUserMask(experiment, user)); + CHECK((mask.GetMask()[345] & (1u << PixelMask::BeamStopPixelBit)) != 0); + CHECK((mask.GetMask()[345] & (1u << PixelMask::UserMaskedPixelBit)) != 0); + + // Loading a mask of a different size is refused. + std::vector wrong_size(experiment.GetPixelsNumConv() + 1, 0); + CHECK_THROWS(mask.LoadBeamStopMask(experiment, wrong_size)); +} diff --git a/viewer/JFJochImageReadingWorker.cpp b/viewer/JFJochImageReadingWorker.cpp index ae7b9cbf..e132d3e5 100644 --- a/viewer/JFJochImageReadingWorker.cpp +++ b/viewer/JFJochImageReadingWorker.cpp @@ -1221,7 +1221,7 @@ QImage JFJochImageReadingWorker::RenderThumbnail_i(int64_t image_number, bool sh for (int yy = ty * by; yy < y1; ++yy) for (int xx = tx * bx; xx < x1; ++xx) { const int32_t v = px[yy * W + xx]; - if (v == GAP_PXL_VALUE || v == ERROR_PXL_VALUE) continue; + if (v == GAP_PXL_VALUE || v == ERROR_PXL_VALUE || v == BEAM_STOP_PXL_VALUE) continue; if (v == SATURATED_PXL_VALUE) { best = static_cast(fg); any = true; continue; } if (!any || v > best) { best = v; any = true; } } diff --git a/viewer/JFJochViewerSidePanel.cpp b/viewer/JFJochViewerSidePanel.cpp index 36445c94..b3c60a86 100644 --- a/viewer/JFJochViewerSidePanel.cpp +++ b/viewer/JFJochViewerSidePanel.cpp @@ -67,6 +67,12 @@ JFJochViewerSidePanel::JFJochViewerSidePanel(QWidget *parent) : QWidget(parent) connect(saturatedPixelsCheckBox, &QCheckBox::toggled, this, &JFJochViewerSidePanel::saturatedPixelsToggled); + auto beamStopCheckBox = new QCheckBox("Show beam stop", this); + beamStopCheckBox->setChecked(true); + beamStopCheckBox->setToolTip("Draw the detected beam-stop shadow (pixel mask bit 9) in coral. " + "Unchecked, it is drawn like any other masked pixel."); + connect(beamStopCheckBox, &QCheckBox::toggled, this, &JFJochViewerSidePanel::showBeamStop); + auto colorSelectButton = new QPushButton("Select feature color", this); connect(colorSelectButton, &QPushButton::clicked, this, [this]() { @@ -95,6 +101,7 @@ JFJochViewerSidePanel::JFJochViewerSidePanel(QWidget *parent) : QWidget(parent) image_feature_grid->addWidget(spotToggleCheckBox, 0, 0); image_feature_grid->addWidget(highlightIceRingToggleCheckBox, 0, 1); image_feature_grid->addWidget(predictionsToggleCheckBox, 1, 0); + image_feature_grid->addWidget(beamStopCheckBox, 1, 1); image_feature_grid->addWidget(saturatedPixelsCheckBox, 2, 0); image_feature_grid->addWidget(highestPixelsComboBox, 2, 1); diff --git a/viewer/JFJochViewerSidePanel.h b/viewer/JFJochViewerSidePanel.h index c067f0f5..1215c6e1 100644 --- a/viewer/JFJochViewerSidePanel.h +++ b/viewer/JFJochViewerSidePanel.h @@ -32,6 +32,7 @@ signals: void setSpotColor(QColor input); void showHighestPixels(int32_t v); void showSaturatedPixels(bool input); + void showBeamStop(bool input); void showPredictions(bool input); void highlightIceRings(bool input); void showROILabels(bool input); diff --git a/viewer/JFJochViewerWindow.cpp b/viewer/JFJochViewerWindow.cpp index 3501c328..939ffafa 100644 --- a/viewer/JFJochViewerWindow.cpp +++ b/viewer/JFJochViewerWindow.cpp @@ -340,6 +340,8 @@ JFJochViewerWindow::JFJochViewerWindow(QWidget *parent, bool dbus, const QString viewer, &JFJochDiffractionImage::showHighestPixels); connect(side_panel, &JFJochViewerSidePanel::showSaturatedPixels, viewer, &JFJochDiffractionImage::showSaturation); + connect(side_panel, &JFJochViewerSidePanel::showBeamStop, + viewer, &JFJochDiffractionImage::showBeamStop); connect(viewer, &JFJochDiffractionImage::writeStatusBar, statusbar, &JFJochViewerStatusBar::display); diff --git a/viewer/image_viewer/JFJochDiffractionImage.cpp b/viewer/image_viewer/JFJochDiffractionImage.cpp index 790602b8..05b99efb 100644 --- a/viewer/image_viewer/JFJochDiffractionImage.cpp +++ b/viewer/image_viewer/JFJochDiffractionImage.cpp @@ -97,6 +97,8 @@ void JFJochDiffractionImage::mouseHover(const QPointF &coord, Qt::KeyboardModifi intensity_str = " Gap "; else if (intensity == ERROR_PXL_VALUE) intensity_str = " Bad pxl "; + else if (intensity == BEAM_STOP_PXL_VALUE) + intensity_str = " Beam stop "; emit writeStatusBar(QString("x=%1 y=%2 %3 d=%4 Å") .arg(coord.x(), 0, 'f', 1) @@ -134,13 +136,15 @@ void JFJochDiffractionImage::ColorRow(size_t y, const PixelColorMap &map, QRgb * for (size_t x = 0; x < W; ++x) { const int32_t v = row[x]; - // The three sentinels are the extremes of the int32 range, so one range test - // separates them from every real pixel value + // The markers occupy the extremes of the int32 range, so one range test separates them + // from every real pixel value (MIN_REAL_PXL_VALUE moves when a marker is added) rgb c; - if (v > GAP_PXL_VALUE && v < SATURATED_PXL_VALUE) + if (v >= MIN_REAL_PXL_VALUE && v < SATURATED_PXL_VALUE) c = map.Apply(static_cast(v)); else if (v == GAP_PXL_VALUE) c = map.gap; + else if (v == BEAM_STOP_PXL_VALUE) + c = map.beam_stop; else c = (v == ERROR_PXL_VALUE) ? map.bad : map.saturated; @@ -975,6 +979,12 @@ void JFJochDiffractionImage::showSaturation(bool input) { updateOverlay(); } +void JFJochDiffractionImage::showBeamStop(bool input) { + show_beam_stop = input; + RenderImage(); + updateOverlay(); +} + void JFJochDiffractionImage::highlightIceRings(bool input) { highlight_ice_rings = input; updateOverlay(); @@ -1008,6 +1018,8 @@ QString JFJochDiffractionImage::PixelLabel(int x, int y) const { return QStringLiteral("Gap"); if (v == ERROR_PXL_VALUE) return QStringLiteral("Err"); + if (v == BEAM_STOP_PXL_VALUE) + return QStringLiteral("Stop"); if (v == SATURATED_PXL_VALUE) return QStringLiteral("Sat"); return QString::number(v); diff --git a/viewer/image_viewer/JFJochDiffractionImage.h b/viewer/image_viewer/JFJochDiffractionImage.h index f56fac9d..c94a7149 100644 --- a/viewer/image_viewer/JFJochDiffractionImage.h +++ b/viewer/image_viewer/JFJochDiffractionImage.h @@ -136,6 +136,7 @@ public slots: void showHighestPixels(int32_t v); void showSaturation(bool input); + void showBeamStop(bool input); void highlightIceRings(bool input); void setHDRMode(bool input); diff --git a/viewer/image_viewer/JFJochFollowerImage.cpp b/viewer/image_viewer/JFJochFollowerImage.cpp index 5aa201e9..48e6041c 100644 --- a/viewer/image_viewer/JFJochFollowerImage.cpp +++ b/viewer/image_viewer/JFJochFollowerImage.cpp @@ -18,6 +18,8 @@ static QString PixelValueText(int32_t value) { return QStringLiteral("Gap"); if (value == ERROR_PXL_VALUE) return QStringLiteral("Err"); + if (value == BEAM_STOP_PXL_VALUE) + return QStringLiteral("Stop"); if (value == SATURATED_PXL_VALUE) return QStringLiteral("Sat"); return QString::number(value); diff --git a/viewer/image_viewer/JFJochImage.cpp b/viewer/image_viewer/JFJochImage.cpp index 3f43121d..e06eb313 100644 --- a/viewer/image_viewer/JFJochImage.cpp +++ b/viewer/image_viewer/JFJochImage.cpp @@ -752,6 +752,7 @@ PixelColorMap JFJochImage::MakeColorMap() const { .bad = bad_color, // Saturation color .saturated = show_saturation ? bad_color : color_scale.Apply(1.0f), + .beam_stop = show_beam_stop ? color_scale.Apply(ColorScaleSpecial::BeamStop) : bad_color, }; } diff --git a/viewer/image_viewer/JFJochImage.h b/viewer/image_viewer/JFJochImage.h index 571f3cb4..8b43034c 100644 --- a/viewer/image_viewer/JFJochImage.h +++ b/viewer/image_viewer/JFJochImage.h @@ -31,7 +31,7 @@ struct PixelColorMap { float inv_range = 0.0f; float inv_range_log = 0.0f; bool hdr = false; - rgb gap{}, bad{}, saturated{}; + rgb gap{}, bad{}, saturated{}, beam_stop{}; [[nodiscard]] rgb Apply(float v) const { float f; @@ -97,6 +97,8 @@ protected: [[nodiscard]] virtual bool AllowROI() const { return false; } bool show_saturation = false; + // Beam-stop shadow drawn in its own colour; off shows it as an ordinary masked pixel. + bool show_beam_stop = true; bool auto_bg = false; bool auto_fg = false; diff --git a/viewer/widgets/JFJochViewerSettingsDock.cpp b/viewer/widgets/JFJochViewerSettingsDock.cpp index 6cca9526..41c16b31 100644 --- a/viewer/widgets/JFJochViewerSettingsDock.cpp +++ b/viewer/widgets/JFJochViewerSettingsDock.cpp @@ -623,6 +623,13 @@ QWidget *JFJochViewerSettingsDock::BuildScalingSection() { limitRes->setChecked(scaling_.GetHighResolutionLimit_A().has_value()); auto *highRes = new NumberLineEdit(0.3f, 5.0f, scaling_.GetHighResolutionLimit_A().value_or(2.0), 1, "Å", this); highRes->setEnabled(limitRes->isChecked()); + auto *limitLowRes = new QCheckBox("Low-resolution limit", this); + limitLowRes->setChecked(scaling_.GetLowResolutionLimit_A().has_value()); + limitLowRes->setToolTip("Drop reflections coarser than this from scaling and merging. They sit behind " + "or beside the beam stop and are measured on a background it has eaten into. " + "On by default at 50 Å, the value XDS configurations use."); + auto *lowRes = new NumberLineEdit(5.0f, 500.0f, scaling_.GetLowResolutionLimit_A().value_or(50.0), 1, "Å", this); + lowRes->setEnabled(limitLowRes->isChecked()); form->addRow("", friedel); form->addRow("", corrections); @@ -633,6 +640,10 @@ QWidget *JFJochViewerSettingsDock::BuildScalingSection() { resRow->addWidget(limitRes); resRow->addWidget(highRes, 1); form->addRow("", resRow); + auto *lowResRow = new QHBoxLayout(); + lowResRow->addWidget(limitLowRes); + lowResRow->addWidget(lowRes, 1); + form->addRow("", lowResRow); section->setContentLayout(form); section->setExpanded(false); // folded on start (only geometry + unit cell start open) @@ -642,6 +653,8 @@ QWidget *JFJochViewerSettingsDock::BuildScalingSection() { scaling_.StillsPartialityRefine(partRefine->isChecked()); scaling_.HighResolutionLimit_A(limitRes->isChecked() ? std::optional(highRes->value()) : std::nullopt); + scaling_.LowResolutionLimit_A(limitLowRes->isChecked() + ? std::optional(lowRes->value()) : std::nullopt); emit scalingChanged(scaling_); }; connect(friedel, &QCheckBox::toggled, this, [emitScaling] { emitScaling(); }); @@ -650,6 +663,9 @@ QWidget *JFJochViewerSettingsDock::BuildScalingSection() { connect(limitRes, &QCheckBox::toggled, this, [emitScaling, highRes](bool on) { highRes->setEnabled(on); emitScaling(); }); connect(highRes, &NumberLineEdit::newValue, this, [emitScaling] { emitScaling(); }); + connect(limitLowRes, &QCheckBox::toggled, this, [emitScaling, lowRes](bool on) { + lowRes->setEnabled(on); emitScaling(); }); + connect(lowRes, &NumberLineEdit::newValue, this, [emitScaling] { emitScaling(); }); return section; } diff --git a/viewer/windows/JFJochProcessingJobsWindow.cpp b/viewer/windows/JFJochProcessingJobsWindow.cpp index 7598fa6f..a2bb326a 100644 --- a/viewer/windows/JFJochProcessingJobsWindow.cpp +++ b/viewer/windows/JFJochProcessingJobsWindow.cpp @@ -245,6 +245,29 @@ int JFJochProcessingJobsWindow::askJob(const ReprocessingInputs &inputs, JobSpec "axis from the whole sweep, then re-integrates at the refined geometry. The refined pass is the " "canonical output; the header-geometry pass is kept as _01. Default on; a no-op for stills."); + // Beam-stop shadow: a projection of a few frames shows where the stop and its holder shadow the + // detector; those pixels go into the mask (bit 9) so nothing behind them is integrated. Cheap and + // useful in every mode, so it is offered for all of them and defaulted on. + auto *beam_stop = new QCheckBox("Detect beam stop", &dlg); + beam_stop->setChecked(true); + beam_stop->setToolTip( + "Project a few frames, find where the beam stop and its holder shadow the detector, and add them " + "to the pixel mask (bit 9, cleared at the start of every run). Reflections behind the stop are " + "attenuated but not flagged, so they otherwise integrate low with a plausible sigma."); + auto *beam_stop_frames = new QSpinBox(&dlg); + beam_stop_frames->setRange(3, 100000); + beam_stop_frames->setValue(60); + beam_stop_frames->setToolTip("Frames projected to find the shadow. Fewer leaves the background too " + "sparsely counted to tell a shadow from noise."); + connect(beam_stop, &QCheckBox::toggled, beam_stop_frames, &QWidget::setEnabled); + auto *beamStopRow = new QWidget(&dlg); + auto *beamStopRowLayout = new QHBoxLayout(beamStopRow); + beamStopRowLayout->setContentsMargins(0, 0, 0, 0); + beamStopRowLayout->addWidget(beam_stop); + beamStopRowLayout->addWidget(new QLabel("images:", &dlg)); + beamStopRowLayout->addWidget(beam_stop_frames); + beamStopRowLayout->addStretch(); + auto *form = new QFormLayout; form->addRow("Start image", start_image); form->addRow("End image", end_image); @@ -255,6 +278,7 @@ int JFJochProcessingJobsWindow::askJob(const ReprocessingInputs &inputs, JobSpec form->addRow(scaling); form->addRow(refineRow); form->addRow(postrefine); + form->addRow(beamStopRow); // Calibration: the calibrant and method come from the panel's Calib page; state them here (with the // azimuthal sector count the rings method depends on) so the run is not a surprise. @@ -307,6 +331,8 @@ int JFJochProcessingJobsWindow::askJob(const ReprocessingInputs &inputs, JobSpec spec.refine_geometry = refine_geometry->isEnabled() && refine_geometry->isChecked(); spec.refine_geometry_frames = refine_frames->value(); spec.rotation_postrefine = postrefine->isEnabled() && postrefine->isChecked(); + spec.detect_beam_stop = beam_stop->isChecked(); + spec.detect_beam_stop_frames = beam_stop_frames->value(); } return result; } @@ -325,6 +351,8 @@ ProcessConfig JFJochProcessingJobsWindow::buildConfig(const JobSpec &spec, const config.write_process_h5 = spec.save_h5; config.write_merged = spec.save_merged; config.spot_finding = inputs.spot_finding; + if (spec.detect_beam_stop) + config.detect_beam_stop = spec.detect_beam_stop_frames; if (spec.mode == ProcessMode::Calibration) { config.calibration_method = spec.calibration.method; config.calibrant_ring_q = spec.calibration.ring_q; @@ -398,6 +426,7 @@ void JFJochProcessingJobsWindow::newJob(ProcessMode mode, CalibrationSelection c scaling.CorrectionSurfaces(dock.GetCorrectionSurfaces()); scaling.StillsPartialityRefine(dock.GetStillsPartialityRefine()); scaling.HighResolutionLimit_A(dock.GetHighResolutionLimit_A()); + scaling.LowResolutionLimit_A(dock.GetLowResolutionLimit_A()); experiment.ImportScalingSettings(scaling); } diff --git a/viewer/windows/JFJochProcessingJobsWindow.h b/viewer/windows/JFJochProcessingJobsWindow.h index 384f0671..162cb278 100644 --- a/viewer/windows/JFJochProcessingJobsWindow.h +++ b/viewer/windows/JFJochProcessingJobsWindow.h @@ -93,6 +93,8 @@ private: bool refine_geometry = false; // stills-only global geometry bundle-adjust (needs a known cell) int refine_geometry_frames = 200; // strong indexed frames fed to that bundle adjustment bool rotation_postrefine = true; // rotation-only two-pass geometry post-refine (default on) + bool detect_beam_stop = true; // project frames and mask the beam-stop shadow (default on) + int detect_beam_stop_frames = 60; // frames projected to find it CalibrationSelection calibration; // Calibration mode: calibrant rings + rings/spots method };