diff --git a/common/AzimuthalIntegrationMapping.cpp b/common/AzimuthalIntegrationMapping.cpp index d8b4901a..0af849c2 100644 --- a/common/AzimuthalIntegrationMapping.cpp +++ b/common/AzimuthalIntegrationMapping.cpp @@ -10,6 +10,7 @@ #include "JFJochException.h" #include "DiffractionGeometry.h" #include "RawToConvertedGeometry.h" +#include "TableChecksum.h" AzimuthalIntegrationMapping::AzimuthalIntegrationMapping(const DiffractionExperiment &experiment, const PixelMask& mask, @@ -49,6 +50,9 @@ AzimuthalIntegrationMapping::AzimuthalIntegrationMapping(const DiffractionExperi SetupConvGeom(experiment.GetDiffractionGeometry(),mask.GetMask()); UpdateMaxBinNumber(); + + pixel_to_bin_checksum = TableChecksum(pixel_to_bin.data(), pixel_to_bin.size() * sizeof(uint16_t)); + corrections_checksum = TableChecksum(corrections.data(), corrections.size() * sizeof(float)); } void AzimuthalIntegrationMapping::SetupConvGeomRows(const DiffractionGeometry &geom, const std::vector &mask, @@ -207,6 +211,39 @@ const std::vector &AzimuthalIntegrationMapping::Resolution() const { return pixel_resolution; } +uint64_t AzimuthalIntegrationMapping::GetPixelToBinChecksum() const { + return pixel_to_bin_checksum; +} + +uint64_t AzimuthalIntegrationMapping::GetCorrectionsChecksum() const { + return corrections_checksum; +} + +std::shared_ptr> +AzimuthalIntegrationMapping::ResolutionMaskBits(std::optional high_res, + std::optional low_res) const { + const std::lock_guard lock(res_mask_mutex); + if (res_mask_bits && res_mask_high == high_res && res_mask_low == low_res) + return res_mask_bits; + + // An unset limit masks nothing at that end. At the high-resolution end 0 does that on its own - no + // pixel has d < 0, and the detector's own edge is where the pixels stop anyway; at the + // low-resolution end every pixel lies above any finite stand-in, so it takes an infinite one. + const float high = high_res.value_or(0.0f); + const float low = low_res.value_or(INFINITY); + + const size_t npixel = pixel_resolution.size(); + auto bits = std::make_shared>(npixel / 32 + (npixel % 32 != 0 ? 1 : 0), 0); + for (size_t i = 0; i < npixel; i++) + if (pixel_resolution[i] > low || pixel_resolution[i] < high) + (*bits)[i / 32] |= 1u << (i % 32); + + res_mask_high = high_res; + res_mask_low = low_res; + res_mask_bits = bits; + return bits; +} + const AzimuthalIntegrationSettings &AzimuthalIntegrationMapping::Settings() const { return settings; } diff --git a/common/AzimuthalIntegrationMapping.h b/common/AzimuthalIntegrationMapping.h index 14337547..716c7c50 100644 --- a/common/AzimuthalIntegrationMapping.h +++ b/common/AzimuthalIntegrationMapping.h @@ -3,6 +3,8 @@ #pragma once +#include +#include #include #include "DiffractionExperiment.h" #include "PixelMask.h" @@ -24,6 +26,23 @@ protected: std::optional polarization_factor; + // Checksums of the two tables the GPU engines upload, taken once here. They are part of the key + // the shared device-table cache looks them up by (CudaSharedTables.h), so an engine that hands + // its own checksum in does not have to hash tens of megabytes on its way to a cache hit - and + // there is one engine per worker per pass. Both vectors are written in the constructor and never + // touched again, so a checksum taken there stays true for the mapping's whole life. + uint64_t pixel_to_bin_checksum = 0; + uint64_t corrections_checksum = 0; + + // Bit-packed "this pixel is outside the resolution limits" mask, memoised for the limits it was + // last built for. Every worker's spot finder wants the identical mask and building it walks every + // pixel, so it is built once and shared. The mapping is const and read from all the worker threads + // at once, hence the mutex; and the mask is handed out as a shared_ptr, so a caller keeps the one + // it was given even if another thread later replaces the cached one. + mutable std::mutex res_mask_mutex; + mutable std::optional res_mask_high, res_mask_low; + mutable std::shared_ptr> res_mask_bits; + size_t nthreads; void UpdateMaxBinNumber(); @@ -49,6 +68,12 @@ public: [[nodiscard]] uint16_t QToBin(float q) const; [[nodiscard]] const std::vector &Corrections() const; [[nodiscard]] const std::vector &Resolution() const; + [[nodiscard]] uint64_t GetPixelToBinChecksum() const; + [[nodiscard]] uint64_t GetCorrectionsChecksum() const; + // Pixels the spot finders must ignore because their resolution falls outside the limits, packed + // 32 pixels to a word (bit set = ignore), in the finders' own layout. + [[nodiscard]] std::shared_ptr> + ResolutionMaskBits(std::optional high_res, std::optional low_res) const; [[nodiscard]] size_t GetWidth() const; [[nodiscard]] size_t GetHeight() const; [[nodiscard]] const AzimuthalIntegrationSettings& Settings() const; diff --git a/common/PixelMask.cpp b/common/PixelMask.cpp index 8b34b2ec..c59ef2dc 100644 --- a/common/PixelMask.cpp +++ b/common/PixelMask.cpp @@ -3,13 +3,16 @@ #include "PixelMask.h" #include "RawToConvertedGeometry.h" +#include "TableChecksum.h" #include "JFJochException.h" #include "JFJochCompressor.h" PixelMask::PixelMask() = default; PixelMask::PixelMask(size_t width, size_t height) - : mask(width*height, 0) {} + : mask(width*height, 0) { + UpdateBinaryMask(); +} PixelMask::PixelMask(const DiffractionExperiment &experiment) : PixelMask(experiment.GetXPixelsNumConv(), @@ -17,7 +20,9 @@ PixelMask::PixelMask(const DiffractionExperiment &experiment) CalcEdgePixels(experiment); } -PixelMask::PixelMask(const std::vector &in_mask) : mask(in_mask) {} +PixelMask::PixelMask(const std::vector &in_mask) : mask(in_mask) { + UpdateBinaryMask(); +} uint32_t PixelMask::LoadMask(const std::vector &input_mask, uint8_t bit) { uint32_t ret = 0; @@ -35,7 +40,7 @@ uint32_t PixelMask::LoadMask(const std::vector &input_mask, uint8_t bi return ret; } -void PixelMask::UpdateRawMask(const DiffractionExperiment &experiment) { +void PixelMask::UpdateDerived(const DiffractionExperiment &experiment) { switch (experiment.GetDetectorType()) { case DetectorType::JUNGFRAU: case DetectorType::EIGER: @@ -46,6 +51,14 @@ void PixelMask::UpdateRawMask(const DiffractionExperiment &experiment) { raw_mask.clear(); break; } + UpdateBinaryMask(); +} + +void PixelMask::UpdateBinaryMask() { + binary_mask.resize(mask.size()); + for (size_t i = 0; i < mask.size(); i++) + binary_mask[i] = (mask[i] != 0); + binary_mask_checksum = TableChecksum(binary_mask.data(), binary_mask.size()); } void PixelMask::CalcEdgePixels_i(const DiffractionExperiment &experiment) { @@ -96,7 +109,7 @@ void PixelMask::CalcEdgePixels_i(const DiffractionExperiment &experiment) { void PixelMask::CalcEdgePixels(const DiffractionExperiment &experiment) { CalcEdgePixels_i(experiment); - UpdateRawMask(experiment); + UpdateDerived(experiment); } const std::vector &PixelMask::GetMaskRaw() const { @@ -110,6 +123,14 @@ const std::vector &PixelMask::GetMask() const { return mask; } +const std::vector &PixelMask::GetBinaryMask() const { + return binary_mask; +} + +uint64_t PixelMask::GetBinaryMaskChecksum() const { + return binary_mask_checksum; +} + const std::vector &PixelMask::GetMask(const DiffractionExperiment& experiment) const { if (experiment.IsGeometryTransformed()) return GetMask(); @@ -180,7 +201,7 @@ void PixelMask::LoadDetectorBadPixelMask(const DiffractionExperiment &experiment CalcEdgePixels_i(experiment); - UpdateRawMask(experiment); + UpdateDerived(experiment); } PixelMaskStatistics PixelMask::GetStatistics() const { @@ -207,12 +228,12 @@ PixelMaskStatistics PixelMask::GetStatistics() const { void PixelMask::LoadUserMask(const DiffractionExperiment& experiment, const std::vector &in_mask) { if (in_mask.size() == mask.size()) { LoadMask(in_mask, UserMaskedPixelBit); - UpdateRawMask(experiment); + UpdateDerived(experiment); } else if (in_mask.size() == experiment.GetModulesNum() * RAW_MODULE_SIZE) { std::vector tmp(experiment.GetPixelsNumConv(), 0); RawToConvertedGeometry(experiment, tmp.data(), in_mask. data()); LoadMask(tmp, UserMaskedPixelBit); - UpdateRawMask(experiment); + UpdateDerived(experiment); } else throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Size of input user mask invalid"); @@ -223,13 +244,13 @@ void PixelMask::LoadBeamStopMask(const DiffractionExperiment& experiment, const throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Size of input beam stop mask invalid"); LoadMask(in_mask, BeamStopPixelBit); - UpdateRawMask(experiment); + UpdateDerived(experiment); } void PixelMask::ClearBeamStopMask(const DiffractionExperiment& experiment) { for (auto &i: mask) i &= ~(1u << BeamStopPixelBit); - UpdateRawMask(experiment); + UpdateDerived(experiment); } void PixelMask::LoadUserMask(const DiffractionExperiment& experiment, const CompressedImage& image) { @@ -307,6 +328,7 @@ void PixelMask::LoadDECTRISBadPixelMask(const std::vector &input_mask) } } raw_mask = {}; // For DECTRIS - there is no raw mask + UpdateBinaryMask(); } void PixelMask::LoadDarkBadPixelMask(const DiffractionExperiment& experiment, const std::vector &input_mask) { @@ -325,5 +347,5 @@ void PixelMask::LoadDarkBadPixelMask(const DiffractionExperiment& experiment, co mask[i] &= ~(1 << NoisyPixelBit); } } - UpdateRawMask(experiment); + UpdateDerived(experiment); } diff --git a/common/PixelMask.h b/common/PixelMask.h index 90bf8bcd..eb7b2534 100644 --- a/common/PixelMask.h +++ b/common/PixelMask.h @@ -20,9 +20,17 @@ struct PixelMaskStatistics { class PixelMask { std::vector mask; std::vector raw_mask; + // One byte per pixel, 1 where the pixel is masked at all - the form the GPU image preprocessor + // uploads - and the checksum the shared device-table cache keys on (CudaSharedTables.h). Both are + // pure functions of `mask`, and an analysis engine is built per worker per pass, so they are + // derived here once instead of in every one of those engines. + std::vector binary_mask; + uint64_t binary_mask_checksum = 0; uint32_t LoadMask(const std::vector& mask, uint8_t bit); - void UpdateRawMask(const DiffractionExperiment& experiment); + // Everything that follows from the mask, recomputed wherever the mask changes. + void UpdateDerived(const DiffractionExperiment& experiment); + void UpdateBinaryMask(); void CalcEdgePixels_i(const DiffractionExperiment& experiment); @@ -55,6 +63,8 @@ public: [[nodiscard]] const std::vector &GetMaskRaw() const; [[nodiscard]] const std::vector &GetMask(const DiffractionExperiment& experiment) const; [[nodiscard]] const std::vector &GetMask() const; + [[nodiscard]] const std::vector &GetBinaryMask() const; + [[nodiscard]] uint64_t GetBinaryMaskChecksum() const; [[nodiscard]] std::vector GetUserMask(const DiffractionExperiment& experiment) const; [[nodiscard]] std::vector GetUserMask() const; diff --git a/common/TableChecksum.h b/common/TableChecksum.h new file mode 100644 index 00000000..aad6a8b6 --- /dev/null +++ b/common/TableChecksum.h @@ -0,0 +1,40 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +// FNV-1a over a block of bytes. Not a cryptographic hash and does not need to be - it exists to +// notice that the bytes behind a reused address changed, not to resist anyone. Its user is the +// shared device-table cache (image_analysis/indexing/CudaSharedTables.h), which keys on it; it lives +// here so the classes that own those tables can compute their own checksum once and hand it in. +// +// Run over eight interleaved lanes and fold them at the end. FNV's multiply is a loop-carried +// dependency, so one lane retires a byte every few cycles however much memory bandwidth is going +// spare; eight independent chains fill that latency. +inline uint64_t TableChecksum(const void *data, size_t bytes) { + constexpr uint64_t PRIME = 1099511628211ULL; + constexpr size_t LANES = 8; + const auto *p = static_cast(data); + uint64_t h[LANES]; + for (size_t l = 0; l < LANES; l++) + h[l] = 1469598103934665603ULL + l; + const size_t n = bytes / LANES * LANES; + for (size_t i = 0; i < n; i += LANES) + for (size_t l = 0; l < LANES; l++) { + h[l] ^= p[i + l]; + h[l] *= PRIME; + } + uint64_t out = 1469598103934665603ULL; + for (size_t l = 0; l < LANES; l++) { + out ^= h[l]; + out *= PRIME; + } + for (size_t i = n; i < bytes; i++) { + out ^= p[i]; + out *= PRIME; + } + return out; +} diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 8ad97ca5..e430407a 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -43,7 +43,6 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp indexer(in_indexer), prediction(CreateBraggPrediction(experiment.IsRotationIndexing())), mask(in_mask), - mask_resolution(experiment.GetPixelsNum(), false), mask_high_res(-1), mask_low_res(-1) { #ifdef JFJOCH_USE_CUDA @@ -301,8 +300,8 @@ ImageSpotFinder &MXAnalysisWithoutFPGA::FixedThresholdFinder() { experiment.GetYPixelsNum()); // It missed every mask update that happened before it existed, so it takes the current one // now. Without this it would find spots outside the resolution limits. - if (mask_high_res.has_value() || mask_low_res.has_value()) - spotFinder->SetResolutionMask(mask_resolution); + if (mask_resolution) + spotFinder->SetResolutionMaskBits(*mask_resolution); } return *spotFinder; } @@ -361,18 +360,14 @@ void MXAnalysisWithoutFPGA::RunROIOnly(DataMessage &output) { void MXAnalysisWithoutFPGA::UpdateMaskResolution(const SpotFindingSettings &settings) { mask_low_res = settings.low_resolution_limit; mask_high_res = settings.high_resolution_limit; - // An unset limit masks nothing at that end. At the high-resolution end 0 does that on its own - no - // pixel has d < 0, and the detector's own edge is where the pixels stop anyway; at the low-resolution - // end every pixel lies above any finite stand-in, so it takes an infinite one. - const float high_res = mask_high_res.value_or(0.0f); - const float low_res = mask_low_res.value_or(INFINITY); - auto const &resolution_map = integration.Resolution(); - for (int i = 0; i < mask_resolution.size(); i++) - mask_resolution[i] = (resolution_map[i] > low_res) || (resolution_map[i] < high_res); + // The mask is a pure function of the resolution map and the two limits, so the mapping builds it - + // once for all the workers, which otherwise each walked every pixel of the detector to arrive at + // the same bits. + mask_resolution = integration.ResolutionMaskBits(mask_high_res, mask_low_res); // The finders keep their own copy (the GPU ones a bit-packed device copy), so the mask is handed // over here - when the limits change - rather than with every image. if (spotFinder) - spotFinder->SetResolutionMask(mask_resolution); - adaptiveSpotFinder->SetResolutionMask(mask_resolution); + spotFinder->SetResolutionMaskBits(*mask_resolution); + adaptiveSpotFinder->SetResolutionMaskBits(*mask_resolution); } diff --git a/image_analysis/MXAnalysisWithoutFPGA.h b/image_analysis/MXAnalysisWithoutFPGA.h index 7af8249b..e77955f1 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.h +++ b/image_analysis/MXAnalysisWithoutFPGA.h @@ -64,7 +64,9 @@ class MXAnalysisWithoutFPGA { // not compressed) and return where it landed. const uint8_t *Decompress(const CompressedImage &image); - std::vector mask_resolution; + // Pixels outside the resolution limits, bit-packed. Built by the integration mapping, which is + // shared by every worker's engine and hands out the same mask to all of them. + std::shared_ptr> mask_resolution; // The limits mask_resolution was built for. Kept as the OPTIONAL the caller passed, so an unset // high-resolution limit compares equal to itself and the mask is not rebuilt on every image. std::optional mask_high_res; diff --git a/image_analysis/azint/AzIntEngineGPU.cu b/image_analysis/azint/AzIntEngineGPU.cu index dce1d43c..9ea95342 100644 --- a/image_analysis/azint/AzIntEngineGPU.cu +++ b/image_analysis/azint/AzIntEngineGPU.cu @@ -111,9 +111,11 @@ AzIntEngineGPU::AzIntEngineGPU(const AzimuthalIntegrationMapping &integration, s // Geometry-only, so shared per GPU: the first engine on this device uploads them, the rest reuse // them. Keyed by the mapping's own vectors, which outlive every engine built from it. gpu_azint_correction = SharedDeviceTable(integration.Corrections().data(), npixel, - integration.Corrections().data(), *stream); + integration.Corrections().data(), + integration.GetCorrectionsChecksum(), *stream); gpu_pixel_to_bin = SharedDeviceTable(integration.GetPixelToBin().data(), npixel, - integration.GetPixelToBin().data(), *stream); + integration.GetPixelToBin().data(), + integration.GetPixelToBinChecksum(), *stream); } void AzIntEngineGPU::Run(const ImagePreprocessorBuffer &image, AzimuthalIntegrationProfile &profile) { diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp index 3136cd17..c45df8fc 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp @@ -5,8 +5,11 @@ #include #include +#include +#include #include #include +#include #include "../../common/JFJochMath.h" // PI (M_PI is not standard, and MSVC does not define it) @@ -32,6 +35,24 @@ double parallax_var_px2(const std::string &material, double thickness_um, double return var / (pixel_um * pixel_um); } +// The radial-offset kernels below are a pure function of these six numbers, and one engine is built +// per worker per pass - 96 of them on a two-pass run - so the table was built 96 times over from the +// same inputs. Build it once and let the rest copy it; it is a few hundred floats. Two workers can +// still race to build the same table, which costs nothing but the second build: the values are +// identical, and emplace keeps whichever arrived first. +struct RadialKernelKey { + float r1_sq, r2, r3; + int n_kern, k_off, k_len; + + bool operator<(const RadialKernelKey &o) const { + return std::tie(r1_sq, r2, r3, n_kern, k_off, k_len) + < std::tie(o.r1_sq, o.r2, o.r3, o.n_kern, o.k_off, o.k_len); + } +}; + +std::mutex radial_kernel_mutex; +std::map> radial_kernel_cache; + } // namespace BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &experiment) @@ -128,10 +149,19 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe // above grow_max. k_off = static_cast(std::ceil(r3 + std::max(grow_max, n_kern - 1))) + 1; k_len = 2 * k_off + 1; - k_diff.clear(); - k_diff.reserve(static_cast(n_kern) * k_len); - for (int j = 0; j < n_kern; ++j) - BuildRadialKernel(static_cast(j)); + const RadialKernelKey kernel_key{r1_sq, r2, r3, n_kern, k_off, k_len}; + { + const std::lock_guard lock(radial_kernel_mutex); + if (const auto it = radial_kernel_cache.find(kernel_key); it != radial_kernel_cache.end()) + k_diff = it->second; + } + if (k_diff.empty()) { + k_diff.reserve(static_cast(n_kern) * k_len); + for (int j = 0; j < n_kern; ++j) + BuildRadialKernel(static_cast(j)); + const std::lock_guard lock(radial_kernel_mutex); + radial_kernel_cache.emplace(kernel_key, k_diff); + } polarization = experiment.GetPolarizationFactor(); } diff --git a/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu b/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu index 1dfb2052..19b86924 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu +++ b/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu @@ -235,17 +235,12 @@ ImagePreprocessorGPU::ImagePreprocessorGPU(const DiffractionExperiment &experime gpu_stats(1), cpu_stats(1), cpu_stats_reg(cpu_stats) { - // Setup mask. The same for every worker, so it is uploaded once per GPU and shared; keyed on the - // PixelMask's own vector, which the derived table is a pure function of. - // Hoist the accessor and index without the bounds check: this runs once per worker over every - // pixel of the detector - 18 million times per engine on a 16 Mpx one, and an engine is built per - // worker per pass - and .at() on each of them stops the loop vectorising for a bound the loop - // itself already respects. - const std::vector &mask_raw = mask.GetMask(); - std::vector mask_vec(npixels); - for (int i = 0; i < npixels; i++) - mask_vec[i] = (mask_raw[i] != 0); - gpu_mask = SharedDeviceTable(mask.GetMask().data(), npixels, mask_vec.data(), *stream); + // Setup mask. The byte-per-pixel form and its checksum come from the PixelMask, which derives them + // whenever the mask changes: they are the same for every worker, and deriving them walks every + // pixel of the detector - 18 million of them on a 16 Mpx one, per engine, with an engine built per + // worker per pass. The table is then uploaded once per GPU and shared by the engines on it. + gpu_mask = SharedDeviceTable(mask.GetBinaryMask().data(), npixels, mask.GetBinaryMask().data(), + mask.GetBinaryMaskChecksum(), *stream); // Setup GPU settings. The current device, not device 0: workers are pinned round-robin across GPUs, // so device 0's SM count can belong to a different card than the one these kernels launch on. diff --git a/image_analysis/indexing/CudaSharedTables.h b/image_analysis/indexing/CudaSharedTables.h index 3c19d5fb..d62d227b 100644 --- a/image_analysis/indexing/CudaSharedTables.h +++ b/image_analysis/indexing/CudaSharedTables.h @@ -10,6 +10,7 @@ #include #include "CUDAMemHelpers.h" +#include "../../common/TableChecksum.h" // Read-only lookup tables that depend only on the detector geometry (pixel -> azimuthal bin, the // per-pixel correction factors, the pixel mask). One analysis engine is built per worker thread, so @@ -26,8 +27,16 @@ // A bare address is not enough on its own to say "same table", though: a host buffer can be mutated // in place, or freed and a new one allocated at the same address, and either would hand the caller a // device copy of something else - silently, since the data is only ever read. So the byte length and -// a checksum of the bytes actually uploaded are part of the key too. Both are computed once per -// engine construction, against an upload of the same buffer, so they cost nothing measurable. +// a checksum of the bytes actually uploaded are part of the key too. +// +// The checksum is part of the KEY, so it is computed before the lookup and a cache hit pays for it as +// well: the tables are tens of megabytes and one engine is built per worker per pass, which put the +// hashing alone at several percent of all CPU samples. The overload below therefore takes a checksum +// the caller already has. That does not weaken anything, because it is the buffer's OWNER that +// computes it - AzimuthalIntegrationMapping and PixelMask hash their tables whenever they write them +// and hand out the result - so the checksum still describes the bytes as they are now. What must not +// be done instead is to memoise the checksum on the address here: that is precisely the reused-address +// hazard the checksum exists to catch. // // Entries are held weakly, so the tables are released once the last engine using them is gone. @@ -40,39 +49,6 @@ namespace jfjoch_cuda_shared_tables { std::map> tables; }; - // FNV-1a. Not a cryptographic hash and does not need to be - it exists to notice that the bytes - // behind a reused address changed, not to resist anyone. - // - // Run over eight interleaved lanes and fold them at the end. FNV's multiply is a loop-carried - // dependency, so one lane retires a byte every few cycles however much memory bandwidth is - // going spare; eight independent chains fill that latency. The tables are tens of megabytes and - // every worker thread hashes them as it builds its engines, which put this at ~5% of all CPU - // samples on a 16M-pixel detector. - inline uint64_t checksum(const void *data, size_t bytes) { - constexpr uint64_t PRIME = 1099511628211ULL; - constexpr size_t LANES = 8; - const auto *p = static_cast(data); - uint64_t h[LANES]; - for (size_t l = 0; l < LANES; l++) - h[l] = 1469598103934665603ULL + l; - const size_t n = bytes / LANES * LANES; - for (size_t i = 0; i < n; i += LANES) - for (size_t l = 0; l < LANES; l++) { - h[l] ^= p[i + l]; - h[l] *= PRIME; - } - uint64_t out = 1469598103934665603ULL; - for (size_t l = 0; l < LANES; l++) { - out ^= h[l]; - out *= PRIME; - } - for (size_t i = n; i < bytes; i++) { - out ^= p[i]; - out *= PRIME; - } - return out; - } - inline Registry ®istry() { static Registry r; return r; @@ -87,16 +63,16 @@ namespace jfjoch_cuda_shared_tables { } // Return the device-resident copy of `host` (`count` elements) for the calling thread's GPU, -// uploading it on `stream` the first time it is asked for. +// uploading it on `stream` the first time it is asked for. `host_checksum` is TableChecksum() of +// those `count * sizeof(T)` bytes, computed by whoever owns the buffer. template std::shared_ptr> SharedDeviceTable(const void *key, size_t count, const T *host, - cudaStream_t stream) { + uint64_t host_checksum, cudaStream_t stream) { int device = 0; jfjoch_cuda_shared_tables::check(cudaGetDevice(&device)); const size_t bytes = count * sizeof(T); - const jfjoch_cuda_shared_tables::TableKey table_key{ - device, key, bytes, jfjoch_cuda_shared_tables::checksum(host, bytes)}; + const jfjoch_cuda_shared_tables::TableKey table_key{device, key, bytes, host_checksum}; auto ® = jfjoch_cuda_shared_tables::registry(); // The upload happens while the lock is held: another worker must not obtain the pointer before @@ -127,3 +103,10 @@ std::shared_ptr> SharedDeviceTable(const void *key, size_t coun reg.tables[table_key] = std::shared_ptr(table); return table; } + +// Same, for a caller with no checksum of its own to hand over. +template +std::shared_ptr> SharedDeviceTable(const void *key, size_t count, const T *host, + cudaStream_t stream) { + return SharedDeviceTable(key, count, host, TableChecksum(host, count * sizeof(T)), stream); +} diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu index 6ba4a8f2..65ed5123 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu @@ -347,9 +347,11 @@ AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping & // Both tables are functions of the detector geometry alone, so they are uploaded once per GPU and // shared: the azimuthal-integration engine in the same worker reads the very same two arrays. gpu_pixel_to_bin = SharedDeviceTable(mapping.GetPixelToBin().data(), npix, - mapping.GetPixelToBin().data(), *stream); + mapping.GetPixelToBin().data(), + mapping.GetPixelToBinChecksum(), *stream); gpu_corrections = SharedDeviceTable(mapping.Corrections().data(), npix, - mapping.Corrections().data(), *stream); + mapping.Corrections().data(), + mapping.GetCorrectionsChecksum(), *stream); } void AdaptiveSpotFinderGPU::ReducePass(const ImagePreprocessorBuffer &image, float clip_k, @@ -476,8 +478,8 @@ void AdaptiveSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image, // so the ordering already guarantees flag_strong has finished. Waiting here only idled the host. } -void AdaptiveSpotFinderGPU::SetResolutionMask(const std::vector &mask) { - ImageSpotFinder::SetResolutionMask(mask); +void AdaptiveSpotFinderGPU::SetResolutionMaskBits(const std::vector &packed_mask) { + ImageSpotFinder::SetResolutionMaskBits(packed_mask); extractor.SetResolutionMask(res_mask_bits); } diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h index f6358505..24e8cb74 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.h @@ -114,7 +114,7 @@ public: AdaptiveSpotFinderGPU &operator=(const AdaptiveSpotFinderGPU &) = delete; void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override; - void SetResolutionMask(const std::vector &mask) override; + void SetResolutionMaskBits(const std::vector &packed_mask) override; const std::vector &ExtractComponents(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override; diff --git a/image_analysis/spot_finding/ImageSpotFinder.cpp b/image_analysis/spot_finding/ImageSpotFinder.cpp index d50063f6..47a44d3f 100644 --- a/image_analysis/spot_finding/ImageSpotFinder.cpp +++ b/image_analysis/spot_finding/ImageSpotFinder.cpp @@ -1,7 +1,6 @@ // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only -#include #include #include "../../common/JFJochException.h" @@ -33,10 +32,19 @@ void ImageSpotFinder::SetResolutionMask(const std::vector &mask) { if (mask.size() != npixel) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "ImageSpotFinder::SetResolutionMask: mask size mismatch"); - std::fill(res_mask_bits.begin(), res_mask_bits.end(), 0); + std::vector packed(OutputSize(), 0); for (size_t i = 0; i < npixel; i++) if (mask[i]) - res_mask_bits[i / 32] |= 1u << (i % 32); + packed[i / 32] |= 1u << (i % 32); + SetResolutionMaskBits(packed); +} + +void ImageSpotFinder::SetResolutionMaskBits(const std::vector &packed_mask) { + if (packed_mask.size() != OutputSize()) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "ImageSpotFinder::SetResolutionMaskBits: mask size mismatch"); + res_mask_bits = packed_mask; + const size_t npixel = static_cast(width) * height; if (npixel % 32 != 0) res_mask_bits.back() |= ~((1u << (npixel % 32)) - 1u); } diff --git a/image_analysis/spot_finding/ImageSpotFinder.h b/image_analysis/spot_finding/ImageSpotFinder.h index 781939a2..6566a84d 100644 --- a/image_analysis/spot_finding/ImageSpotFinder.h +++ b/image_analysis/spot_finding/ImageSpotFinder.h @@ -53,7 +53,11 @@ public: // Pixels to ignore, one bool per pixel (true = ignore). Set when the resolution limits change, // not per image: the GPU finders keep a bit-packed device copy of it, and re-uploading that for // every image would cost more than the extraction it feeds. - virtual void SetResolutionMask(const std::vector &mask); + void SetResolutionMask(const std::vector &mask); + // The same mask already packed 32 pixels to a word, which is how the finders keep it. Every + // worker's finder is given the identical mask, so the packing is done once by whoever owns the + // resolution map rather than by each of them (AzimuthalIntegrationMapping::ResolutionMaskBits). + virtual void SetResolutionMaskBits(const std::vector &packed_mask); // Every connected component of the last Detect() with at most max-pix pixels. min-pix is NOT // applied here on purpose - it is the only spot setting that changes between the passes of the diff --git a/image_analysis/spot_finding/ImageSpotFinderGPU.cu b/image_analysis/spot_finding/ImageSpotFinderGPU.cu index 27f3f82a..4270d165 100644 --- a/image_analysis/spot_finding/ImageSpotFinderGPU.cu +++ b/image_analysis/spot_finding/ImageSpotFinderGPU.cu @@ -258,8 +258,8 @@ ImageSpotFinderGPU::ImageSpotFinderGPU(int32_t in_width, int32_t in_height, gpu_out_1 = CudaDevicePtr(OutputSize()); } -void ImageSpotFinderGPU::SetResolutionMask(const std::vector &mask) { - ImageSpotFinder::SetResolutionMask(mask); +void ImageSpotFinderGPU::SetResolutionMaskBits(const std::vector &packed_mask) { + ImageSpotFinder::SetResolutionMaskBits(packed_mask); extractor.SetResolutionMask(res_mask_bits); } diff --git a/image_analysis/spot_finding/ImageSpotFinderGPU.h b/image_analysis/spot_finding/ImageSpotFinderGPU.h index 5f13f238..df9e7662 100644 --- a/image_analysis/spot_finding/ImageSpotFinderGPU.h +++ b/image_analysis/spot_finding/ImageSpotFinderGPU.h @@ -25,7 +25,7 @@ public: ~ImageSpotFinderGPU() override = default; void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override; - void SetResolutionMask(const std::vector &mask) override; + void SetResolutionMaskBits(const std::vector &packed_mask) override; const std::vector &ExtractComponents(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override; };