diff --git a/common/BraggIntegrationSettings.cpp b/common/BraggIntegrationSettings.cpp index 3fe0a861..0528084e 100644 --- a/common/BraggIntegrationSettings.cpp +++ b/common/BraggIntegrationSettings.cpp @@ -87,3 +87,12 @@ 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 60ab4059..2c15921d 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,4 +28,5 @@ 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 98b1aed6..459af31e 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() && mapping && profile; + && !pixel_reference_.empty(); 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(), nrefl, msg.number); + i_outcome.reflections = BraggIntegrate2D(outcome.experiment, image, prediction.GetReflections(), mapping, profile, 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 7674392e..d73fc2f1 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 = nullptr, - const AzimuthalIntegrationProfile *profile = nullptr); + const AzimuthalIntegrationMapping &mapping, + const AzimuthalIntegrationProfile &profile); 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 a5347236..96d8430d 100644 --- a/image_analysis/MXAnalysisAfterFPGA.cpp +++ b/image_analysis/MXAnalysisAfterFPGA.cpp @@ -25,10 +25,12 @@ double stddev(const std::vector &v) { } -MXAnalysisAfterFPGA::MXAnalysisAfterFPGA(const DiffractionExperiment &in_experiment, IndexAndRefine &indexer) +MXAnalysisAfterFPGA::MXAnalysisAfterFPGA(const DiffractionExperiment &in_experiment, IndexAndRefine &indexer, + AzimuthalIntegrationMapping &mapping) : experiment(in_experiment), indexer(indexer), - prediction(CreateBraggPrediction(experiment.IsRotationIndexing())) { + prediction(CreateBraggPrediction(experiment.IsRotationIndexing())), + mapping(mapping) { if (experiment.IsSpotFindingEnabled()) find_spots = true; } @@ -103,7 +105,8 @@ void MXAnalysisAfterFPGA::ReadFromCPU(DeviceOutput *output, const SpotFindingSet } } -void MXAnalysisAfterFPGA::Process(DataMessage &message, const SpotFindingSettings& spot_finding_settings) { +void MXAnalysisAfterFPGA::Process(DataMessage &message, const SpotFindingSettings &spot_finding_settings, + const AzimuthalIntegrationProfile &profile) { if (find_spots && (state == State::Enabled)) { const auto t0 = std::chrono::steady_clock::now(); SpotAnalyze(experiment, spot_finding_settings, spots, message); @@ -111,7 +114,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); + indexer.ProcessImage(message, spot_finding_settings, message.image, *prediction, mapping, profile); } if (spot_finding_timing_active) { diff --git a/image_analysis/MXAnalysisAfterFPGA.h b/image_analysis/MXAnalysisAfterFPGA.h index ccddf87a..e3b708ac 100644 --- a/image_analysis/MXAnalysisAfterFPGA.h +++ b/image_analysis/MXAnalysisAfterFPGA.h @@ -28,9 +28,10 @@ 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); + MXAnalysisAfterFPGA(const DiffractionExperiment& experiment, IndexAndRefine &indexer, AzimuthalIntegrationMapping &mapping); void ReadFromFPGA(const DeviceOutput* output, const SpotFindingSettings& settings, @@ -40,5 +41,5 @@ public: const SpotFindingSettings &settings, size_t module_number); - void Process(DataMessage &message, const SpotFindingSettings& settings); + void Process(DataMessage &message, const SpotFindingSettings& settings, const AzimuthalIntegrationProfile &profile); }; diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 6a638213..62b6379f 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 b616ecc1..35ec4849 100644 --- a/image_analysis/bragg_integration/BraggIntegrate2D.cpp +++ b/image_analysis/bragg_integration/BraggIntegrate2D.cpp @@ -8,65 +8,62 @@ #include namespace { + template + float Median(std::vector &values) { + if (values.empty()) + return 0.0f; -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()); - const size_t middle = values.size() / 2; - std::nth_element(values.begin(), values.begin() + middle, values.end()); + if (values.size() % 2 == 1) + return static_cast(values[middle]); - if (values.size() % 2 == 1) - return static_cast(values[middle]); + const T upper = values[middle]; + std::nth_element(values.begin(), values.begin() + middle - 1, values.begin() + middle); + const T lower = values[middle - 1]; - const T upper = values[middle]; - std::nth_element(values.begin(), values.begin() + middle - 1, values.begin() + middle); - const T lower = values[middle - 1]; + return 0.5f * static_cast(lower + upper); + } - return 0.5f * static_cast(lower + upper); -} + 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); -void MarkReflectionMask(std::vector &mask, - size_t xpixel, size_t ypixel, - const Reflection &r, float r_2, float r_2_sq) { + 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; - 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; + 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); - - return mask; -} + for (size_t i = 0; i < npredicted; i++) + MarkReflectionMask(mask, xpixel, ypixel, predicted.at(i), r_2, r_2_sq); + return mask; + } } // namespace template @@ -75,7 +72,6 @@ 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); @@ -164,24 +239,36 @@ std::vector IntegrateInternal(const DiffractionExperiment &experimen std::vector buffer; auto ptr = reinterpret_cast(image.GetUncompressedPtr(buffer)); - const float r_3 = settings.GetR3(); - const float r_1_sq = settings.GetR1() * settings.GetR1(); + const float r_1 = 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 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 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); - 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 (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); + } 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); @@ -193,6 +280,8 @@ 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()) @@ -200,17 +289,23 @@ std::vector BraggIntegrate2D(const DiffractionExperiment &experiment switch (image.GetMode()) { case CompressedImageMode::Int8: - return IntegrateInternal(experiment, image, predicted, npredicted, INT8_MIN, INT8_MAX, image_number); + return IntegrateInternal(experiment, image, mapping, profile, predicted, npredicted, + INT8_MIN, INT8_MAX, image_number); case CompressedImageMode::Int16: - return IntegrateInternal(experiment, image, predicted, npredicted, INT16_MIN, INT16_MAX, image_number); + return IntegrateInternal(experiment, image, mapping, profile, predicted, npredicted, + INT16_MIN, INT16_MAX, image_number); case CompressedImageMode::Int32: - return IntegrateInternal(experiment, image, predicted, npredicted, INT32_MIN, INT32_MAX, image_number); + return IntegrateInternal(experiment, image, mapping, profile, predicted, npredicted, + INT32_MIN, INT32_MAX, image_number); case CompressedImageMode::Uint8: - return IntegrateInternal(experiment, image, predicted, npredicted, UINT8_MAX, UINT8_MAX, image_number); + return IntegrateInternal(experiment, image, mapping, profile, predicted, npredicted, + UINT8_MAX, UINT8_MAX, image_number); case CompressedImageMode::Uint16: - return IntegrateInternal(experiment, image, predicted, npredicted, UINT16_MAX, UINT16_MAX, image_number); + return IntegrateInternal(experiment, image, mapping, profile, predicted, npredicted, + UINT16_MAX, UINT16_MAX, image_number); case CompressedImageMode::Uint32: - return IntegrateInternal(experiment, image, predicted, npredicted, UINT32_MAX, UINT32_MAX, image_number); + return IntegrateInternal(experiment, image, mapping, profile, 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 353757b5..3045ff7d 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 -#ifndef JFJOCH_BRAGGINTEGRATE2D_H -#define JFJOCH_BRAGGINTEGRATE2D_H +#pragma once #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); - -#endif //JFJOCH_BRAGGINTEGRATE2D_H + int64_t image_number); \ No newline at end of file diff --git a/receiver/JFJochReceiver.cpp b/receiver/JFJochReceiver.cpp index 64e82e33..1aa8ce1e 100644 --- a/receiver/JFJochReceiver.cpp +++ b/receiver/JFJochReceiver.cpp @@ -31,7 +31,8 @@ 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) { + indexer(experiment, indexing_thread_pool), + az_int_mapping(experiment, pixel_mask) { logger.Info("Initializing receiver"); // Ensure there is nothing running for now if (!image_buffer.Finalize(std::chrono::seconds(1))) @@ -43,13 +44,7 @@ JFJochReceiver::JFJochReceiver(const DiffractionExperiment &in_experiment, current_status.SetEfficiency({}); current_status.SetStatus(JFJochReceiverStatus{}); // GetStatus() is virtual function and cannot be called yet! - 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); + plots.Setup(experiment, az_int_mapping); push_images_to_writer = (experiment.GetImageNum() > 0) && (!experiment.GetFilePrefix().empty()); } @@ -111,13 +106,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 83918987..79be2802 100644 --- a/receiver/JFJochReceiver.h +++ b/receiver/JFJochReceiver.h @@ -80,7 +80,7 @@ protected: std::vector> adu_histogram_module; PixelMask pixel_mask; - std::unique_ptr az_int_mapping; + AzimuthalIntegrationMapping az_int_mapping; std::optional max_delay; std::mutex max_delay_mutex; diff --git a/receiver/JFJochReceiverFPGA.cpp b/receiver/JFJochReceiverFPGA.cpp index 802b7faf..a3ac6384 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); + analyzer = std::make_unique(experiment, indexer, az_int_mapping); } catch (const JFJochException &e) { frame_transformation_ready.count_down(); logger.Error("Thread setup error {}", e.what()); @@ -309,9 +309,6 @@ 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 { @@ -338,7 +335,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(); @@ -428,7 +425,7 @@ void JFJochReceiverFPGA::FrameTransformationThread(uint32_t threadid) { experiment.GetYPixelsNum(), experiment.GetImageMode(), CompressionAlgorithm::NO_COMPRESSION); - analyzer->Process(message, local_spot_finding_settings); + analyzer->Process(message, local_spot_finding_settings, az_int_profile_image); auto status = image_buffer.GetStatus(); message.receiver_buf_available = status.available_slots; @@ -681,5 +678,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 e986956e..8c04b1cf 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 75ba2a42..149515c0 100644 --- a/tests/BraggIntegrate2DTest.cpp +++ b/tests/BraggIntegrate2DTest.cpp @@ -59,6 +59,10 @@ 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]; @@ -74,7 +78,7 @@ TEST_CASE("BraggIntegrate2D_RejectsReflectionsFromBackgroundUsingR2MaskAndMedian CompressedImage image(image_data, width, height); std::vector predicted = {r1, r2}; - auto integrated = BraggIntegrate2D(experiment, image, predicted, predicted.size(), 17); + auto integrated = BraggIntegrate2D(experiment, image, predicted, mapping, profile, predicted.size(), 17); REQUIRE(integrated.size() == 2); diff --git a/tools/jfjoch_scale.cpp b/tools/jfjoch_scale.cpp index 7f169940..e6da83f2 100644 --- a/tools/jfjoch_scale.cpp +++ b/tools/jfjoch_scale.cpp @@ -148,8 +148,6 @@ 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();