image_analysis: share the read-only GPU lookup tables per device

One analysis engine is built per worker thread, and each uploaded its own copy of
tables that are pure functions of the detector geometry: the pixel -> azimuthal bin
map and the per-pixel corrections (both in AzIntEngineGPU AND again in
AdaptiveSpotFinderGPU, from the same mapping), plus the pixel mask. On an 18 Mpx
detector that is ~224 MB per worker; with 32 workers ~7 GB of device memory held 32
identical copies.

Upload each table once per GPU instead and hand every engine on that device a shared
pointer to it. The cache is keyed by (device, source-vector address) because workers
are pinned round-robin across GPUs, so on a multi-GPU node each device keeps its own
copy - a kernel may only read memory resident on the device it runs on - and the
table is freed on the device that allocated it. Entries are held weakly, so a table
goes away with the last engine using it.

Measured on an 18 Mpx detector, 32 worker threads, 16 GB card: the stills path went
from exhausting the card (OOM in de-novo indexing) to 8.6 GB peak, and a normal
rotation run from 14.6 GB to 7.4 GB - it had been running within 1.6 GB of the limit,
so any larger detector or second GPU consumer would have tipped it over. Per-worker
footprint drops 403 -> 173 MB. Merge statistics are unchanged on a six-crystal
regression subset, including two-pass runs where the second pass rebuilds the mapping
on refined geometry, and wall time is unchanged (13.5-13.8 s vs 13.8-14.1 s).

Also take the launch configuration from the current device rather than device 0 in
AzIntEngineGPU and ImagePreprocessorGPU: with round-robin pinning, device 0's SM count
and shared-memory size can belong to a different card than the one the kernels use.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

image_analysis/azint/AzIntEngineGPU.cu

