Build the detector's lookup tables once, not once per worker

The image loop gives every worker its own analysis engine, so a run builds ninety-six of them. Each
one derived, from scratch, tables that are the same in all of them: the byte-per-pixel mask, the
resolution mask, the radial kernel, and the checksum that names the shared device tables.

The checksum was the worst of it, because it is part of the cache KEY and so is computed before the
lookup - a hit still hashed the whole table. On a 16 Mpx detector that is the bin table, the
corrections and the mask, 126 MB an engine, about twelve gigabytes over a run, to answer a question
whose answer had not changed. The header said it cost nothing measurable; a profile says otherwise,
and says it is worst exactly during the ramp when the machine has nothing else to do.

It cannot simply be remembered against the address, which is what it exists to catch: a buffer can
be freed and another allocated where it was, and the cache would then hand back a device copy of
something else. So the owner of the bytes computes it instead. The azimuthal mapping writes its two
tables in its constructor and never again. The pixel mask re-derives its binary form and its
checksum on every path that changes the mask, and all of those paths are now private to the class.
The key therefore still describes the bytes as they are at the moment of the lookup.

The resolution mask was two passes over every pixel - a float comparison into a vector<bool>, then a
bit-by-bit repack - in each of the ninety-six. It is one pass now, writing the packed form directly,
built once for the limits asked for and handed out as a shared pointer so a worker keeps the mask it
was given. The radial kernel is cached on the six numbers it is derived from.

