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; };