Jungfraujoch/image_analysis/beam_stop/ShadowFinder.cpp

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#include "ShadowFinder.h"

#include <algorithm>
#include <cmath>
#include <limits>
#include <atomic>
#include <future>
#include <thread>
#include <queue>
#include <type_traits>

#include <spdlog/spdlog.h>

#include "../../common/CUDAWrapper.h"
#include "../../common/ParallelFor.h"
#include "../../common/JFJochException.h"

// 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. Its area is how many pixels back
// a ring's countability test, which decides where an azimuthal comparison is possible at all.
constexpr int    POOL_PX = 5;
constexpr double MEAN_POOLED_PIXELS = POOL_PX * POOL_PX;

// A ring with fewer valid pixels than this says nothing about whether it was counted.
constexpr int MIN_RING_PIXELS = 32;

// 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, i.e. every pixel within Chebyshev distance r of a set one.
//
// This was a multi-source BFS, which is what the distance is defined by - but on a full rectangle
// with no obstacles the 8-connected graph distance IS the Chebyshev distance (a path stepping
// towards the target never has to leave the frame), so the result is a dilation by the (2r+1) square
// clipped to the frame. A square is the product of a horizontal and a vertical segment and max is
// associative, so it separates into one pass along x and one along y - O(1) per pixel whatever r is,
// no queue, and no 4-bytes-per-pixel distance array. The values are 0/1, so the running "is any set"
// is just a count of set pixels in the window.
std::vector<char> dilate(const std::vector<char> &in, int W, int H, int r, size_t nthreads) {
    if (r <= 0)
        return in;
    std::vector<char> tmp(in.size()), out(in.size());

    ParallelChunks(H, nthreads, [&](int ylo, int yhi) {
        for (int y = ylo; y < yhi; y++) {
            const char *src = in.data() + static_cast<size_t>(y) * W;
            char *dst = tmp.data() + static_cast<size_t>(y) * W;
            int count = 0;
            for (int x = 0; x <= std::min(r, W - 1); x++)
                count += src[x];
            for (int x = 0; x < W; x++) {
                dst[x] = count > 0;
                if (x + r + 1 < W) count += src[x + r + 1];
                if (x - r >= 0)    count -= src[x - r];
            }
        }
    });

    // The vertical pass walks a column, which strides by a whole row. Take a strip of columns at a
    // time so each row it touches is read contiguously instead of one byte per cache line.
    constexpr int STRIP = 64;
    ParallelChunks((W + STRIP - 1) / STRIP, nthreads, [&](int slo, int shi) {
        std::vector<int> count(STRIP);
        for (int st = slo; st < shi; st++) {
            const int x0 = st * STRIP, xn = std::min(STRIP, W - x0);
            std::fill(count.begin(), count.begin() + xn, 0);
            for (int y = 0; y <= std::min(r, H - 1); y++)
                for (int i = 0; i < xn; i++)
                    count[i] += tmp[static_cast<size_t>(y) * W + x0 + i];
            for (int y = 0; y < H; y++) {
                for (int i = 0; i < xn; i++)
                    out[static_cast<size_t>(y) * W + x0 + i] = count[i] > 0;
                if (y + r + 1 < H)
                    for (int i = 0; i < xn; i++)
                        count[i] += tmp[static_cast<size_t>(y + r + 1) * W + x0 + i];
                if (y - r >= 0)
                    for (int i = 0; i < xn; i++)
                        count[i] -= tmp[static_cast<size_t>(y - r) * W + x0 + i];
            }
        }
    });
    return out;
}

// Erosion by `r` = dilation of the complement. The dilation is clipped to the frame and cannot seed
// outside it, so outside the frame contributes nothing - a pixel within r of the edge is eroded only
// by what the frame actually holds.
std::vector<char> erode(const std::vector<char> &in, int W, int H, int r, size_t nthreads) {
    std::vector<char> comp(in.size());
    ParallelChunks(static_cast<int>(in.size()), nthreads, [&](int lo, int hi) {
        for (int i = lo; i < hi; i++) comp[i] = !in[i];
    });
    const auto grown = dilate(comp, W, H, r, nthreads);
    std::vector<char> out(in.size());
    ParallelChunks(static_cast<int>(out.size()), nthreads, [&](int lo, int hi) {
        for (int i = lo; i < hi; i++) out[i] = !grown[i];
    });
    return out;
}

