diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index cc2489d2..394dbaed 100644 --- a/image_analysis/CMakeLists.txt +++ b/image_analysis/CMakeLists.txt @@ -30,9 +30,7 @@ ADD_LIBRARY(JFJochImageAnalysis STATIC RotationParameters.cpp RotationParameters.h WriteMmcif.cpp - WriteMmcif.h - azint/AzIntCPU.cpp - azint/AzIntCPU.h) + WriteMmcif.h) FIND_PACKAGE(Eigen3 3.4 REQUIRED NO_MODULE) # provides Eigen3::Eigen @@ -44,5 +42,6 @@ ADD_SUBDIRECTORY(geom_refinement) ADD_SUBDIRECTORY(lattice_search) ADD_SUBDIRECTORY(scale_merge) ADD_SUBDIRECTORY(image_preprocessing) +ADD_SUBDIRECTORY(azint) -TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochImagePreprocessing JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge gemmi) +TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochAzIntEngine JFJochImagePreprocessing JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge gemmi) diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 3623250e..5a319f82 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -11,6 +11,15 @@ #include "bragg_prediction/BraggPredictionFactory.h" #include "image_preprocessing/ImagePreprocessor.h" +#include "azint/AzIntEngineCPU.h" +#include "spot_finding/ImageSpotFinderCPU.h" +#ifdef JFJOCH_USE_CUDA +#include "azint/AzIntEngineGPU.h" +#include "spot_finding/ImageSpotFinderGPU.h" +#include "../common/CUDAWrapper.h" +#endif + + MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_experiment, const AzimuthalIntegration &in_integration, const PixelMask &in_mask, @@ -19,7 +28,6 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp integration(in_integration), npixels(experiment.GetPixelsNum()), xpixels(experiment.GetXPixelsNum()), - spotFinder(CreateImageSpotFinder(experiment.GetXPixelsNum(), experiment.GetYPixelsNum())), indexer(in_indexer), prediction(CreateBraggPrediction(experiment.IsRotationIndexing())), mask(in_mask), @@ -27,7 +35,18 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp mask_high_res(-1), mask_low_res(-1) { preprocessor = std::make_unique(in_experiment, in_mask); - azint = std::make_unique(integration); + +#ifdef JFJOCH_USE_CUDA + if (get_gpu_count() == 0) { +#endif + spotFinder = std::make_unique(experiment.GetXPixelsNum(), experiment.GetYPixelsNum()); + azint = std::make_unique(integration); +#ifdef JFJOCH_USE_CUDA + } else { + spotFinder = std::make_unique(experiment.GetXPixelsNum(), experiment.GetYPixelsNum()); + azint = std::make_unique(integration); + } +#endif } void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, diff --git a/image_analysis/MXAnalysisWithoutFPGA.h b/image_analysis/MXAnalysisWithoutFPGA.h index 33b449be..4f0f10ea 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.h +++ b/image_analysis/MXAnalysisWithoutFPGA.h @@ -14,7 +14,7 @@ #include "bragg_prediction/BraggPrediction.h" #include "spot_finding/ImageSpotFinder.h" #include "indexing/IndexerThreadPool.h" -#include "azint/AzIntCPU.h" +#include "azint/AzIntEngine.h" #include "IndexAndRefine.h" #include "image_preprocessing/ImagePreprocessor.h" @@ -30,7 +30,7 @@ class MXAnalysisWithoutFPGA { size_t npixels; size_t xpixels; - std::unique_ptr azint; + std::unique_ptr azint; std::unique_ptr spotFinder; IndexAndRefine &indexer; std::unique_ptr prediction; diff --git a/image_analysis/azint/AzIntCPU.h b/image_analysis/azint/AzIntCPU.h deleted file mode 100644 index 18bbfc3f..00000000 --- a/image_analysis/azint/AzIntCPU.h +++ /dev/null @@ -1,16 +0,0 @@ -// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#pragma once - -#include - -#include "../../common/AzimuthalIntegration.h" -#include "../../common/AzimuthalIntegrationProfile.h" - -class AzIntCPU { - const AzimuthalIntegration& integration; -public: - AzIntCPU(const AzimuthalIntegration& integration); - void Run(const std::vector &image, AzimuthalIntegrationProfile &profile) const; -}; \ No newline at end of file diff --git a/image_analysis/azint/AzIntEngine.cpp b/image_analysis/azint/AzIntEngine.cpp new file mode 100644 index 00000000..d5b28cea --- /dev/null +++ b/image_analysis/azint/AzIntEngine.cpp @@ -0,0 +1,11 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "AzIntEngine.h" + +AzIntEngine::AzIntEngine(const AzimuthalIntegration &integration) +: integration(integration), +npixel(integration.GetPixelToBin().size()), +azint_sum(integration.GetBinNumber(), 0.0f), +azint_count(integration.GetBinNumber(), 0.0f), +azint_bins(integration.GetBinNumber()){} diff --git a/image_analysis/azint/AzIntEngine.h b/image_analysis/azint/AzIntEngine.h new file mode 100644 index 00000000..97761aa3 --- /dev/null +++ b/image_analysis/azint/AzIntEngine.h @@ -0,0 +1,21 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include "../../common/AzimuthalIntegration.h" +#include "../../common/AzimuthalIntegrationProfile.h" + +class AzIntEngine { +protected: + const uint16_t azint_bins; + const size_t npixel; + const AzimuthalIntegration& integration; + std::vector azint_sum; + std::vector azint_sum2; + std::vector azint_count; +public: + AzIntEngine(const AzimuthalIntegration& integration); + virtual ~AzIntEngine() = default; + virtual void Run(const std::vector &image, AzimuthalIntegrationProfile &profile) = 0; +}; diff --git a/image_analysis/azint/AzIntCPU.cpp b/image_analysis/azint/AzIntEngineCPU.cpp similarity index 66% rename from image_analysis/azint/AzIntCPU.cpp rename to image_analysis/azint/AzIntEngineCPU.cpp index 77187e63..7c3ede74 100644 --- a/image_analysis/azint/AzIntCPU.cpp +++ b/image_analysis/azint/AzIntEngineCPU.cpp @@ -1,23 +1,20 @@ // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only -#include "AzIntCPU.h" +#include "AzIntEngineCPU.h" -AzIntCPU::AzIntCPU(const AzimuthalIntegration &integration) - : integration(integration) {} +AzIntEngineCPU::AzIntEngineCPU(const AzimuthalIntegration &integration) + : AzIntEngine(integration) {} -void AzIntCPU::Run(const std::vector &image, AzimuthalIntegrationProfile &profile) const { - const auto azint_bins = integration.GetBinNumber(); - - std::vector azint_sum(azint_bins); - std::vector azint_count(azint_bins); +void AzIntEngineCPU::Run(const std::vector &image, AzimuthalIntegrationProfile &profile){ for (int i = 0; i < azint_count.size(); i++) { azint_sum[i] = 0.0f; + azint_sum2[i] = 0.0f; azint_count[i] = 0; } - if (image.size() != integration.GetPixelToBin().size()) + if (image.size() != npixel) throw std::runtime_error("ImageSpotFinder::AzimIntegration: Mismatch in size"); const uint16_t *pixel_to_bin = integration.GetPixelToBin().data(); @@ -28,9 +25,11 @@ void AzIntCPU::Run(const std::vector &image, AzimuthalIntegrationProfil if (bin < azint_bins) { float val = static_cast(image[i]) * corrections[i]; azint_sum[bin] += val; + azint_sum2[bin] += val * val; ++azint_count[bin]; } } + profile.Clear(integration); profile.Add(azint_sum, azint_count); } diff --git a/image_analysis/azint/AzIntEngineCPU.h b/image_analysis/azint/AzIntEngineCPU.h new file mode 100644 index 00000000..f05ca5d5 --- /dev/null +++ b/image_analysis/azint/AzIntEngineCPU.h @@ -0,0 +1,12 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include "AzIntEngine.h" + +class AzIntEngineCPU : public AzIntEngine { +public: + AzIntEngineCPU(const AzimuthalIntegration& integration); + void Run(const std::vector &image, AzimuthalIntegrationProfile &profile) override; +}; \ No newline at end of file diff --git a/image_analysis/azint/AzIntEngineGPU.cu b/image_analysis/azint/AzIntEngineGPU.cu new file mode 100644 index 00000000..57653d1f --- /dev/null +++ b/image_analysis/azint/AzIntEngineGPU.cu @@ -0,0 +1,156 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "AzIntEngineGPU.h" + +inline void cuda_err(cudaError_t val) { + if (val != cudaSuccess) + throw JFJochException(JFJochExceptionCategory::GPUCUDAError, cudaGetErrorString(val)); +} + +__global__ +void gpu_azim_shared( + const uint16_t *__restrict__ pixel_to_bin, + const float *__restrict__ corrections, + const int32_t *__restrict__ input_buffer, + float *__restrict__ azint_sum, + float *__restrict__ azint_sum2, + uint32_t *__restrict__ azint_count, + size_t num_pixels, + int azint_bins) { + extern __shared__ float shared[]; + + float *s_sum = shared; + float *s_sum2 = &s_sum[azint_bins]; + uint32_t *s_count = (uint32_t *) &s_sum2[azint_bins]; + + // Initialize shared memory + for (int i = threadIdx.x; i < azint_bins; i += blockDim.x) { + s_sum[i] = 0.0f; + s_sum2[i] = 0.0f; + s_count[i] = 0; + } + + __syncthreads(); + + for (size_t idx = blockIdx.x * blockDim.x + threadIdx.x; + idx < num_pixels; + idx += blockDim.x * gridDim.x) { + uint16_t bin = pixel_to_bin[idx]; + + int32_t v = input_buffer[idx]; + bool valid = (v != INT32_MIN) & (v != INT32_MAX); + + if (bin < azint_bins && valid) { + const float val = static_cast(v) * corrections[idx]; + const float val2 = val * val; + atomicAdd(&s_sum[bin], val); + atomicAdd(&s_sum2[bin], val2); + atomicAdd(&s_count[bin], 1); + } + } + + __syncthreads(); + + // Merge to global memory + for (unsigned int i = threadIdx.x; i < azint_bins; i += blockDim.x) { + atomicAdd(&azint_sum[i], s_sum[i]); + atomicAdd(&azint_sum2[i], s_sum2[i]); + atomicAdd(&azint_count[i], s_count[i]); + } +} + +__global__ +void gpu_azim( + const uint16_t *__restrict__ pixel_to_bin, + const float *__restrict__ corrections, + const int32_t *__restrict__ input_buffer, + float *__restrict__ azint_sum, + float *__restrict__ azint_sum2, + uint32_t *__restrict__ azint_count, + size_t num_pixels, + int azint_bins) { + for (size_t idx = blockIdx.x * blockDim.x + threadIdx.x; + idx < num_pixels; + idx += blockDim.x * gridDim.x) { + uint16_t bin = pixel_to_bin[idx]; + + int32_t v = input_buffer[idx]; + bool valid = (v != INT32_MIN) & (v != INT32_MAX); + + if (bin < azint_bins && valid) { + const float val = static_cast(v) * corrections[idx]; + const float val2 = val * val; + atomicAdd(&azint_sum[bin], val); + atomicAdd(&azint_sum2[bin], val2); + atomicAdd(&azint_count[bin], 1); + } + } +} + +AzIntEngineGPU::AzIntEngineGPU(const AzimuthalIntegration &integration) + : AzIntEngine(integration), + gpu_azint_correction(npixel), + gpu_pixel_to_bin(npixel), + gpu_sum(azint_bins), + gpu_sum2(azint_bins), + gpu_count(azint_bins), + cpu_sum_reg(azint_sum), + cpu_sum2_reg(azint_sum2), + cpu_count_reg(azint_count), + gpu_image(npixel) { + + cudaDeviceProp prop; + cudaGetDeviceProperties(&prop, 0); + + threads = 128; + blocks = 4 * prop.multiProcessorCount; + shared_size = prop.sharedMemPerBlock; + shared_needed = azint_bins * sizeof(float) + azint_bins * sizeof(float) + azint_bins * sizeof(uint32_t); + + cudaMemcpy(gpu_azint_correction, integration.Corrections().data(), sizeof(float) * npixel, + cudaMemcpyHostToDevice); + cudaMemcpy(gpu_pixel_to_bin, integration.GetPixelToBin().data(), sizeof(uint16_t) * npixel, + cudaMemcpyHostToDevice); +} + +void AzIntEngineGPU::Run(const std::vector &image, AzimuthalIntegrationProfile &profile) { + if (image.size() != integration.GetPixelToBin().size()) + throw std::runtime_error("ImageSpotFinder::AzimIntegration: Mismatch in size"); + cuda_err(cudaMemsetAsync(gpu_sum, 0, sizeof(float) * azint_bins, stream)); + cuda_err(cudaMemsetAsync(gpu_sum2, 0, sizeof(float) * azint_bins, stream)); + cuda_err(cudaMemsetAsync(gpu_count, 0, sizeof(uint32_t) * azint_bins, stream)); + cudaMemcpyAsync(gpu_image, image.data(), sizeof(int32_t) * npixel, cudaMemcpyHostToDevice, stream); + + if (shared_needed < shared_size) { + gpu_azim_shared<<>>( + gpu_pixel_to_bin, + gpu_azint_correction, + gpu_image, + gpu_sum, + gpu_sum2, + gpu_count, + npixel, + azint_bins + ); + } else { + gpu_azim<<>>( + gpu_pixel_to_bin, + gpu_azint_correction, + gpu_image, + gpu_sum, + gpu_sum2, + gpu_count, + npixel, + azint_bins + ); + } + + cudaMemcpyAsync(azint_sum.data(), gpu_sum, sizeof(float) * azint_bins, cudaMemcpyDeviceToHost, stream); + cudaMemcpyAsync(azint_sum2.data(), gpu_sum2, sizeof(float) * azint_bins, cudaMemcpyDeviceToHost, stream); + cudaMemcpyAsync(azint_count.data(), gpu_count, sizeof(uint32_t) * azint_bins, cudaMemcpyDeviceToHost, stream); + cuda_err(cudaStreamSynchronize(stream)); + + profile.Clear(integration); + profile.Add(azint_sum, azint_count); +} diff --git a/image_analysis/azint/AzIntEngineGPU.h b/image_analysis/azint/AzIntEngineGPU.h new file mode 100644 index 00000000..e0778d79 --- /dev/null +++ b/image_analysis/azint/AzIntEngineGPU.h @@ -0,0 +1,30 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include "AzIntEngine.h" +#include "../indexing/CUDAMemHelpers.h" + +class AzIntEngineGPU : public AzIntEngine { + CudaStream stream; + int threads; + int blocks; + int shared_needed; + int shared_size; + + CudaDevicePtr gpu_azint_correction; + CudaDevicePtr gpu_pixel_to_bin; + + CudaDevicePtr gpu_sum; + CudaDevicePtr gpu_sum2; + CudaDevicePtr gpu_count; + CudaRegisteredVector cpu_sum_reg; + CudaRegisteredVector cpu_sum2_reg; + CudaRegisteredVector cpu_count_reg; + + CudaDevicePtr gpu_image; +public: + AzIntEngineGPU(const AzimuthalIntegration& integration); + void Run(const std::vector &image, AzimuthalIntegrationProfile &profile) override; +}; diff --git a/image_analysis/azint/CMakeLists.txt b/image_analysis/azint/CMakeLists.txt new file mode 100644 index 00000000..006e914a --- /dev/null +++ b/image_analysis/azint/CMakeLists.txt @@ -0,0 +1,5 @@ +ADD_LIBRARY(JFJochAzIntEngine STATIC AzIntEngine.cpp AzIntEngine.h AzIntEngineCPU.cpp AzIntEngineCPU.h) +TARGET_LINK_LIBRARIES(JFJochAzIntEngine JFJochCommon) +IF (JFJOCH_CUDA_AVAILABLE) + TARGET_SOURCES(JFJochAzIntEngine PRIVATE AzIntEngineGPU.cu AzIntEngineGPU.h) +ENDIF() \ No newline at end of file