mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 68 Wiki Activity
Files
6368c001739e74dca22cb3996a0fa43e06a084cb
Jungfraujoch/image_analysis/beam_stop/ShadowFinder.cpp
T
jungfrauandClaude Opus 5 6368c00173 Decode and accumulate the beam-stop projection on the GPU 
The pre-scan decompressed its frames on the host and folded them into a
per-pixel projection there. On a 16M-pixel detector that is 60 frames of 72 MB
to decompress and 20 bytes per pixel to read and write back per frame - about
40 GB of memory traffic - and it was the whole cost of the phase once the mask
was no longer the bottleneck.

Only the compressed chunk crosses PCIe now. BSLZ4DecoderGPU already exposes the
raw decoded bytes (Decode(), the path its own tests use), which is what this
needs: the projection is defined on the RAW STORED COUNTS with the pixel type's
sentinel skipped, not on the preprocessed image, so nothing here goes through
the preprocessor. Sums, maxima and counts are integers, so the device result is
identical to the host's rather than merely close.

Frames are folded in batches of four. The fold reads and writes the whole
accumulator whatever the batch holds, so per frame it was spending most of the
bandwidth on the accumulator rather than on the data; four is where that stops
mattering, and every frame beyond it is another full frame of device memory,
which costs more in cudaMalloc - device-synchronizing - than it saves.

The accumulator is built on a thread of its own. It allocates and clears
several hundred megabytes, and doing that in the constructor stalled the caller
before it had read its first frame.

Frames the device cannot take - anything but bitshuffle+LZ4 - still go to a host
shard, so a run mixing compressions needs no second code path, and a build
without CUDA is unchanged.

RotationScaleMergeGPU set the CUDA device in its constructor and never put it
back. CUDA's current device is per-thread, so that silently re-pinned the
calling thread for the rest of its life, and the destructor freed several
gigabytes against whatever device happened to be current by then - CudaDevicePtr
records no device of its own. Every entry point now sets the device on entry and
restores it on exit.

ParallelFor/ParallelChunks moved to common/ParallelFor.h; two files had copies
and a third wants them.

Measured on a 16M-pixel rotation dataset: pre-scan 4.78 s -> 2.37 s -> ~2.0 s,
shadow unchanged at 139126 pixels (22143 on a 2M-pixel dataset). Full 24-crystal
battery: same space group on all 24, none failed, 15m32s -> 14m49s.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
2026-08-15 18:38:26 -04:00

712 lines
30 KiB
C++
Raw Blame History

// 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
}
void ShadowFinder::SetShardCount(size_t n) {
const size_t npixels = static_cast<size_t>(width) * height;
shards.clear();
shards.resize(std::max<size_t>(1, n));
for (auto &p : shards) {
p.max_value.assign(npixels, 0);
p.sum_value.assign(npixels, 0);
p.valid_count.assign(npixels, 0);
p.frames = 0;
}
}
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
if (shards.size() == 1
#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)
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
if (ShadowAccumulatorGPU::Supports(data.image) && Gpu()) {
// 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.
ShadowAccumulatorGPU *acc = Gpu();
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];
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;
}
Reference in New Issue View Git Blame Copy Permalink