// Pixels of `passable` reachable from any of `seeds` (8-connected flood).
std::vector<char> flood(const std::vector<char> &passable, int W, int H, const std::vector<int> &seeds) {
    std::vector<char> visited(passable.size(), 0);
    std::queue<int> q;
    for (const int s : seeds)
        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;
        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.
//
// The flood is run over the bounding box of `region` grown by one, not the whole detector. Outside
// that box every pixel is background and the box's surrounding ring is background too, so the whole
// outside is one border-connected component: a background pixel inside the box is border-connected
// exactly when it reaches the ring. The beam stop occupies a small part of a detector, so this is
// the same answer over a fraction of the pixels.
std::vector<char> fill_holes(const std::vector<char> &region, int W, int H) {
    int x0 = W, x1 = -1, y0 = H, y1 = -1;
    for (int y = 0; y < H; y++)
        for (int x = 0; x < W; x++)
            if (region[static_cast<size_t>(y) * W + x]) {
                x0 = std::min(x0, x); x1 = std::max(x1, x);
                y0 = std::min(y0, y); y1 = std::max(y1, y);
            }
    if (x1 < 0)
        return region;   // nothing to enclose
    x0 = std::max(0, x0 - 1); x1 = std::min(W - 1, x1 + 1);
    y0 = std::max(0, y0 - 1); y1 = std::min(H - 1, y1 + 1);

    const int BW = x1 - x0 + 1, BH = y1 - y0 + 1;
    std::vector<char> bg_visited(static_cast<size_t>(BW) * BH, 0);
    std::queue<int> q;   // indices into the box
    auto push = [&](int bx, int by) {
        const int j = by * BW + bx;
        if (!region[static_cast<size_t>(by + y0) * W + bx + x0] && !bg_visited[j]) {
            bg_visited[j] = 1; q.push(j);
        }
    };
    for (int bx = 0; bx < BW; bx++) { push(bx, 0); push(bx, BH - 1); }
    for (int by = 0; by < BH; by++) { push(0, by); push(BW - 1, by); }
    while (!q.empty()) {
        const int i = q.front(); q.pop();
        const int by = i / BW, bx = i % BW;
        for (int dy = -1; dy <= 1; dy++)
            for (int dx = -1; dx <= 1; dx++) {
                const int yy = by + dy, xx = bx + dx;
                if (yy < 0 || yy >= BH || xx < 0 || xx >= BW)
                    continue;
                const int j = yy * BW + xx;
                if (!region[static_cast<size_t>(yy + y0) * W + xx + x0] && !bg_visited[j]) {
                    bg_visited[j] = 1; q.push(j);
                }
            }
    }

    std::vector<char> out = region;
    for (int by = 0; by < BH; by++)
        for (int bx = 0; bx < BW; bx++) {
            const size_t i = static_cast<size_t>(by + y0) * W + bx + x0;
            if (!region[i] && !bg_visited[by * BW + bx])
                out[i] = 1;
        }
    return out;
}

// Sum of `in` over the k x k box centred on each pixel, zero outside the frame.
//
// Each row's and each column's running sum keeps exactly the terms it had, in the order it had them,
// so the floating-point form rounds identically - only the traversal changes. Walking one column at
// a time strides a whole row per step and misses on every access, so the vertical pass takes a strip
// of columns together and reads each row it touches contiguously. Rows, and strips, are independent.
template <typename T>
std::vector<T> box_sum(const std::vector<T> &in, int W, int H, int k, size_t nthreads) {
    const int half = k / 2;
    std::vector<T> row(in.size()), out(in.size());

    ParallelChunks(H, nthreads, [&](int ylo, int yhi) {
        for (int y = ylo; y < yhi; y++) {
            const T *src = in.data() + static_cast<size_t>(y) * W;
            T *dst = row.data() + static_cast<size_t>(y) * W;
            T s = 0;
            for (int x = 0; x <= std::min(half, W - 1); x++)
                s += src[x];
            for (int x = 0; x < W; x++) {
                dst[x] = s;
                if (x + half + 1 < W) s += src[x + half + 1];
                if (x - half >= 0)    s -= src[x - half];
            }
        }
    });

    constexpr int STRIP = 64;
    ParallelChunks((W + STRIP - 1) / STRIP, nthreads, [&](int slo, int shi) {
        std::vector<T> s(STRIP);
        for (int st = slo; st < shi; st++) {
            const int x0 = st * STRIP, xn = std::min(STRIP, W - x0);
            std::fill(s.begin(), s.begin() + xn, T{0});
            for (int y = 0; y <= std::min(half, H - 1); y++)
                for (int i = 0; i < xn; i++)
                    s[i] += row[static_cast<size_t>(y) * W + x0 + i];
            for (int y = 0; y < H; y++) {
                for (int i = 0; i < xn; i++)
                    out[static_cast<size_t>(y) * W + x0 + i] = s[i];
                if (y + half + 1 < H)
                    for (int i = 0; i < xn; i++)
                        s[i] += row[static_cast<size_t>(y + half + 1) * W + x0 + i];
                if (y - half >= 0)
                    for (int i = 0; i < xn; i++)
                        s[i] -= row[static_cast<size_t>(y - half) * W + x0 + i];
            }
        }
    });
    return out;
}

// The values of each ring, laid end to end, with the ring's slice given by offset[r]..offset[r+1].
// radius and pooled do not change over the three baseline iterations, so this is built once and
// every iteration is an order statistic of the same, already-sorted, ring.
struct RingValues {
    std::vector<float> values;
    std::vector<int>   offset;
};

RingValues bin_by_ring(const std::vector<float> &values, const std::vector<char> &valid,
                       const std::vector<int> &radius, int max_radius, size_t nthreads) {
    RingValues rv;
    rv.offset.assign(max_radius + 2, 0);
    for (size_t i = 0; i < values.size(); i++)
        if (valid[i])
            rv.offset[radius[i] + 1]++;
    for (int r = 0; r <= max_radius; r++)
        rv.offset[r + 1] += rv.offset[r];

    rv.values.resize(rv.offset[max_radius + 1]);
    std::vector<int> cursor(rv.offset.begin(), rv.offset.end() - 1);
    for (size_t i = 0; i < values.size(); i++)
        if (valid[i])
            rv.values[cursor[radius[i]]++] = values[i];

    // Sorted once; the three iterations then only pick a rank and count a prefix.
    ParallelFor(max_radius + 1, nthreads, [&](int r) {
        std::sort(rv.values.begin() + rv.offset[r], rv.values.begin() + rv.offset[r + 1]);
    });
    return rv;
}

}  // namespace

