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