diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 63a55f58..8ad97ca5 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -50,8 +50,6 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp if (get_gpu_count() == 0) { #endif preprocessor_buffer = std::make_unique(experiment.GetPixelsNum()); - spotFinder = std::make_unique(experiment.GetXPixelsNum(), experiment.GetYPixelsNum()); - azint = std::make_unique(integration); preprocessor = std::make_unique(in_experiment, in_mask); bragg_engine = std::make_unique(in_experiment); if (experiment.ROI().size() >= 1) @@ -70,8 +68,6 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp // enable_fused_adaptive_gpu = true, so on the GPU path the copy is off in practice. preprocessor = std::make_unique(in_experiment, in_mask, stream, /*copy_image_to_host=*/!enable_fused_adaptive_gpu); - spotFinder = std::make_unique(experiment.GetXPixelsNum(), experiment.GetYPixelsNum(), stream); - azint = std::make_unique(integration, stream); bragg_engine = std::make_unique(in_experiment, stream); if (experiment.ROI().size() >= 1) roi = std::make_unique(experiment, stream); @@ -144,7 +140,7 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, if (!fused) { const auto azint_start_time = std::chrono::steady_clock::now(); - azint->Run(*preprocessor_buffer, profile); + AzInt().Run(*preprocessor_buffer, profile); const auto azint_end_time = std::chrono::steady_clock::now(); output.azint_time_s = std::chrono::duration(azint_end_time - azint_start_time).count(); } @@ -160,7 +156,7 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, ImageSpotFinder &finder = spot_finding_settings.adaptive_threshold ? static_cast(*adaptiveSpotFinder) - : *spotFinder; + : FixedThresholdFinder(); const auto integrate_fn = [this](const std::vector &predicted, size_t npredicted, int64_t image_number) { return bragg_engine->Run(*preprocessor_buffer, predicted, npredicted, image_number); @@ -293,6 +289,36 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, integration.Settings(), spot_finding_settings.ice_ring_width_Q_recipA); } +ImageSpotFinder &MXAnalysisWithoutFPGA::FixedThresholdFinder() { + if (!spotFinder) { +#ifdef JFJOCH_USE_CUDA + if (stream) + spotFinder = std::make_unique(experiment.GetXPixelsNum(), + experiment.GetYPixelsNum(), stream); + else +#endif + spotFinder = std::make_unique(experiment.GetXPixelsNum(), + 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); + } + return *spotFinder; +} + +AzIntEngine &MXAnalysisWithoutFPGA::AzInt() { + if (!azint) { +#ifdef JFJOCH_USE_CUDA + if (stream) + azint = std::make_unique(integration, stream); + else +#endif + azint = std::make_unique(integration); + } + return *azint; +} + void MXAnalysisWithoutFPGA::RebuildROI() { if (experiment.ROI().empty()) { roi.reset(); @@ -346,6 +372,7 @@ void MXAnalysisWithoutFPGA::UpdateMaskResolution(const SpotFindingSettings &sett // 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. - spotFinder->SetResolutionMask(mask_resolution); + if (spotFinder) + spotFinder->SetResolutionMask(mask_resolution); adaptiveSpotFinder->SetResolutionMask(mask_resolution); } diff --git a/image_analysis/MXAnalysisWithoutFPGA.h b/image_analysis/MXAnalysisWithoutFPGA.h index 58512091..7af8249b 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.h +++ b/image_analysis/MXAnalysisWithoutFPGA.h @@ -36,9 +36,16 @@ class MXAnalysisWithoutFPGA { size_t npixels; size_t xpixels; + // Built on first use: the fused adaptive finder produces the azimuthal profile as a by-product, + // so on the rugnux path this engine is constructed and then never run. std::unique_ptr azint; + AzIntEngine &AzInt(); std::unique_ptr roi; + // Built on first use. Which finder an image takes arrives with its SpotFindingSettings, and + // with adaptive detection on - the default everywhere but the broker - this one is never asked + // for; on the GPU it is ~14 MB and 15 device allocations per worker. std::unique_ptr spotFinder; + ImageSpotFinder &FixedThresholdFinder(); // Self-calibrating finder, used when spot settings request adaptive detection. Kept alongside the // default finder because the choice arrives with the per-image settings, not at construction. It is // an AdaptiveSpotFinderCPU by default; on the GPU path, when the fused engine is enabled (rugnux diff --git a/image_analysis/image_preprocessing/BSLZ4DecoderGPU.cu b/image_analysis/image_preprocessing/BSLZ4DecoderGPU.cu index 8068328d..29cd70d0 100644 --- a/image_analysis/image_preprocessing/BSLZ4DecoderGPU.cu +++ b/image_analysis/image_preprocessing/BSLZ4DecoderGPU.cu @@ -231,6 +231,18 @@ BSLZ4DecoderGPU::BSLZ4DecoderGPU(size_t in_max_uncompressed_bytes, std::shared_p // tens; sizing this from the UNCOMPRESSED size cost ~73 MB per worker to hold ~4 MB. } +// The stored pixel depth is fixed within a dataset, so in practice this runs once - but it is not +// promised anywhere, and sizing for the widest type instead would hold twice the memory a 16-bit +// detector needs. Grown with slack because cudaMalloc and cudaFree synchronise the whole device. +void BSLZ4DecoderGPU::EnsureUncompressedCapacity(size_t bytes) { + if (bytes <= max_uncompressed_bytes) + return; + const size_t want = std::max(bytes, max_uncompressed_bytes + max_uncompressed_bytes / 2); + cuda_err(cudaStreamSynchronize(*stream)); + gpu_shuffled = CudaDevicePtr(want); + max_uncompressed_bytes = want; +} + void BSLZ4DecoderGPU::EnsureCompressedCapacity(size_t bytes) { if (bytes <= compressed_capacity) return; @@ -259,8 +271,7 @@ BSLZ4ShuffledImage BSLZ4DecoderGPU::DecodeShuffled(const CompressedImage &image) if (clen < 12) throw JFJochException(JFJochExceptionCategory::Compression, "bslz4 chunk shorter than its header"); - if (total_bytes > max_uncompressed_bytes) - throw JFJochException(JFJochExceptionCategory::Compression, "bslz4 image larger than the decoder was sized for"); + EnsureUncompressedCapacity(total_bytes); if (be64(src) != total_bytes) throw JFJochException(JFJochExceptionCategory::Compression, "bslz4 header size does not match the image"); diff --git a/image_analysis/image_preprocessing/BSLZ4DecoderGPU.h b/image_analysis/image_preprocessing/BSLZ4DecoderGPU.h index bdc81485..98da7b5f 100644 --- a/image_analysis/image_preprocessing/BSLZ4DecoderGPU.h +++ b/image_analysis/image_preprocessing/BSLZ4DecoderGPU.h @@ -73,6 +73,7 @@ class BSLZ4DecoderGPU { size_t max_blocks = 0; void EnsureCompressedCapacity(size_t bytes); + void EnsureUncompressedCapacity(size_t bytes); void EnsureBlockCapacity(size_t nblocks); public: diff --git a/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu b/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu index 426393f0..37c772af 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu +++ b/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu @@ -295,7 +295,9 @@ bool ImagePreprocessorGPU::AnalyzeCompressed(ImagePreprocessorBuffer &processed_ return false; if (!bslz4_decoder) - bslz4_decoder = std::make_unique(npixels * sizeof(uint32_t), stream); + // Sized for the depth this image actually has, not for the widest one there could be: on + // 16-bit data the difference is half of a full frame per worker thread. + bslz4_decoder = std::make_unique(image.GetUncompressedSize(), stream); // LZ4 on the device, then ONE kernel that un-transposes the bitshuffle blocks and preprocesses // them as it goes. The decompressed image is never materialised: the fused kernel reads the diff --git a/image_analysis/indexing/CudaSharedTables.h b/image_analysis/indexing/CudaSharedTables.h index d7c30d95..3c19d5fb 100644 --- a/image_analysis/indexing/CudaSharedTables.h +++ b/image_analysis/indexing/CudaSharedTables.h @@ -42,14 +42,35 @@ namespace jfjoch_cuda_shared_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(data); - uint64_t h = 1469598103934665603ULL; - for (size_t i = 0; i < bytes; i++) { - h ^= p[i]; - h *= 1099511628211ULL; + 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; } - return h; + for (size_t i = n; i < bytes; i++) { + out ^= p[i]; + out *= PRIME; + } + return out; } inline Registry ®istry() {