ShadowFinder::ShadowFinder(const DiffractionExperiment &experiment, const PixelMask &mask)
    : width(static_cast<int>(experiment.GetXPixelsNumConv())),
      height(static_cast<int>(experiment.GetYPixelsNumConv())),
      beam_x(experiment.GetBeamX_pxl()),
      beam_y(experiment.GetBeamY_pxl()),
      pixel_mask(mask.GetMask(experiment)) {
    if (pixel_mask.size() != static_cast<size_t>(width) * height)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "ShadowFinder: pixel mask does not match the detector");
    SetShardCount(1);
#ifdef JFJOCH_USE_CUDA
    if (get_gpu_count() > 0) {
        const size_t npixels = static_cast<size_t>(width) * height;
        gpu_pending = std::async(std::launch::async, [npixels] {
            return std::make_unique<ShadowAccumulatorGPU>(npixels);
        });
    }
#endif
}

// A shard's accumulators are allocated when a frame is first added to it, not here: with a GPU they
// are never used at all, and on a 16 Mpx detector eight of them are 2.9 GB to allocate and clear -
// which measured 0.8 s of the pre-scan, all of it wasted.
void ShadowFinder::SetShardCount(size_t n) {
    shards.clear();
    shards.resize(std::max<size_t>(1, n));
}

