From 59281f6330fbed591052960b16e34dbabe05a398 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Wed, 10 Jun 2026 16:36:37 +0200 Subject: [PATCH] Revert "Bragg integration: option to use azimuthal integration profile" This reverts commit e4230bc14ee9a0b5a127d24d1cb71ca486820de4. --- common/BraggIntegrationSettings.cpp | 9 - common/BraggIntegrationSettings.h | 5 +- image_analysis/IndexAndRefine.cpp | 14 +- image_analysis/IndexAndRefine.h | 8 +- image_analysis/MXAnalysisAfterFPGA.cpp | 11 +- image_analysis/MXAnalysisAfterFPGA.h | 5 +- image_analysis/MXAnalysisWithoutFPGA.cpp | 2 +- .../bragg_integration/BraggIntegrate2D.cpp | 219 +++++------------- .../bragg_integration/BraggIntegrate2D.h | 10 +- receiver/JFJochReceiver.cpp | 25 +- receiver/JFJochReceiver.h | 2 +- receiver/JFJochReceiverFPGA.cpp | 11 +- receiver/JFJochReceiverLite.cpp | 4 +- tests/BraggIntegrate2DTest.cpp | 6 +- tools/jfjoch_scale.cpp | 2 + 15 files changed, 115 insertions(+), 218 deletions(-) diff --git a/common/BraggIntegrationSettings.cpp b/common/BraggIntegrationSettings.cpp index 0528084e..3fe0a861 100644 --- a/common/BraggIntegrationSettings.cpp +++ b/common/BraggIntegrationSettings.cpp @@ -87,12 +87,3 @@ float BraggIntegrationSettings::GetDMinLimit_A() const { float BraggIntegrationSettings::GetMinimumSigmaInRegardsToI() const { return minimum_sigma_in_regards_to_i; } - -bool BraggIntegrationSettings::IsUseAzimProfile() const { - return use_azim_profile; -} - -BraggIntegrationSettings &BraggIntegrationSettings::UseAzimProfile(bool input) { - use_azim_profile = input; - return *this; -} diff --git a/common/BraggIntegrationSettings.h b/common/BraggIntegrationSettings.h index 2c15921d..60ab4059 100644 --- a/common/BraggIntegrationSettings.h +++ b/common/BraggIntegrationSettings.h @@ -12,14 +12,14 @@ class BraggIntegrationSettings { float d_min_limit_A = 1.0; std::optional fixed_profile_radius; float minimum_sigma_in_regards_to_i = 0.02; - bool use_azim_profile = false; + public: BraggIntegrationSettings& R1(float input); BraggIntegrationSettings& R2(float input); BraggIntegrationSettings& R3(float input); BraggIntegrationSettings& DMinLimit_A(float input); BraggIntegrationSettings& FixedProfileRadius_recipA(std::optional input); - BraggIntegrationSettings& UseAzimProfile(bool input); + [[nodiscard]] float GetR1() const; [[nodiscard]] float GetR2() const; @@ -28,5 +28,4 @@ public: [[nodiscard]] float GetDMinLimit_A() const; [[nodiscard]] float GetMinimumSigmaInRegardsToI() const; - [[nodiscard]] bool IsUseAzimProfile() const; }; diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 459af31e..98b1aed6 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -239,8 +239,8 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, const CompressedImage &image, BraggPrediction &prediction, const IndexAndRefine::IndexingOutcome &outcome, - const AzimuthalIntegrationMapping &mapping, - const AzimuthalIntegrationProfile &profile) { + const AzimuthalIntegrationMapping *mapping, + const AzimuthalIntegrationProfile *profile) { if (!outcome.lattice_candidate) return; @@ -299,11 +299,11 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, // integrated reflections that flow into the normal save/merge). const bool use_pixel_refine = experiment.GetIndexingSettings().GetGeomRefinementAlgorithm() == GeomRefinementAlgorithmEnum::PixelRefine - && !pixel_reference_.empty(); + && !pixel_reference_.empty() && mapping && profile; if (use_pixel_refine) { auto integration_start_time = std::chrono::steady_clock::now(); - PixelRefineIntegrate(msg, image, prediction, outcome, mapping, profile, i_outcome); + PixelRefineIntegrate(msg, image, prediction, outcome, *mapping, *profile, i_outcome); msg.integrated_reflections = i_outcome.reflections.size(); auto integration_end_time = std::chrono::steady_clock::now(); msg.integration_time_s = std::chrono::duration(integration_end_time - integration_start_time).count(); @@ -314,7 +314,7 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, msg.bragg_prediction_time_s = std::chrono::duration(pred_end_time - pred_start_time).count(); auto integration_start_time = std::chrono::steady_clock::now(); - i_outcome.reflections = BraggIntegrate2D(outcome.experiment, image, prediction.GetReflections(), mapping, profile, nrefl, msg.number); + i_outcome.reflections = BraggIntegrate2D(outcome.experiment, image, prediction.GetReflections(), nrefl, msg.number); msg.integrated_reflections = i_outcome.reflections.size(); auto integration_end_time = std::chrono::steady_clock::now(); msg.integration_time_s = std::chrono::duration(integration_end_time - integration_start_time).count(); @@ -357,8 +357,8 @@ void IndexAndRefine::ProcessImage(DataMessage &msg, const SpotFindingSettings &spot_finding_settings, const CompressedImage &image, BraggPrediction &prediction, - const AzimuthalIntegrationMapping &mapping, - const AzimuthalIntegrationProfile &profile) { + const AzimuthalIntegrationMapping *mapping, + const AzimuthalIntegrationProfile *profile) { if (!indexer_ || !spot_finding_settings.indexing) return; diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index d73fc2f1..7674392e 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -66,8 +66,8 @@ class IndexAndRefine { const CompressedImage &image, BraggPrediction &prediction, const IndexingOutcome &outcome, - const AzimuthalIntegrationMapping &mapping, - const AzimuthalIntegrationProfile &profile); + const AzimuthalIntegrationMapping *mapping, + const AzimuthalIntegrationProfile *profile); std::unique_ptr scaling_engine; void ScaleImage(DataMessage &msg, IntegrationOutcome& outcome); @@ -93,8 +93,8 @@ public: void ForceRotationIndexerLattice(const CrystalLattice& lattice); void ProcessImage(DataMessage &msg, const SpotFindingSettings &settings, const CompressedImage &image, BraggPrediction &prediction, - const AzimuthalIntegrationMapping &mapping, - const AzimuthalIntegrationProfile &profile); + const AzimuthalIntegrationMapping *mapping = nullptr, + const AzimuthalIntegrationProfile *profile = nullptr); IndexAndRefine& ReferenceIntensities(std::vector &reference); ScalingResult ScaleAllImages(const std::vector &reference, size_t nthreads = 0); diff --git a/image_analysis/MXAnalysisAfterFPGA.cpp b/image_analysis/MXAnalysisAfterFPGA.cpp index 96d8430d..a5347236 100644 --- a/image_analysis/MXAnalysisAfterFPGA.cpp +++ b/image_analysis/MXAnalysisAfterFPGA.cpp @@ -25,12 +25,10 @@ double stddev(const std::vector &v) { } -MXAnalysisAfterFPGA::MXAnalysisAfterFPGA(const DiffractionExperiment &in_experiment, IndexAndRefine &indexer, - AzimuthalIntegrationMapping &mapping) +MXAnalysisAfterFPGA::MXAnalysisAfterFPGA(const DiffractionExperiment &in_experiment, IndexAndRefine &indexer) : experiment(in_experiment), indexer(indexer), - prediction(CreateBraggPrediction(experiment.IsRotationIndexing())), - mapping(mapping) { + prediction(CreateBraggPrediction(experiment.IsRotationIndexing())) { if (experiment.IsSpotFindingEnabled()) find_spots = true; } @@ -105,8 +103,7 @@ void MXAnalysisAfterFPGA::ReadFromCPU(DeviceOutput *output, const SpotFindingSet } } -void MXAnalysisAfterFPGA::Process(DataMessage &message, const SpotFindingSettings &spot_finding_settings, - const AzimuthalIntegrationProfile &profile) { +void MXAnalysisAfterFPGA::Process(DataMessage &message, const SpotFindingSettings& spot_finding_settings) { if (find_spots && (state == State::Enabled)) { const auto t0 = std::chrono::steady_clock::now(); SpotAnalyze(experiment, spot_finding_settings, spots, message); @@ -114,7 +111,7 @@ void MXAnalysisAfterFPGA::Process(DataMessage &message, const SpotFindingSetting spot_finding_time_total += (t1 - t0); if (spot_finding_settings.indexing) - indexer.ProcessImage(message, spot_finding_settings, message.image, *prediction, mapping, profile); + indexer.ProcessImage(message, spot_finding_settings, message.image, *prediction); } if (spot_finding_timing_active) { diff --git a/image_analysis/MXAnalysisAfterFPGA.h b/image_analysis/MXAnalysisAfterFPGA.h index e3b708ac..ccddf87a 100644 --- a/image_analysis/MXAnalysisAfterFPGA.h +++ b/image_analysis/MXAnalysisAfterFPGA.h @@ -28,10 +28,9 @@ class MXAnalysisAfterFPGA { std::chrono::duration spot_finding_time_total{0.0}; bool spot_finding_timing_active = false; - const AzimuthalIntegrationMapping &mapping; public: - MXAnalysisAfterFPGA(const DiffractionExperiment& experiment, IndexAndRefine &indexer, AzimuthalIntegrationMapping &mapping); + MXAnalysisAfterFPGA(const DiffractionExperiment& experiment, IndexAndRefine &indexer); void ReadFromFPGA(const DeviceOutput* output, const SpotFindingSettings& settings, @@ -41,5 +40,5 @@ public: const SpotFindingSettings &settings, size_t module_number); - void Process(DataMessage &message, const SpotFindingSettings& settings, const AzimuthalIntegrationProfile &profile); + void Process(DataMessage &message, const SpotFindingSettings& settings); }; diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 62b6379f..6a638213 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -92,7 +92,7 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, if (spot_finding_settings.indexing) indexer.ProcessImage(output, spot_finding_settings, CompressedImage(preprocessor_buffer->getBuffer(), experiment.GetXPixelsNum(), experiment.GetYPixelsNum()), - *prediction, integration, profile); + *prediction, &integration, &profile); } output.max_viable_pixel_value = ret.max_value; diff --git a/image_analysis/bragg_integration/BraggIntegrate2D.cpp b/image_analysis/bragg_integration/BraggIntegrate2D.cpp index 35ec4849..b616ecc1 100644 --- a/image_analysis/bragg_integration/BraggIntegrate2D.cpp +++ b/image_analysis/bragg_integration/BraggIntegrate2D.cpp @@ -8,62 +8,65 @@ #include namespace { - template - float Median(std::vector &values) { - if (values.empty()) - return 0.0f; - const size_t middle = values.size() / 2; - std::nth_element(values.begin(), values.begin() + middle, values.end()); +template +float Median(std::vector &values) { + if (values.empty()) + return 0.0f; - if (values.size() % 2 == 1) - return static_cast(values[middle]); + const size_t middle = values.size() / 2; + std::nth_element(values.begin(), values.begin() + middle, values.end()); - const T upper = values[middle]; - std::nth_element(values.begin(), values.begin() + middle - 1, values.begin() + middle); - const T lower = values[middle - 1]; + if (values.size() % 2 == 1) + return static_cast(values[middle]); - return 0.5f * static_cast(lower + upper); - } + const T upper = values[middle]; + std::nth_element(values.begin(), values.begin() + middle - 1, values.begin() + middle); + const T lower = values[middle - 1]; - void MarkReflectionMask(std::vector &mask, - size_t xpixel, size_t ypixel, - const Reflection &r, float r_2, float r_2_sq) { - int64_t x0 = std::floor(r.predicted_x - r_2 - 1.0f); - int64_t x1 = std::ceil(r.predicted_x + r_2 + 1.0f); - int64_t y0 = std::floor(r.predicted_y - r_2 - 1.0f); - int64_t y1 = std::ceil(r.predicted_y + r_2 + 1.0f); + return 0.5f * static_cast(lower + upper); +} - if (x0 < 0) - x0 = 0; - if (y0 < 0) - y0 = 0; - if (x1 >= static_cast(xpixel)) - x1 = static_cast(xpixel) - 1; - if (y1 >= static_cast(ypixel)) - y1 = static_cast(ypixel) - 1; +void MarkReflectionMask(std::vector &mask, + size_t xpixel, size_t ypixel, + const Reflection &r, float r_2, float r_2_sq) { - for (int64_t y = y0; y <= y1; y++) { - for (int64_t x = x0; x <= x1; x++) { - const float dist_sq = (x - r.predicted_x) * (x - r.predicted_x) - + (y - r.predicted_y) * (y - r.predicted_y); - if (dist_sq < r_2_sq) - mask[y * xpixel + x] = 1; - } + int64_t x0 = std::floor(r.predicted_x - r_2 - 1.0f); + int64_t x1 = std::ceil(r.predicted_x + r_2 + 1.0f); + int64_t y0 = std::floor(r.predicted_y - r_2 - 1.0f); + int64_t y1 = std::ceil(r.predicted_y + r_2 + 1.0f); + + if (x0 < 0) + x0 = 0; + if (y0 < 0) + y0 = 0; + if (x1 >= static_cast(xpixel)) + x1 = static_cast(xpixel) - 1; + if (y1 >= static_cast(ypixel)) + y1 = static_cast(ypixel) - 1; + + for (int64_t y = y0; y <= y1; y++) { + for (int64_t x = x0; x <= x1; x++) { + const float dist_sq = (x - r.predicted_x) * (x - r.predicted_x) + + (y - r.predicted_y) * (y - r.predicted_y); + if (dist_sq < r_2_sq) + mask[y * xpixel + x] = 1; } } +} - std::vector BuildReflectionMask(const std::vector &predicted, - size_t npredicted, - size_t xpixel, size_t ypixel, - float r_2, float r_2_sq) { - std::vector mask(xpixel * ypixel, 0); +std::vector BuildReflectionMask(const std::vector &predicted, + size_t npredicted, + size_t xpixel, size_t ypixel, + float r_2, float r_2_sq) { + std::vector mask(xpixel * ypixel, 0); - for (size_t i = 0; i < npredicted; i++) - MarkReflectionMask(mask, xpixel, ypixel, predicted.at(i), r_2, r_2_sq); + for (size_t i = 0; i < npredicted; i++) + MarkReflectionMask(mask, xpixel, ypixel, predicted.at(i), r_2, r_2_sq); + + return mask; +} - return mask; - } } // namespace template @@ -72,6 +75,7 @@ void IntegrateReflection(Reflection &r, const T *image, const std::vector -void IntegrateReflectionAzim(Reflection &r, const T *image, - size_t xpixel, size_t ypixel, - int64_t special_value, int64_t saturation, - float r_1, float r_1_sq, - const std::vector &pixel_to_bin, - const std::vector &bkg_mean, - float minimum_sigma_in_regards_to_i) { - int64_t x0 = std::floor(r.predicted_x - r_1 - 1.0); - int64_t x1 = std::ceil(r.predicted_x + r_1 + 1.0); - int64_t y0 = std::floor(r.predicted_y - r_1 - 1.0); - int64_t y1 = std::ceil(r.predicted_y + r_1 + 1.0); - x0 = std::max(0L, x0); - x1 = std::min(static_cast(xpixel - 1), x1); - y0 = std::max(0L, y0); - y1 = std::min(static_cast(ypixel - 1), y1); - - int64_t I_sum = 0; - int64_t I_npixel_inner = 0; - int64_t I_npixel_integrated = 0; - int64_t I_sum_x = 0; - int64_t I_sum_y = 0; - - std::vector bkg_values; - bkg_values.reserve(static_cast((x1 - x0 + 1) * (y1 - y0 + 1))); - - float bkg = 0.0; - - for (int64_t y = y0; y <= y1; y++) { - for (int64_t x = x0; x <= x1; x++) { - const float dist_sq = (x - r.predicted_x) * (x - r.predicted_x) - + (y - r.predicted_y) * (y - r.predicted_y); - const auto pixel = image[y * xpixel + x]; - auto azim_bint = pixel_to_bin[y * xpixel + x]; - - if (dist_sq < r_1_sq) - I_npixel_inner++; - - if (azim_bint >= bkg_mean.size()) - continue; - - if (pixel == special_value || pixel == saturation) - continue; - - if (dist_sq < r_1_sq) { - bkg += bkg_values.at(azim_bint); - I_sum += pixel; - I_sum_x += x * pixel; - I_sum_y += y * pixel; - I_npixel_integrated++; - } - } - } - - if ((I_npixel_integrated == I_npixel_inner) && (bkg_values.size() > 5)) { - r.bkg = bkg / static_cast(I_npixel_integrated); - r.I = static_cast(I_sum) - bkg; - if (I_sum > 0) { - r.observed_x = static_cast(I_sum_x) / static_cast(I_sum); - r.observed_y = static_cast(I_sum_y) / static_cast(I_sum); - } - - // sigma is max of the: - // - 1 (for zero photons) - // - Poisson noise (sqrt(I_sum)) (for in between) - // - minimum_sigma_in_regards_to_i of Intensity (for very large numbers) - r.sigma = std::max(1.0f, r.I * minimum_sigma_in_regards_to_i); - if (I_sum > 0) - r.sigma = std::max(r.sigma, std::sqrt(static_cast(I_sum))); - r.observed = true; - } else { - r.I = 0; - r.bkg = 0; - r.sigma = NAN; - r.observed = false; - } -} - template std::vector IntegrateInternal(const DiffractionExperiment &experiment, const CompressedImage &image, - const AzimuthalIntegrationMapping &mapping, - const AzimuthalIntegrationProfile &profile, const std::vector &predicted, size_t npredicted, int64_t special_value, int64_t saturation, int64_t image_number) { + std::vector ret; ret.reserve(npredicted); @@ -239,36 +164,24 @@ std::vector IntegrateInternal(const DiffractionExperiment &experimen std::vector buffer; auto ptr = reinterpret_cast(image.GetUncompressedPtr(buffer)); - const float r_1 = settings.GetR1(); - const float r_2 = settings.GetR2(); const float r_3 = settings.GetR3(); - - const float r_1_sq = r_1 * r_1; - const float r_2_sq = r_2 * r_2; - const float r_3_sq = r_3 * r_3; + const float r_1_sq = settings.GetR1() * settings.GetR1(); + const float r_2 = settings.GetR2(); + const float r_2_sq = settings.GetR2() * settings.GetR2(); + const float r_3_sq = settings.GetR3() * settings.GetR3(); const float minimum_sigma_in_regards_to_i = settings.GetMinimumSigmaInRegardsToI(); const auto reflection_mask = BuildReflectionMask(predicted, npredicted, image.GetWidth(), image.GetHeight(), r_2, r_2_sq); - const auto &pixel_to_bin = mapping.GetPixelToBin(); - const auto p = profile.GetResult(); - for (int i = 0; i < npredicted; i++) { auto r = predicted.at(i); - if (settings.IsUseAzimProfile()) { - IntegrateReflectionAzim(r, ptr, image.GetWidth(), image.GetHeight(), - special_value, saturation, - r_1, r_1_sq, pixel_to_bin, p, minimum_sigma_in_regards_to_i); - } else { - IntegrateReflection(r, ptr, reflection_mask, image.GetWidth(), image.GetHeight(), special_value, saturation, - r_3, r_1_sq, r_2_sq, r_3_sq, minimum_sigma_in_regards_to_i); - } + IntegrateReflection(r, ptr, reflection_mask, image.GetWidth(), image.GetHeight(), special_value, saturation, + r_3, r_1_sq, r_2_sq, r_3_sq, minimum_sigma_in_regards_to_i); if (r.observed) { if (experiment.GetPolarizationFactor()) - r.rlp /= geom.CalcAzIntPolarizationCorr(r.predicted_x, r.predicted_y, - experiment.GetPolarizationFactor().value()); + r.rlp /= geom.CalcAzIntPolarizationCorr(r.predicted_x, r.predicted_y, experiment.GetPolarizationFactor().value()); r.image_scale_corr = r.rlp / r.partiality; r.image_number = static_cast(image_number); ret.emplace_back(r); @@ -280,8 +193,6 @@ std::vector IntegrateInternal(const DiffractionExperiment &experimen std::vector BraggIntegrate2D(const DiffractionExperiment &experiment, const CompressedImage &image, const std::vector &predicted, - const AzimuthalIntegrationMapping &mapping, - const AzimuthalIntegrationProfile &profile, size_t npredicted, int64_t image_number) { if (image.GetCompressedSize() == 0 || predicted.empty()) @@ -289,23 +200,17 @@ std::vector BraggIntegrate2D(const DiffractionExperiment &experiment switch (image.GetMode()) { case CompressedImageMode::Int8: - return IntegrateInternal(experiment, image, mapping, profile, predicted, npredicted, - INT8_MIN, INT8_MAX, image_number); + return IntegrateInternal(experiment, image, predicted, npredicted, INT8_MIN, INT8_MAX, image_number); case CompressedImageMode::Int16: - return IntegrateInternal(experiment, image, mapping, profile, predicted, npredicted, - INT16_MIN, INT16_MAX, image_number); + return IntegrateInternal(experiment, image, predicted, npredicted, INT16_MIN, INT16_MAX, image_number); case CompressedImageMode::Int32: - return IntegrateInternal(experiment, image, mapping, profile, predicted, npredicted, - INT32_MIN, INT32_MAX, image_number); + return IntegrateInternal(experiment, image, predicted, npredicted, INT32_MIN, INT32_MAX, image_number); case CompressedImageMode::Uint8: - return IntegrateInternal(experiment, image, mapping, profile, predicted, npredicted, - UINT8_MAX, UINT8_MAX, image_number); + return IntegrateInternal(experiment, image, predicted, npredicted, UINT8_MAX, UINT8_MAX, image_number); case CompressedImageMode::Uint16: - return IntegrateInternal(experiment, image, mapping, profile, predicted, npredicted, - UINT16_MAX, UINT16_MAX, image_number); + return IntegrateInternal(experiment, image, predicted, npredicted, UINT16_MAX, UINT16_MAX, image_number); case CompressedImageMode::Uint32: - return IntegrateInternal(experiment, image, mapping, profile, predicted, npredicted, - UINT32_MAX, UINT32_MAX, image_number); + return IntegrateInternal(experiment, image, predicted, npredicted, UINT32_MAX, UINT32_MAX, image_number); default: throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Image mode not supported"); diff --git a/image_analysis/bragg_integration/BraggIntegrate2D.h b/image_analysis/bragg_integration/BraggIntegrate2D.h index 3045ff7d..353757b5 100644 --- a/image_analysis/bragg_integration/BraggIntegrate2D.h +++ b/image_analysis/bragg_integration/BraggIntegrate2D.h @@ -1,17 +1,17 @@ // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only -#pragma once +#ifndef JFJOCH_BRAGGINTEGRATE2D_H +#define JFJOCH_BRAGGINTEGRATE2D_H #include #include "../../common/DiffractionExperiment.h" #include "../../common/Reflection.h" -#include "../../common/AzimuthalIntegrationProfile.h" std::vector BraggIntegrate2D(const DiffractionExperiment &experiment, const CompressedImage &image, const std::vector &predicted, - const AzimuthalIntegrationMapping &mapping, - const AzimuthalIntegrationProfile &profile, size_t npredicted, - int64_t image_number); \ No newline at end of file + int64_t image_number); + +#endif //JFJOCH_BRAGGINTEGRATE2D_H diff --git a/receiver/JFJochReceiver.cpp b/receiver/JFJochReceiver.cpp index 1aa8ce1e..64e82e33 100644 --- a/receiver/JFJochReceiver.cpp +++ b/receiver/JFJochReceiver.cpp @@ -31,8 +31,7 @@ JFJochReceiver::JFJochReceiver(const DiffractionExperiment &in_experiment, serialmx_filter(in_experiment), numa_policy(in_numa_policy), pixel_mask(in_pixel_mask), - indexer(experiment, indexing_thread_pool), - az_int_mapping(experiment, pixel_mask) { + indexer(experiment, indexing_thread_pool) { logger.Info("Initializing receiver"); // Ensure there is nothing running for now if (!image_buffer.Finalize(std::chrono::seconds(1))) @@ -44,7 +43,13 @@ JFJochReceiver::JFJochReceiver(const DiffractionExperiment &in_experiment, current_status.SetEfficiency({}); current_status.SetStatus(JFJochReceiverStatus{}); // GetStatus() is virtual function and cannot be called yet! - plots.Setup(experiment, az_int_mapping); + auto start_time_point = std::chrono::steady_clock::now(); + az_int_mapping = std::make_unique(experiment, pixel_mask); + auto end_time_point = std::chrono::steady_clock::now(); + auto duration = std::chrono::duration(end_time_point - start_time_point); + logger.Info("Azimuthal integration mapping done in {:.5f} s with {} threads", duration.count(), az_int_mapping->GetNThreads()); + + plots.Setup(experiment, *az_int_mapping); push_images_to_writer = (experiment.GetImageNum() > 0) && (!experiment.GetFilePrefix().empty()); } @@ -106,13 +111,13 @@ void JFJochReceiver::SendStartMessage() { StartMessage message{}; experiment.FillMessage(message); message.arm_date = time_UTC(std::chrono::system_clock::now()); - message.az_int_q_bin_count = az_int_mapping.GetQBinCount(); - message.az_int_bin_to_q = az_int_mapping.GetBinToQ(); - message.az_int_bin_to_two_theta = az_int_mapping.GetBinToTwoTheta(); - message.az_int_phi_bin_count = az_int_mapping.GetAzimuthalBinCount(); - if (az_int_mapping.GetAzimuthalBinCount() > 1) { - message.az_int_bin_to_phi = az_int_mapping.GetBinToPhi(); - message.az_int_map = az_int_mapping.GetPixelToBin(); + message.az_int_q_bin_count = az_int_mapping->GetQBinCount(); + message.az_int_bin_to_q = az_int_mapping->GetBinToQ(); + message.az_int_bin_to_two_theta = az_int_mapping->GetBinToTwoTheta(); + message.az_int_phi_bin_count = az_int_mapping->GetAzimuthalBinCount(); + if (az_int_mapping->GetAzimuthalBinCount() > 1) { + message.az_int_bin_to_phi = az_int_mapping->GetBinToPhi(); + message.az_int_map = az_int_mapping->GetPixelToBin(); } message.writer_notification_zmq_addr = image_pusher.GetWriterNotificationSocketAddress(); message.rois = experiment.ROI().ExportMetadata(); diff --git a/receiver/JFJochReceiver.h b/receiver/JFJochReceiver.h index 79be2802..83918987 100644 --- a/receiver/JFJochReceiver.h +++ b/receiver/JFJochReceiver.h @@ -80,7 +80,7 @@ protected: std::vector> adu_histogram_module; PixelMask pixel_mask; - AzimuthalIntegrationMapping az_int_mapping; + std::unique_ptr az_int_mapping; std::optional max_delay; std::mutex max_delay_mutex; diff --git a/receiver/JFJochReceiverFPGA.cpp b/receiver/JFJochReceiverFPGA.cpp index a3ac6384..802b7faf 100644 --- a/receiver/JFJochReceiverFPGA.cpp +++ b/receiver/JFJochReceiverFPGA.cpp @@ -297,7 +297,7 @@ void JFJochReceiverFPGA::FrameTransformationThread(uint32_t threadid) { try { numa_policy.Bind(threadid); - analyzer = std::make_unique(experiment, indexer, az_int_mapping); + analyzer = std::make_unique(experiment, indexer); } catch (const JFJochException &e) { frame_transformation_ready.count_down(); logger.Error("Thread setup error {}", e.what()); @@ -309,6 +309,9 @@ void JFJochReceiverFPGA::FrameTransformationThread(uint32_t threadid) { frame_transformation_ready.count_down(); + uint16_t az_int_min_bin = std::floor(az_int_mapping->QToBin(experiment.GetLowQForBkgEstimate_recipA())); + uint16_t az_int_max_bin = std::ceil(az_int_mapping->QToBin(experiment.GetHighQForBkgEstimate_recipA())); + int64_t image_number; while (images_to_go.Get(image_number) != 0) { try { @@ -335,7 +338,7 @@ void JFJochReceiverFPGA::FrameTransformationThread(uint32_t threadid) { ImageMetadata metadata(experiment); - AzimuthalIntegrationProfile az_int_profile_image(az_int_mapping); + AzimuthalIntegrationProfile az_int_profile_image(*az_int_mapping); auto local_spot_finding_settings = GetSpotFindingSettings(); @@ -425,7 +428,7 @@ void JFJochReceiverFPGA::FrameTransformationThread(uint32_t threadid) { experiment.GetYPixelsNum(), experiment.GetImageMode(), CompressionAlgorithm::NO_COMPRESSION); - analyzer->Process(message, local_spot_finding_settings, az_int_profile_image); + analyzer->Process(message, local_spot_finding_settings); auto status = image_buffer.GetStatus(); message.receiver_buf_available = status.available_slots; @@ -678,5 +681,5 @@ void JFJochReceiverFPGA::LoadCalibrationToFPGA(uint16_t data_stream) { acquisition_device[data_stream].InitializeROIMap(experiment, roi_map); // Initialize data processing - acquisition_device[data_stream].InitializeDataProcessing(experiment, az_int_mapping); + acquisition_device[data_stream].InitializeDataProcessing(experiment, *az_int_mapping); } diff --git a/receiver/JFJochReceiverLite.cpp b/receiver/JFJochReceiverLite.cpp index 8c04b1cf..e986956e 100644 --- a/receiver/JFJochReceiverLite.cpp +++ b/receiver/JFJochReceiverLite.cpp @@ -242,7 +242,7 @@ void JFJochReceiverLite::DataAnalysisThread(uint32_t id) { measurement_started.wait(); try { - analysis = std::make_unique(experiment, az_int_mapping, pixel_mask, indexer); + analysis = std::make_unique(experiment, *az_int_mapping, pixel_mask, indexer); } catch (const JFJochException &e) { Cancel(e); return; @@ -279,7 +279,7 @@ void JFJochReceiverLite::DataAnalysisThread(uint32_t id) { auto image_start_time = std::chrono::high_resolution_clock::now(); - AzimuthalIntegrationProfile profile(az_int_mapping); + AzimuthalIntegrationProfile profile(*az_int_mapping); analysis->Analyze(data_msg, profile, GetSpotFindingSettings()); auto image_end_time = std::chrono::high_resolution_clock::now(); diff --git a/tests/BraggIntegrate2DTest.cpp b/tests/BraggIntegrate2DTest.cpp index 149515c0..75ba2a42 100644 --- a/tests/BraggIntegrate2DTest.cpp +++ b/tests/BraggIntegrate2DTest.cpp @@ -59,10 +59,6 @@ TEST_CASE("BraggIntegrate2D_RejectsReflectionsFromBackgroundUsingR2MaskAndMedian Reflection r1 = MakeReflection(10.0f, 10.0f); Reflection r2 = MakeReflection(16.0f, 10.0f); - PixelMask mask(experiment); - AzimuthalIntegrationMapping mapping(experiment, mask); - AzimuthalIntegrationProfile profile(mapping); - for (int y = 0; y < static_cast(height); y++) { for (int x = 0; x < static_cast(width); x++) { auto &pixel = image_data[y * width + x]; @@ -78,7 +74,7 @@ TEST_CASE("BraggIntegrate2D_RejectsReflectionsFromBackgroundUsingR2MaskAndMedian CompressedImage image(image_data, width, height); std::vector predicted = {r1, r2}; - auto integrated = BraggIntegrate2D(experiment, image, predicted, mapping, profile, predicted.size(), 17); + auto integrated = BraggIntegrate2D(experiment, image, predicted, predicted.size(), 17); REQUIRE(integrated.size() == 2); diff --git a/tools/jfjoch_scale.cpp b/tools/jfjoch_scale.cpp index e6da83f2..7f169940 100644 --- a/tools/jfjoch_scale.cpp +++ b/tools/jfjoch_scale.cpp @@ -148,6 +148,8 @@ int main(int argc, char **argv) { partiality_model = PartialityModel::Fixed; else if (strcmp(optarg, "rot") == 0) partiality_model = PartialityModel::Rotation; + else if (strcmp(optarg, "postref") == 0) + partiality_model = PartialityModel::Postrefinement; else { logger.Error("Invalid partiality mode: {}", optarg); print_usage();