From 70e11d9ab8fea02fcfb6bf337586cf57da56c36b Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 12 Oct 2026 08:07:10 +0200 Subject: [PATCH] rugnux pre-scan: projection and defective pixels per card, the sample read once, beam stop memoised The beam-stop projection (ShadowAccumulatorGPU) and the defective-pixel sums (HotPixelFinderGPU) are now kept one per card, fed by workers pinned t % gpus like the other pre-scan workers, and folded into the first card when they are read (cudaMemcpyPeer + an integer add / max-with-count kernel). All of them are integer sums or a maximum with the first-value rule, so the fold is order-free and the result is identical to one card's. With one card nothing is folded. The pre-scan sample is read once and kept (compressed as stored; up to 2 GiB, a frame past that is read again where it is wanted) for the spot measurement and the defective pixels, which read every frame again each before - a CBF/marCCD/SMV frame was decoded on the host three times. A beam-stop detection that finds no shadow leaves the mask as it was, so every probe copy of the run repeated it on the same inputs. Its answer (shadow, ring centre, mean projection) is now kept in a memo shared with the copies, keyed on the sampled frames, the geometry the rings are drawn in, the polarization and the binary mask; up to 4 entries. Debug TIMING lines for the pre-scan's parts (read and projection, beam stop, defective pixels, spot measurement). New tests: two accumulators/finders folded match one fed every frame. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- .../beam_stop/ShadowAccumulatorGPU.cu | 38 +- .../beam_stop/ShadowAccumulatorGPU.h | 7 + image_analysis/beam_stop/ShadowFinder.cpp | 81 +++- image_analysis/beam_stop/ShadowFinder.h | 26 +- image_analysis/hot_pixels/HotPixels.cpp | 53 ++- image_analysis/hot_pixels/HotPixels.h | 16 +- image_analysis/hot_pixels/HotPixelsGPU.cu | 32 ++ image_analysis/hot_pixels/HotPixelsGPU.h | 8 + rugnux/Rugnux.h | 28 +- rugnux/RugnuxPreScan.cpp | 416 +++++++++++------- tests/HotPixelFinderTest.cpp | 73 +++ tests/ShadowFinderTest.cpp | 46 ++ 12 files changed, 626 insertions(+), 198 deletions(-) diff --git a/image_analysis/beam_stop/ShadowAccumulatorGPU.cu b/image_analysis/beam_stop/ShadowAccumulatorGPU.cu index 1499f1fec..dff15e593 100644 --- a/image_analysis/beam_stop/ShadowAccumulatorGPU.cu +++ b/image_analysis/beam_stop/ShadowAccumulatorGPU.cu @@ -49,6 +49,22 @@ __global__ void accumulate_kernel(const T *__restrict__ frames, int nframes, siz } } +// Another card's projection folded into this one, with the update Add makes: the maximum is only taken +// from a projection that counted the pixel. +__global__ void fold_kernel(const int64_t *__restrict__ other_max, const int64_t *__restrict__ other_sum, + const uint32_t *__restrict__ other_count, int64_t *__restrict__ max_value, + int64_t *__restrict__ sum_value, uint32_t *__restrict__ valid_count, size_t npixels) { + for (size_t i = blockIdx.x * static_cast(blockDim.x) + threadIdx.x; i < npixels; + i += static_cast(blockDim.x) * gridDim.x) { + if (other_count[i] == 0) + continue; + if (valid_count[i] == 0 || other_max[i] > max_value[i]) + max_value[i] = other_max[i]; + sum_value[i] += other_sum[i]; + valid_count[i] += other_count[i]; + } +} + template void launch(const uint8_t *raw, int nframes, size_t frame_bytes, int64_t *max_value, int64_t *sum_value, uint32_t *valid_count, size_t npixels, int blocks, @@ -76,7 +92,6 @@ ShadowAccumulatorGPU::ShadowAccumulatorGPU(size_t npixels) cuda_err(cudaMemsetAsync(gpu_sum, 0, sizeof(int64_t) * npixels, *stream)); cuda_err(cudaMemsetAsync(gpu_count, 0, sizeof(uint32_t) * npixels, *stream)); - int device = 0; cuda_err(cudaGetDevice(&device)); cudaDeviceProp prop{}; cuda_err(cudaGetDeviceProperties(&prop, device)); @@ -176,6 +191,27 @@ void ShadowAccumulatorGPU::Add(const CompressedImage &image) { FoldPending(); } +void ShadowAccumulatorGPU::Fold(ShadowAccumulatorGPU &other) { + { + // The other card's last batch is folded on its own stream, with its card current. + ScopedGpuPin pin(other.device); + other.FoldPending(); + cuda_err(cudaStreamSynchronize(*other.stream)); + } + FoldPending(); + CudaDevicePtr other_max(npixels), other_sum(npixels); + CudaDevicePtr other_count(npixels); + cuda_err(cudaMemcpyPeerAsync(other_max, device, other.gpu_max, other.device, sizeof(int64_t) * npixels, *stream)); + cuda_err(cudaMemcpyPeerAsync(other_sum, device, other.gpu_sum, other.device, sizeof(int64_t) * npixels, *stream)); + cuda_err(cudaMemcpyPeerAsync(other_count, device, other.gpu_count, other.device, sizeof(uint32_t) * npixels, + *stream)); + fold_kernel<<>>(other_max, other_sum, other_count, gpu_max, gpu_sum, gpu_count, npixels); + cuda_err(cudaGetLastError()); + // The temporaries are pooled: they go once the stream has used them (see CudaDevicePtr). + cuda_err(cudaStreamSynchronize(*stream)); + frames += other.frames; +} + void ShadowAccumulatorGPU::Download(std::vector &max_value, std::vector &sum_value, std::vector &valid_count) { FoldPending(); diff --git a/image_analysis/beam_stop/ShadowAccumulatorGPU.h b/image_analysis/beam_stop/ShadowAccumulatorGPU.h index 4abeaf14f..af387557d 100644 --- a/image_analysis/beam_stop/ShadowAccumulatorGPU.h +++ b/image_analysis/beam_stop/ShadowAccumulatorGPU.h @@ -24,6 +24,7 @@ // Only bitshuffle+LZ4 images are handled; anything else is left to the host path. class ShadowAccumulatorGPU { const size_t npixels; + int device = 0; // the card it was built on, where everything below lives std::shared_ptr stream; std::unique_ptr decoder; @@ -61,6 +62,12 @@ public: void Add(const CompressedImage &image); [[nodiscard]] uint32_t GetFrameCount() const { return frames; } + [[nodiscard]] int GetDevice() const { return device; } + + // Add the projection another card accumulated to this one. Integer sums and a maximum taken as + // Add takes it, so the result does not depend on how the frames were split between the cards or + // on the order the cards are folded in. Called with this accumulator's card current. + void Fold(ShadowAccumulatorGPU &other); // The beam-stop mask and the mean projection, made where the projection is (ShadowMaskGPU.h), so // that it does not have to come back at all. diff --git a/image_analysis/beam_stop/ShadowFinder.cpp b/image_analysis/beam_stop/ShadowFinder.cpp index d1e0170c1..bd4a365cb 100644 --- a/image_analysis/beam_stop/ShadowFinder.cpp +++ b/image_analysis/beam_stop/ShadowFinder.cpp @@ -412,15 +412,30 @@ ShadowFinder::ShadowFinder(const DiffractionExperiment &experiment, const PixelM throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "ShadowFinder: pixel mask does not match the detector"); #ifdef JFJOCH_USE_CUDA - if (get_gpu_count() > 0) { + const int32_t cards = get_gpu_count(); + card_mutex = std::make_unique(std::max(cards, 1)); + for (int32_t c = 0; c < cards; c++) { const size_t npixels = static_cast(width) * height; - gpu_pending = std::async(std::launch::async, [npixels] { + gpu_pending.push_back(std::async(std::launch::async, [npixels, c] { + set_gpu(c); return std::make_unique(npixels); - }); + })); } #endif } +ShadowFinder::~ShadowFinder() { +#ifdef JFJOCH_USE_CUDA + // Each accumulator goes with its own card current: its pooled buffers are freed there. + Gpu(0); + for (auto &g : gpus) + if (g) { + ScopedGpuPin pin(g->GetDevice()); + g.reset(); + } +#endif +} + void ShadowFinder::BeamCenter(float x, float y) { std::unique_lock ul(m); beam_x = x; @@ -432,7 +447,8 @@ ShadowFinder::Projection ShadowFinder::Reduce() const { // The device holds its own projection. Bring it back; when every frame went to the GPU it is the // whole answer. Projection out; - if (Gpu() && gpu->GetFrameCount() > 0) { + if (ShadowAccumulatorGPU *gpu = Gpu(); gpu && gpu->GetFrameCount() > 0) { + ScopedGpuPin pin(gpu->GetDevice()); gpu->Download(out.max_value, out.sum_value, out.valid_count); out.frames = gpu->GetFrameCount(); } @@ -487,10 +503,13 @@ void ShadowFinder::AddImage(const DataMessage &data, std::vector &buffe "ShadowFinder: image size does not match the detector"); #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); + // On the card of the calling thread. The frames of one card queue here - but each is only a chunk + // upload plus two kernels, and nothing was decompressed on the host to get this far. + int device = 0; + if (get_gpu_count() > 0) + cudaGetDevice(&device); + if (ShadowAccumulatorGPU *acc = ShadowAccumulatorGPU::Supports(data.image) ? Gpu(device) : nullptr) { + std::unique_lock ul(card_mutex[device]); try { acc->Add(data.image); return; @@ -536,20 +555,43 @@ void ShadowFinder::AddImage(const DataMessage &data, std::vector &buffe } #ifdef JFJOCH_USE_CUDA -ShadowAccumulatorGPU *ShadowFinder::Gpu() const { +ShadowAccumulatorGPU *ShadowFinder::Gpu(int device) const { std::unique_lock ul(gpu_mutex); - if (gpu_pending.valid()) { + for (auto &pending : gpu_pending) { try { - gpu = gpu_pending.get(); + gpus.push_back(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(); + gpus.emplace_back(); } } - return gpu.get(); + gpu_pending.clear(); + return device < static_cast(gpus.size()) ? gpus[device].get() : nullptr; +} + +ShadowAccumulatorGPU *ShadowFinder::Gpu() const { + Gpu(0); + std::unique_lock ul(gpu_mutex); + ShadowAccumulatorGPU *into = nullptr; + for (auto &g : gpus) { + if (!g || g->GetFrameCount() == 0) + continue; + if (!into) { + into = g.get(); + continue; + } + { + ScopedGpuPin pin(into->GetDevice()); + into->Fold(*g); + } + // Folded in, so it must not be counted again - and its memory is no longer needed. + ScopedGpuPin pin(g->GetDevice()); + g.reset(); + } + return into; } #endif @@ -557,14 +599,16 @@ uint32_t ShadowFinder::GetFrameCount() const { std::unique_lock ul(m); uint32_t frames = host.frames; #ifdef JFJOCH_USE_CUDA - if (gpu) frames += gpu->GetFrameCount(); + std::unique_lock gl(gpu_mutex); + for (const auto &g : gpus) + if (g) frames += g->GetFrameCount(); #endif return frames; } const ShadowFinder::Projection &ShadowFinder::Reduced() const { #ifdef JFJOCH_USE_CUDA - if (gpu && gpu->GetFrameCount() > 0) { + if (ShadowAccumulatorGPU *gpu = Gpu(); gpu && gpu->GetFrameCount() > 0) { if (!reduced) reduced = Reduce(); return *reduced; @@ -583,8 +627,10 @@ void ShadowFinder::ReleaseProjection() { std::vector ShadowFinder::GetMeanProjection() const { std::unique_lock ul(m); #ifdef JFJOCH_USE_CUDA - if (Gpu() && gpu->GetFrameCount() > 0 && host.frames == 0) + if (ShadowAccumulatorGPU *gpu = Gpu(); gpu && gpu->GetFrameCount() > 0 && host.frames == 0) { + ScopedGpuPin pin(gpu->GetDevice()); return gpu->MeanProjection(pixel_mask); + } #endif const Projection &p = Reduced(); const auto &sum_value = p.sum_value; @@ -605,7 +651,7 @@ std::vector ShadowFinder::GetMask(size_t nthreads) const { nthreads = AvailableCpus(); #ifdef JFJOCH_USE_CUDA // Where every frame went to the device the mask is made there too, from the projection as it lies. - if (Gpu() && gpu->GetFrameCount() > 0 && host.frames == 0) { + if (ShadowAccumulatorGPU *gpu = Gpu(); gpu && gpu->GetFrameCount() > 0 && host.frames == 0) { const float diag = std::hypot(static_cast(width), static_cast(height)); if (!std::isfinite(beam_x) || !std::isfinite(beam_y) || std::fabs(beam_x - width * 0.5f) > 4.0f * diag || std::fabs(beam_y - height * 0.5f) > 4.0f * diag) @@ -622,6 +668,7 @@ std::vector ShadowFinder::GetMask(size_t nthreads) const { setup.distance_mm = geometry.GetDetectorDistance_mm(); setup.has_polarization = polarization.has_value(); setup.polarization = polarization.value_or(0.0f); + ScopedGpuPin pin(gpu->GetDevice()); return gpu->Mask(setup, pixel_mask); } #endif diff --git a/image_analysis/beam_stop/ShadowFinder.h b/image_analysis/beam_stop/ShadowFinder.h index 866e3df4d..1b1af40ad 100644 --- a/image_analysis/beam_stop/ShadowFinder.h +++ b/image_analysis/beam_stop/ShadowFinder.h @@ -80,19 +80,24 @@ class ShadowFinder { 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 - // 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 + // One per card when GPUs are available. Frames a card can decode are accumulated on the card of the + // thread adding them instead of on the host - only the compressed chunk crosses PCIe - and the + // cards' projections are folded into one, with the host projection, when the mask is read. Frames + // a card cannot take (anything but bitshuffle+LZ4) still go to the host, so a run mixing + // compressions is handled without a second code path. + // Each is 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 + // stall the caller before it has read its first frame. The first AddImage waits for them, by // which time the reads have been running for a while. - mutable std::future> gpu_pending; - mutable std::unique_ptr gpu; - mutable std::mutex gpu_mutex; + mutable std::vector>> gpu_pending; + mutable std::vector> gpus; // by card; null where none + mutable std::mutex gpu_mutex; // guards the two above + std::unique_ptr card_mutex; // a card's frames are added one at a time - // The accumulator once its construction has finished, or null if there is none. + // The accumulator of card `device` once the constructions have finished, or null if it has none. + ShadowAccumulatorGPU *Gpu(int device) const; + // The one accumulator every other card's projection has been folded into, or null if there is + // none. What the projection is read from; frames added after the first call are not seen. [[nodiscard]] ShadowAccumulatorGPU *Gpu() const; #endif @@ -119,6 +124,7 @@ public: static constexpr uint32_t TRANSMITTING = 2; ShadowFinder(const DiffractionExperiment &experiment, const PixelMask &mask); + ~ShadowFinder(); // The centre the rings are drawn about. It starts as the file's, which is the only one there // is when the finder is built; a caller that has measured one replaces it before reading the diff --git a/image_analysis/hot_pixels/HotPixels.cpp b/image_analysis/hot_pixels/HotPixels.cpp index a5c33257e..139e3caca 100644 --- a/image_analysis/hot_pixels/HotPixels.cpp +++ b/image_analysis/hot_pixels/HotPixels.cpp @@ -9,9 +9,12 @@ #include #include +#include "../../common/CUDAWrapper.h" #include "../../common/JFJochMath.h" #include "../../common/ParallelFor.h" +#include + namespace { // The lower median of v[0..n), reordering it. @@ -240,21 +243,57 @@ void HotPixelFinder::AddLevels(const std::vector §or_level, const s } } +HotPixelFinder::~HotPixelFinder() { +#ifdef JFJOCH_USE_CUDA + // Each device half goes with its own card current, where its buffers are. + for (auto &g : gpus) + if (g) { + ScopedGpuPin pin(g->GetDevice()); + g.reset(); + } +#endif +} + #ifdef JFJOCH_USE_CUDA void HotPixelFinder::PrepareDevice() { std::lock_guard lock(m); - if (!gpu) - gpu = std::make_unique(key.get(), width * height, key_begin, nrings, SECTORS, - HotPixelLevelRules{MIN_SECTOR_PIXELS, MIN_RING_PIXELS, - LIT_NSIGMA, LIT_OFFSET}); + if (!gpus.empty()) + return; + std::vector>> built; + for (int32_t c = 0; c < get_gpu_count(); c++) + built.push_back(std::async(std::launch::async, [this, c] { + set_gpu(c); + return std::make_unique(key.get(), width * height, key_begin, nrings, SECTORS, + HotPixelLevelRules{MIN_SECTOR_PIXELS, MIN_RING_PIXELS, + LIT_NSIGMA, LIT_OFFSET}); + })); + for (auto &b : built) + gpus.push_back(b.get()); } void HotPixelFinder::AddDeviceImage(const int32_t *device_image, HotPixelFinderGPU::Frame &frame) { PrepareDevice(); - gpu->Add(device_image, frame); + int device = 0; + cudaGetDevice(&device); + gpus[device]->Add(device_image, frame); std::lock_guard lock(m); frames++; } + +HotPixelFinderGPU *HotPixelFinder::FoldedGpu() { + if (gpus.empty()) + return nullptr; + for (size_t c = 1; c < gpus.size(); c++) { + { + ScopedGpuPin pin(gpus[0]->GetDevice()); + gpus[0]->Fold(*gpus[c]); + } + ScopedGpuPin pin(gpus[c]->GetDevice()); + gpus[c].reset(); + } + gpus.resize(1); + return gpus[0].get(); +} #endif HotPixelFinder::Result HotPixelFinder::GetMask(double oscillation_deg, double spacing_deg, size_t nthreads) { @@ -269,6 +308,10 @@ HotPixelFinder::Result HotPixelFinder::GetMask(double oscillation_deg, double sp // Counted in integers by blocks of rows in parallel, so the totals do not depend on the split. std::vector lit(nrings, 0.0), seen(nrings, 0.0); #ifdef JFJOCH_USE_CUDA + HotPixelFinderGPU *gpu = FoldedGpu(); + std::optional pin; + if (gpu) + pin.emplace(gpu->GetDevice()); if (gpu) { std::vector device_lit, device_seen; gpu->ChanceCounts(device_lit, device_seen); diff --git a/image_analysis/hot_pixels/HotPixels.h b/image_analysis/hot_pixels/HotPixels.h index 0e3b221d4..489d708f9 100644 --- a/image_analysis/hot_pixels/HotPixels.h +++ b/image_analysis/hot_pixels/HotPixels.h @@ -91,6 +91,7 @@ public: }; HotPixelFinder(const DiffractionExperiment &experiment, const PixelMask &mask, size_t nthreads); + ~HotPixelFinder(); // One preprocessed frame (INT32_MIN = masked or error value, INT32_MAX = saturated). Thread // safe; `scratch` is the caller's own, reused across its calls. @@ -98,12 +99,13 @@ public: #ifdef JFJOCH_USE_CUDA // The same for a frame preprocessed on the device, on `frame`'s stream - the caller's own, as - // `scratch` above. The per-pixel sums are then kept on the device and replace the host's when the - // mask is read, so a finder is fed one way or the other, not both. Thread safe. + // `scratch` above, on the calling thread's card. The per-pixel sums are then kept on the cards and + // replace the host's when the mask is read, so a finder is fed one way or the other, not both. + // Thread safe. void AddDeviceImage(const int32_t *device_image, HotPixelFinderGPU::Frame &frame); - // Build the device half now rather than with the first device frame, so the workers do not wait on - // it one behind the other. + // Build the device halves, one per card, now rather than with the first device frame, so the + // workers do not wait on them one behind the other. void PrepareDevice(); #endif @@ -147,7 +149,11 @@ private: std::unique_ptr n_error_ring_ok; std::unique_ptr error_level_sum; #ifdef JFJOCH_USE_CUDA - std::unique_ptr gpu; // built by PrepareDevice or the first device frame + // By card, built by PrepareDevice or the first device frame. Each card's frames add to its own sums, + // and the cards are folded into the first when the mask is read. + std::vector> gpus; + // The first card's device half with every other card's sums added to it, or null without one. + HotPixelFinderGPU *FoldedGpu(); #endif // Each ring-sector's level and lit threshold, from the frame's order statistics, and the frame's diff --git a/image_analysis/hot_pixels/HotPixelsGPU.cu b/image_analysis/hot_pixels/HotPixelsGPU.cu index 3f54f8c89..a5725533b 100644 --- a/image_analysis/hot_pixels/HotPixelsGPU.cu +++ b/image_analysis/hot_pixels/HotPixelsGPU.cu @@ -222,6 +222,12 @@ __global__ void candidate_kernel(size_t npixels, int sectors, uint32_t frames, i if (out) out[slot] = static_cast(i); } +template +__global__ void add_kernel(size_t n, const T *__restrict__ in, T *__restrict__ out) { + const size_t i = blockIdx.x * static_cast(blockDim.x) + threadIdx.x; + if (i < n) out[i] += in[i]; +} + __global__ void iota_kernel(size_t n, uint32_t *__restrict__ out) { const size_t i = blockIdx.x * static_cast(blockDim.x) + threadIdx.x; if (i < n) out[i] = static_cast(i); @@ -254,6 +260,7 @@ HotPixelFinderGPU::HotPixelFinderGPU(const int32_t *host_key, size_t npixels, sum_value(npixels, ALLOC), error_level_sum(npixels, ALLOC), key_frames(std::max(nkeys, 1), ALLOC), key_level_sum(std::max(nkeys, 1), ALLOC), key_local_level_sum(std::max(nkeys, 1), ALLOC) { + cuda_err(cudaGetDevice(&device)); cuda_err(cudaEventCreateWithFlags(&last_accumulate, cudaEventDisableTiming)); CudaStream stream; cuda_err(cudaMemcpyAsync(key, host_key, npixels * sizeof(int32_t), cudaMemcpyHostToDevice, stream)); @@ -325,6 +332,31 @@ void HotPixelFinderGPU::Add(const int32_t *device_image, Frame &frame) { cuda_err(cudaEventRecord(last_accumulate, stream)); } +void HotPixelFinderGPU::Fold(HotPixelFinderGPU &other) { + CudaStream stream; + cuda_err(cudaStreamWaitEvent(stream, last_accumulate, 0)); + cuda_err(cudaStreamWaitEvent(stream, other.last_accumulate, 0)); + // One scratch buffer for every array in turn: the copies and additions are ordered on one stream. + CudaDevicePtr scratch(std::max(npixels, std::max(nkeys, 1)) * sizeof(int64_t), ALLOC); + const auto add = [&](T *mine, const T *theirs, size_t n) { + if (n == 0) return; + cuda_err(cudaMemcpyPeerAsync(scratch, device, theirs, other.device, n * sizeof(T), stream)); + add_kernel<<((n + THREADS - 1) / THREADS), THREADS, 0, stream>>>( + n, reinterpret_cast(scratch.get()), mine); + cuda_err(cudaGetLastError()); + }; + add(n_lit.get(), other.n_lit.get(), npixels); + add(n_error.get(), other.n_error.get(), npixels); + add(n_error_ring_ok.get(), other.n_error_ring_ok.get(), npixels); + add(sum_value.get(), other.sum_value.get(), npixels); + add(error_level_sum.get(), other.error_level_sum.get(), npixels); + add(key_frames.get(), other.key_frames.get(), nkeys); + add(key_level_sum.get(), other.key_level_sum.get(), nkeys); + add(key_local_level_sum.get(), other.key_local_level_sum.get(), nkeys); + cuda_err(cudaEventRecord(last_accumulate, stream)); + cuda_err(cudaStreamSynchronize(stream)); +} + void HotPixelFinderGPU::ChanceCounts(std::vector &lit, std::vector &seen) { CudaStream stream; cuda_err(cudaStreamWaitEvent(stream, last_accumulate, 0)); diff --git a/image_analysis/hot_pixels/HotPixelsGPU.h b/image_analysis/hot_pixels/HotPixelsGPU.h index c8bcd8132..7d9cad18b 100644 --- a/image_analysis/hot_pixels/HotPixelsGPU.h +++ b/image_analysis/hot_pixels/HotPixelsGPU.h @@ -28,6 +28,7 @@ struct HotPixelLevelRules { // only what can decide it comes back: per-ring counts for the chance rate, then the few pixels that can // be persistent or carry the error value. class HotPixelFinderGPU { + int device = 0; // the card it was built on, where everything below lives const size_t npixels; const size_t nkeys; const int nrings; @@ -74,6 +75,13 @@ public: HotPixelFinderGPU(const HotPixelFinderGPU &) = delete; HotPixelFinderGPU &operator=(const HotPixelFinderGPU &) = delete; + [[nodiscard]] int GetDevice() const { return device; } + + // Add the sums another card accumulated to this one's. Every sum is an integer, so the result does + // not depend on how the frames were split between the cards or on the order they are folded in. + // Called with this finder's card current. + void Fold(HotPixelFinderGPU &other); + // Add one frame, preprocessed on frame.stream, as HotPixelFinder::AddImage does. Queued on that // stream; nothing is waited for on the host. void Add(const int32_t *device_image, Frame &frame); diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 0e3af02f4..df0d1fa21 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -45,6 +45,7 @@ class JFJochReader; class JFJochHDF5Reader; +struct JFJochReaderRawImage; struct PipelineLocals; // RugnuxPipeline.h, private to the implementation // Offline reprocessing of a stored Jungfraujoch HDF5 dataset, shared by rugnux (all of its --mode @@ -1119,9 +1120,32 @@ class Rugnux { // first config_.prepass_fraction of the sweep, but at least PREPASS_MIN_ROTATION_DEG (Rugnux.cpp). [[nodiscard]] int PrepassEnd(int start_image, int end_image) const; // The pre-scan's defective-pixel mask (HotPixelFinder), read on `sample` with the rings drawn - // about `ring_centre` where the background measured one. + // about `ring_centre` where the background measured one. frame(i, scratch) is sample frame i, as + // the pre-scan kept it or read into scratch; null where it cannot be read. void MaskDefectivePixels(int start_image, const std::vector &sample, - const std::optional> &ring_centre); + const std::optional> &ring_centre, + const std::function &frame); + // The beam-stop detection of the pre-scan, kept under everything it reads (the key in PreScan): a + // detection that finds no shadow leaves the mask as it was, and every copy of the run that + // pre-scans after it would detect again on the same frames, geometry and mask. Shared with the + // copies the run makes of itself. + struct BeamStopMemo { + // A handful covers the geometries a run's probes pre-scan at; each holds the mean projection, + // 4 bytes a pixel, which the beam-centre capture and the diagnostic picture read. + static constexpr size_t MAX_ENTRIES = 4; + struct Entry { + std::vector key; + std::vector> shadow; // the mask's non-zero pixels + std::shared_ptr> projection; + std::optional> ring_centre; + std::string centre_line; // what the detection said about the centre + size_t frames_read = 0; + size_t shadow_frames = 0; + }; + std::mutex m; + std::vector entries; + }; + std::shared_ptr beam_stop_memo_ = std::make_shared(); // Powder calibration, second pass: re-integrate the same images into a profile binned with the // geometry the first fit produced, and fit that. Nothing but the binning changes, which is the // point - which pixel went into which (q, phi) bin was decided by the geometry the run started diff --git a/rugnux/RugnuxPreScan.cpp b/rugnux/RugnuxPreScan.cpp index e965f3ca7..7447ba72c 100644 --- a/rugnux/RugnuxPreScan.cpp +++ b/rugnux/RugnuxPreScan.cpp @@ -149,6 +149,31 @@ namespace { // Measured over the rotation battery: the radius chosen is the whole-sample answer on 38 of 38 // crystals, from an average of 32 frames of the 60 instead of all of them. constexpr std::array WIDTH_TIER_STRIDE = {8, 4, 2, 1}; + + // The sample is read once and kept for everything that reads it after the projection - the spot + // measurement and the defective pixels used to read it again each, and a CBF, SMV or marCCD frame + // is decoded on the host every time it is read. Frames are kept as the reader gives them, + // compressed where the file is: a few MB a frame on a 16 Mpx EIGER, 4 bytes a pixel where the + // reader decodes. A frame that would take the kept frames past this is let go and read again + // where it is wanted, so a sample of large decoded frames costs at most this much memory. + constexpr size_t PRESCAN_KEPT_BYTES = size_t{2} << 30; + + // fn(t) for t < n on threads of their own, thread t on card t % gpus as an image-loop worker is: an + // engine built on a worker's card must be used, and freed, with that card current. Which card a + // frame goes to changes nothing it gives: the spot finders are exact, and the projection and the + // defective-pixel sums are integers, kept per card and added together when they are read. + template + void OnWorkers(size_t n, Fn &&fn) { + std::vector> futures; + futures.reserve(n); + for (size_t t = 0; t < n; t++) + futures.emplace_back(std::async(std::launch::async, [&fn, t] { + if (const int32_t gpus = get_gpu_count(); gpus > 0) + pin_gpu(static_cast(t) % gpus); + fn(t); + })); + for (auto &f : futures) f.get(); + } } void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, RugnuxObserver *observer) { @@ -193,7 +218,37 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru observer->OnPhase(want_shadow ? "Beam stop" : "Spot width"); Step("Pre-scan"); - ShadowFinder finder(experiment_, pixel_mask_); + // What the beam-stop detection reads: the frames, the geometry the rings are drawn in and the mask. + // A detection that finds no shadow leaves the mask as it was, so every copy of the run that + // pre-scans after it - a probe, a re-check - would detect again on the same inputs; it takes the + // answer kept under them instead (BeamStopMemo). + std::vector shadow_key; + std::optional shadow_memo; + if (want_shadow) { + const auto &g = experiment_.GetDiffractionGeometry(); + const uint64_t mask_checksum = pixel_mask_.GetBinaryMaskChecksum(); + shadow_key = {static_cast(start_image), static_cast(config_.stride), + static_cast(experiment_.GetXPixelsNumConv()), + static_cast(experiment_.GetYPixelsNumConv()), + static_cast(experiment_.IsGeometryTransformed()), + static_cast(mask_checksum >> 32), static_cast(mask_checksum & 0xffffffffu), + g.GetBeamX_pxl(), g.GetBeamY_pxl(), g.GetDetectorDistance_mm(), g.GetPixelSize_mm(), + g.GetWavelength_A(), g.GetPoniRot1_rad(), g.GetPoniRot2_rad(), g.GetPoniRot3_rad(), + static_cast(experiment_.GetPolarizationFactor().has_value()), + experiment_.GetPolarizationFactor().value_or(0.0f)}; + for (const float v : g.GetDetectorMatrix().arr()) + shadow_key.push_back(v); + shadow_key.insert(shadow_key.end(), shadow_sample.begin(), shadow_sample.end()); + std::lock_guard lock(beam_stop_memo_->m); + for (const auto &e : beam_stop_memo_->entries) + if (e.key == shadow_key) + shadow_memo = e; + } + const bool find_shadow = want_shadow && !shadow_memo; + + std::optional finder; + if (find_shadow) + finder.emplace(experiment_, pixel_mask_); // The spot width is read in the same pass as the projection rather than in a pre-pass of its own. // The image is preprocessed a second time for it: the shadow is accumulated from raw counts and @@ -253,7 +308,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru }; // The spots of one image, and what is taken off them. - const auto find_spots = [&](PreScanWorker &w, CompressedImage &image, int image_idx, + const auto find_spots = [&](PreScanWorker &w, const CompressedImage &image, int image_idx, bool for_width, std::vector &curves, std::vector &spot_q, std::vector &footprint) { try { @@ -341,65 +396,101 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // Per frame, in sample order, so the gain does not depend on how the workers interleaved. std::vector gain_levels(want_gain ? ordinals.size() : 0); const size_t gain_frame_stride = (ordinals.size() + background_gain::GAIN_MAX_FRAMES - 1) / background_gain::GAIN_MAX_FRAMES; - { + + // The sample as read here, by position, for the spot measurement and the defective pixels below + // (PRESCAN_KEPT_BYTES). + const bool keep_frames = want_spots || want_hot; + std::vector kept(keep_frames ? ordinals.size() : 0); + std::vector is_kept(ordinals.size(), 0); + std::atomic kept_bytes{0}; + // Frame i of the sample: the one kept, or read into `scratch`. Null where it cannot be read. + const auto frame = [&](size_t i, JFJochReaderRawImage &scratch) -> const CompressedImage * { + if (is_kept[i]) + return &kept[i].image; + const int image_idx = start_image + ordinals[i] * config_.stride; + try { + if (reader_.ReadRawImage(image_idx, scratch)) + return &scratch.image; + } catch (const std::exception &e) { + if (IsFatalResourceError(e)) throw; + logger.Warning("Pre-scan: failed to load image {}: {}", image_idx, e.what()); + } + return nullptr; + }; + + if (find_shadow || want_gain || keep_frames) { + TimingMark mark(logger, "pre-scan read and projection"); 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, [&] { - std::vector shadow_buffer; - std::vector gain_pixels; - JFJochReaderRawImage raw_image; - for (size_t i = next.fetch_add(1); i < ordinals.size(); i = next.fetch_add(1)) { - const int ordinal = ordinals[i]; - const int image_idx = start_image + ordinal * config_.stride; - bool read = false; - try { - read = reader_.ReadRawImage(image_idx, raw_image); - } catch (const std::exception &e) { - if (IsFatalResourceError(e)) throw; - logger.Warning("Pre-scan: failed to load image {}: {}", image_idx, e.what()); - continue; + OnWorkers(nworkers, [&](size_t) { + std::vector shadow_buffer; + std::vector gain_pixels; + JFJochReaderRawImage scratch; + for (size_t i = next.fetch_add(1); i < ordinals.size(); i = next.fetch_add(1)) { + const int ordinal = ordinals[i]; + const int image_idx = start_image + ordinal * config_.stride; + JFJochReaderRawImage &raw_image = keep_frames ? kept[i] : scratch; + bool read = false; + try { + read = reader_.ReadRawImage(image_idx, raw_image); + } catch (const std::exception &e) { + if (IsFatalResourceError(e)) throw; + logger.Warning("Pre-scan: failed to load image {}: {}", image_idx, e.what()); + continue; + } + if (!read) continue; + frames_read.fetch_add(1, std::memory_order_relaxed); + // Gated on find_shadow, not on the frame set alone: AddImage projects the whole + // image, and with --detect-beam-stop=off (or a shadow already in the mask) that is + // the pass's largest per-frame cost paid for a result nobody reads. The set still + // decides WHICH frames feed it, so the mask is built from the frames it always was. + if (find_shadow && shadow_set.contains(ordinal)) { + DataMessage msg{}; + msg.image = raw_image.image; + msg.number = ordinal; + msg.original_number = image_idx; + finder->AddImage(msg, shadow_buffer); + } + if (want_gain && i % gain_frame_stride == 0) { + const auto &img = raw_image.image; + const size_t n = img.GetWidth() * img.GetHeight(); + const uint8_t *ptr = img.GetUncompressedPtr(shadow_buffer); + gain_pixels.resize(n); + bool measurable = true; + switch (img.GetMode()) { + case CompressedImageMode::Int32: + std::memcpy(gain_pixels.data(), ptr, n * sizeof(int32_t)); + break; + case CompressedImageMode::Uint16: + std::copy_n(reinterpret_cast(ptr), n, gain_pixels.begin()); + break; + case CompressedImageMode::Int16: + std::copy_n(reinterpret_cast(ptr), n, gain_pixels.begin()); + break; + default: + measurable = false; // the CCD readers give 16- or 32-bit integers } - if (!read) continue; - frames_read.fetch_add(1, std::memory_order_relaxed); - // Gated on want_shadow, not on the frame set alone: AddImage projects the whole - // image, and with --detect-beam-stop=off (or a shadow already in the mask) that is - // the pass's largest per-frame cost paid for a result nobody reads. The set still - // decides WHICH frames feed it, so the mask is built from the frames it always was. - if (want_shadow && shadow_set.contains(ordinal)) { - DataMessage msg{}; - msg.image = raw_image.image; - msg.number = ordinal; - msg.original_number = image_idx; - finder.AddImage(msg, shadow_buffer); - } - if (want_gain && i % gain_frame_stride == 0) { - const auto &img = raw_image.image; - const size_t n = img.GetWidth() * img.GetHeight(); - const uint8_t *ptr = img.GetUncompressedPtr(shadow_buffer); - gain_pixels.resize(n); - switch (img.GetMode()) { - case CompressedImageMode::Int32: - std::memcpy(gain_pixels.data(), ptr, n * sizeof(int32_t)); - break; - case CompressedImageMode::Uint16: - std::copy_n(reinterpret_cast(ptr), n, gain_pixels.begin()); - break; - case CompressedImageMode::Int16: - std::copy_n(reinterpret_cast(ptr), n, gain_pixels.begin()); - break; - default: - continue; // the CCD readers give 16- or 32-bit integers - } + if (measurable) gain_levels[i] = background_gain::MeasureFrame(gain_pixels, img.GetWidth(), img.GetHeight(), pixel_mask_.GetMask(), experiment_.GetSaturationLimit()); + } + if (keep_frames) { + // The reader's own scratch is not part of the image (JFJochReaderRawImage). + raw_image.read_buffer = {}; + const size_t bytes = raw_image.image_buffer.capacity(); + if (kept_bytes.fetch_add(bytes) + bytes <= PRESCAN_KEPT_BYTES) + is_kept[i] = 1; + else { + kept_bytes.fetch_sub(bytes); + raw_image = {}; } } - })); - for (auto &f : futures) f.get(); - } + } + }); + logger.Debug("Pre-scan: {} of {} frames kept for the steps after the projection, {:.0f} MB", + std::count(is_kept.begin(), is_kept.end(), 1), ordinals.size(), kept_bytes.load() / 1e6); + } else + frames_read = shadow_memo->frames_read; // Applied before anything below copies the integration settings (bragg_before_adaptive_), so every // pass and every fallback integrates with it. At or below 1 - an ADSC or a marCCD, whose pixel @@ -425,6 +516,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru } const auto spot_measurement = [&] { + TimingMark mark(logger, "pre-scan spot measurement"); prescan_mapping = std::make_unique(prescan_x, prescan_mask, *azint_geometry_); std::vector> curves_of(ordinals.size()); @@ -459,20 +551,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru std::vector workers(nworkers); std::vector worker_built(nworkers, 0); - // fn(t) on threads of their own, thread t on card t % gpus as an image-loop worker is: worker - // t's engine is built on that card and must be used, and freed, with it current - and each - // tier runs on fresh threads. Which card finds a frame's spots does not change them. - const auto on_workers = [&](auto &&fn) { - std::vector> futures; - futures.reserve(nworkers); - for (size_t t = 0; t < nworkers; t++) - futures.emplace_back(std::async(std::launch::async, [&, t] { - if (const int32_t gpus = get_gpu_count(); gpus > 0) - pin_gpu(static_cast(t) % gpus); - fn(t); - })); - for (auto &f : futures) f.get(); - }; + // Each tier runs on fresh threads (OnWorkers), worker t on the card its engine was built on. bool width_settled = !want_width; // The bandwidth reads the same spots, but its slope needs more of them than the width does: @@ -482,7 +561,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru size_t phase_begin = 0; for (const size_t phase_end : tier_end) { std::atomic next{phase_begin}; - on_workers([&](size_t t) { + OnWorkers(nworkers, [&](size_t t) { if (!worker_built[t]) { workers[t] = make_worker(); worker_built[t] = 1; @@ -496,19 +575,11 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru const bool for_shape = width_settled && !shapes_settled; if (!(for_width || for_shape)) continue; // nothing to measure on this frame - bool read = false; - try { - read = reader_.ReadRawImage(image_idx, w.raw_image); - } catch (const std::exception &e) { - if (IsFatalResourceError(e)) throw; - logger.Warning("Pre-scan: failed to load image {}: {}", image_idx, - e.what()); - continue; - } - if (!read) continue; + const CompressedImage *image = frame(i, w.raw_image); + if (!image) continue; if (for_shape) shape_only[i] = 1; - find_spots(w, w.raw_image.image, image_idx, for_width || for_shape, curves_of[i], + find_spots(w, *image, image_idx, for_width || for_shape, curves_of[i], spot_q_of[i], footprint_of[i]); } }); @@ -533,7 +604,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru shape_curves.insert(shape_curves.end(), c.begin(), c.end()); shapes_settled = width_settled && shape_curves.size() >= spot_width::BANDWIDTH_POOL_SPOTS; } - on_workers([&](size_t t) { PreScanWorker released = std::move(workers[t]); }); + OnWorkers(nworkers, [&](size_t t) { PreScanWorker released = std::move(workers[t]); }); for (size_t i = 0; i < ordinals.size(); i++) if (!shape_only[i]) @@ -717,23 +788,52 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // frame of its own, and an unmasked shadow biases it by a few pixels at most, which is // three orders below what this comparison would notice. The centre the run reports and // consumes is the one fitted below instead, with the shadow out of the way. - projection = std::make_shared>(finder.GetMeanProjection()); - if (const auto ring_center = FindBeamCenterFromBackground(experiment_, pixel_mask_, - *projection)) { - finder.BeamCenter(ring_center->beam_x_pxl, ring_center->beam_y_pxl); - ring_centre = std::make_pair(ring_center->beam_x_pxl, ring_center->beam_y_pxl); - logger.Info("Beam stop: comparing each pixel against its ring about the measured " - "centre ({:.2f},{:.2f}) +- {:.2f} px, {:.2f} px from the file's", - ring_center->beam_x_pxl, ring_center->beam_y_pxl, ring_center->sigma_pxl, - std::hypot(ring_center->beam_x_pxl - experiment_.GetBeamX_pxl(), - ring_center->beam_y_pxl - experiment_.GetBeamY_pxl())); + std::vector shadow; + size_t shadow_frames = 0; + if (shadow_memo) { + projection = shadow_memo->projection; + ring_centre = shadow_memo->ring_centre; + shadow.assign(experiment_.GetXPixelsNumConv() * experiment_.GetYPixelsNumConv(), 0); + for (const auto &[i, v] : shadow_memo->shadow) + shadow[i] = v; + shadow_frames = shadow_memo->shadow_frames; + logger.Info("{} - kept from an earlier pre-scan of this run on the same frames, geometry and mask", + shadow_memo->centre_line); } else { - logger.Info("Beam stop: the background does not place a centre on this run, so the " - "rings are drawn about the file's ({:.2f},{:.2f})", - experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl()); + TimingMark mark(logger, "pre-scan beam stop"); + std::string centre_line; + projection = std::make_shared>(finder->GetMeanProjection()); + if (const auto ring_center = FindBeamCenterFromBackground(experiment_, pixel_mask_, + *projection)) { + finder->BeamCenter(ring_center->beam_x_pxl, ring_center->beam_y_pxl); + ring_centre = std::make_pair(ring_center->beam_x_pxl, ring_center->beam_y_pxl); + centre_line = fmt::format("Beam stop: comparing each pixel against its ring about the measured " + "centre ({:.2f},{:.2f}) +- {:.2f} px, {:.2f} px from the file's", + ring_center->beam_x_pxl, ring_center->beam_y_pxl, ring_center->sigma_pxl, + std::hypot(ring_center->beam_x_pxl - experiment_.GetBeamX_pxl(), + ring_center->beam_y_pxl - experiment_.GetBeamY_pxl())); + } else { + centre_line = fmt::format("Beam stop: the background does not place a centre on this run, so the " + "rings are drawn about the file's ({:.2f},{:.2f})", + experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl()); + } + logger.Info("{}", centre_line); + shadow = finder->GetMask(config_.nthreads); + shadow_frames = finder->GetFrameCount(); + finder.reset(); + + BeamStopMemo::Entry entry{shadow_key, {}, projection, ring_centre, centre_line, frames_read.load(), + shadow_frames}; + for (size_t i = 0; i < shadow.size(); i++) + if (shadow[i] != 0) + entry.shadow.emplace_back(static_cast(i), shadow[i]); + std::lock_guard lock(beam_stop_memo_->m); + auto &entries = beam_stop_memo_->entries; + // A handful covers the geometries a run's probes pre-scan at; the oldest go first. + if (entries.size() == BeamStopMemo::MAX_ENTRIES) + entries.erase(entries.begin()); + entries.push_back(std::move(entry)); } - const auto shadow = finder.GetMask(config_.nthreads); - finder.ReleaseProjection(); const auto shadowed = std::count_if(shadow.begin(), shadow.end(), [](uint32_t v) { return v != 0; }); pixel_mask_.LoadBeamStopMask(experiment_, shadow); opaque_shadow = shadow; @@ -747,11 +847,11 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru logger.Warning("Beam stop: no shadow found in {} images. The comparison is made about " "the centre named above, so this may point at that centre rather than " "at the experiment", - finder.GetFrameCount()); + shadow_frames); else { logger.Info("Beam stop shadow: {} pixels ({:.2f}% of the detector) found in {} images", shadowed, 100.0 * static_cast(shadowed) / static_cast(shadow.size()), - finder.GetFrameCount()); + shadow_frames); Note(fmt::format("beam stop masked: {:.2f}% of the detector", 100.0 * static_cast(shadowed) / static_cast(shadow.size()))); } @@ -785,7 +885,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // second pass of a two-pass run the shadow is already in the mask, so this does not run again - // by then the centre has been post-refined and the file's value is no longer what is in it. if (want_hot) - MaskDefectivePixels(start_image, shadow_sample, ring_centre); + MaskDefectivePixels(start_image, shadow_sample, ring_centre, frame); // The centre is fitted with the shadow masked but NOT the transmitting pieces the shadow finder // adds after it. Their dimming is multiplicative in azimuth, which the fit's per-sector amplitude @@ -888,15 +988,17 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru } void Rugnux::MaskDefectivePixels(int start_image, const std::vector &sample, - const std::optional> &ring_centre) { + const std::optional> &ring_centre, + const std::function &frame) { Logger logger("Rugnux"); + TimingMark mark(logger, "pre-scan defective pixels"); hot_pixels_measured_ = true; - // The frames are read a second time, because the rings each pixel is compared against have to - // be drawn about the beam, and the centre the file states can be far enough out that a ring - // crosses the background's own radial fall-off - the pixels on its bright side then read high on - // every frame. The centre the beam-stop projection has just measured does not exist until every - // frame of the first read is in. + // The rings each pixel is compared against have to be drawn about the beam, and the centre the + // file states can be far enough out that a ring crosses the background's own radial fall-off - the + // pixels on its bright side then read high on every frame. The centre the beam-stop projection + // measures does not exist until every frame of the first read is in, so the frames are gone + // through a second time - from memory, as the pre-scan kept them, where it could. DiffractionExperiment about = experiment_; if (ring_centre) about.BeamX_pxl(ring_centre->first).BeamY_pxl(ring_centre->second); @@ -909,65 +1011,63 @@ void Rugnux::MaskDefectivePixels(int start_image, const std::vector &sample finder.PrepareDevice(); #endif std::atomic next{0}; - std::vector> futures; const size_t nworkers = std::min(nthreads, std::min(PRESCAN_MAX_WORKERS, sample.size())); - for (size_t t = 0; t < nworkers; t++) - futures.emplace_back(std::async(std::launch::async, [&] { - // With a GPU the frame is decoded, preprocessed and analysed on the device - only the - // compressed chunk crosses PCIe - and the result is the same to the bit (HotPixelsGPU.h). - std::unique_ptr preprocessor; - std::unique_ptr image; + OnWorkers(nworkers, [&](size_t) { + // With a GPU the frame is decoded, preprocessed and analysed on the device - only the + // compressed chunk crosses PCIe - and the result is the same to the bit (HotPixelsGPU.h). + std::unique_ptr preprocessor; + std::unique_ptr image; #ifdef JFJOCH_USE_CUDA - std::optional frame; - if (get_gpu_count() > 0) { - auto stream = std::make_shared(); - preprocessor = std::make_unique(experiment_, pixel_mask_, stream, - /*copy_image_to_host=*/false); - image = std::make_unique(experiment_.GetPixelsNum(), - /*host_mirror=*/false); - frame.emplace(stream); - } + std::optional device_frame; + if (get_gpu_count() > 0) { + auto stream = std::make_shared(); + preprocessor = std::make_unique(experiment_, pixel_mask_, stream, + /*copy_image_to_host=*/false); + image = std::make_unique(experiment_.GetPixelsNum(), + /*host_mirror=*/false); + device_frame.emplace(stream); + } #endif - if (!preprocessor) { - preprocessor = std::make_unique(experiment_, pixel_mask_); - image = std::make_unique(experiment_.GetPixelsNum()); - } - JFJochReaderRawImage raw; - std::vector buffer; - std::vector scratch; - for (size_t i = next.fetch_add(1); i < sample.size(); i = next.fetch_add(1)) { - const int image_idx = start_image + sample[i] * config_.stride; + if (!preprocessor) { + preprocessor = std::make_unique(experiment_, pixel_mask_); + image = std::make_unique(experiment_.GetPixelsNum()); + } + JFJochReaderRawImage raw; + std::vector buffer; + std::vector scratch; + for (size_t i = next.fetch_add(1); i < sample.size(); i = next.fetch_add(1)) { + const int image_idx = start_image + sample[i] * config_.stride; + const CompressedImage *compressed = frame(i, raw); + if (!compressed) continue; + try { + // As in the image loops (MXAnalysisWithoutFPGA::Analyze): a frame the device cannot + // decode goes to the host decoder, so a GPU-side problem costs speed, not the frame. + ImageStatistics stats; + bool decoded_on_device = false; try { - if (!reader_.ReadRawImage(image_idx, raw)) continue; - // As in the image loops (MXAnalysisWithoutFPGA::Analyze): a frame the device cannot - // decode goes to the host decoder, so a GPU-side problem costs speed, not the frame. - ImageStatistics stats; - bool decoded_on_device = false; - try { - decoded_on_device = preprocessor->AnalyzeCompressed(*image, raw.image, stats); - } catch (const JFJochException &e) { - cuda_throw_if_context_lost(); - logger.Warning("Defective pixels: device decoding of image {} failed ({}), " - "decompressing it on the host", image_idx, e.what()); - cuda_clear_error(); - } - if (!decoded_on_device) - preprocessor->Analyze(*image, raw.image.GetUncompressedPtr(buffer), raw.image.GetMode()); - } catch (const std::exception &e) { - if (IsFatalResourceError(e)) throw; - logger.Warning("Defective pixels: failed to read image {}: {}", image_idx, e.what()); - continue; + decoded_on_device = preprocessor->AnalyzeCompressed(*image, *compressed, stats); + } catch (const JFJochException &e) { + cuda_throw_if_context_lost(); + logger.Warning("Defective pixels: device decoding of image {} failed ({}), " + "decompressing it on the host", image_idx, e.what()); + cuda_clear_error(); } -#ifdef JFJOCH_USE_CUDA - if (frame) { - finder.AddDeviceImage(image->getGPUBuffer(), *frame); - continue; - } -#endif - finder.AddImage(image->data(), scratch); + if (!decoded_on_device) + preprocessor->Analyze(*image, compressed->GetUncompressedPtr(buffer), compressed->GetMode()); + } catch (const std::exception &e) { + if (IsFatalResourceError(e)) throw; + logger.Warning("Defective pixels: failed to read image {}: {}", image_idx, e.what()); + continue; } - })); - for (auto &f : futures) f.get(); +#ifdef JFJOCH_USE_CUDA + if (device_frame) { + finder.AddDeviceImage(image->getGPUBuffer(), *device_frame); + continue; + } +#endif + finder.AddImage(image->data(), scratch); + } + }); const auto goniometer = experiment_.GetGoniometer(); const double oscillation = std::abs(goniometer->GetIncrement_deg()); diff --git a/tests/HotPixelFinderTest.cpp b/tests/HotPixelFinderTest.cpp index dc897ffa5..a4015e8b8 100644 --- a/tests/HotPixelFinderTest.cpp +++ b/tests/HotPixelFinderTest.cpp @@ -4,7 +4,9 @@ #include #include +#include #include +#include #include #include "../common/CUDAWrapper.h" @@ -161,4 +163,75 @@ TEST_CASE("HotPixelFinder_DeviceMatchesHost", "[HotPixelFinder]") { CHECK(result.error == expected.error); CHECK(result.mask == expected.mask); } + +// The pre-scan keeps one device half per card and folds the cards together when the mask is read. Two +// halves, each given part of the frames, folded into one must give exactly the chance counts and the +// candidates of one half given every frame. +TEST_CASE("HotPixelFinderGPU_FoldedCardsMatchOneCard", "[HotPixelFinder]") { + if (get_gpu_count() == 0) { + WARN("No CUDA GPU present. Skipping HotPixelFinderGPU_FoldedCardsMatchOneCard"); + return; + } + // Four rings of two sectors: a key by band of rows and half of the detector, one column masked. + constexpr int NRINGS = 4, SECTORS = 2; + const size_t npixels = static_cast(W) * H; + std::vector key(npixels); + std::vector key_begin(NRINGS * SECTORS + 1, 0); + for (int y = 0; y < H; y++) + for (int x_ = 0; x_ < W; x_++) { + const int k = x_ == 7 ? -1 : (y * NRINGS / H) * SECTORS + x_ * SECTORS / W; + key[I(x_, y)] = k; + if (k >= 0) key_begin[k + 1]++; + } + for (int k = 0; k < NRINGS * SECTORS; k++) + key_begin[k + 1] += key_begin[k]; + const HotPixelLevelRules rules{HotPixelFinder::MIN_SECTOR_PIXELS, HotPixelFinder::MIN_RING_PIXELS, + HotPixelFinder::LIT_NSIGMA, HotPixelFinder::LIT_OFFSET}; + HotPixelFinderGPU one(key.data(), npixels, key_begin, NRINGS, SECTORS, rules); + HotPixelFinderGPU first(key.data(), npixels, key_begin, NRINGS, SECTORS, rules); + HotPixelFinderGPU second(key.data(), npixels, key_begin, NRINGS, SECTORS, rules); + HotPixelFinderGPU::Frame frame(std::make_shared()); + CudaDevicePtr device_image(npixels); + + std::mt19937 rng(11); + std::uniform_real_distribution u(0.0, 1.0); + std::vector image(npixels); + for (int f = 0; f < NFRAMES; f++) { + for (size_t i = 0; i < npixels; i++) { + const double d = u(rng); + image[i] = d < 0.01 ? INT32_MIN : d < 0.015 ? INT32_MAX : std::poisson_distribution(8.0)(rng); + } + for (int p = 0; p < 40; p++) + image[I(20 + p, 30 + p)] += f % 3 == 0 ? 0 : 40; // lit on most frames + REQUIRE(cudaMemcpy(device_image, image.data(), npixels * sizeof(int32_t), cudaMemcpyHostToDevice) + == cudaSuccess); + one.Add(device_image, frame); + (f % 3 == 0 ? second : first).Add(device_image, frame); + REQUIRE(cudaStreamSynchronize(*frame.stream) == cudaSuccess); + } + first.Fold(second); + + std::vector lit_one, seen_one, lit_folded, seen_folded; + one.ChanceCounts(lit_one, seen_one); + first.ChanceCounts(lit_folded, seen_folded); + CHECK(lit_folded == lit_one); + CHECK(seen_folded == seen_one); + const std::vector k_chance(NRINGS, 5); + const auto expected = one.GetCandidates(NFRAMES, 10, true, k_chance); + const auto result = first.GetCandidates(NFRAMES, 10, true, k_chance); + CHECK(expected.index.size() >= 40); + // The candidates are written by atomic slot, so their order is the device's; compare them by pixel. + const auto by_pixel = [](const HotPixelFinderGPU::Candidates &c) { + std::vector> v; + for (size_t j = 0; j < c.index.size(); j++) + v.emplace_back(c.index[j], c.n_lit[j], c.n_error[j], c.n_error_ring_ok[j], c.sum_value[j], + c.error_level_sum[j]); + std::sort(v.begin(), v.end()); + return v; + }; + CHECK(by_pixel(result) == by_pixel(expected)); + CHECK(result.key_frames == expected.key_frames); + CHECK(result.key_level_sum == expected.key_level_sum); + CHECK(result.key_local_level_sum == expected.key_local_level_sum); +} #endif diff --git a/tests/ShadowFinderTest.cpp b/tests/ShadowFinderTest.cpp index f01e2c0b4..2fdbc7691 100644 --- a/tests/ShadowFinderTest.cpp +++ b/tests/ShadowFinderTest.cpp @@ -558,4 +558,50 @@ TEST_CASE("ShadowFinder_DeviceMaskMatchesHost", "[ShadowFinder]") { CHECK(r.device_mask == r.host_mask); } } + +// The pre-scan keeps one accumulator per card and folds them together when the mask is read. Two +// accumulators, each given half of the frames, folded into one must hold exactly what one given every +// frame holds - including a pixel whose only valid values are negative and fell all to the second half, +// which the maximum must take although the first half never counted the pixel. +TEST_CASE("ShadowAccumulatorGPU_FoldedCardsMatchOneCard", "[ShadowFinder]") { + if (get_gpu_count() == 0) + SKIP("no GPU"); + constexpr size_t N = static_cast(W) * H; + std::vector> frames; + uint32_t state = 12345; + const auto next = [&] { state = state * 1664525u + 1013904223u; return state >> 8; }; + for (int f = 0; f < 2 * NFRAMES; f++) { + frames.emplace_back(N); + for (auto &v : frames.back()) + v = next() % 50 == 0 ? INT32_MIN : static_cast(next() % 1100) - 100; + frames.back()[I(3, 3)] = f % 2 == 0 ? INT32_MIN : -7 - f; + frames.back()[I(5, 5)] = INT32_MIN; + } + JFJochBitShuffleCompressor compressor(CompressionAlgorithm::BSHUF_LZ4); + std::vector> compressed; + for (const auto &frame : frames) + compressed.push_back(compressor.Compress(frame)); + const auto image = [&](size_t f) { + return CompressedImage(compressed[f].data(), compressed[f].size(), W, H, CompressedImageMode::Int32, + CompressionAlgorithm::BSHUF_LZ4); + }; + + ShadowAccumulatorGPU one(N), first(N), second(N); + for (size_t f = 0; f < frames.size(); f++) { + one.Add(image(f)); + (f % 2 == 0 ? first : second).Add(image(f)); + } + first.Fold(second); + + std::vector max_one, sum_one, max_folded, sum_folded; + std::vector count_one, count_folded; + one.Download(max_one, sum_one, count_one); + first.Download(max_folded, sum_folded, count_folded); + CHECK(first.GetFrameCount() == one.GetFrameCount()); + CHECK(max_folded[I(3, 3)] == -8); + CHECK(count_folded[I(5, 5)] == 0); + CHECK(max_folded == max_one); + CHECK(sum_folded == sum_one); + CHECK(count_folded == count_one); +} #endif