@@ -91,8 +91,6 @@ void gpu_azim(
 AzIntEngineGPU::AzIntEngineGPU(const AzimuthalIntegrationMapping &integration, std::shared_ptr<CudaStream> stream)
     : AzIntEngine(integration),
       stream(stream),
-      gpu_azint_correction(npixel),
-      gpu_pixel_to_bin(npixel),
       gpu_sum(azint_bins),
       gpu_sum2(azint_bins),
       gpu_count(azint_bins),
@@ -100,22 +98,22 @@ AzIntEngineGPU::AzIntEngineGPU(const AzimuthalIntegrationMapping &integration, s
       cpu_sum2_reg(azint_sum2),
       cpu_count_reg(azint_count) {
 
+    int device = 0;
+    cuda_err(cudaGetDevice(&device));       // this worker's GPU, not necessarily 0
     cudaDeviceProp prop{};
-    cuda_err(cudaGetDeviceProperties(&prop, 0));
+    cuda_err(cudaGetDeviceProperties(&prop, device));
 
     threads = 128;
     blocks = 4 * prop.multiProcessorCount;
     shared_size = prop.sharedMemPerBlock;
     shared_needed = azint_bins * (2 * sizeof(float) + sizeof(uint32_t));
 
-    // On this engine's stream, like every other operation it issues: the streams are non-blocking, so a
-    // NULL-stream copy is no longer ordered against the kernels that read what it uploads. The one-time
-    // synchronise leaves the constructor with the uploads settled rather than in flight.
-    cuda_err(cudaMemcpyAsync(gpu_azint_correction, integration.Corrections().data(), sizeof(float) * npixel,
-                             cudaMemcpyHostToDevice, *stream));
-    cuda_err(cudaMemcpyAsync(gpu_pixel_to_bin, integration.GetPixelToBin().data(), sizeof(uint16_t) * npixel,
-                             cudaMemcpyHostToDevice, *stream));
-    cuda_err(cudaStreamSynchronize(*stream));
+    // 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);
+    gpu_pixel_to_bin = SharedDeviceTable(integration.GetPixelToBin().data(), npixel,
+                                         integration.GetPixelToBin().data(), *stream);
 }
 
 void AzIntEngineGPU::Run(const ImagePreprocessorBuffer &image, AzimuthalIntegrationProfile &profile) {
@@ -127,12 +125,12 @@ void AzIntEngineGPU::Run(const ImagePreprocessorBuffer &image, AzimuthalIntegrat
 
     if (shared_needed < shared_size) {
         gpu_azim_shared<<<blocks, threads, shared_needed, *stream>>>(
-        gpu_pixel_to_bin,gpu_azint_correction,image.getGPUBuffer(), gpu_sum, gpu_sum2,
+        gpu_pixel_to_bin->get(),gpu_azint_correction->get(),image.getGPUBuffer(), gpu_sum, gpu_sum2,
         gpu_count, npixel, azint_bins
         );
     } else {
         gpu_azim<<<blocks, threads, 0, *stream>>>(
-            gpu_pixel_to_bin,gpu_azint_correction,image.getGPUBuffer(), gpu_sum, gpu_sum2,
+            gpu_pixel_to_bin->get(),gpu_azint_correction->get(),image.getGPUBuffer(), gpu_sum, gpu_sum2,
             gpu_count, npixel, azint_bins
         );
     }

image_analysis/azint/AzIntEngineGPU.h

@@ -5,6 +5,7 @@
 
 #include "AzIntEngine.h"
 #include "../indexing/CUDAMemHelpers.h"
+#include "../indexing/CudaSharedTables.h"
 
 class AzIntEngineGPU : public AzIntEngine {
     std::shared_ptr<CudaStream> stream;
@@ -13,8 +14,10 @@ class AzIntEngineGPU : public AzIntEngine {
     size_t shared_needed;
     size_t shared_size;
 
-    CudaDevicePtr<float> gpu_azint_correction;
-    CudaDevicePtr<uint16_t> gpu_pixel_to_bin;
+    // Geometry-only tables: one copy per GPU, shared with every other engine on it (see
+    // CudaSharedTables.h) instead of one copy per worker thread.
+    std::shared_ptr<CudaDevicePtr<float>> gpu_azint_correction;
+    std::shared_ptr<CudaDevicePtr<uint16_t>> gpu_pixel_to_bin;
 
     CudaDevicePtr<float> gpu_sum;
     CudaDevicePtr<float> gpu_sum2;

image_analysis/image_preprocessing/ImagePreprocessorGPU.cu

@@ -93,25 +93,23 @@ ImagePreprocessorGPU::ImagePreprocessorGPU(const DiffractionExperiment &experime
                                            std::shared_ptr<CudaStream> stream)
     : ImagePreprocessor(experiment),
       stream(stream),
-      gpu_mask(npixels),
       gpu_decompressed_image(npixels * sizeof(uint32_t)),  // Overshoot - if input image is 1- or 2-byte, then it is still fine, while memory loss is minimal
       gpu_stats(1),
       cpu_stats(1),
       cpu_stats_reg(cpu_stats) {
-    // Setup mask
+    // 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.
     std::vector<uint8_t> mask_vec(npixels);
     for (int i = 0; i < npixels; i++)
         mask_vec[i] = (mask.GetMask().at(i) != 0);
+    gpu_mask = SharedDeviceTable(mask.GetMask().data(), npixels, mask_vec.data(), *stream);
 
-    // On this engine's stream, like every other operation it issues: the streams are non-blocking, so a
-    // NULL-stream copy is no longer ordered against the kernels that read the mask. Synchronise before
-    // leaving the constructor - mask_vec is a local and the copy must not outlive it.
-    cudaMemcpyAsync(gpu_mask, mask_vec.data(), npixels, cudaMemcpyHostToDevice, *stream);
-    cudaStreamSynchronize(*stream);
-
-    // Setup GPU settings
+    // 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.
+    int device = 0;
+    cudaGetDevice(&device);
     cudaDeviceProp prop{};
-    cudaGetDeviceProperties(&prop, 0);
+    cudaGetDeviceProperties(&prop, device);
 
     threads = 128;
     blocks = 4 * prop.multiProcessorCount;
@@ -154,7 +152,7 @@ ImageStatistics ImagePreprocessorGPU::Analyze(ImagePreprocessorBuffer &processed
     cudaMemcpyAsync(gpu_stats, cpu_stats.data(), sizeof(ImageStatistics), cudaMemcpyHostToDevice, *stream);
     preprocess_kernel<T> <<< blocks, threads, 0, *stream >>>(
         reinterpret_cast<const T *>(gpu_decompressed_image.get()),
-        gpu_mask,
+        gpu_mask->get(),
         processed_image.getGPUBuffer(),
         gpu_stats,
         sat_value,

image_analysis/image_preprocessing/ImagePreprocessorGPU.h

@@ -5,12 +5,14 @@
 
 #include "ImagePreprocessor.h"
 #include "../indexing/CUDAMemHelpers.h"
+#include "../indexing/CudaSharedTables.h"
 
 class ImagePreprocessorGPU : public ImagePreprocessor {
     std::shared_ptr<CudaStream> stream;
     int threads;
     int blocks;
-    CudaDevicePtr<uint8_t> gpu_mask;
+    // Geometry-only, so one copy per GPU shared with every other engine on it (CudaSharedTables.h).
+    std::shared_ptr<CudaDevicePtr<uint8_t>> gpu_mask;
     CudaDevicePtr<uint8_t> gpu_decompressed_image;
     CudaDevicePtr<ImageStatistics> gpu_stats;
 

image_analysis/indexing/CudaSharedTables.h

@@ -0,0 +1,77 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <map>
+#include <memory>
+#include <mutex>
+#include <utility>
+
+#include "CUDAMemHelpers.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
+// each of those used to upload its own copy: on a 18 Mpx detector that is ~220 MB per thread, and
+// 32 threads spent ~7 GB of device memory on identical data.
+//
+// Upload once per GPU instead and hand every engine on that GPU a shared pointer to the same table.
+// The cache is keyed by (device, key) because a worker thread is pinned round-robin to a device
+// (pin_gpu()), so on a multi-GPU node each device keeps its own copy - a kernel may only read memory
+// resident on the device it runs on. `key` identifies the table's source data; use the address of the
+// host vector that produced it, which lives in the experiment / integration mapping and therefore
+// outlives every engine.
+//
+// Entries are held weakly, so the tables are released once the last engine using them is gone.
+
+namespace jfjoch_cuda_shared_tables {
+    struct Registry {
+        std::mutex m;
+        std::map<std::pair<int, const void *>, std::weak_ptr<void>> tables;
+    };
+
+    inline Registry &registry() {
+        static Registry r;
+        return r;
+    }
+
+    // Not called cuda_err: the .cu files that include this header define their own such helper in an
+    // anonymous namespace, and a second one at global scope would make every call ambiguous.
+    inline void check(cudaError_t val) {
+        if (val != cudaSuccess)
+            throw JFJochException(JFJochExceptionCategory::GPUCUDAError, cudaGetErrorString(val));
+    }
+}
+
+// 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.
+template <typename T>
+std::shared_ptr<CudaDevicePtr<T>> SharedDeviceTable(const void *key, size_t count, const T *host,
+                                                    cudaStream_t stream) {
+    int device = 0;
+    jfjoch_cuda_shared_tables::check(cudaGetDevice(&device));
+
+    auto &reg = jfjoch_cuda_shared_tables::registry();
+    // The upload happens while the lock is held: another worker must not obtain the pointer before
+    // its content is on the device.
+    std::lock_guard lock(reg.m);
+    auto &slot = reg.tables[{device, key}];
+    if (auto cached = slot.lock())
+        return std::static_pointer_cast<CudaDevicePtr<T>>(cached);
+
+    // Free on the device that allocated it - the last engine to drop the table may well be a worker
+    // pinned to a different GPU.
+    std::shared_ptr<CudaDevicePtr<T>> table(new CudaDevicePtr<T>(count), [device](CudaDevicePtr<T> *p) {
+        int current = 0;
+        cudaGetDevice(&current);
+        cudaSetDevice(device);
+        delete p;
+        cudaSetDevice(current);
+    });
+    jfjoch_cuda_shared_tables::check(
+        cudaMemcpyAsync(table->get(), host, count * sizeof(T), cudaMemcpyHostToDevice, stream));
+    jfjoch_cuda_shared_tables::check(cudaStreamSynchronize(stream));
+
+    slot = std::shared_ptr<void>(table);
+    return table;
+}

image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu

@@ -170,8 +170,6 @@ AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping &
       stream(std::move(in_stream)),
       nbins(in_mapping.GetBinNumber()),
       npix(in_mapping.GetPixelToBin().size()),
-      gpu_pixel_to_bin(npix),
-      gpu_corrections(npix),
       gpu_sum(nbins),
       gpu_sum2(nbins),
       gpu_count(nbins),
@@ -203,14 +201,12 @@ AdaptiveSpotFinderGPU::AdaptiveSpotFinderGPU(const AzimuthalIntegrationMapping &
     shared_clip = static_cast<size_t>(nbins) * (2 * sizeof(float) + sizeof(uint32_t));
     use_shared = (shared_plain < prop.sharedMemPerBlock);
 
-    // On this engine's stream, like every other operation it issues: the streams are non-blocking, so a
-    // NULL-stream copy is no longer ordered against the kernels that read what it uploads. The one-time
-    // synchronise leaves the constructor with the uploads settled rather than in flight.
-    cuda_err(cudaMemcpyAsync(gpu_pixel_to_bin, mapping.GetPixelToBin().data(), sizeof(uint16_t) * npix,
-                             cudaMemcpyHostToDevice, *stream));
-    cuda_err(cudaMemcpyAsync(gpu_corrections, mapping.Corrections().data(), sizeof(float) * npix,
-                             cudaMemcpyHostToDevice, *stream));
-    cuda_err(cudaStreamSynchronize(*stream));
+    // 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);
+    gpu_corrections = SharedDeviceTable(mapping.Corrections().data(), npix,
+                                        mapping.Corrections().data(), *stream);
 }
 
 void AdaptiveSpotFinderGPU::ReducePass(const ImagePreprocessorBuffer &image, float clip_k,
@@ -218,12 +214,12 @@ void AdaptiveSpotFinderGPU::ReducePass(const ImagePreprocessorBuffer &image, flo
     if (use_shared) {
         const size_t shared = accumulate_corrected ? shared_plain : shared_clip;
         reduce_rings_shared<<<reduce_blocks, reduce_threads, shared, *stream>>>(
-            gpu_pixel_to_bin, gpu_corrections, image.getGPUBuffer(), gpu_mean, gpu_sigma,
+            gpu_pixel_to_bin->get(), gpu_corrections->get(), image.getGPUBuffer(), gpu_mean, gpu_sigma,
             clip_k, accumulate_corrected, gpu_sum, gpu_sum2, gpu_count, gpu_sum_corr, gpu_sum2_corr,
             npix, nbins);
     } else {
         reduce_rings_global<<<reduce_blocks, reduce_threads, 0, *stream>>>(
-            gpu_pixel_to_bin, gpu_corrections, image.getGPUBuffer(), gpu_mean, gpu_sigma,
+            gpu_pixel_to_bin->get(), gpu_corrections->get(), image.getGPUBuffer(), gpu_mean, gpu_sigma,
             clip_k, accumulate_corrected, gpu_sum, gpu_sum2, gpu_count, gpu_sum_corr, gpu_sum2_corr,
             npix, nbins);
     }
@@ -327,7 +323,7 @@ void AdaptiveSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image,
     cuda_err(cudaMemcpyAsync(gpu_thr, host_thr.data(), sizeof(float) * nbins, cudaMemcpyHostToDevice, *stream));
     cuda_err(cudaMemsetAsync(gpu_strong, 0, OutputByteSize(), *stream));
     flag_strong<<<flag_blocks, flag_threads, 0, *stream>>>(
-        image.getGPUBuffer(), gpu_pixel_to_bin, gpu_thr, gpu_strong, npix, nbins);
+        image.getGPUBuffer(), gpu_pixel_to_bin->get(), gpu_thr, gpu_strong, npix, nbins);
     cuda_err(cudaMemcpyAsync(output_buffer.data(), gpu_strong, OutputByteSize(), cudaMemcpyDeviceToHost, *stream));
     cuda_err(cudaStreamSynchronize(*stream));
 }

image_analysis/spot_finding/AdaptiveSpotFinderGPU.h

@@ -30,6 +30,7 @@
 #include "../../common/AzimuthalIntegrationProfile.h"
 #include "../../common/AzimuthalIntegrationMapping.h"
 #include "../indexing/CUDAMemHelpers.h"
+#include "../indexing/CudaSharedTables.h"
 
 class AdaptiveSpotFinderGPU : public ImageSpotFinder {
     const AzimuthalIntegrationMapping &mapping;
@@ -46,9 +47,10 @@ class AdaptiveSpotFinderGPU : public ImageSpotFinder {
     size_t shared_clip = 0;    // per-block shared bytes for a sigma-clip pass (raw rings only)
     bool use_shared = true;    // false -> nbins too large for shared memory, use the global-atomics kernel
 
-    // Static mapping inputs (uploaded once).
-    CudaDevicePtr<uint16_t> gpu_pixel_to_bin;
-    CudaDevicePtr<float>    gpu_corrections;
+    // Static mapping inputs: geometry-only, so one copy per GPU shared with every other engine on it
+    // (see CudaSharedTables.h) rather than one copy per worker thread.
+    std::shared_ptr<CudaDevicePtr<uint16_t>> gpu_pixel_to_bin;
+    std::shared_ptr<CudaDevicePtr<float>>    gpu_corrections;
 
     // Raw per-ring accumulators (re-zeroed each pass) + derived stats used to clip and threshold.
     // double, like the CPU engine's ring accumulators: the ring sigma is the cancelling difference