From 3621d77093f585d9341d30cbdd60e8bfda7b8024 Mon Sep 17 00:00:00 2001 From: kferjaoui Date: Wed, 21 Oct 2026 16:33:59 +0200 Subject: [PATCH] PixelHistogramCUDA: take frames as InputType, keep a floating-point axis The third template parameter is now the frame type copied to the device, so raw uint16 frames can be histogrammed without a host-side conversion. The axis (range, edges, bin arithmetic) is AxisType: InputType for floating-point input, float for integers of at most 16 bits; wider integers are rejected at compile time. The kernel converts each value. Test: uint16 frames on 1.2-ADU bins match PixelHistogram fed the same values as float. --- include/aare/hist/PixelHistogramCUDA.hpp | 30 +++++++++++++++------ include/aare/hist/pixelhistogram_kernel.cuh | 9 ++++--- src/hist/PixelHistogramCUDA.test.cpp | 21 +++++++++++++++ src/hist/PixelHistogramCUDA_testimpl.cu | 29 ++++++++++---------- src/hist/PixelHistogramCUDA_testimpl.hpp | 9 ++++--- 5 files changed, 69 insertions(+), 29 deletions(-) diff --git a/include/aare/hist/PixelHistogramCUDA.hpp b/include/aare/hist/PixelHistogramCUDA.hpp index b332196b..ac7f9de8 100644 --- a/include/aare/hist/PixelHistogramCUDA.hpp +++ b/include/aare/hist/PixelHistogramCUDA.hpp @@ -21,7 +21,7 @@ namespace aare::cuda { template class PixelHistogramCUDA { #ifdef AARE_TESTS @@ -30,6 +30,20 @@ class PixelHistogramCUDA { FRIEND_TEST(known_pattern_transpose_and_cap) #endif + // Frames are copied to the device as InputType, as the caller provides + // them. The histogram axis (range, bin edges, the bin arithmetic) must be + // floating point: an integral input gets a float axis, converted in the + // kernel. That conversion is exact up to 16-bit integers. + static_assert(std::is_floating_point_v || + (std::is_integral_v && sizeof(InputType) <= 2), + "PixelHistogramCUDA: InputType must be floating point or an " + "integer of at most 16 bits"); + + public: + using AxisType = std::conditional_t, + InputType, float>; + + private: // nrows is declared first: its initializer validates all constructor // arguments, so it must run before anything below allocates. size_t nrows; @@ -51,7 +65,7 @@ class PixelHistogramCUDA { // Device buffers for frames and histogram DeviceBuffer d_histogram; - std::array, NUM_SLOTS> d_frames_batch; + std::array, NUM_SLOTS> d_frames_batch; // Events for synchronization std::array copy_done; @@ -120,7 +134,7 @@ class PixelHistogramCUDA { // Initialize device buffers for frames for (int i = 0; i < NUM_SLOTS; ++i) { d_frames_batch[i] = - DeviceBuffer(nrows * ncols * m_batch_size); + DeviceBuffer(nrows * ncols * m_batch_size); } } @@ -141,7 +155,7 @@ class PixelHistogramCUDA { m_registered_input.unregister(); } - void submit_batch(NDView frames) { + void submit_batch(NDView frames) { auto frame_shape = frames.shape(); if (static_cast(frame_shape[0]) > m_batch_size || static_cast(frame_shape[1]) != nrows || @@ -156,7 +170,7 @@ class PixelHistogramCUDA { const int slot = m_next_batch++ % NUM_SLOTS; // pointer to the first frame in the batch - const AxisType *h_frames_batch = frames.data(); + const InputType *h_frames_batch = frames.data(); // Make sure `copy_stream` waits on any kernel (using `slot`) to finish // execution @@ -166,7 +180,7 @@ class PixelHistogramCUDA { // Copy frames to device const size_t npix = nrows * ncols; CUDA_CHECK(cudaMemcpyAsync(d_frames_batch[slot].get(), h_frames_batch, - npix * batch_size * sizeof(AxisType), + npix * batch_size * sizeof(InputType), cudaMemcpyHostToDevice, copy_stream.get())); // Record the copy_done event (in the copy_stream) right after the copy @@ -191,11 +205,11 @@ class PixelHistogramCUDA { cudaEventRecord(kernel_done[slot].get(), compute_stream.get())); } - void fill_histogram(NDView frames) { + void fill_histogram(NDView frames) { auto N_frames = static_cast(frames.shape(0)); for (size_t i = 0; i < N_frames; i += m_batch_size) { size_t batch_end = std::min(i + m_batch_size, N_frames); - NDView batch = frames.sub_view(i, batch_end); + NDView batch = frames.sub_view(i, batch_end); submit_batch(batch); } } diff --git a/include/aare/hist/pixelhistogram_kernel.cuh b/include/aare/hist/pixelhistogram_kernel.cuh index f86d248b..bfeb5d2b 100644 --- a/include/aare/hist/pixelhistogram_kernel.cuh +++ b/include/aare/hist/pixelhistogram_kernel.cuh @@ -33,8 +33,10 @@ bin_to_pixel_major(const DeviceStorageType *__restrict__ d_histogram_binMajor, count < saturation_value ? count : saturation_value); } -template -__global__ void fill_bin_major(const AxisType *__restrict__ d_frames_batch, +// Frames arrive as InputType; each value is converted to the floating-point +// AxisType before binning. +template +__global__ void fill_bin_major(const InputType *__restrict__ d_frames_batch, DeviceStorageType *__restrict__ d_histogram, const size_t npix, const size_t nbins, AxisType xmin, AxisType xmax, AxisType scale, @@ -46,7 +48,8 @@ __global__ void fill_bin_major(const AxisType *__restrict__ d_frames_batch, for (size_t fid = 0; fid < batch_size; ++fid) { - AxisType val_pix_at_fid = d_frames_batch[fid * npix + pix_id]; + const AxisType val_pix_at_fid = + static_cast(d_frames_batch[fid * npix + pix_id]); if (!(val_pix_at_fid >= xmin && val_pix_at_fid < xmax)) // also skips NaN continue; diff --git a/src/hist/PixelHistogramCUDA.test.cpp b/src/hist/PixelHistogramCUDA.test.cpp index 1e53ad78..feb0a858 100644 --- a/src/hist/PixelHistogramCUDA.test.cpp +++ b/src/hist/PixelHistogramCUDA.test.cpp @@ -160,6 +160,27 @@ TEST_CASE("PixelHistogramCUDA result does not depend on the batch size", REQUIRE((cpu.values() == gpu.values)); } +TEST_CASE("PixelHistogramCUDA bins uint16 frames on a float axis", "[cuda]") { + // 1000 bins over 1200 ADU: bins 1.2 wide, so an integral axis (scale + // 1000 / 1200 == 0 in integers) would put everything in bin 0. + const size_t nrows = 7, ncols = 13, nframes = 100, nbins = 1000; + const float xmin = 3900.0f, xmax = 5100.0f; + auto frames = make_synth(nrows, ncols, nframes); // integer-valued floats + + NDArray raw(frames.shape()); + for (ssize_t i = 0; i < frames.size(); ++i) + raw.data()[i] = static_cast(frames.data()[i]); + + // The CPU reference gets the same values as float. + PixelHistogram cpu(nrows, ncols, nbins, xmin, xmax); + fill_cpu(cpu, frames.view()); + + auto gpu = test_cuda::run_gpu(raw.view(), nbins, xmin, xmax, + /*batch_size=*/8); + REQUIRE((cpu.bin_edges() == gpu.edges)); + REQUIRE((cpu.values() == gpu.values)); +} + TEST_CASE("PixelHistogramCUDA accumulates over several fills with values() " "in between", "[cuda]") { diff --git a/src/hist/PixelHistogramCUDA_testimpl.cu b/src/hist/PixelHistogramCUDA_testimpl.cu index 0194cb14..0b649a32 100644 --- a/src/hist/PixelHistogramCUDA_testimpl.cu +++ b/src/hist/PixelHistogramCUDA_testimpl.cu @@ -67,28 +67,29 @@ namespace aare::test_cuda { void check_known_pattern() { cuda::known_pattern_transpose_and_cap_impl(); } -template -DeviceResult run_gpu(NDView frames, size_t nbins, float xmin, +template +DeviceResult run_gpu(NDView frames, size_t nbins, float xmin, float xmax, size_t batch_size) { using D = std::conditional_t; - cuda::PixelHistogramCUDA hist(static_cast(frames.shape(1)), - static_cast(frames.shape(2)), - nbins, xmin, xmax, batch_size); + cuda::PixelHistogramCUDA hist( + static_cast(frames.shape(1)), + static_cast(frames.shape(2)), nbins, xmin, xmax, batch_size); hist.register_input_buffer(frames.data(), frames.total_bytes()); hist.fill_histogram(frames); hist.unregister_input_buffer(); return {hist.bin_centers(), hist.bin_edges(), hist.values()}; } -#define INSTANTIATE(H) \ - template DeviceResult run_gpu(NDView, size_t, float, \ - float, size_t); -INSTANTIATE(uint8_t) -INSTANTIATE(uint16_t) -INSTANTIATE(uint32_t) -INSTANTIATE(uint64_t) -INSTANTIATE(float) -INSTANTIATE(double) +#define INSTANTIATE(H, In) \ + template DeviceResult run_gpu(NDView, size_t, float, \ + float, size_t); +INSTANTIATE(uint8_t, float) +INSTANTIATE(uint16_t, float) +INSTANTIATE(uint32_t, float) +INSTANTIATE(uint64_t, float) +INSTANTIATE(float, float) +INSTANTIATE(double, float) +INSTANTIATE(uint32_t, uint16_t) #undef INSTANTIATE std::vector> diff --git a/src/hist/PixelHistogramCUDA_testimpl.hpp b/src/hist/PixelHistogramCUDA_testimpl.hpp index 6f29fb6c..5433c175 100644 --- a/src/hist/PixelHistogramCUDA_testimpl.hpp +++ b/src/hist/PixelHistogramCUDA_testimpl.hpp @@ -14,10 +14,11 @@ template struct DeviceResult { NDArray values; }; -// Histogram `frames` (nframes, nrows, ncols) on the GPU. Explicitly -// instantiated in the .cu for the storage types under test. -template -DeviceResult run_gpu(NDView frames, size_t nbins, float xmin, +// Histogram `frames` (nframes, nrows, ncols) on the GPU, with frames of type +// In. Explicitly instantiated in the .cu for the storage and input types +// under test; the axis is float for all of them. +template +DeviceResult run_gpu(NDView frames, size_t nbins, float xmin, float xmax, size_t batch_size); // Fill one histogram from consecutive parts of `frames`, part i ending at