From 24e36ae740eb7900a9d920350ec7ad0f9984458a Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 3 Oct 2026 10:46:49 +0200 Subject: [PATCH] Pre-scan: parallel beam-stop mask, leaner background beam-centre fit Exact: p.mtz and the pre-scan products (shadow mask, mean projection, defective-pixel mask, ring and capture centres, compared as hashes and hex floats) are bit-identical to rc174 on three in-house rotation sets, GPU and CPU builds. - ShadowFinder::GetMask: the serial parts run in parallel - connected components by row band joined with union-find (both the shadow and the transmitting-arm searches, and the hole fill), ring binning and the harmonic sector gather by blocks, gap bridging by line; ring pixel counts read off the ring offsets. Mean projection filled in parallel. - ShadowFinder host accumulation: one band-locked projection instead of a 20 B/px shard per pre-scan worker (2.7 GB zeroed and folded on a 16M detector); SetShardCount and the shard argument are gone. - FindBeamCenterFromBackground: the usable-pixel test is made once, the in-band pixels are kept in pixel order so the clipping rounds no longer sweep the whole detector, the 67 MB cell map is gone and the per-iteration block fold runs in parallel - same sums, same order. - HotPixelFinder::GetMask: the chance-rate counts in parallel (integers). Measured on a loaded box (load ~25 from other jobs), pre-scan window: GPU 5.9-6.5 s -> 3.2-3.4 s, CPU 8.4-9.0 s -> 6.1-7.4 s. The GPU-build pre-scan now ends with its background spot measurement (CPU spot finder on ~120 frames, ~13 core-s on 8 workers). Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB --- image_analysis/beam_stop/ShadowFinder.cpp | 533 ++++++++++-------- image_analysis/beam_stop/ShadowFinder.h | 59 +- .../BeamCenterFromBackground.cpp | 82 ++- rugnux/HotPixels.cpp | 22 +- rugnux/Rugnux.cpp | 19 +- tests/ShadowFinderTest.cpp | 38 +- 6 files changed, 427 insertions(+), 326 deletions(-) diff --git a/image_analysis/beam_stop/ShadowFinder.cpp b/image_analysis/beam_stop/ShadowFinder.cpp index 715ad1198..95846c533 100644 --- a/image_analysis/beam_stop/ShadowFinder.cpp +++ b/image_analysis/beam_stop/ShadowFinder.cpp @@ -10,7 +10,6 @@ #include #include #include -#include #include #include @@ -79,7 +78,7 @@ constexpr int MIN_RING_PIXELS = 32; constexpr float BLOCKED_RING_RATIO = 0.35f; // 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). +// at GetMask() time, and all are O(pixels) rather than O(pixels * radius). namespace { // A per-pixel array of GetMask(). A std::vector zeroes what it allocates on the thread that makes it, @@ -186,10 +185,12 @@ Plane erode(const Plane &in, int W, int H, int r, size_t nthreads) { // continues on both sides of it is one shadow - but a gap can be wider than BRIDGE_PX reaches (17 px // between the rows of PILATUS modules), and an arm crossing one fell apart into pieces each too small // to be believed. -Plane bridge_gaps(const Plane ®ion, const Plane &valid, int W, int H) { +Plane bridge_gaps(const Plane ®ion, const Plane &valid, int W, int H, size_t nthreads) { Plane out = region; + // The lines of one direction are independent: each reads `region` and only ever sets its own pixels. auto walk = [&](int n_lines, int len, auto index) { - for (int line = 0; line < n_lines; line++) { + ParallelChunks(n_lines, nthreads, [&](int lo, int hi) { + for (int line = lo; line < hi; line++) { int k = 0; while (k < len) { if (valid[index(line, k)]) { k++; continue; } @@ -199,12 +200,99 @@ Plane bridge_gaps(const Plane ®ion, const Plane &valid, int for (int j = start; j < k; j++) out[index(line, j)] = 1; } } + }); }; walk(H, W, [W](int y, int x) { return static_cast(y) * W + x; }); walk(W, H, [W](int x, int y) { return static_cast(y) * W + x; }); return out; } +// The 8-connected components of `member`: each member pixel gets the index of its component, dense +// from 0 and in no particular order, and every other pixel -1. +// +// Labelled in parallel. Each band of rows is flooded on its own, then the pieces that touch across a +// band boundary are joined. Which pixels share a component is all a caller reads, and that does not +// depend on how the rows were split. +struct Components { + Plane id; + int count = 0; +}; + +Components label_components(const Plane &member, int W, int H, size_t nthreads) { + const int bands = std::min(64, H); // never more bands than rows, so none is empty + std::vector band_row(bands + 1); + for (int b = 0; b <= bands; b++) + band_row[b] = static_cast(static_cast(b) * H / bands); + + Components out; + out.id = Plane(member.size()); + std::vector band_pieces(bands, 0); + ParallelFor(bands, nthreads, [&](int b) { + const size_t lo = static_cast(band_row[b]) * W, hi = static_cast(band_row[b + 1]) * W; + std::fill(out.id.begin() + lo, out.id.begin() + hi, -1); + std::vector stack; + int pieces = 0; + for (size_t start = lo; start < hi; start++) { + if (!member[start] || out.id[start] >= 0) + continue; + out.id[start] = pieces; + stack.push_back(start); + while (!stack.empty()) { + const size_t i = stack.back(); stack.pop_back(); + const int y = static_cast(i / W), x = static_cast(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 < band_row[b] || yy >= band_row[b + 1] || xx < 0 || xx >= W) + continue; + const size_t j = static_cast(yy) * W + xx; + if (member[j] && out.id[j] < 0) { out.id[j] = pieces; stack.push_back(j); } + } + } + pieces++; + } + band_pieces[b] = pieces; + }); + + // A piece is named by its band's first index plus its number in the band, and the pieces are + // joined across each boundary row by union-find. + std::vector first(bands + 1, 0); + for (int b = 0; b < bands; b++) + first[b + 1] = first[b] + band_pieces[b]; + std::vector parent(first[bands]); + for (size_t k = 0; k < parent.size(); k++) + parent[k] = static_cast(k); + const auto find = [&](int k) { + while (parent[k] != k) { parent[k] = parent[parent[k]]; k = parent[k]; } + return k; + }; + for (int b = 0; b + 1 < bands; b++) { + const size_t above = static_cast(band_row[b + 1] - 1) * W, below = above + W; + for (int x = 0; x < W; x++) { + if (!member[above + x]) + continue; + for (int xx = std::max(0, x - 1); xx <= std::min(W - 1, x + 1); xx++) + if (member[below + xx]) { + const int ra = find(first[b] + out.id[above + x]); + const int rb = find(first[b + 1] + out.id[below + xx]); + if (ra != rb) parent[std::max(ra, rb)] = std::min(ra, rb); + } + } + } + std::vector dense(parent.size(), -1), component(parent.size()); + for (size_t k = 0; k < parent.size(); k++) { + const int root = find(static_cast(k)); + if (dense[root] < 0) dense[root] = out.count++; + component[k] = dense[root]; + } + ParallelFor(bands, nthreads, [&](int b) { + const size_t lo = static_cast(band_row[b]) * W, hi = static_cast(band_row[b + 1]) * W; + for (size_t i = lo; i < hi; i++) + if (out.id[i] >= 0) out.id[i] = component[first[b] + out.id[i]]; + }); + return out; +} + // 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 @@ -212,52 +300,53 @@ Plane bridge_gaps(const Plane ®ion, const Plane &valid, int // 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. -Plane fill_holes(const Plane ®ion, int W, int H) { +Plane fill_holes(const Plane ®ion, int W, int H, size_t nthreads) { + std::vector row_x0(H, W), row_x1(H, -1); + ParallelChunks(H, nthreads, [&](int ylo, int yhi) { + for (int y = ylo; y < yhi; y++) + for (int x = 0; x < W; x++) + if (region[static_cast(y) * W + x]) { + row_x0[y] = std::min(row_x0[y], x); + row_x1[y] = x; + } + }); 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(y) * W + x]) { - x0 = std::min(x0, x); x1 = std::max(x1, x); - y0 = std::min(y0, y); y1 = std::max(y1, y); - } + if (row_x1[y] >= 0) { + x0 = std::min(x0, row_x0[y]); x1 = std::max(x1, row_x1[y]); + y0 = std::min(y0, y); 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); + // The background of the box, in components; one that reaches the box's edge is outside. const int BW = x1 - x0 + 1, BH = y1 - y0 + 1; - std::vector bg_visited(static_cast(BW) * BH, 0); - std::queue q; // indices into the box - auto push = [&](int bx, int by) { - const int j = by * BW + bx; - if (!region[static_cast(by + y0) * W + bx + x0] && !bg_visited[j]) { - bg_visited[j] = 1; q.push(j); - } + Plane background(static_cast(BW) * BH); + ParallelChunks(BH, nthreads, [&](int lo, int hi) { + for (int by = lo; by < hi; by++) + for (int bx = 0; bx < BW; bx++) + background[static_cast(by) * BW + bx] = !region[static_cast(by + y0) * W + bx + x0]; + }); + const auto pieces = label_components(background, BW, BH, nthreads); + std::vector outside(pieces.count, 0); + const auto edge = [&](int bx, int by) { + const int c = pieces.id[static_cast(by) * BW + bx]; + if (c >= 0) outside[c] = 1; }; - 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(yy + y0) * W + xx + x0] && !bg_visited[j]) { - bg_visited[j] = 1; q.push(j); - } - } - } + for (int bx = 0; bx < BW; bx++) { edge(bx, 0); edge(bx, BH - 1); } + for (int by = 0; by < BH; by++) { edge(0, by); edge(BW - 1, by); } Plane out = region; - for (int by = 0; by < BH; by++) - for (int bx = 0; bx < BW; bx++) { - const size_t i = static_cast(by + y0) * W + bx + x0; - if (!region[i] && !bg_visited[by * BW + bx]) - out[i] = 1; - } + ParallelChunks(BH, nthreads, [&](int lo, int hi) { + for (int by = lo; by < hi; by++) + for (int bx = 0; bx < BW; bx++) { + const int c = pieces.id[static_cast(by) * BW + bx]; + if (c >= 0 && !outside[c]) + out[static_cast(by + y0) * W + bx + x0] = 1; + } + }); return out; } @@ -321,19 +410,40 @@ struct RingValues { RingValues bin_by_ring(const Plane &values, const Plane &valid, const Plane &radius, int max_radius, size_t nthreads) { + // Counted and scattered by blocks of pixels in parallel: each block writes its values of a ring + // after those of the blocks before it, so every ring holds its values in pixel order, as a single + // pass would leave them - and they are sorted below in any case. + constexpr int BLOCKS = 64; + const size_t n = values.size(); + const auto block_begin = [n](int b) { return n * b / BLOCKS; }; + const size_t rings = static_cast(max_radius) + 1; + std::vector cursor(BLOCKS * rings, 0); + ParallelFor(BLOCKS, nthreads, [&](int b) { + int *count = cursor.data() + b * rings; + for (size_t i = block_begin(b); i < block_begin(b + 1); i++) + if (valid[i]) + count[radius[i]]++; + }); + 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]; + for (size_t r = 0; r < rings; r++) { + int at = rv.offset[r]; + for (int b = 0; b < BLOCKS; b++) { + const int count = cursor[b * rings + r]; + cursor[b * rings + r] = at; + at += count; + } + rv.offset[r + 1] = at; + } rv.values.resize(rv.offset[max_radius + 1]); - std::vector 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]; + ParallelFor(BLOCKS, nthreads, [&](int b) { + int *next = cursor.data() + b * rings; + for (size_t i = block_begin(b); i < block_begin(b + 1); i++) + if (valid[i]) + rv.values[next[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) { @@ -355,7 +465,6 @@ ShadowFinder::ShadowFinder(const DiffractionExperiment &experiment, const PixelM if (pixel_mask.size() != static_cast(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(width) * height; @@ -372,98 +481,37 @@ void ShadowFinder::BeamCenter(float x, float y) { beam_y = y; } -// 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(1, n)); -} - +#ifdef JFJOCH_USE_CUDA 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]; - + // The device holds its own projection. Bring it back; when every frame went to the GPU it is the + // whole answer. Projection out; - const size_t npixels = static_cast(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(1, std::min(std::thread::hardware_concurrency(), - shards.size() * 2)); - const size_t chunk = (npixels + nthreads - 1) / nthreads; - std::vector> 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]; - } - } - })); + if (Gpu() && gpu->GetFrameCount() > 0) { + gpu->Download(out.max_value, out.sum_value, out.valid_count); + out.frames = gpu->GetFrameCount(); } - for (auto &f : futures) f.get(); + if (host.frames == 0) + return out; + + // A pixel is touched by one worker only and the sums and counts are integers, so the result is + // the same as folding on one thread. + out.frames += host.frames; + ParallelChunks(static_cast(out.max_value.size()), std::thread::hardware_concurrency(), [&](int lo, int hi) { + for (int i = lo; i < hi; i++) { + if (host.valid_count[i] == 0) + continue; + if (out.valid_count[i] == 0 || host.max_value[i] > out.max_value[i]) + out.max_value[i] = host.max_value[i]; + out.sum_value[i] += host.sum_value[i]; + out.valid_count[i] += host.valid_count[i]; + } + }); return out; } +#endif template -void ShadowFinder::Add(const T *ptr, Projection &p) { +void ShadowFinder::Add(const T *ptr, size_t begin, size_t end) { // 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 @@ -474,27 +522,23 @@ void ShadowFinder::Add(const T *ptr, Projection &p) { else masked = std::numeric_limits::max(); - for (size_t i = 0; i < p.max_value.size(); i++) { + for (size_t i = begin; i < end; i++) { const T v = ptr[i]; if (v == masked) continue; const int64_t vi = static_cast(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]++; + if (host.valid_count[i] == 0 || vi > host.max_value[i]) + host.max_value[i] = vi; + host.sum_value[i] += vi; + host.valid_count[i]++; } - p.frames++; } -void ShadowFinder::AddImage(const DataMessage &data, std::vector &buffer, size_t shard) { +void ShadowFinder::AddImage(const DataMessage &data, std::vector &buffer) { if (static_cast(data.image.GetWidth()) * data.image.GetHeight() != static_cast(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 @@ -511,25 +555,37 @@ void ShadowFinder::AddImage(const DataMessage &data, std::vector &buffe } #endif - Projection &p = shards[shard]; - if (p.max_value.empty()) { - const size_t npixels = static_cast(width) * height; - p.max_value.assign(npixels, 0); - p.sum_value.assign(npixels, 0); - p.valid_count.assign(npixels, 0); + const size_t npixels = static_cast(width) * height; + { + std::unique_lock ul(host_mutex); + if (host.max_value.empty()) { + host.max_value.resize(npixels); + host.sum_value.resize(npixels); + host.valid_count.resize(npixels); + } } const auto ptr = data.image.GetUncompressedPtr(buffer); - switch (data.image.GetMode()) { - case CompressedImageMode::Int8: Add(reinterpret_cast(ptr), p); break; - case CompressedImageMode::Uint8: Add(reinterpret_cast(ptr), p); break; - case CompressedImageMode::Int16: Add(reinterpret_cast(ptr), p); break; - case CompressedImageMode::Uint16: Add(reinterpret_cast(ptr), p); break; - case CompressedImageMode::Int32: Add(reinterpret_cast(ptr), p); break; - case CompressedImageMode::Uint32: Add(reinterpret_cast(ptr), p); break; - default: - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - "ShadowFinder: unsupported image mode"); + const size_t rows_per_band = (static_cast(height) + BANDS - 1) / BANDS; + const size_t first = next_band.fetch_add(1); + for (size_t b = 0; b < BANDS; b++) { + const size_t band = (first + b) % BANDS; + const size_t begin = std::min(npixels, band * rows_per_band * width); + const size_t end = std::min(npixels, (band + 1) * rows_per_band * width); + std::lock_guard lock(band_mutex[band]); + switch (data.image.GetMode()) { + case CompressedImageMode::Int8: Add(reinterpret_cast(ptr), begin, end); break; + case CompressedImageMode::Uint8: Add(reinterpret_cast(ptr), begin, end); break; + case CompressedImageMode::Int16: Add(reinterpret_cast(ptr), begin, end); break; + case CompressedImageMode::Uint16: Add(reinterpret_cast(ptr), begin, end); break; + case CompressedImageMode::Int32: Add(reinterpret_cast(ptr), begin, end); break; + case CompressedImageMode::Uint32: Add(reinterpret_cast(ptr), begin, end); break; + default: + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "ShadowFinder: unsupported image mode"); + } } + std::unique_lock ul(host_mutex); + host.frames++; } #ifdef JFJOCH_USE_CUDA @@ -552,8 +608,7 @@ ShadowAccumulatorGPU *ShadowFinder::Gpu() const { uint32_t ShadowFinder::GetFrameCount() const { std::unique_lock ul(m); - uint32_t frames = 0; - for (const auto &p : shards) frames += p.frames; + uint32_t frames = host.frames; #ifdef JFJOCH_USE_CUDA if (gpu) frames += gpu->GetFrameCount(); #endif @@ -561,14 +616,21 @@ uint32_t ShadowFinder::GetFrameCount() const { } const ShadowFinder::Projection &ShadowFinder::Reduced() const { - if (!reduced) - reduced = Reduce(); - return *reduced; +#ifdef JFJOCH_USE_CUDA + if (gpu && gpu->GetFrameCount() > 0) { + if (!reduced) + reduced = Reduce(); + return *reduced; + } +#endif + return host; } void ShadowFinder::ReleaseProjection() { +#ifdef JFJOCH_USE_CUDA std::unique_lock ul(m); reduced.reset(); +#endif } std::vector ShadowFinder::GetMeanProjection() const { @@ -577,10 +639,12 @@ std::vector ShadowFinder::GetMeanProjection() const { const auto &sum_value = p.sum_value; const auto &valid_count = p.valid_count; - std::vector mean(static_cast(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(static_cast(sum_value[i]) / valid_count[i]); + std::vector mean(static_cast(width) * height); + ParallelChunks(static_cast(mean.size()), std::thread::hardware_concurrency(), [&](int lo, int hi) { + for (int i = lo; i < hi; i++) + mean[i] = valid_count[i] > 0 && pixel_mask[i] == 0 + ? static_cast(static_cast(sum_value[i]) / valid_count[i]) : NAN; + }); return mean; } @@ -727,9 +791,8 @@ std::vector ShadowFinder::GetMask(size_t nthreads) const { // 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 ring_pixels(max_radius + 1, 0); - for (int i = 0; i < n_pixels; i++) - if (valid[i]) - ring_pixels[radius[i]]++; + for (int rad = 0; rad <= max_radius; rad++) + ring_pixels[rad] = rings.offset[rad + 1] - rings.offset[rad]; // A ring lies inside the stop when its background is a fraction of what this detector's // background typically is. Counting statistics cannot decide this: on a bright dataset the @@ -789,36 +852,20 @@ std::vector ShadowFinder::GetMask(size_t nthreads) const { // and the stop. What keeps the test specific instead is size, since the background wanders by a // pixel or two at a time and hardware does not. const Plane bridged = dilate(low, W, H, BRIDGE_PX, nthreads); - Plane region = filled_plane(n_pixels, 0, nthreads); + Plane region(n_pixels); { - Plane seen = filled_plane(n_pixels, 0, nthreads); - std::vector component; - std::queue q; - for (int start = 0; start < n_pixels; start++) { - if (!bridged[start] || seen[start]) - continue; - component.clear(); - int n_low = 0; - seen[start] = 1; - q.push(start); - while (!q.empty()) { - const int i = q.front(); q.pop(); - component.push_back(i); - n_low += low[i]; - 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 (bridged[j] && !seen[j]) { seen[j] = 1; q.push(j); } - } - } - if (n_low >= MIN_SHADOW_PIXELS) - for (const int i : component) - region[i] = low[i]; - } + const auto pieces = label_components(bridged, W, H, nthreads); + std::vector> n_low(pieces.count); + ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) { + for (int i = lo; i < hi; i++) + if (pieces.id[i] >= 0 && low[i]) + n_low[pieces.id[i]].fetch_add(1, std::memory_order_relaxed); + }); + ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) { + for (int i = lo; i < hi; i++) + region[i] = pieces.id[i] >= 0 && n_low[pieces.id[i]].load(std::memory_order_relaxed) >= MIN_SHADOW_PIXELS + ? low[i] : 0; + }); } // The rings that lie wholly inside the stop are decided by the ring walk above rather than by @@ -866,7 +913,7 @@ std::vector ShadowFinder::GetMask(size_t nthreads) const { region = erode(dilate(region, W, H, 2, nthreads), W, H, 2, nthreads); - region = fill_holes(region, W, H); + region = fill_holes(region, W, H, nthreads); // Expose recorded reflections - done last, with no fill afterwards, so a spot the shadow @@ -903,15 +950,31 @@ std::vector ShadowFinder::GetMask(size_t nthreads) const { // allowed to explain a dim sector away - where it would ask for more than the median, the median // stands - so the step can only drop what it found before, never find something new. const int n_bands = max_radius / HARMONIC_BAND_PX + 1; - std::vector> sector_values(static_cast(n_bands) * HARMONIC_SECTORS); - for (int i = 0; i < n_pixels; i++) { - if (!valid[i] || region[i]) - continue; - const float dx = static_cast(i % W) - beam_x, dy = static_cast(i / W) - beam_y; - const double phi = std::atan2(dy, dx) + std::numbers::pi; - const int sector = std::min(HARMONIC_SECTORS - 1, static_cast(phi / (2.0 * std::numbers::pi) * HARMONIC_SECTORS)); - sector_values[static_cast(radius[i] / HARMONIC_BAND_PX) * HARMONIC_SECTORS + sector].push_back(ratio[i]); - } + const size_t n_sectors = static_cast(n_bands) * HARMONIC_SECTORS; + // Gathered by blocks of rows in parallel and joined in block order. Only the median of each + // sector is read, and that is the same whatever order its values were gathered in. + constexpr int SECTOR_BLOCKS = 64; + std::vector>> block_values(SECTOR_BLOCKS); + ParallelFor(SECTOR_BLOCKS, nthreads, [&](int b) { + auto &values = block_values[b]; + values.resize(n_sectors); + const int lo = static_cast(static_cast(n_pixels) * b / SECTOR_BLOCKS); + const int hi = static_cast(static_cast(n_pixels) * (b + 1) / SECTOR_BLOCKS); + for (int i = lo; i < hi; i++) { + if (!valid[i] || region[i]) + continue; + const float dx = static_cast(i % W) - beam_x, dy = static_cast(i / W) - beam_y; + const double phi = std::atan2(dy, dx) + std::numbers::pi; + const int sector = std::min(HARMONIC_SECTORS - 1, static_cast(phi / (2.0 * std::numbers::pi) * HARMONIC_SECTORS)); + values[static_cast(radius[i] / HARMONIC_BAND_PX) * HARMONIC_SECTORS + sector].push_back(ratio[i]); + } + }); + std::vector> sector_values(n_sectors); + ParallelFor(static_cast(n_sectors), nthreads, [&](int k) { + for (const auto &values : block_values) + sector_values[k].insert(sector_values[k].end(), values[k].begin(), values[k].end()); + }); + block_values.clear(); std::vector harm_c(n_bands, 0.0f), harm_s(n_bands, 0.0f); ParallelFor(n_bands, nthreads, [&](int band) { std::vector med(HARMONIC_SECTORS, -1.0), c(HARMONIC_SECTORS), s(HARMONIC_SECTORS); @@ -978,36 +1041,24 @@ std::vector ShadowFinder::GetMask(size_t nthreads) const { && poisson_deficit_sigma(pooled[i] * counted, baseline[radius[i]] * model * counted) > MIN_DEFICIT_SIGMA; } }); - const Plane joined = bridge_gaps(dilate(dim, W, H, BRIDGE_PX, nthreads), valid, W, H); - Plane seen = filled_plane(n_pixels, 0, nthreads); - std::vector component; - std::queue q; - for (int start = 0; start < n_pixels; start++) { - if (!joined[start] || seen[start]) - continue; - component.clear(); - int n_dim = 0, n_low = 0; - seen[start] = 1; - q.push(start); - while (!q.empty()) { - const int i = q.front(); q.pop(); - component.push_back(i); - n_dim += dim[i]; - n_low += dim[i] && low[i]; - 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 (joined[j] && !seen[j]) { seen[j] = 1; q.push(j); } - } + const Plane joined = bridge_gaps(dilate(dim, W, H, BRIDGE_PX, nthreads), valid, W, H, nthreads); + const auto pieces = label_components(joined, W, H, nthreads); + std::vector> n_dim(pieces.count), n_low(pieces.count); + ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) { + for (int i = lo; i < hi; i++) + if (pieces.id[i] >= 0 && dim[i]) { + n_dim[pieces.id[i]].fetch_add(1, std::memory_order_relaxed); + if (low[i]) + n_low[pieces.id[i]].fetch_add(1, std::memory_order_relaxed); + } + }); + ParallelChunks(n_pixels, nthreads, [&](int lo, int hi) { + for (int i = lo; i < hi; i++) { + const int c = pieces.id[i]; + if (c >= 0 && dim[i] && n_dim[c].load(std::memory_order_relaxed) >= MIN_SHADOW_PIXELS + && n_low[c].load(std::memory_order_relaxed) >= MIN_CORE_PIXELS) + mask[i] = TRANSMITTING; } - if (n_dim >= MIN_SHADOW_PIXELS && n_low >= MIN_CORE_PIXELS) - for (const int i : component) - if (dim[i]) - mask[i] = TRANSMITTING; - } + }); return mask; } diff --git a/image_analysis/beam_stop/ShadowFinder.h b/image_analysis/beam_stop/ShadowFinder.h index 1b7e66ff2..866e3df4d 100644 --- a/image_analysis/beam_stop/ShadowFinder.h +++ b/image_analysis/beam_stop/ShadowFinder.h @@ -4,6 +4,7 @@ #pragma once #include +#include #include #include #include @@ -36,9 +37,9 @@ // // 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: workers call AddImage concurrently, each naming a shard of its own (see -// SetShardCount) - so no two threads touch the same accumulator and nothing is locked while -// an image is added. The shards are summed when the projection is read. +// Thread-safe: workers call AddImage concurrently. The projection is split into bands of rows, each +// with a lock of its own, and a worker adding a frame starts at a different band from the one before +// it, so workers meet only when they reach the same band. class ShadowFinder { mutable std::mutex m; @@ -62,22 +63,27 @@ class ShadowFinder { std::vector pixel_mask; // pixels already masked carry no background to test - // Per-pixel projection over the frames added so far (converted geometry). One set per shard: - // the sums and counts are integers, so summing the shards is exact and the result does not - // depend on how the frames were spread over them. + // Per-pixel projection over the frames added so far (converted geometry). The sums and counts are + // integers and the maximum is a maximum, so the result does not depend on the order the frames + // arrive in. struct Projection { std::vector max_value; std::vector sum_value; std::vector valid_count; uint32_t frames = 0; }; - std::vector shards; + // The frames added on the host. Allocated by the first of them: with a GPU there are usually none. + Projection host; + std::mutex host_mutex; // guards the allocation and the frame count, not the sums + static constexpr size_t BANDS = 64; + std::mutex band_mutex[BANDS]; + std::atomic next_band{0}; #ifdef JFJOCH_USE_CUDA // Present when a GPU is available. Frames it can decode are accumulated there instead of on the // host - only the compressed chunk crosses PCIe - and its projection is folded in with the - // shards when the mask is read. Frames it cannot take (anything but bitshuffle+LZ4) still go to - // a host shard, so a run mixing compressions is handled without a second code path. + // host projection when the mask is read. Frames it cannot take (anything but bitshuffle+LZ4) still + // go to the host, so a run mixing compressions is handled without a second code path. // Built on a thread of its own: it allocates and clears several hundred megabytes of device // memory, and cudaMalloc synchronises the whole device, so doing it in the constructor would // stall the caller before it has read its first frame. The first AddImage waits for it, by @@ -90,17 +96,22 @@ class ShadowFinder { [[nodiscard]] ShadowAccumulatorGPU *Gpu() const; #endif - template void Add(const T *ptr, Projection &p); + // Add the pixels [begin, end) of one frame to the host projection. + template void Add(const T *ptr, size_t begin, size_t end); - // Sum the shards into one projection. max_value is only taken from a shard that actually - // counted the pixel - a shard that never saw it holds 0, which would beat a genuinely +#ifdef JFJOCH_USE_CUDA + // The device's projection with the host's folded in. max_value is only taken from a projection + // that actually counted the pixel - one that never saw it holds 0, which would beat a genuinely // negative maximum. [[nodiscard]] Projection Reduce() const; - // The projection, reduced on its first read and kept. The ring-centre fit, the mask and the - // beam-centre capture all read the same one, and on a 16 Mpx detector each reduction is 360 MB - // brought back from the device into fresh memory. Called with `m` held. + // That projection, made on its first read and kept. The ring-centre fit, the mask and the + // beam-centre capture all read the same one, and on a 16 Mpx detector each is 360 MB brought back + // from the device into fresh memory. mutable std::optional reduced; +#endif + // The projection the frames added so far make: the host's, or the one above where the device + // took frames. Called with `m` held. [[nodiscard]] const Projection &Reduced() const; public: @@ -116,20 +127,16 @@ public: // hardware. The projection is not centred on anything, so this may be set after the frames. void BeamCenter(float x, float y); - // Give each worker a shard to accumulate into. Must be called before the first AddImage, - // and costs 20 bytes per pixel per shard. - void SetShardCount(size_t n); - - // Accumulate one full converted-geometry image into shard `shard`. Gap / masked pixels - // (the pixel type's sentinel extreme) are skipped. `buffer` is scratch space for - // decompression, reused across the calls of one worker. - void AddImage(const DataMessage &data, std::vector &buffer, size_t shard = 0); + // Accumulate one full converted-geometry image. Gap / masked pixels (the pixel type's sentinel + // extreme) are skipped. `buffer` is scratch space for decompression, reused across the calls of + // one worker. + void AddImage(const DataMessage &data, std::vector &buffer); // Compute the shadow mask (SHADOW, TRANSMITTING or 0 = keep), of the converted pixel count. // TRANSMITTING marks the pieces of hardware that let part of the beam through, added after the // shadow proper; both are masked, and a consumer that must not see those pieces can tell them apart. - // Recomputed on each call from the projection, which is summed over the shards on the first read - // of it (GetMask or GetMeanProjection): frames added after that are not seen. + // Recomputed on each call from the projection, which is put together on the first read of it + // (GetMask or GetMeanProjection): frames added after that are not seen. // nthreads = 0 asks for all hardware threads. The per-pixel passes over a 16M-pixel detector // dominate this, and they are all exactly parallel. [[nodiscard]] std::vector GetMask(size_t nthreads = 0) const; @@ -140,7 +147,7 @@ public: [[nodiscard]] std::vector GetMeanProjection() const; // Let go of the projection the two above read, once the caller has what it wants of it: on a - // 16 Mpx detector it is 360 MB. A later read sums the shards again. + // 16 Mpx detector it is 360 MB. A later read puts it together again. void ReleaseProjection(); [[nodiscard]] uint32_t GetFrameCount() const; diff --git a/image_analysis/geom_refinement/BeamCenterFromBackground.cpp b/image_analysis/geom_refinement/BeamCenterFromBackground.cpp index b2c306b9e..19f468a65 100644 --- a/image_analysis/geom_refinement/BeamCenterFromBackground.cpp +++ b/image_analysis/geom_refinement/BeamCenterFromBackground.cpp @@ -7,6 +7,7 @@ #include #include +#include "../../common/CompressedImage.h" #include "../../common/JFJochMath.h" #include "../../common/ParallelFor.h" @@ -111,7 +112,6 @@ FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const Pixe float beam_y = start ? start->second : geom.GetBeamY_pxl(); constexpr int n_cells = RADIAL_BINS * SECTORS; - std::vector cell_of(n_pixels); std::vector sum(n_cells), sum_sq(n_cells), sum_jx(n_cells), sum_jy(n_cells); std::vector count(n_cells), count_all(n_cells); std::vector profile(RADIAL_BINS), d_profile(RADIAL_BINS), clip_limit(n_cells); @@ -127,6 +127,38 @@ FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const Pixe std::vector block_jy(static_cast(BLOCKS) * n_cells); std::vector block_count(static_cast(BLOCKS) * n_cells); + // Whether a pixel can take part at all, which does not depend on the centre. + std::vector> usable(n_pixels); + ParallelFor(BLOCKS, nthreads, [&](int b) { + for (size_t i = static_cast(block_row[b]) * W; i < static_cast(block_row[b + 1]) * W; i++) + usable[i] = pixel_mask[i] == 0 && std::isfinite(mean[i]); + }); + // The pixels of each block that fall in the band at the current centre, with their cell and value, + // in pixel order from the block's first pixel on. The clipping rounds read these instead of the + // whole detector, in the same order. + std::vector> band_cell(n_pixels); + std::vector> band_value(n_pixels); + std::vector band_pixels(BLOCKS); + + // The blocks' cells folded in block order. Each cell is folded on its own, so the cells are split + // over the threads and every cell is still summed in the same order. + const auto fold = [&](bool with_jacobian) { + ParallelChunks(n_cells, nthreads, [&](int c0, int c1) { + for (int c = c0; c < c1; c++) { + double s = 0, ss = 0, jx = 0, jy = 0; + int32_t n = 0; + for (int b = 0; b < BLOCKS; b++) { + const size_t k = static_cast(b) * n_cells + c; + s += block_sum[k]; ss += block_sum_sq[k]; + if (with_jacobian) { jx += block_jx[k]; jy += block_jy[k]; } + n += block_count[k]; + } + sum[c] = s; sum_sq[c] = ss; count[c] = n; + if (with_jacobian) { sum_jx[c] = jx; sum_jy[c] = jy; } + } + }); + }; + // Most pixels lie outside the band. Those clearly outside it in tan(2theta) = rho / lz - by a // margin far above float rounding - skip the square root and both atan2 below; every pixel the exact // test would keep still reaches it. @@ -148,11 +180,13 @@ FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const Pixe std::fill(b_jx, b_jx + n_cells, 0.0); std::fill(b_jy, b_jy + n_cells, 0.0); std::fill(b_count, b_count + n_cells, 0); + int32_t *cells = band_cell.data() + static_cast(block_row[b]) * W; + float *values = band_value.data() + static_cast(block_row[b]) * W; + size_t n_band = 0; for (int y = block_row[b]; y < block_row[b + 1]; y++) { for (int x = 0; x < W; x++) { const size_t i = static_cast(y) * W + x; - cell_of[i] = -1; - if (pixel_mask[i] != 0 || !std::isfinite(mean[i])) + if (!usable[i]) continue; const float u = (x - beam_x) * pixel_size; const float v = (y - beam_y) * pixel_size; @@ -183,7 +217,9 @@ FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const Pixe const float g_x = lz * lx / (rho * denominator); const float g_y = lz * ly / (rho * denominator); const float g_z = -rho / denominator; - cell_of[i] = cell; + cells[n_band] = cell; + values[n_band] = mean[i]; + n_band++; b_count[cell]++; b_sum[cell] += mean[i]; b_sum_sq[cell] += static_cast(mean[i]) * mean[i]; @@ -191,18 +227,9 @@ FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const Pixe b_jy[cell] += -pixel_size * (g_x * rot[1] + g_y * rot[4] + g_z * rot[7]); } } + band_pixels[b] = n_band; }); - for (int c = 0; c < n_cells; c++) { - double s = 0, ss = 0, jx = 0, jy = 0; - int32_t n = 0; - for (int b = 0; b < BLOCKS; b++) { - const size_t k = static_cast(b) * n_cells + c; - s += block_sum[k]; ss += block_sum_sq[k]; - jx += block_jx[k]; jy += block_jy[k]; - n += block_count[k]; - } - sum[c] = s; sum_sq[c] = ss; sum_jx[c] = jx; sum_jy[c] = jy; count[c] = n; - } + fold(true); count_all = count; // the Jacobian sums belong to the unclipped pixel set @@ -220,26 +247,19 @@ FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const Pixe std::fill(b_sum, b_sum + n_cells, 0.0); std::fill(b_sum_sq, b_sum_sq + n_cells, 0.0); std::fill(b_count, b_count + n_cells, 0); - const size_t lo = static_cast(block_row[b]) * W; - const size_t hi = static_cast(block_row[b + 1]) * W; - for (size_t i = lo; i < hi; i++) { - const int32_t c = cell_of[i]; - if (c < 0 || clip_limit[c] < 0.0f || mean[i] > clip_limit[c]) + const int32_t *cells = band_cell.data() + static_cast(block_row[b]) * W; + const float *values = band_value.data() + static_cast(block_row[b]) * W; + for (size_t j = 0; j < band_pixels[b]; j++) { + const int32_t c = cells[j]; + const float value = values[j]; + if (clip_limit[c] < 0.0f || value > clip_limit[c]) continue; b_count[c]++; - b_sum[c] += mean[i]; - b_sum_sq[c] += static_cast(mean[i]) * mean[i]; + b_sum[c] += value; + b_sum_sq[c] += static_cast(value) * value; } }); - for (int c = 0; c < n_cells; c++) { - double s = 0, ss = 0; - int32_t n = 0; - for (int b = 0; b < BLOCKS; b++) { - const size_t k = static_cast(b) * n_cells + c; - s += block_sum[k]; ss += block_sum_sq[k]; n += block_count[k]; - } - sum[c] = s; sum_sq[c] = ss; count[c] = n; - } + fold(false); } // Radial profile: the median over the sectors that have a mean, on rings that are diff --git a/rugnux/HotPixels.cpp b/rugnux/HotPixels.cpp index 0af35f712..5b898f8ca 100644 --- a/rugnux/HotPixels.cpp +++ b/rugnux/HotPixels.cpp @@ -281,11 +281,25 @@ HotPixelFinder::Result HotPixelFinder::GetMask(double oscillation_deg, double sp // The chance rate per ring, from the pixels lit on no more than half of their frames: whatever // lights those - reflections, zingers, noise above the bound - lights a defect-free pixel too. + // Counted in integers by blocks of rows in parallel, so the totals do not depend on the split. + std::vector> block_lit(BANDS), block_seen(BANDS); + const size_t rows_per_band = (height + BANDS - 1) / BANDS; + ParallelFor(static_cast(BANDS), nthreads, [&](int b) { + block_lit[b].assign(nrings, 0); + block_seen[b].assign(nrings, 0); + const size_t begin = std::min(width * height, b * rows_per_band * width); + const size_t end = std::min(width * height, (b + 1) * rows_per_band * width); + for (size_t i = begin; i < end; i++) + if (key[i] >= 0 && n_valid(i) > 0 && 2 * n_lit[i] <= n_valid(i)) { + block_lit[b][key[i] / SECTORS] += n_lit[i]; + block_seen[b][key[i] / SECTORS] += n_valid(i); + } + }); std::vector lit(nrings, 0.0), seen(nrings, 0.0); - for (size_t i = 0; i < width * height; i++) - if (key[i] >= 0 && n_valid(i) > 0 && 2 * n_lit[i] <= n_valid(i)) { - lit[key[i] / SECTORS] += n_lit[i]; - seen[key[i] / SECTORS] += n_valid(i); + for (int r = 0; r < nrings; r++) + for (size_t b = 0; b < BANDS; b++) { + lit[r] += static_cast(block_lit[b][r]); + seen[r] += static_cast(block_seen[b][r]); } std::vector k_chance(nrings, n + 1); for (int r = 0; r < nrings; r++) diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index ca1dab038..d99d8d8e1 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -244,10 +244,10 @@ namespace { // as signal. The margin is capped at a tenth of the sweep so a short run still has a sample. constexpr int PRESCAN_END_MARGIN_IMAGES = 5; - // Workers reading the pre-scan sample. Each owns a shard of the beam-stop projection so no two - // threads touch the same accumulator, and a shard costs 20 bytes per pixel - 362 MB on a 16M - // detector - so this is capped well below the worker count of the run proper. The accumulation - // is memory-bound rather than compute-bound, so a handful of workers already saturates it. + // Workers reading the pre-scan sample. Each holds detector-sized buffers of its own, and pages of + // fresh memory are slow to fault in when many threads do it at once, so this is capped well below + // the worker count of the run proper. The beam-stop projection is memory-bound rather than + // compute-bound, so a handful of workers already saturates it. constexpr size_t PRESCAN_MAX_WORKERS = 8; // The spot width is measured on a GROWING share of the pre-scan sample: every eighth frame of @@ -1014,9 +1014,9 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // Read the sample on several workers. The reader serialises on the HDF5 lock, but the // decompression, the projection and the spot finding - which is all of the cost on a large - // detector - run in parallel. Each worker accumulates into a shard of its own, so nothing is - // locked while an image is added, and the per-frame results are stitched together in sample - // order below so the beam centre sees the same input however the workers interleaved. + // detector - run in parallel. The projection is integer sums, so the order the workers add their + // frames in does not reach it, and the per-frame results are stitched together in sample order + // below so the beam centre sees the same input however the workers interleaved. // // Two passes over the sample. The first builds the projection, which is all the shadow, the // defective pixels and the beam-centre capture below read; the second finds the spots, for the @@ -1026,13 +1026,12 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru const std::vector ordinals(sample.begin(), sample.end()); const size_t nworkers = std::min(std::max(config_.nthreads, 1), std::min(PRESCAN_MAX_WORKERS, ordinals.size())); - finder.SetShardCount(nworkers); { std::atomic next{0}; std::vector> futures; futures.reserve(nworkers); for (size_t t = 0; t < nworkers; t++) - futures.emplace_back(std::async(std::launch::async, [&, t] { + futures.emplace_back(std::async(std::launch::async, [&] { std::vector shadow_buffer; JFJochReaderRawImage raw_image; for (size_t i = next.fetch_add(1); i < ordinals.size(); i = next.fetch_add(1)) { @@ -1057,7 +1056,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru msg.image = raw_image.image; msg.number = ordinal; msg.original_number = image_idx; - finder.AddImage(msg, shadow_buffer, t); + finder.AddImage(msg, shadow_buffer); } } })); diff --git a/tests/ShadowFinderTest.cpp b/tests/ShadowFinderTest.cpp index 99473260e..cc89345f4 100644 --- a/tests/ShadowFinderTest.cpp +++ b/tests/ShadowFinderTest.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include "../common/DetectorSetup.h" @@ -168,16 +169,15 @@ TEST_CASE("ShadowFinder_MaskDoesNotDependOnTheThreadCount", "[ShadowFinder]") { CHECK(finder.GetMask(8) == one); } -// Workers accumulate into shards of their own and the shards are summed when the projection is read, -// so which worker saw which frame must not reach the answer - including the maximum, which only one -// shard holds when the reflection is on a single frame. -TEST_CASE("ShadowFinder_ShardingDoesNotChangeTheProjection", "[ShadowFinder]") { +// Workers add their frames concurrently, each starting at a different band of the projection, so +// which worker added which frame, and in what order, must not reach the answer - including the +// maximum, which one frame alone holds when the reflection is on a single frame. +TEST_CASE("ShadowFinder_ConcurrentWorkersDoNotChangeTheProjection", "[ShadowFinder]") { const DiffractionExperiment x = TestExperiment(); const PixelMask pixel_mask(x); ShadowFinder serial(x, pixel_mask); - ShadowFinder sharded(x, pixel_mask); - sharded.SetShardCount(4); + ShadowFinder concurrent(x, pixel_mask); std::vector> frames; std::vector buffer; @@ -185,22 +185,32 @@ TEST_CASE("ShadowFinder_ShardingDoesNotChangeTheProjection", "[ShadowFinder]") { frames.push_back(Scene(/*cross=*/false, /*reflection=*/f == 0)); DataMessage msg{}; msg.image = CompressedImage(frames.back(), W, H); - serial.AddImage(msg, buffer, 0); - sharded.AddImage(msg, buffer, static_cast(f) % 4); + serial.AddImage(msg, buffer); } + std::vector workers; + for (int t = 0; t < 4; t++) + workers.emplace_back([&, t] { + std::vector worker_buffer; + for (int f = NFRAMES - 1 - t; f >= 0; f -= 4) { + DataMessage msg{}; + msg.image = CompressedImage(frames[f], W, H); + concurrent.AddImage(msg, worker_buffer); + } + }); + for (auto &w : workers) w.join(); - CHECK(serial.GetFrameCount() == sharded.GetFrameCount()); + CHECK(serial.GetFrameCount() == concurrent.GetFrameCount()); const auto a = serial.GetMeanProjection(); - const auto b = sharded.GetMeanProjection(); + const auto b = concurrent.GetMeanProjection(); REQUIRE(a.size() == b.size()); // NAN marks a pixel nothing counted, and NAN != NAN, so compare the bits rather than the values. CHECK(memcmp(a.data(), b.data(), a.size() * sizeof(float)) == 0); - // The reflection is on one frame, so its maximum lives in a single shard. If the fold lost it, - // the mask would swallow the reflection instead of giving it back. - CHECK(serial.GetMask(1) == sharded.GetMask(1)); - CHECK(sharded.GetMask(1)[I(C - 14, C - 2)] == 0); + // The reflection is on one frame, so only that frame's maximum sees it. If it were lost, the mask + // would swallow the reflection instead of giving it back. + CHECK(serial.GetMask(1) == concurrent.GetMask(1)); + CHECK(concurrent.GetMask(1)[I(C - 14, C - 2)] == 0); } // Four opaque arms and a centred disk: the scene is invariant under a quarter turn, so the mask must