From ac202a55a161b57db39aecb6c15744e4a32bfdda Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 28 Aug 2026 10:38:10 +0200 Subject: [PATCH] Spot finding: the strong-pixel limit follows the detector An image with 65535 or more strong pixels was given up on and reported ZERO spots - silently, no log line, indistinguishable from a frame that did not diffract. 65535 is one pixel in 64 of the JUNGFRAU 4M the number was written for; left fixed while the detectors grew it became one in 276 of an 18-megapixel EIGER, which a strongly diffracting crystal passes on its best frames. On the strong rotation set just added to the battery it cost 767 of 1800 images: peakCountUnfiltered 0 and resolutionEstimate NaN across two blocks of the sweep, the two where the crystal diffracts hardest. Make the bar one pixel in 64 everywhere, and never below the value that stood here, so no smaller detector loses ground. It lived in three places - the host extractor, StrongPixelSet, and SpotExtractorGPU's buffer capacity - now one function. The bar was there for a reason and raising it alone would not have been safe. sparseccl walks a sliding window of the last two lines and tests every pixel in it, which is quadratic in how many strong pixels a line pair holds: a handful for the silicon-tracker hits upstream wrote it for, four thousand for a lit detector line, and 76 seconds for a fully lit frame. But the pixels arrive in raster order, so the window need not be walked at all - a pixel's earlier 8-neighbours are the one to its left and the at most three above it, which is what the GPU extractor already finds by binary search. Keeping the previous line's range and a forward-only cursor gives the same edge set and the same unions in the same order, so the labels are identical, and the fully lit frame now takes 0.16 s. Verified bit-identical on real frames, on fully dense frames, across occupancy 1e-5 to 5e-2, and on 4000 randomised images including ones with blank lines; SpotExtractorGPU's host-vs-device parity test passes untouched. ImagePreprocessorBufferGPU's gather staging was sized to the old constant, with a comment tying it to the caller's give-up. Raising that give-up without it would have run the gather off the end of the device buffer, so it follows the same limit now. Byte-identical .hkl on three battery crystals that never reach the bar. On the strong set, with symmetry, cell and geometry pinned so only the spot list moves: better in every resolution shell, CC1/2 97.8 -> 98.5%, R_meas 30.5 -> 28.4%, ISa 3.36 -> 3.58, indexing rate 0.772 -> 0.824. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01FBumeJVx4oeXxiBRpkrE5H --- docs/CHANGELOG.md | 2 + .../ImagePreprocessorBufferGPU.cu | 5 +- .../ImagePreprocessorBufferGPU.h | 9 ++- .../spot_finding/ImageSpotFinder.cpp | 6 +- .../spot_finding/SpotExtractorGPU.cu | 29 ++++---- .../spot_finding/SpotExtractorGPU.h | 7 +- .../spot_finding/SpotFindingSettings.h | 16 +++++ .../spot_finding/StrongPixelSet.cpp | 66 +++++++++++-------- image_analysis/spot_finding/StrongPixelSet.h | 1 + 9 files changed, 89 insertions(+), 52 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 14505b792..a1dac33e0 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,6 +1,8 @@ # Changelog ## 1.0.0 ### 1.0.0-rc.166 +* Spot finding no longer throws away a whole image when it holds many strong pixels: the limit follows the detector (one pixel in 64) instead of standing at the 65535 that suited a 4-megapixel detector, which a strongly diffracting crystal on an 18-megapixel one passes on its best frames. +* The connected-component search is linear in the strong pixels rather than quadratic in how many of them a detector line holds; a fully lit image is labelled in 0.16 s instead of 76 s, and the spots it finds are unchanged. * In `jfjoch_viewer`, Alt and the mouse wheel step through the dataset one image at a time. ### 1.0.0-rc.165 diff --git a/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu b/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu index 54265f5c3..e45f97dd3 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu +++ b/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu @@ -16,8 +16,9 @@ ImagePreprocessorBufferGPU::ImagePreprocessorBufferGPU(size_t npixel, bool host_ gpu_image(npixel), // A no-op when the mirror was not allocated: CudaRegisteredVector skips an empty vector. buffer_reg(buffer), - gpu_gather_index(MAX_GATHER), - gpu_gather_value(MAX_GATHER) { + max_gather(StrongPixelLimit(npixel)), + gpu_gather_index(max_gather), + gpu_gather_value(max_gather) { } int32_t *ImagePreprocessorBufferGPU::getGPUBuffer() { diff --git a/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h b/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h index dd223be15..45fbb7bc6 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h +++ b/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h @@ -7,15 +7,18 @@ #include "ImagePreprocessorBuffer.h" #include "../indexing/CUDAMemHelpers.h" +#include "../spot_finding/SpotFindingSettings.h" class ImagePreprocessorBufferGPU : public ImagePreprocessorBuffer { CudaDevicePtr gpu_image; CudaRegisteredVector buffer_reg; // Staging for Gather(). Its only caller is ImageSpotFinder::ExtractSpots, which gives up on a frame - // with UINT16_MAX or more strong pixels (the connected-component search rejects it anyway), so that - // is the largest gather that can be asked for. - static constexpr size_t MAX_GATHER = UINT16_MAX; + // with StrongPixelLimit or more strong pixels (the connected-component search rejects it anyway), + // so that is the largest gather that can be asked for. It follows the detector, so this has to as + // well - sized to a constant while the caller's bar was raised, the gather would run off the end + // of the buffer rather than merely lose the frame. + const size_t max_gather; CudaDevicePtr gpu_gather_index; CudaDevicePtr gpu_gather_value; // Own stream: every analysis engine synchronises its own stream before it returns, so the device diff --git a/image_analysis/spot_finding/ImageSpotFinder.cpp b/image_analysis/spot_finding/ImageSpotFinder.cpp index 47a44d3f7..8fe18f98a 100644 --- a/image_analysis/spot_finding/ImageSpotFinder.cpp +++ b/image_analysis/spot_finding/ImageSpotFinder.cpp @@ -71,9 +71,9 @@ void ImageSpotFinder::ExtractComponentsHost(const ImagePreprocessorBuffer &image } components.clear(); - // The connected-component search rejects a frame with this many strong pixels, so their values are - // of no use. - if (strong_pixel.size() >= UINT16_MAX) + // The connected-component search gives up on a frame with this many strong pixels, so their values + // are of no use - not even worth gathering off the device. + if (strong_pixel.size() >= StrongPixelLimit(static_cast(width) * height)) return; image.Gather(strong_pixel, strong_pixel_value); diff --git a/image_analysis/spot_finding/SpotExtractorGPU.cu b/image_analysis/spot_finding/SpotExtractorGPU.cu index 2a5e9c38c..2795de1b9 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.cu +++ b/image_analysis/spot_finding/SpotExtractorGPU.cu @@ -185,7 +185,7 @@ __global__ void finish_components(const uint32_t *__restrict__ index, const int3 const int n = static_cast(min(*nstrong, capacity)); if (threadIdx.x == 0) *nout = 0; __syncthreads(); - // Same give-up as StrongPixelSet::FindComponentsImage - except that here the count is known + // The same give-up the host makes at StrongPixelLimit - except that here the count is known // before a single pixel has been written anywhere, so the frame costs nothing to reject. if (n == 0 || static_cast(n) >= capacity) return; @@ -275,16 +275,17 @@ SpotExtractorGPU::SpotExtractorGPU(int32_t in_width, int32_t in_height, std::sha : stream(std::move(in_stream)), width(in_width), nwords((static_cast(in_width) * in_height + 31) / 32), + max_strong(StrongPixelLimit(static_cast(in_width) * in_height)), gpu_res_mask(nwords), gpu_nstrong(1), - gpu_index(MAX_STRONG), - gpu_value(MAX_STRONG), - gpu_parent(MAX_STRONG), - gpu_root(MAX_STRONG), - gpu_label(MAX_STRONG), - gpu_count(MAX_STRONG), - gpu_spot(MAX_STRONG), - gpu_spot_out(MAX_STRONG), + gpu_index(max_strong), + gpu_value(max_strong), + gpu_parent(max_strong), + gpu_root(max_strong), + gpu_label(max_strong), + gpu_count(max_strong), + gpu_spot(max_strong), + gpu_spot_out(max_strong), gpu_nspot(1), host_nspot(1), host_spot(SPOT_PREFIX) { @@ -326,15 +327,15 @@ void SpotExtractorGPU::Extract(const uint32_t *gpu_strong, const int32_t *gpu_im scan_block_counts<<<1, FINISH_THREADS, 0, *stream>>>(gpu_block_count, gpu_block_offset, gpu_nstrong, compact_blocks); scatter_bits<<>>(gpu_strong, gpu_res_mask, gpu_block_offset, gpu_image, - gpu_index, gpu_value, nwords, MAX_STRONG); + gpu_index, gpu_value, nwords, max_strong); // Fixed grids reading the strong-pixel count from device memory: the host never learns it, so it // never has to synchronise in the middle of the frame. - init_parent<<<512, THREADS, 0, *stream>>>(gpu_parent, gpu_nstrong, MAX_STRONG); - union_neighbours<<<512, THREADS, 0, *stream>>>(gpu_index, gpu_parent, gpu_nstrong, MAX_STRONG, width); - resolve_roots<<<512, THREADS, 0, *stream>>>(gpu_parent, gpu_root, gpu_nstrong, MAX_STRONG); + init_parent<<<512, THREADS, 0, *stream>>>(gpu_parent, gpu_nstrong, max_strong); + union_neighbours<<<512, THREADS, 0, *stream>>>(gpu_index, gpu_parent, gpu_nstrong, max_strong, width); + resolve_roots<<<512, THREADS, 0, *stream>>>(gpu_parent, gpu_root, gpu_nstrong, max_strong); finish_components<<<1, FINISH_THREADS, 0, *stream>>>(gpu_index, gpu_value, gpu_root, gpu_label, gpu_count, gpu_spot, gpu_spot_out, gpu_nspot, - gpu_nstrong, MAX_STRONG, width, max_pix); + gpu_nstrong, max_strong, width, max_pix); cuda_err(cudaMemcpyAsync(host_nspot, gpu_nspot, sizeof(uint32_t), cudaMemcpyDeviceToHost, *stream)); cuda_err(cudaMemcpyAsync(host_spot, gpu_spot_out, SPOT_PREFIX * sizeof(SpotExtractorGPUSpot), cudaMemcpyDeviceToHost, *stream)); diff --git a/image_analysis/spot_finding/SpotExtractorGPU.h b/image_analysis/spot_finding/SpotExtractorGPU.h index e00a86cc6..4cbc93735 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.h +++ b/image_analysis/spot_finding/SpotExtractorGPU.h @@ -51,9 +51,10 @@ class SpotExtractorGPU { const int32_t width; const size_t nwords; - // The connected-component search gives up above this many strong pixels (see - // StrongPixelSet::FindComponentsImage), so nothing larger is ever built. - static constexpr uint32_t MAX_STRONG = UINT16_MAX; + // Strong pixels this engine's buffers hold, and above which the extraction gives up on the frame - + // StrongPixelLimit, so it follows the detector rather than standing at a constant. 104 bytes of + // device memory apiece, 29 MB on an 18-megapixel detector. + const uint32_t max_strong; // Spots copied back together with their count in one transfer. A frame with more than this many // surviving spots - far past anything indexable - simply takes a second copy. static constexpr uint32_t SPOT_PREFIX = 4096; diff --git a/image_analysis/spot_finding/SpotFindingSettings.h b/image_analysis/spot_finding/SpotFindingSettings.h index 5ecc300d7..a2bdde060 100644 --- a/image_analysis/spot_finding/SpotFindingSettings.h +++ b/image_analysis/spot_finding/SpotFindingSettings.h @@ -3,9 +3,25 @@ #pragma once +#include #include #include +// Strong pixels above which the connected-component search gives up on a frame, unlabelled: an image +// with this much of the detector over threshold is not a diffraction pattern, and it is what the +// spot extractor's buffers are sized to. It is no longer a time limit - the search is linear in the +// strong pixels either way, and a fully lit 18-megapixel frame labels in 0.16 s. +// +// The bar has to be a FRACTION of the detector. It stood at a fixed 65535, which is one pixel in 64 +// of the JUNGFRAU 4M it was written for; left fixed while the detectors grew it became one pixel in +// 276 of an 18-megapixel EIGER - a bar a strongly diffracting crystal clears on its best frames, +// which were then dropped whole, and in silence. One in 64 everywhere, and never below the value +// that used to stand here, so no smaller detector loses ground. +constexpr uint32_t StrongPixelLimit(size_t pixel_count) { + const auto limit = static_cast(pixel_count / 64); + return limit > UINT16_MAX ? limit : UINT16_MAX; +} + struct SpotFindingSettings { bool enable = true; float signal_to_noise_threshold = 4.0; // STRONG_PIXEL in XDS diff --git a/image_analysis/spot_finding/StrongPixelSet.cpp b/image_analysis/spot_finding/StrongPixelSet.cpp index 7af8acfba..6c631236a 100644 --- a/image_analysis/spot_finding/StrongPixelSet.cpp +++ b/image_analysis/spot_finding/StrongPixelSet.cpp @@ -4,6 +4,9 @@ // SparseCCL code taken from https://github.com/acts-project/traccc/blob/main/core/include/traccc/clusterization/detail/sparse_ccl.hpp // (c) 2021-2022 CERN for the benefit of the ACTS project // Mozilla Public License Version 2.0 +// +// The union-find and the two-scan structure are theirs. How a pixel's earlier neighbours are FOUND +// is not: see sparseccl below. #include @@ -18,16 +21,6 @@ void StrongPixelSet::AddStrongPixel(uint16_t col, uint16_t line, int32_t photons ++strong_pixel_count; } -bool is_far_enough(strong_pixel pixel0, strong_pixel pixel1) { - return (pixel1.line - pixel0.line) > 1; -} - -bool is_adjacent(strong_pixel pixel0, strong_pixel pixel1) { - auto line_diff = pixel0.line - pixel1.line; - auto col_diff = pixel0.col - pixel1.col; - return line_diff <= 1 && line_diff >= -1 && col_diff <= 1 && col_diff >= -1; -} - uint32_t StrongPixelSet::find_root(uint32_t e) { uint32_t r = e; while (L[r] != r) @@ -52,18 +45,38 @@ std::vector StrongPixelSet::sparseccl() { unsigned int labels = 0; - // first scan: pixel association - uint32_t start_j = 0; + // First scan: pixel association. The pixels arrive in raster order - line ascending, column + // ascending within a line - which upstream uses to walk a sliding window of the last two lines, + // testing every pixel in it for adjacency. That is quadratic in how many strong pixels a line + // pair holds: fine for the silicon-tracker hits it was written for, but a flooded detector line + // holds four thousand of them, and labelling a fully lit frame took 76 seconds. + // + // Since the columns ascend, the window need not be walked. A pixel's earlier 8-neighbours are + // exactly the one to its left and the at most three above it, so keep the previous line's range + // and a cursor into it that only ever moves forward - the same four neighbours the GPU extractor + // finds by binary search. Same edge set, same unions in the same order, therefore the same + // labels; the flooded frame now takes 0.16 s. + uint32_t line_begin = 0; // first pixel of the line being scanned + uint32_t prev_begin = 0, prev_end = 0; // the pixels of the line above it + uint32_t up = 0; // cursor into [prev_begin, prev_end) for (uint32_t i = 0; i < pixels.size(); ++i) { L[i] = i; - uint32_t ai = i; - for (uint32_t j = start_j; j < i; ++j) { - if (is_adjacent(pixels[i], pixels[j])) { - ai = make_union(ai, find_root(j)); - } else if (is_far_enough(pixels[j], pixels[i])) { - ++start_j; - } + if (i > 0 && pixels[i].line != pixels[i - 1].line) { + // The line above is the previous one only if it really is the line above: a line with no + // strong pixel at all leaves nothing to join to. + prev_begin = (pixels[i].line == pixels[i - 1].line + 1) ? line_begin : i; + prev_end = i; + line_begin = i; + up = prev_begin; } + uint32_t ai = i; + while (up < prev_end && pixels[up].col + 1 < pixels[i].col) + ++up; + for (uint32_t j = up; j < prev_end && pixels[j].col <= pixels[i].col + 1; ++j) + ai = make_union(ai, find_root(j)); + // The pixel to the left comes last, as it did when the window was walked in order. + if (i > line_begin && pixels[i - 1].col + 1 == pixels[i].col) + ai = make_union(ai, find_root(i - 1)); } // second scan: transitive closure @@ -85,19 +98,18 @@ std::vector StrongPixelSet::sparseccl() { void StrongPixelSet::FindComponentsImage(const SpotFindingSettings &settings, std::vector &spots) { - // Avoid spot finding, when more than 65536 strong pixel count (will be super slow) - if (!pixels.empty() && (strong_pixel_count < UINT16_MAX)) { - for (const auto &spot: sparseccl()) { - if (spot.PixelCount() <= settings.max_pix_per_spot) - spots.push_back(spot); - } + // No StrongPixelLimit test here: the caller knows how big the image is and has already applied it. + for (const auto &spot: sparseccl()) { + if (spot.PixelCount() <= settings.max_pix_per_spot) + spots.push_back(spot); } } void StrongPixelSet::FindSpots(const DiffractionExperiment &experiment, const SpotFindingSettings &settings, std::vector &spots, uint16_t module_number) { - // Avoid spot finding, when more than 65536 strong pixel count (will be super slow) - if (!pixels.empty() && (strong_pixel_count < UINT16_MAX)) { + // Per module, so the bar is the module's own - and ReadFPGAOutput has already refused anything + // past max_strong_pixel_per_module, far below it. + if (!pixels.empty() && (strong_pixel_count < StrongPixelLimit(RAW_MODULE_SIZE))) { for (const auto &spot: sparseccl()) { if ((spot.PixelCount() <= settings.max_pix_per_spot) && (spot.PixelCount() >= settings.min_pix_per_spot.value_or(2))) { diff --git a/image_analysis/spot_finding/StrongPixelSet.h b/image_analysis/spot_finding/StrongPixelSet.h index 2e53c28c2..5026122d3 100644 --- a/image_analysis/spot_finding/StrongPixelSet.h +++ b/image_analysis/spot_finding/StrongPixelSet.h @@ -38,6 +38,7 @@ public: // Every connected component of at most max-pix pixels. min-pix is deliberately NOT applied: it is // the only spot setting that changes between the passes of the per-image min-pix search, so one // connected-component search serves all of them (see ImageSpotFinder::ExtractComponents). + // The caller is responsible for not handing over more than StrongPixelLimit pixels. void FindComponentsImage(const SpotFindingSettings &settings, std::vector &spots); uint32_t GetStrongPixelCount() const; };