Nothing computes a different value; only who computes it changes. Byte-identical merged output on a
16 Mpx set and on a small one.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011n8riB6X59oRjkrSHzNPAU
This commit is contained in:
jungfrau
2026-08-23 12:59:58 -04:00
co-authored by Claude Opus 5
parent 27020d27e9
commit 5ee0f22a61
17 changed files with 248 additions and 93 deletions
+37
View File
@@ -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<uint32_t> &mask,
@@ -207,6 +211,39 @@ const std::vector<float> &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<const std::vector<uint32_t>>
AzimuthalIntegrationMapping::ResolutionMaskBits(std::optional<float> high_res,
std::optional<float> 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<std::vector<uint32_t>>(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;
}
+25
View File
@@ -3,6 +3,8 @@
#pragma once
#include <memory>
#include <mutex>
#include <optional>
#include "DiffractionExperiment.h"
#include "PixelMask.h"
@@ -24,6 +26,23 @@ protected:
std::optional<float> 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<float> res_mask_high, res_mask_low;
mutable std::shared_ptr<const std::vector<uint32_t>> res_mask_bits;
size_t nthreads;
void UpdateMaxBinNumber();
@@ -49,6 +68,12 @@ public:
[[nodiscard]] uint16_t QToBin(float q) const;
[[nodiscard]] const std::vector<float> &Corrections() const;
[[nodiscard]] const std::vector<float> &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<const std::vector<uint32_t>>
ResolutionMaskBits(std::optional<float> high_res, std::optional<float> low_res) const;
[[nodiscard]] size_t GetWidth() const;
[[nodiscard]] size_t GetHeight() const;
[[nodiscard]] const AzimuthalIntegrationSettings& Settings() const;
+32 -10
View File
@@ -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<uint32_t> &in_mask) : mask(in_mask) {}
PixelMask::PixelMask(const std::vector<uint32_t> &in_mask) : mask(in_mask) {
UpdateBinaryMask();
}
uint32_t PixelMask::LoadMask(const std::vector<uint32_t> &input_mask, uint8_t bit) {
uint32_t ret = 0;
@@ -35,7 +40,7 @@ uint32_t PixelMask::LoadMask(const std::vector<uint32_t> &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<uint32_t> &PixelMask::GetMaskRaw() const {
@@ -110,6 +123,14 @@ const std::vector<uint32_t> &PixelMask::GetMask() const {
return mask;
}
const std::vector<uint8_t> &PixelMask::GetBinaryMask() const {
return binary_mask;
}
uint64_t PixelMask::GetBinaryMaskChecksum() const {
return binary_mask_checksum;
}
const std::vector<uint32_t> &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<uint32_t> &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<uint32_t> 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<uint32_t> &input_mask)
}
}
raw_mask = {}; // For DECTRIS - there is no raw mask
UpdateBinaryMask();
}
void PixelMask::LoadDarkBadPixelMask(const DiffractionExperiment& experiment, const std::vector<uint32_t> &input_mask) {
@@ -325,5 +347,5 @@ void PixelMask::LoadDarkBadPixelMask(const DiffractionExperiment& experiment, co
mask[i] &= ~(1 << NoisyPixelBit);
}
}
UpdateRawMask(experiment);
UpdateDerived(experiment);
}
+11 -1
View File
@@ -20,9 +20,17 @@ struct PixelMaskStatistics {
class PixelMask {
std::vector<uint32_t> mask;
std::vector<uint32_t> 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<uint8_t> binary_mask;
uint64_t binary_mask_checksum = 0;
uint32_t LoadMask(const std::vector<uint32_t>& 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<uint32_t> &GetMaskRaw() const;
[[nodiscard]] const std::vector<uint32_t> &GetMask(const DiffractionExperiment& experiment) const;
[[nodiscard]] const std::vector<uint32_t> &GetMask() const;
[[nodiscard]] const std::vector<uint8_t> &GetBinaryMask() const;
[[nodiscard]] uint64_t GetBinaryMaskChecksum() const;
[[nodiscard]] std::vector<uint32_t> GetUserMask(const DiffractionExperiment& experiment) const;
[[nodiscard]] std::vector<uint32_t> GetUserMask() const;
+40
View File
@@ -0,0 +1,40 @@
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#pragma once
#include <cstddef>
#include <cstdint>
// 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<const unsigned char *>(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;
}
+8 -13
View File
@@ -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);
}
+3 -1
View File
@@ -64,7 +64,9 @@ class MXAnalysisWithoutFPGA {
// not compressed) and return where it landed.
const uint8_t *Decompress(const CompressedImage &image);
std::vector<bool> 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<const std::vector<uint32_t>> 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<float> mask_high_res;
+4 -2
View File
@@ -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) {
@@ -5,8 +5,11 @@
#include <algorithm>
#include <cmath>
#include <map>
#include <mutex>
#include <numeric>
#include <string>
#include <tuple>
#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<RadialKernelKey, std::vector<float>> 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<int>(std::ceil(r3 + std::max<double>(grow_max, n_kern - 1))) + 1;
k_len = 2 * k_off + 1;
k_diff.clear();
k_diff.reserve(static_cast<size_t>(n_kern) * k_len);
for (int j = 0; j < n_kern; ++j)
BuildRadialKernel(static_cast<float>(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<size_t>(n_kern) * k_len);
for (int j = 0; j < n_kern; ++j)
BuildRadialKernel(static_cast<float>(j));
const std::lock_guard lock(radial_kernel_mutex);
radial_kernel_cache.emplace(kernel_key, k_diff);
}
polarization = experiment.GetPolarizationFactor();
}
@@ -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<uint32_t> &mask_raw = mask.GetMask();
std::vector<uint8_t> 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.
+22 -39
View File
@@ -10,6 +10,7 @@
#include <utility>
#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<TableKey, std::weak_ptr<void>> 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<const unsigned char *>(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 &registry() {
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 <typename T>
std::shared_ptr<CudaDevicePtr<T>> 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 &reg = 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<CudaDevicePtr<T>> SharedDeviceTable(const void *key, size_t coun
reg.tables[table_key] = std::shared_ptr<void>(table);
return table;
}
// Same, for a caller with no checksum of its own to hand over.
template <typename T>
std::shared_ptr<CudaDevicePtr<T>> SharedDeviceTable(const void *key, size_t count, const T *host,
cudaStream_t stream) {
return SharedDeviceTable(key, count, host, TableChecksum(host, count * sizeof(T)), stream);
}
@@ -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<bool> &mask) {
ImageSpotFinder::SetResolutionMask(mask);
void AdaptiveSpotFinderGPU::SetResolutionMaskBits(const std::vector<uint32_t> &packed_mask) {
ImageSpotFinder::SetResolutionMaskBits(packed_mask);
extractor.SetResolutionMask(res_mask_bits);
}
@@ -114,7 +114,7 @@ public:
AdaptiveSpotFinderGPU &operator=(const AdaptiveSpotFinderGPU &) = delete;
void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override;
void SetResolutionMask(const std::vector<bool> &mask) override;
void SetResolutionMaskBits(const std::vector<uint32_t> &packed_mask) override;
const std::vector<DiffractionSpot> &ExtractComponents(const ImagePreprocessorBuffer &image,
const SpotFindingSettings &settings) override;
@@ -1,7 +1,6 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include <algorithm>
#include <bit>
#include "../../common/JFJochException.h"
@@ -33,10 +32,19 @@ void ImageSpotFinder::SetResolutionMask(const std::vector<bool> &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<uint32_t> 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<uint32_t> &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<size_t>(width) * height;
if (npixel % 32 != 0)
res_mask_bits.back() |= ~((1u << (npixel % 32)) - 1u);
}
@@ -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<bool> &mask);
void SetResolutionMask(const std::vector<bool> &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<uint32_t> &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
@@ -258,8 +258,8 @@ ImageSpotFinderGPU::ImageSpotFinderGPU(int32_t in_width, int32_t in_height,
gpu_out_1 = CudaDevicePtr<uint32_t>(OutputSize());
}
void ImageSpotFinderGPU::SetResolutionMask(const std::vector<bool> &mask) {
ImageSpotFinder::SetResolutionMask(mask);
void ImageSpotFinderGPU::SetResolutionMaskBits(const std::vector<uint32_t> &packed_mask) {
ImageSpotFinder::SetResolutionMaskBits(packed_mask);
extractor.SetResolutionMask(res_mask_bits);
}
@@ -25,7 +25,7 @@ public:
~ImageSpotFinderGPU() override = default;
void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override;
void SetResolutionMask(const std::vector<bool> &mask) override;
void SetResolutionMaskBits(const std::vector<uint32_t> &packed_mask) override;
const std::vector<DiffractionSpot> &ExtractComponents(const ImagePreprocessorBuffer &image,
const SpotFindingSettings &settings) override;
};