diff --git a/image_analysis/beam_stop/SHADOW_FINDER.md b/image_analysis/beam_stop/SHADOW_FINDER.md new file mode 100644 index 00000000..64f5b37b --- /dev/null +++ b/image_analysis/beam_stop/SHADOW_FINDER.md @@ -0,0 +1,138 @@ +# Beam-stop shadow finder + +`ShadowFinder` detects the beam-stop shadow — the central **disk** plus its **holder arm** — +from a small number of images, and returns a per-pixel mask of the shadowed region. It exists +so the pipeline can mask the beam stop **automatically** instead of relying on a hand-drawn user +mask. + +It is written in the accumulate-then-finalize shape of `DarkMaskAnalysis`: feed frames with +`AddImage()`, then read the mask once with `GetMask()`. **It is not yet wired into any workflow** +— see *Wiring plan* below. The class is self-contained (`common/` only) and unit-testable. + +## API + +```cpp +ShadowFinder finder(experiment); // geometry (beam centre, converted size) +std::vector scratch; +for (int i : first_N_frames) // ~20–50 spread frames is plenty (see below) + finder.AddImage(read_frame(i), scratch); // DataMessage with .image, like DarkMaskAnalysis +std::vector shadow = finder.GetMask(); // 1 = shadow, 0 = keep; converted geometry +``` + +- **Input** (`AddImage`): a full **converted-geometry** image inside a `DataMessage`, dispatched by + pixel type exactly like `DarkMaskAnalysis::AnalyzeImage`. A pixel equal to its type's sentinel + extreme (`INT*_MIN` for signed, `UINT*_MAX` for unsigned) is treated as *no data* (module gap / + masked) and skipped. For signed types the opposite extreme (`INT*_MAX`, a genuine saturation + clamp) is kept, so saturated reflections still register as bright. Thread-safe. +- **Output** (`GetMask`): a `0/1` vector of size `GetPixelsNumConv()`, ready to feed + `PixelMask::LoadMask(..., bit)`. Recomputed from the accumulators on each call — call it once at + the end. +- **Settings** (`ShadowFinderSettings`): plain struct with defaults (below). When wired, this can + move onto `DiffractionExperiment` the way `DarkMaskSettings` does. + +## Algorithm + +The shadow is a **stationary, azimuthal anomaly**: air scatter makes the background fall off +smoothly with radius but stay roughly constant around each ring, so the beam stop is a localized +dip at a fixed place. Detection therefore compares each pixel to the typical background **at its +own radius**, not to a global threshold. + +1. **Projection.** Per pixel accumulate the **max** and the **mean** (sum + count) over the added + frames, skipping the no-data sentinel. The stop is stationary, so it stays dark while spots and + the air-scatter halo build up around it. +2. **Robust radial baseline.** Take the **median of the mean-projection in each 1-pixel ring** + about the beam centre. The median ignores the shadow (a minority of each ring) and the spots. + Iterate 3×, excluding the shadow it finds, so the baseline is not pulled down by the shadow. +3. **Ratio.** `ratio = mean / baseline(r)`. Shadowed pixels sit far below 1; the smooth halo and + symmetric ice rings cancel to ≈1 and are **not** flagged. +4. **Core + connectivity.** `core = ratio < shadow_ratio`, then keep only the part **connected to + the beam centre** (flood from a small central seed through a `bridge_px`-dilated copy, so the + arm crosses module gaps). The centre anchor is what makes it specific to the beam stop rather + than any random low patch. +5. **Reflection guard.** A pixel whose **max**-projection reached `min_reflection` recorded a real + reflection — a beam stop cannot block a reflection that was measured — so it is **never masked** + (a small-cluster test drops single-pixel zingers). +6. **Central disk (low-res guard).** From the azimuthal **blocked fraction** about the beam centre, + mask a full disk out to the radius where nearly every azimuth is blocked (a disk blocks all + directions; a thin arm/gap does not). Returns nothing when the centre is not behind a disk, and + is **capped just inside the innermost reflection** so a real edge spot is never swallowed. This + guards the low-resolution zone the stop makes un-measurable, including a disk the raw data has + already voided. +7. **Soft boundary.** Grow the core+disk into partially-shadowed pixels (`ratio < penumbra_ratio`) + within `penumbra_max_px`, then a morphological close rounds the contour and hole-filling closes + the disk interior. +8. **Carve reflections last.** Remove any lit pixel from the final mask (no fill afterwards) so a + spot the geometry still covered is given back rather than re-enclosed. + +Three findings from the prototype are baked in: the **mean** ratio (not the max) does the +detecting — the max projection is a great *picture* but a noisier detector; the **centre anchor** +is what distinguishes the stop from stray low patches; and **a lit pixel is never masked**. +Detection converges in **~2–3 frames** (the shadow cuts the halo present in every frame); 20–50 is +just safe headroom. + +## Parameters (`ShadowFinderSettings`) + +| field | default | meaning | +|---|---|---| +| `shadow_ratio` | 0.35 | core threshold: mean below this fraction of the radial baseline | +| `penumbra_ratio` | 0.72 | soft boundary grows out to this ratio … | +| `penumbra_max_px` | 14 | … but no further than this from the core | +| `bridge_px` | 6 | bridge module gaps / breaks the arm crosses | +| `min_reflection` | 25 | max-projection value at/above which a pixel is a protected reflection | + +## Complexity & memory + +Per `AddImage`: `O(pixels)`. `GetMask`: a handful of `O(pixels)` sweeps (three ring-median passes +and a few BFS-based morphology ops), i.e. a couple of seconds on a full ~18 MP detector — a +one-time finalize. Accumulators are `int32 + int64 + uint32` per pixel (~0.3 GB at 18 MP); the +finalize adds a few transient `float`/`char` buffers. Fine for the offline (Rugnux/viewer) +context it targets; not intended for the online per-frame hot path. + +## Wiring plan (deferred — implement later) + +The detector is deliberately reader-agnostic and bit-agnostic: one algorithm, fed frames by each +caller, with the caller deciding where the resulting mask lands. Two intended consumers: + +**1. Build.** Add `beam_stop/ShadowFinder.cpp` + `.h` to the `JFJochImageAnalysis` source list in +`image_analysis/CMakeLists.txt` (mirroring `dark_mask_analysis/DarkMaskAnalysis.*` at lines 27–28). + +**2. Mask bit + loader** (`common/PixelMask.{h,cpp}`): add `BeamStopShadowBit = 9` (bit 9 is free; +`PixelMask.h:31-36`), a `LoadBeamStopMask` mirroring `LoadDarkBadPixelMask` (`.cpp:298-315`), a +`beam_stop_pixel` counter in `PixelMaskStatistics` / `GetStatistics` (`.cpp:186-205`), and +optionally `GetBeamStopMask()` (mirroring `GetUserMask`). Leave `GetUserMask` untouched (keeps the +user-mask channel clean). Application is then automatic: the CPU/GPU preprocessor masks any pixel +with a non-zero bit (`ImagePreprocessorCPU.cpp:10`, `ImagePreprocessorGPU.cu:104`), and HDF5 +write/read + CBOR already carry the full bitfield, so bit 9 persists to `_process.h5` and reloads +with no extra code. + +**3. Rugnux offline processing → bit 9** (dataset-specific, per run). Add a `detect_beam_stop` +flag to `ProcessConfig` (`rugnux/Rugnux.h`) and a `--detect-beam-stop` CLI flag mirroring +`--detect-ice-rings` in `rugnux/rugnux_cli.cpp`. In `Rugnux::Run()`, after geometry/goniometer +setup and **before** the azimuthal mapping (`Rugnux.cpp:281`) and the `start_message.pixel_mask` +persist (`:296`), run a `ShadowFinder` pre-pass over a spread sample (`select_equally_spaced_image_ordinals` ++ `reader_.GetRawImage`, as `RefineStillsGeometry` does at `:96-157`) and +`pixel_mask_.LoadBeamStopMask(...)`. Everything downstream then applies and persists it for free. + +**4. Viewer standalone → user mask (bit 8)** and a viewer preview. A "Find beam stop" action reads +the first N frames via the worker's `file_reader`, runs `ShadowFinder`, and holds the result as +bit 9 for a distinct-coloured, toggleable overlay; a "Transfer to user mask" action ORs it into the +user mask through the existing chokepoint `JFJochImageReadingWorker::UpdateUserMask_i` +(`.cpp:783`, same idiom as `LoadUserMaskTIFF` add-mode `:851-858`). The distinct colour needs a +`ColorScaleSpecial::BeamStop` + sentinel branch in `JFJochReaderImage.cpp:117-124`; the on/off +toggle mirrors "Show saturated pixels". + +**5. Docs.** Add a "Bit 9 – beam-stop shadow" row to the bit table in `docs/PIXEL_MASK.md:9-21`. + +Not needed: the online broker / receiver / FPGA / REST-API / React-frontend paths — the auto +beam-stop mask is an **offline** (Rugnux + viewer) feature and never goes to the detector. + +## Known limitations / future refinements + +- **Boundary smoothing** uses a morphological close; the prototype used a Gaussian level-set for a + slightly smoother edge. Upgrade if the edge looks blocky. +- **Off-centre disk**: the low-res disk guard keys on the *stored* beam centre. If that is off the + physical stop (seen on one stills detector), a small already-dead disk gets no extra circular + guard — a data-driven disk centre would close it. +- **Unsigned saturation**: for unsigned images a detector-saturating reflection collides with the + no-data sentinel and is skipped; harmless in practice (such reflections do not sit under the + stop), and absent entirely on the signed/preprocessed path. diff --git a/image_analysis/beam_stop/ShadowFinder.cpp b/image_analysis/beam_stop/ShadowFinder.cpp new file mode 100644 index 00000000..c2044337 --- /dev/null +++ b/image_analysis/beam_stop/ShadowFinder.cpp @@ -0,0 +1,338 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "ShadowFinder.h" + +#include +#include +#include +#include +#include + +#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. +// --------------------------------------------------------------------------------- +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) { + if (r <= 0) + return in; + std::vector dist(in.size(), -1); + std::queue q; + for (size_t i = 0; i < in.size(); i++) + if (in[i]) { dist[i] = 0; q.push(static_cast(i)); } + while (!q.empty()) { + const int i = q.front(); q.pop(); + if (dist[i] >= r) + continue; + const int y = i / W, x = i % W; + for (int dy = -1; dy <= 1; dy++) + for (int dx = -1; dx <= 1; dx++) { + const int yy = y + dy, xx = x + dx; + if (yy < 0 || yy >= H || xx < 0 || xx >= W) + continue; + const int j = yy * W + xx; + if (dist[j] < 0) { dist[j] = dist[i] + 1; q.push(j); } + } + } + std::vector out(in.size()); + for (size_t i = 0; i < out.size(); i++) + out[i] = (dist[i] >= 0) ? 1 : 0; + 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) { + 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); + std::vector out(in.size()); + for (size_t i = 0; i < out.size(); i++) + out[i] = !grown[i]; + return out; +} + +// 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 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); + } + while (!q.empty()) { + const int i = q.front(); q.pop(); + const int y = i / W, x = i % W; + for (int dy = -1; dy <= 1; dy++) + for (int dx = -1; dx <= 1; dx++) { + const int yy = y + dy, xx = x + dx; + if (yy < 0 || yy >= H || xx < 0 || xx >= W) + continue; + const int j = yy * W + xx; + if (passable[j] && !visited[j]) { visited[j] = 1; q.push(j); } + } + } + return visited; +} + +// Fill holes: background not reachable from the image border becomes region. +std::vector FillHoles(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); } }; + for (int x = 0; x < W; x++) { push(x); push((H - 1) * W + x); } + for (int y = 0; y < H; y++) { push(y * W); push(y * W + W - 1); } + while (!q.empty()) { + const int i = q.front(); q.pop(); + const int y = i / W, x = i % W; + for (int dy = -1; dy <= 1; dy++) + for (int dx = -1; dx <= 1; dx++) { + const int yy = y + dy, xx = x + dx; + if (yy < 0 || yy >= H || xx < 0 || xx >= W) + continue; + const int j = yy * W + xx; + if (!region[j] && !bg_visited[j]) { bg_visited[j] = 1; q.push(j); } + } + } + std::vector out = region; + for (size_t i = 0; i < out.size(); i++) + if (!region[i] && !bg_visited[i]) + out[i] = 1; + 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> bins(max_radius + 1); + for (size_t i = 0; i < values.size(); i++) + if (use[i]) + bins[radius[i]].push_back(values[i]); + std::vector median(max_radius + 1, 0.0f); + for (int r = 0; r <= max_radius; r++) { + auto &b = bins[r]; + if (!b.empty()) { + const size_t k = b.size() / 2; + std::nth_element(b.begin(), b.begin() + k, b.end()); + median[r] = b[k]; + } + } + 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]]++; + } + 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; +} + +} // namespace + +// --------------------------------------------------------------------------------- + +ShadowFinder::ShadowFinder(const DiffractionExperiment &experiment, ShadowFinderSettings in_settings) + : width(static_cast(experiment.GetXPixelsNumConv())), + height(static_cast(experiment.GetYPixelsNumConv())), + beam_x(experiment.GetBeamX_pxl()), + beam_y(experiment.GetBeamY_pxl()), + settings(in_settings), + max_value(static_cast(width) * height, 0), + sum_value(static_cast(width) * height, 0), + valid_count(static_cast(width) * height, 0) {} + +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. + T masked; + if constexpr (std::is_signed_v) + masked = std::numeric_limits::min(); + else + masked = std::numeric_limits::max(); + + std::unique_lock ul(m); + for (size_t i = 0; i < max_value.size(); i++) { + const T v = ptr[i]; + if (v == masked) + continue; + const int32_t vi = static_cast(v); + if (valid_count[i] == 0 || vi > max_value[i]) + max_value[i] = vi; + sum_value[i] += vi; + valid_count[i]++; + } + frames++; +} + +void ShadowFinder::AddImage(const DataMessage &data, std::vector buffer) { + if (static_cast(data.image.GetWidth()) * data.image.GetHeight() != max_value.size()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "ShadowFinder: image size does not match the detector"); + + const auto ptr = data.image.GetUncompressedPtr(buffer); + switch (data.image.GetMode()) { + case CompressedImageMode::Int8: Add(reinterpret_cast(ptr)); break; + case CompressedImageMode::Uint8: Add(reinterpret_cast(ptr)); break; + case CompressedImageMode::Int16: Add(reinterpret_cast(ptr)); break; + case CompressedImageMode::Uint16: Add(reinterpret_cast(ptr)); break; + case CompressedImageMode::Int32: Add(reinterpret_cast(ptr)); break; + case CompressedImageMode::Uint32: Add(reinterpret_cast(ptr)); break; + default: + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "ShadowFinder: unsupported image mode"); + } +} + +uint32_t ShadowFinder::GetFrameCount() const { + std::unique_lock ul(m); + return frames; +} + +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; + + std::vector mask(N, 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); + 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) { + 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; + } + + // --- robust radial baseline; iterate to keep the shadow out of its own baseline --- + std::vector ratio(N, 1.0f); + std::vector excluded(N, 0); + for (int iter = 0; iter < 3; iter++) { + std::vector use(N); + for (int i = 0; i < N; i++) + use[i] = valid[i] && !excluded[i]; + const auto baseline = RingMedian(mean, use, radius, max_radius); + for (int i = 0; i < N; 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; + } + + // --- 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; + + 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]; + + // --- 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); + for (int y = 0; y < H; y++) + for (int x = 0; x < W; x++) { + const int i = y * W + x; + if (!lit[i]) continue; + int neighbours = 0; + for (int dy = -1; dy <= 1; dy++) + for (int dx = -1; dx <= 1; dx++) { + const int yy = y + dy, xx = x + dx; + if ((dx || dy) && yy >= 0 && yy < H && xx >= 0 && xx < W && lit[yy * W + xx]) + neighbours++; + } + 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) + region[i] = 1; + + region = Erode(Dilate(region, W, H, 2), W, H, 2); // close: round the boundary + region = FillHoles(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++) + if (reflection_grown[i]) + region[i] = 0; + + for (int i = 0; i < N; 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 new file mode 100644 index 00000000..14ee4adc --- /dev/null +++ b/image_analysis/beam_stop/ShadowFinder.h @@ -0,0 +1,76 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include + +#include "../../common/CompressedImage.h" +#include "../../common/DiffractionExperiment.h" +#include "../../common/JFJochMessages.h" // DataMessage + +// 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. +// +// 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. +// 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; + + // 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 + uint32_t frames = 0; + + template void Add(const T *ptr); + +public: + ShadowFinder(const DiffractionExperiment &experiment, ShadowFinderSettings settings = {}); + + // 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). + 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). + [[nodiscard]] std::vector GetMask() const; + + [[nodiscard]] uint32_t GetFrameCount() const; +};