ShadowFinder::Projection ShadowFinder::Reduce() const {
#ifdef JFJOCH_USE_CUDA
    // The device holds its own projection. Bring it back and let it take part in the fold below as
    // one more shard; when every frame went to the GPU it is the whole answer.
    Projection device;
    if (Gpu() && gpu->GetFrameCount() > 0) {
        gpu->Download(device.max_value, device.sum_value, device.valid_count);
        device.frames = gpu->GetFrameCount();
        bool host_empty = true;
        for (const auto &p : shards)
            host_empty = host_empty && (p.frames == 0);
        if (host_empty)
            return device;
    }
#endif

    // Only when that shard actually holds something: its accumulators are allocated on first use, so
    // an unused shard is empty rather than zeroed, and returning it would hand the callers below a
    // projection they index by pixel.
    if (shards.size() == 1 && shards[0].frames > 0
#ifdef JFJOCH_USE_CUDA
        && !(gpu && gpu->GetFrameCount() > 0)
#endif
    )
        return shards[0];

    Projection out;
    const size_t npixels = static_cast<size_t>(width) * height;
    out.max_value.assign(npixels, 0);
    out.sum_value.assign(npixels, 0);
    out.valid_count.assign(npixels, 0);
    for (const auto &p : shards)
        out.frames += p.frames;
#ifdef JFJOCH_USE_CUDA
    out.frames += device.frames;
#endif

    // Each worker owns a slice of the pixels and folds every shard into it. The sums and counts are
    // integers and a pixel is touched by one worker only, so the result is the same as folding them
    // one shard at a time on one thread - this is several hundred megabytes per shard and is limited
    // by memory rather than by arithmetic.
    const size_t nthreads = std::max<size_t>(1, std::min<size_t>(std::thread::hardware_concurrency(),
                                                                 shards.size() * 2));
    const size_t chunk = (npixels + nthreads - 1) / nthreads;
    std::vector<std::future<void>> futures;
    futures.reserve(nthreads);
    for (size_t t = 0; t < nthreads; t++) {
        const size_t lo = t * chunk, hi = std::min(npixels, lo + chunk);
        if (lo >= hi) break;
        futures.emplace_back(std::async(std::launch::async, [&, lo, hi] {
#ifdef JFJOCH_USE_CUDA
            const Projection *extra[1] = {&device};
            for (const auto *pp : extra) {
                const auto &p = *pp;
                if (p.frames == 0) continue;
                for (size_t i = lo; i < hi; i++) {
                    if (p.valid_count[i] == 0)
                        continue;
                    if (out.valid_count[i] == 0 || p.max_value[i] > out.max_value[i])
                        out.max_value[i] = p.max_value[i];
                    out.sum_value[i] += p.sum_value[i];
                    out.valid_count[i] += p.valid_count[i];
                }
            }
#endif
            for (const auto &p : shards) {
                if (p.frames == 0) continue;   // never used, and its accumulators were never allocated
                for (size_t i = lo; i < hi; i++) {
                    if (p.valid_count[i] == 0)
                        continue;
                    if (out.valid_count[i] == 0 || p.max_value[i] > out.max_value[i])
                        out.max_value[i] = p.max_value[i];
                    out.sum_value[i] += p.sum_value[i];
                    out.valid_count[i] += p.valid_count[i];
                }
            }
        }));
    }
    for (auto &f : futures) f.get();
    return out;
}

template<class T>
void ShadowFinder::Add(const T *ptr, Projection &p) {
    // 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 is a genuine saturated value and is kept, so a saturated
    // reflection still registers as bright.
    T masked;
    if constexpr (std::is_signed_v<T>)
        masked = std::numeric_limits<T>::min();
    else
        masked = std::numeric_limits<T>::max();

    for (size_t i = 0; i < p.max_value.size(); i++) {
        const T v = ptr[i];
        if (v == masked)
            continue;
        const int64_t vi = static_cast<int64_t>(v);
        if (p.valid_count[i] == 0 || vi > p.max_value[i])
            p.max_value[i] = vi;
        p.sum_value[i] += vi;
        p.valid_count[i]++;
    }
    p.frames++;
}

