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