diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index c66d2156..7d3e7e6d 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -767,6 +767,8 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) { h_var_bkg.resize(new_capacity); h_obs_x.resize(new_capacity); h_obs_y.resize(new_capacity); h_ok.resize(new_capacity); h_has_obs.resize(new_capacity); + + capacity = new_capacity; } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h index f8e71f60..1a86f9b6 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h @@ -60,6 +60,7 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { std::vector h_I, h_sigma, h_bkg, h_var_bkg, h_obs_x, h_obs_y; std::vector h_ok, h_has_obs; + void EnsureCapacity(size_t n); public: diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index c6180928..d92ae0b6 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -203,16 +203,24 @@ RotationScaleMerge::RotationScaleMerge(const DiffractionExperiment &experiment, void RotationScaleMerge::DivideOutIncidentFlux() { frame_flux.assign(n_frames, 1.0); - std::vector mean_bkg(n_frames, NAN), finite; - for (int f = 0; f < n_frames; ++f) { + // A frame at a time, not an observation at a time: each frame's mean is a running sum over its + // own contiguous range, so giving a thread whole frames keeps every sum in one thread and in the + // order it had. Splitting by observation would cut a frame across two threads and the sum would + // have to be recombined, which is a different sequence of roundings. + std::vector mean_bkg(n_frames, NAN); + ParallelFor(n_frames, nthreads, [&](int f) { double sum = 0.0; int n = 0; for (int i = frame_start[f]; i < frame_start[f] + frame_count[f]; ++i) if (std::isfinite(partials[i].bkg) && partials[i].bkg > 0.0f) { sum += partials[i].bkg; ++n; } - if (n == 0) continue; - mean_bkg[f] = sum / n; - finite.push_back(mean_bkg[f]); - } + if (n > 0) mean_bkg[f] = sum / n; + }); + // Collected in frame order, serially: it is one value per frame, and the median only depends on + // the set, but keeping the order removes the question. + std::vector finite; + finite.reserve(n_frames); + for (int f = 0; f < n_frames; ++f) + if (std::isfinite(mean_bkg[f])) finite.push_back(mean_bkg[f]); if (finite.empty()) return; @@ -224,8 +232,12 @@ void RotationScaleMerge::DivideOutIncidentFlux() { for (int f = 0; f < n_frames; ++f) if (std::isfinite(mean_bkg[f])) frame_flux[f] = mean_bkg[f] / median; - for (auto &o : partials) - o.rlp = static_cast(o.rlp / frame_flux[o.frame]); + ParallelChunks(static_cast(partials.size()), nthreads, [&](int lo, int hi) { + for (int i = lo; i < hi; ++i) { + auto &o = partials[i]; + o.rlp = static_cast(o.rlp / frame_flux[o.frame]); + } + }); } void RotationScaleMerge::Ingest() { diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu index 192f5e92..8e1d144a 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu +++ b/image_analysis/spot_finding/AdaptiveSpotFinderGPU.cu @@ -441,8 +441,8 @@ void AdaptiveSpotFinderGPU::Detect(const ImagePreprocessorBuffer &image, cuda_err(cudaMemsetAsync(gpu_strong, 0, OutputByteSize(), *stream)); flag_strong<<>>( image.getGPUBuffer(), gpu_pixel_to_bin->get(), gpu_thr, gpu_strong, npix, nbins); - // The bit buffer stays on the device - ExtractComponents reads it there. - cuda_err(cudaStreamSynchronize(*stream)); + // The bit buffer stays on the device and ExtractComponents reads it there, on this same stream, + // so the ordering already guarantees flag_strong has finished. Waiting here only idled the host. } void AdaptiveSpotFinderGPU::SetResolutionMask(const std::vector &mask) {