void ShadowFinder::AddImage(const DataMessage &data, std::vector<uint8_t> &buffer, size_t shard) {
    if (static_cast<size_t>(data.image.GetWidth()) * data.image.GetHeight()
        != static_cast<size_t>(width) * height)
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "ShadowFinder: image size does not match the detector");
    if (shard >= shards.size())
        throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                              "ShadowFinder: shard out of range");

#ifdef JFJOCH_USE_CUDA
    // One device, so the frames queue here - but each is only a chunk upload plus two kernels, and
    // nothing was decompressed on the host to get this far.
    if (ShadowAccumulatorGPU *acc = ShadowAccumulatorGPU::Supports(data.image) ? Gpu() : nullptr) {
        std::unique_lock ul(gpu_mutex);
        try {
            acc->Add(data.image);
            return;
        } catch (const std::exception &e) {
            spdlog::warn("Beam stop: GPU accumulate failed ({}), falling back to the host", e.what());
        }
    }
#endif

    Projection &p = shards[shard];
    if (p.max_value.empty()) {
        const size_t npixels = static_cast<size_t>(width) * height;
        p.max_value.assign(npixels, 0);
        p.sum_value.assign(npixels, 0);
        p.valid_count.assign(npixels, 0);
    }
    const auto ptr = data.image.GetUncompressedPtr(buffer);
    switch (data.image.GetMode()) {
        case CompressedImageMode::Int8:   Add(reinterpret_cast<const int8_t   *>(ptr), p); break;
        case CompressedImageMode::Uint8:  Add(reinterpret_cast<const uint8_t  *>(ptr), p); break;
        case CompressedImageMode::Int16:  Add(reinterpret_cast<const int16_t  *>(ptr), p); break;
        case CompressedImageMode::Uint16: Add(reinterpret_cast<const uint16_t *>(ptr), p); break;
        case CompressedImageMode::Int32:  Add(reinterpret_cast<const int32_t  *>(ptr), p); break;
        case CompressedImageMode::Uint32: Add(reinterpret_cast<const uint32_t *>(ptr), p); break;
        default:
            throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
                                  "ShadowFinder: unsupported image mode");
    }
}

#ifdef JFJOCH_USE_CUDA
ShadowAccumulatorGPU *ShadowFinder::Gpu() const {
    std::unique_lock ul(gpu_mutex);
    if (gpu_pending.valid()) {
        try {
            gpu = gpu_pending.get();
        } catch (const std::exception &e) {
            // A GPU that cannot hold the projection is not a reason to fail: the host path gives
            // the same answer, only slower.
            spdlog::warn("Beam stop: GPU projection unavailable ({}), accumulating on the host",
                         e.what());
            gpu.reset();
        }
    }
    return gpu.get();
}
#endif

uint32_t ShadowFinder::GetFrameCount() const {
    std::unique_lock ul(m);
    uint32_t frames = 0;
    for (const auto &p : shards) frames += p.frames;
#ifdef JFJOCH_USE_CUDA
    if (gpu) frames += gpu->GetFrameCount();
#endif
    return frames;
}

std::vector<float> ShadowFinder::GetMeanProjection() const {
    std::unique_lock ul(m);
    const Projection p = Reduce();
    const auto &sum_value = p.sum_value;
    const auto &valid_count = p.valid_count;

    std::vector<float> mean(static_cast<size_t>(width) * height, NAN);
    for (size_t i = 0; i < mean.size(); i++)
        if (valid_count[i] > 0 && pixel_mask[i] == 0)
            mean[i] = static_cast<float>(static_cast<double>(sum_value[i]) / valid_count[i]);
    return mean;
}

