diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index ed9dca95..cc2489d2 100644 --- a/image_analysis/CMakeLists.txt +++ b/image_analysis/CMakeLists.txt @@ -30,7 +30,9 @@ ADD_LIBRARY(JFJochImageAnalysis STATIC RotationParameters.cpp RotationParameters.h WriteMmcif.cpp - WriteMmcif.h) + WriteMmcif.h + azint/AzIntCPU.cpp + azint/AzIntCPU.h) FIND_PACKAGE(Eigen3 3.4 REQUIRED NO_MODULE) # provides Eigen3::Eigen diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 922ce31b..3623250e 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -26,7 +26,8 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp mask_resolution(experiment.GetPixelsNum(), false), mask_high_res(-1), mask_low_res(-1) { - preprocessor = std::make_unique(in_experiment, in_mask, spotFinder->GetInputBuffer()); + preprocessor = std::make_unique(in_experiment, in_mask); + azint = std::make_unique(integration); } void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, @@ -44,12 +45,12 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, output.compression_time_s = std::chrono::duration(compression_end_time - compression_start_time).count(); const auto preprocessing_start_time = std::chrono::steady_clock::now(); - auto ret = preprocessor->Analyze(image_ptr, output.image.GetMode()); + auto ret = preprocessor->Analyze(spotFinder->GetInputBuffer(), image_ptr, output.image.GetMode()); const auto preprocessing_end_time = std::chrono::steady_clock::now(); output.preprocessing_time_s = std::chrono::duration(preprocessing_end_time - preprocessing_start_time).count(); const auto azint_start_time = std::chrono::steady_clock::now(); - spotFinder->AzimIntegration(integration, profile); + azint->Run(spotFinder->GetInputBuffer(), 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(); diff --git a/image_analysis/MXAnalysisWithoutFPGA.h b/image_analysis/MXAnalysisWithoutFPGA.h index 21d4d06a..33b449be 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.h +++ b/image_analysis/MXAnalysisWithoutFPGA.h @@ -14,6 +14,7 @@ #include "bragg_prediction/BraggPrediction.h" #include "spot_finding/ImageSpotFinder.h" #include "indexing/IndexerThreadPool.h" +#include "azint/AzIntCPU.h" #include "IndexAndRefine.h" #include "image_preprocessing/ImagePreprocessor.h" @@ -29,6 +30,7 @@ class MXAnalysisWithoutFPGA { size_t npixels; size_t xpixels; + std::unique_ptr azint; std::unique_ptr spotFinder; IndexAndRefine &indexer; std::unique_ptr prediction; diff --git a/image_analysis/azint/AzIntCPU.cpp b/image_analysis/azint/AzIntCPU.cpp new file mode 100644 index 00000000..77187e63 --- /dev/null +++ b/image_analysis/azint/AzIntCPU.cpp @@ -0,0 +1,36 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "AzIntCPU.h" + +AzIntCPU::AzIntCPU(const AzimuthalIntegration &integration) + : integration(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); + + for (int i = 0; i < azint_count.size(); i++) { + azint_sum[i] = 0.0f; + azint_count[i] = 0; + } + + if (image.size() != integration.GetPixelToBin().size()) + throw std::runtime_error("ImageSpotFinder::AzimIntegration: Mismatch in size"); + + const uint16_t *pixel_to_bin = integration.GetPixelToBin().data(); + const float *corrections = integration.Corrections().data(); + + for (int i = 0; i < image.size(); i++) { + const uint16_t bin = pixel_to_bin[i]; + if (bin < azint_bins) { + float val = static_cast(image[i]) * corrections[i]; + azint_sum[bin] += val; + ++azint_count[bin]; + } + } + profile.Clear(integration); + profile.Add(azint_sum, azint_count); +} diff --git a/image_analysis/azint/AzIntCPU.h b/image_analysis/azint/AzIntCPU.h new file mode 100644 index 00000000..18bbfc3f --- /dev/null +++ b/image_analysis/azint/AzIntCPU.h @@ -0,0 +1,16 @@ +// 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/image_preprocessing/ImagePreprocessor.cpp b/image_analysis/image_preprocessing/ImagePreprocessor.cpp index 521b2c06..85273698 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessor.cpp +++ b/image_analysis/image_preprocessing/ImagePreprocessor.cpp @@ -5,42 +5,41 @@ #include "ImagePreprocessor.h" ImagePreprocessor::ImagePreprocessor(const DiffractionExperiment &experiment, - const PixelMask &mask, - std::vector &processed_image) + const PixelMask &mask) : npixels(experiment.GetPixelsNum()), experiment(experiment), - processed_image(processed_image), mask_1bit(npixels, false), saturation_limit(experiment.GetSaturationLimit()) { - if (processed_image.size() != npixels) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Processed image size mismatch"); - for (int i = 0; i < npixels; i++) mask_1bit[i] = (mask.GetMask().at(i) != 0); } -ImageStatistics ImagePreprocessor::Analyze(const uint8_t *image_ptr, CompressedImageMode image_mode) { +ImageStatistics ImagePreprocessor::Analyze(std::vector &processed_image, const uint8_t *image_ptr, CompressedImageMode image_mode) { switch (image_mode) { case CompressedImageMode::Int8: - return Analyze(image_ptr, INT8_MIN, INT8_MAX); + return Analyze(processed_image, image_ptr, INT8_MIN, INT8_MAX); case CompressedImageMode::Int16: - return Analyze(image_ptr, INT16_MIN, INT16_MAX); + return Analyze(processed_image, image_ptr, INT16_MIN, INT16_MAX); case CompressedImageMode::Int32: - return Analyze(image_ptr, INT32_MIN, INT32_MAX); + return Analyze(processed_image, image_ptr, INT32_MIN, INT32_MAX); case CompressedImageMode::Uint8: - return Analyze(image_ptr, UINT8_MAX, UINT8_MAX); + return Analyze(processed_image, image_ptr, UINT8_MAX, UINT8_MAX); case CompressedImageMode::Uint16: - return Analyze(image_ptr, UINT16_MAX, UINT16_MAX); + return Analyze(processed_image, image_ptr, UINT16_MAX, UINT16_MAX); case CompressedImageMode::Uint32: - return Analyze(image_ptr, UINT32_MAX, UINT32_MAX); + return Analyze(processed_image, image_ptr, UINT32_MAX, UINT32_MAX); default: throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "RGB/float mode not supported"); } } template -ImageStatistics ImagePreprocessor::Analyze(const uint8_t *input, T err_pixel_val, T sat_pixel_val) { +ImageStatistics ImagePreprocessor::Analyze(std::vector &processed_image, const uint8_t *input, T err_pixel_val, T sat_pixel_val) { + + if (processed_image.size() != npixels) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Processed image size mismatch"); + auto image = reinterpret_cast(input); ImageStatistics ret{}; diff --git a/image_analysis/image_preprocessing/ImagePreprocessor.h b/image_analysis/image_preprocessing/ImagePreprocessor.h index 5ad3bd5d..8f325673 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessor.h +++ b/image_analysis/image_preprocessing/ImagePreprocessor.h @@ -21,14 +21,11 @@ class ImagePreprocessor { protected: const size_t npixels; const DiffractionExperiment &experiment; - std::vector &processed_image; std::vector mask_1bit; const int64_t saturation_limit; - template ImageStatistics Analyze(const uint8_t *input, T err_value, T sat_value); + template ImageStatistics Analyze(std::vector &processed_image, const uint8_t *input, T err_value, T sat_value); public: - ImagePreprocessor(const DiffractionExperiment &experiment, - const PixelMask &mask, - std::vector &processed_image); + ImagePreprocessor(const DiffractionExperiment &experiment, const PixelMask &mask); - ImageStatistics Analyze(const uint8_t *decompressed_image, CompressedImageMode image_mode); + ImageStatistics Analyze(std::vector &processed_image, const uint8_t *decompressed_image, CompressedImageMode image_mode); }; diff --git a/image_analysis/spot_finding/ImageSpotFinder.cpp b/image_analysis/spot_finding/ImageSpotFinder.cpp index fd561ec9..92719aac 100644 --- a/image_analysis/spot_finding/ImageSpotFinder.cpp +++ b/image_analysis/spot_finding/ImageSpotFinder.cpp @@ -43,35 +43,3 @@ std::vector ImageSpotFinder::ExtractSpots(const SpotFindingSett pixel_set.FindSpotsImage(settings, vec); return vec; } - -void ImageSpotFinder::AzimIntegration(const AzimuthalIntegration &integration, AzimuthalIntegrationProfile &profile) { - const auto azint_bins = integration.GetBinNumber(); - std::vector azint_sum(azint_bins); - std::vector azint_sum2(azint_bins); - std::vector azint_count(azint_bins); - - for (int i = 0; i < azint_count.size(); i++) { - azint_sum[i] = 0.0f; - azint_sum2[i] = 0.0f; - azint_count[i] = 0; - } - - auto &pixel_to_bin = integration.GetPixelToBin(); - auto &corrections = integration.Corrections(); - - if (pixel_to_bin.size() != width * height || corrections.size() != width * height) - throw std::runtime_error("ImageSpotFinder::AzimIntegration: Mismatch in size"); - - for (int i = 0; i < width * height; i++) { - const uint16_t bin = pixel_to_bin[i]; - if (bin < azint_bins) { - float val = static_cast(input_buffer[i]) * corrections[i]; - azint_sum[bin] += val; - ++azint_count[bin]; - } - } - profile.Clear(integration); - profile.Add(azint_sum, azint_count); -} - - diff --git a/image_analysis/spot_finding/ImageSpotFinder.h b/image_analysis/spot_finding/ImageSpotFinder.h index bd6b6833..3c0d3005 100644 --- a/image_analysis/spot_finding/ImageSpotFinder.h +++ b/image_analysis/spot_finding/ImageSpotFinder.h @@ -30,8 +30,6 @@ public: virtual ~ImageSpotFinder() = default; virtual std::vector Run(const SpotFindingSettings &settings, const std::vector &res_mask) = 0; - - void AzimIntegration(const AzimuthalIntegration &integration, AzimuthalIntegrationProfile &profile); };