std::vector<uint32_t> ShadowFinder::GetMask(size_t nthreads) const {
    std::unique_lock ul(m);
    if (nthreads == 0)
        nthreads = std::max(1u, std::thread::hardware_concurrency());
    const Projection p = Reduce();
    const auto &max_value = p.max_value;
    const auto &sum_value = p.sum_value;
    const auto &valid_count = p.valid_count;
    const uint32_t frames = p.frames;

    const int W = width, H = height;
    const int n_pixels = W * H;

    std::vector<uint32_t> mask(n_pixels, 0);
    if (frames == 0)
        return mask;

    // mean projection, usable pixels and radius from the beam centre
    std::vector<float> mean(n_pixels, 0.0f);
    std::vector<char>  valid(n_pixels, 0);
    std::vector<int>   radius(n_pixels, 0);
    std::atomic<int> max_radius_atomic{0};
    ParallelChunks(H, nthreads, [&](int ylo, int yhi) {
        int local_max = 0;
        for (int y = ylo; y < yhi; y++)
            for (int x = 0; x < W; x++) {
                const int i = y * W + x;
                if (valid_count[i] > 0 && pixel_mask[i] == 0) {
                    mean[i] = static_cast<float>(static_cast<double>(sum_value[i]) / valid_count[i]);
                    valid[i] = 1;
                }
                const float dx = x - beam_x, dy = y - beam_y;
                radius[i] = static_cast<int>(std::lround(std::sqrt(dx * dx + dy * dy)));
                local_max = std::max(local_max, radius[i]);
            }
        // max is associative, so folding the per-worker maxima gives the same answer whatever
        // order they finish in.
        int prev = max_radius_atomic.load();
        while (prev < local_max && !max_radius_atomic.compare_exchange_weak(prev, local_max)) {}
    });
    const int max_radius = max_radius_atomic.load();


    // 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.
    // The count is a count: at most 25 pixels, so an integer box sum is exact and costs half the
    // memory of the floating-point one it replaces. The background itself stays in double - its
    // running sum adds and subtracts across a whole row, and in float the rounding of the two would
    // not cancel, which moves pixels across the shadow threshold below.
    std::vector<double>  num(n_pixels);
    std::vector<int32_t> den(n_pixels);
    ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) {
        for (int i = lo; i < hi; i++) {
            num[i] = valid[i] ? mean[i] : 0.0;
            den[i] = valid[i] ? 1 : 0;
        }
    });
    const auto pooled_sum = box_sum(num, W, H, POOL_PX, nthreads);
    const auto pooled_count = box_sum(den, W, H, POOL_PX, nthreads);
    std::vector<float> pooled(n_pixels, 0.0f);
    ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) {
        for (int i = lo; i < hi; i++)
            if (pooled_count[i] > 0)
                pooled[i] = static_cast<float>(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.
    //
    // The iteration only ever excludes pixels whose pooled background is below a cut, and dividing by
    // a positive baseline is monotone - so the excluded pixels of a ring are exactly the lowest ones,
    // and the pixels the next median is taken over are exactly the rest. Each iteration's median is
    // therefore an order statistic of the ring's values, which do not change: bin and sort the rings
    // once, then each iteration picks a rank and counts a prefix. That replaces nine full-image
    // passes (use / ratio / excluded, three times) with one, and three re-binnings with none.
    const RingValues rings = bin_by_ring(pooled, valid, radius, max_radius, nthreads);
    std::vector<float> baseline(max_radius + 1, 0.0f);
    {
        std::vector<int> excluded_in_ring(max_radius + 1, 0);
        for (int iter = 0; iter < 3; iter++)
            ParallelFor(max_radius + 1, nthreads, [&](int r) {
                const int lo = rings.offset[r], hi = rings.offset[r + 1];
                const int n = hi - lo;
                const int m = excluded_in_ring[r];
                const int avail = n - m;
                baseline[r] = (avail <= 0) ? 0.0f : rings.values[lo + m + avail / 2];
                // Counted the same way the per-pixel test below is written, so the two agree bit for
                // bit; the ring is sorted, so this is the length of a prefix.
                const float d = std::max(baseline[r], 1e-6f);
                int excl = 0;
                while (excl < n && rings.values[lo + excl] / d < SHADOW_RATIO)
                    excl++;
                excluded_in_ring[r] = excl;
            });
    }
    std::vector<float> ratio(n_pixels, 1.0f);
    ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) {
        for (int i = lo; i < hi; i++)
            if (valid[i])
                ratio[i] = pooled[i] / std::max(baseline[radius[i]], 1e-6f);
    });


    // A ring whose background was never counted carries no information to test a pixel against.
    // Walking outward, every ring before the first countable one lies wholly inside the stop - a
    // ring fully within the disk has no unshadowed pixel for the median to find, which is exactly
    // where an azimuthal comparison must fail. Those rings are shadow in their entirety.
    // Innermost rings hold only a handful of pixels, too few to judge, so they are stepped over
    // rather than allowed to end the walk.
    std::vector<int> ring_pixels(max_radius + 1, 0);
    for (int i = 0; i < n_pixels; i++)
        if (valid[i])
            ring_pixels[radius[i]]++;

    // A ring lies inside the stop when its background is a fraction of the background further out.
    // Counting statistics cannot decide this: on a bright dataset the shadow is still well counted.
    // The comparison is only ever used to answer "is this whole ring inside the stop", never to
    // judge an individual pixel, so taking the largest background over an outward window is safe
    // here in a way it would not be per pixel. It is taken only over the rings this same walk is
    // willing to judge, though: at the corner of the detector a ring holds a handful of pixels and
    // its median is one pixel's mean, so one recorded reflection out there would otherwise become
    // the background every ring inside it is compared against.
    std::vector<float> outward_max(max_radius + 2, 0.0f);
    for (int rad = max_radius; rad >= 0; rad--)
        outward_max[rad] = std::max(ring_pixels[rad] >= MIN_RING_PIXELS ? baseline[rad] : 0.0f,
                                    outward_max[rad + 1]);

    int blocked_out_to = -1;
    for (int rad = 0; rad <= max_radius; rad++) {
        if (ring_pixels[rad] < MIN_RING_PIXELS)
            continue;
        if (baseline[rad] >= SHADOW_RATIO * outward_max[rad])
            break;
        blocked_out_to = rad;
    }

    std::vector<char> low(n_pixels, 0);
    ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) {
        for (int i = lo; i < hi; i++) {
            if (!valid[i])
                continue;
            if (radius[i] <= blocked_out_to) {
                low[i] = 1;
                continue;
            }
            const double counted = static_cast<double>(frames) * pooled_count[i];
            low[i] = ratio[i] < SHADOW_RATIO && baseline[radius[i]] * counted >= MIN_EXPECTED_COUNTS;
        }
    });


    // The shadow is the low region connected to the beam centre, bridging the gaps it crosses.
    const std::vector<char> bridged = dilate(low, W, H, BRIDGE_PX, nthreads);
    std::vector<int> seeds;
    for (int i = 0; i < n_pixels; i++)
        if (radius[i] < 4)
            seeds.push_back(i);
    const std::vector<char> connected = flood(bridged, W, H, seeds);
    std::vector<char> 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<char> lit(n_pixels, 0);
    ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) {
        for (int i = lo; i < hi; i++)
            lit[i] = (valid_count[i] > 0) && (max_value[i] >= MIN_REFLECTION);
    });
    std::vector<char> reflection(n_pixels, 0);
    ParallelChunks(H, nthreads, [&](int ylo, int yhi) {
        for (int y = ylo; y < yhi; 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);
            }
    });


    // Grow the soft boundary, round it and fill the disk interior.
    const std::vector<char> penumbra = dilate(region, W, H, PENUMBRA_MAX_PX, nthreads);
    for (int i = 0; i < n_pixels; i++)
        if (penumbra[i] && valid[i] && ratio[i] < PENUMBRA_RATIO)
            region[i] = 1;


    region = erode(dilate(region, W, H, 2, nthreads), W, H, 2, nthreads);

    region = fill_holes(region, W, H);


    // 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<char> reflection_grown = dilate(reflection, W, H, 1, nthreads);
    for (int i = 0; i < n_pixels; i++)
        if (reflection_grown[i])
            region[i] = 0;

    for (int i = 0; i < n_pixels; i++)
        mask[i] = region[i] ? 1 : 0;
    return mask;
}