From 953089009f88056ad71c81e4f60f4e534fabd12b Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 2 Jun 2026 15:50:11 +0200 Subject: [PATCH] AzimuthalIntegrationProfile: Add standard deviation + save pixel count + improve handling of special cases (missing image) --- common/AzimuthalIntegrationProfile.cpp | 35 ++++++++- common/AzimuthalIntegrationProfile.h | 7 +- common/JFJochMessages.h | 3 + frame_serialize/CBORStream2Deserializer.cpp | 4 ++ frame_serialize/CBORStream2Serializer.cpp | 4 ++ image_analysis/MXAnalysisWithoutFPGA.cpp | 3 + image_analysis/azint/AzIntEngineCPU.cpp | 2 +- image_analysis/azint/AzIntEngineGPU.cu | 2 +- receiver/JFJochReceiverFPGA.cpp | 1 + tests/AzimuthalIntegrationTest.cpp | 78 ++++++++++++++++++--- tests/CBORTest.cpp | 13 +++- writer/HDF5DataFilePluginAzInt.cpp | 34 ++++++++- writer/HDF5DataFilePluginAzInt.h | 5 ++ 13 files changed, 175 insertions(+), 16 deletions(-) diff --git a/common/AzimuthalIntegrationProfile.cpp b/common/AzimuthalIntegrationProfile.cpp index d1f128c9..4ccfb609 100644 --- a/common/AzimuthalIntegrationProfile.cpp +++ b/common/AzimuthalIntegrationProfile.cpp @@ -7,8 +7,15 @@ inline float sum_to_count(float sum, uint64_t count) { if (count == 0) return NAN; - else - return sum / (static_cast(count)); + return sum / (static_cast(count)); +} + +inline float calc_std(float sum, float sum2, uint64_t count) { + if (count == 0 || count == 1) + return NAN; + const auto fp_count = static_cast(count); + const float variance = (sum2 - sum * sum / fp_count) / (fp_count - 1); + return std::sqrt(variance); } AzimuthalIntegrationProfile::AzimuthalIntegrationProfile(const AzimuthalIntegrationMapping &mapping) @@ -46,13 +53,20 @@ void AzimuthalIntegrationProfile::Add(int64_t bin, int64_t value) { count[bin]++; } -void AzimuthalIntegrationProfile::Add(const std::vector &in_sum, const std::vector &in_count) { +void AzimuthalIntegrationProfile::Add(const std::vector &in_sum, + const std::vector &in_sum2, + const std::vector &in_count) { std::unique_lock ul(m); if ((in_sum.size() == sum.size()) && (in_count.size() == count.size())) { for (int i = 0; i < sum.size(); i++) { sum[i] += in_sum[i]; count[i] += in_count[i]; } + if (in_sum2.size() == sum2.size()) { + for (int i = 0; i < sum.size(); i++) { + sum2[i] += in_sum2[i]; + } + } } else if (!in_sum.empty() && !in_count.empty()) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Mismatch in size of sum/count datasets"); } @@ -66,6 +80,20 @@ std::vector AzimuthalIntegrationProfile::GetResult() const { return rad_int_profile; } +std::vector AzimuthalIntegrationProfile::GetStd() const { + std::unique_lock ul(m); + + std::vector rad_int_profile(sum.size(), 0); + for (int i = 0; i < sum.size(); i++) + rad_int_profile[i] = calc_std(sum[i], sum2[i], count[i]); + return rad_int_profile; +} + +std::vector AzimuthalIntegrationProfile::GetPixelCount() const { + std::unique_lock ul(m); + return count; +} + std::vector AzimuthalIntegrationProfile::GetResult1D() const { std::unique_lock ul(m); @@ -151,6 +179,7 @@ AzimuthalIntegrationProfile &AzimuthalIntegrationProfile::operator+=(const Azimu for (int i = 0; i < sum.size(); i++) { sum[i] += other.sum[i]; + sum2[i] += other.sum2[i]; count[i] += other.count[i]; } diff --git a/common/AzimuthalIntegrationProfile.h b/common/AzimuthalIntegrationProfile.h index fb78fbaa..87d681db 100644 --- a/common/AzimuthalIntegrationProfile.h +++ b/common/AzimuthalIntegrationProfile.h @@ -32,10 +32,15 @@ public: void Clear(const AzimuthalIntegrationMapping &mapping); void SetTitle(const std::string& input); void Add(const DeviceOutput &result); - void Add(const std::vector &sum, const std::vector &count); + void Add(const std::vector &sum, + const std::vector &sum2, + const std::vector &count); void Add(int64_t bin, int64_t value); std::vector GetResult() const; + std::vector GetStd() const; + std::vector GetPixelCount() const; + std::vector GetResult1D() const; float GetMeanValueOfBins(uint16_t min_bin, uint16_t max_bin) const; diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index 44a617d3..de830204 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -107,6 +107,9 @@ struct DataMessage { std::vector spot_plot_one_over_d_square; std::vector az_int_profile; + std::vector az_int_profile_std; + std::vector az_int_profile_count; + std::optional bkg_estimate; std::optional indexing_result; diff --git a/frame_serialize/CBORStream2Deserializer.cpp b/frame_serialize/CBORStream2Deserializer.cpp index d5ff05ff..0b8e0a57 100644 --- a/frame_serialize/CBORStream2Deserializer.cpp +++ b/frame_serialize/CBORStream2Deserializer.cpp @@ -711,6 +711,10 @@ namespace { message.spot_count_ice_rings = GetCBORUInt(value); else if (key == "az_int_profile") GetCBORFloatArray(value, message.az_int_profile); + else if (key == "az_int_profile_std") + GetCBORFloatArray(value, message.az_int_profile_std); + else if (key == "az_int_profile_count") + GetCBORUInt64Array(value, message.az_int_profile_count); else if (key == "indexing_result") message.indexing_result = GetCBORBool(value); else if (key == "indexing_lattice") { diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index 6473c0c4..4966be47 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -742,6 +742,10 @@ void CBORStream2Serializer::SerializeImageInternal(CborEncoder &mapEncoder, cons CBOR_ENC(mapEncoder, "spot_count_low_res", message.spot_count_low_res); CBOR_ENC(mapEncoder, "spot_count_indexed", message.spot_count_indexed); CBOR_ENC(mapEncoder, "az_int_profile", message.az_int_profile); + if (!message.az_int_profile_std.empty()) + CBOR_ENC(mapEncoder, "az_int_profile_std", message.az_int_profile_std); + CBOR_ENC(mapEncoder, "az_int_profile_count", message.az_int_profile_count); + CBOR_ENC(mapEncoder, "indexing_result", message.indexing_result); if (message.indexing_lattice) CBOR_ENC(mapEncoder, "indexing_lattice", message.indexing_lattice->GetVector()); diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index c7e55fb6..7bade492 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -100,6 +100,9 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, output.error_pixel_count = ret.error_pixel_count; output.saturated_pixel_count = ret.saturated_pixel_count; output.az_int_profile = profile.GetResult(); + output.az_int_profile_count = profile.GetPixelCount(); + output.az_int_profile_std = profile.GetStd(); + output.bkg_estimate = profile.GetBkgEstimate(integration.Settings()); } diff --git a/image_analysis/azint/AzIntEngineCPU.cpp b/image_analysis/azint/AzIntEngineCPU.cpp index afaa6ebb..f8dc7321 100644 --- a/image_analysis/azint/AzIntEngineCPU.cpp +++ b/image_analysis/azint/AzIntEngineCPU.cpp @@ -31,5 +31,5 @@ void AzIntEngineCPU::Run(const ImagePreprocessorBuffer &image, AzimuthalIntegrat } profile.Clear(integration); - profile.Add(azint_sum, azint_count); + profile.Add(azint_sum, azint_sum2, azint_count); } diff --git a/image_analysis/azint/AzIntEngineGPU.cu b/image_analysis/azint/AzIntEngineGPU.cu index dfc81076..132f9761 100644 --- a/image_analysis/azint/AzIntEngineGPU.cu +++ b/image_analysis/azint/AzIntEngineGPU.cu @@ -139,5 +139,5 @@ void AzIntEngineGPU::Run(const ImagePreprocessorBuffer &image, AzimuthalIntegrat cuda_err(cudaStreamSynchronize(*stream)); profile.Clear(integration); - profile.Add(azint_sum, azint_count); + profile.Add(azint_sum, azint_sum2, azint_count); } diff --git a/receiver/JFJochReceiverFPGA.cpp b/receiver/JFJochReceiverFPGA.cpp index a237e0d2..18dc2d12 100644 --- a/receiver/JFJochReceiverFPGA.cpp +++ b/receiver/JFJochReceiverFPGA.cpp @@ -435,6 +435,7 @@ void JFJochReceiverFPGA::FrameTransformationThread(uint32_t threadid) { message.receiver_buf_in_preparation = status.preparation_slots; message.receiver_buf_in_sending = status.sending_slots; message.az_int_profile = az_int_profile_image.GetResult(); + message.az_int_profile_count = az_int_profile_image.GetPixelCount(); message.bkg_estimate = az_int_profile_image.GetBkgEstimate(experiment.GetAzimuthalIntegrationSettings()); scan_result.Add(message); diff --git a/tests/AzimuthalIntegrationTest.cpp b/tests/AzimuthalIntegrationTest.cpp index ac95a083..ed49b531 100644 --- a/tests/AzimuthalIntegrationTest.cpp +++ b/tests/AzimuthalIntegrationTest.cpp @@ -171,18 +171,19 @@ TEST_CASE("AzimuthalIntegrationProfile","[AzimuthalIntegration]") { AzimuthalIntegrationProfile profile(mapping); - std::vector sum(mapping.GetBinNumber()); + std::vector sum(mapping.GetBinNumber()), sum2(mapping.GetBinNumber()); std::vector count(mapping.GetBinNumber()); for (int i = 0; i < mapping.GetBinNumber(); i++) { sum[i] = i * i * 4; + sum2[i] = i * i * i * i * 4; count[i] = i; } - REQUIRE_NOTHROW(profile.Add(sum, count)); - REQUIRE_NOTHROW(profile.Add(sum, count)); + REQUIRE_NOTHROW(profile.Add(sum, sum2, count)); + REQUIRE_NOTHROW(profile.Add(sum, sum2, count)); std::vector sum_wr(mapping.GetBinNumber() - 1); - REQUIRE_THROWS(profile.Add(sum_wr, count)); + REQUIRE_THROWS(profile.Add(sum_wr, sum2, count)); auto plot = profile.GetPlot(); REQUIRE(plot.GetPlots().size() == 1); @@ -198,6 +199,63 @@ TEST_CASE("AzimuthalIntegrationProfile","[AzimuthalIntegration]") { } } +TEST_CASE("AzimuthalIntegrationProfile_GetStd","[AzimuthalIntegration]") { + DiffractionExperiment x(DetJF4M()); + x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000); + x.QSpacingForAzimInt_recipA(0.1).QRangeForAzimInt_recipA(0.1, 4); + + PixelMask pixel_mask(x); + AzimuthalIntegrationMapping mapping(x, pixel_mask); + + AzimuthalIntegrationProfile profile(mapping); + + REQUIRE(mapping.GetBinNumber() >= 4); + + std::vector sum(mapping.GetBinNumber()), sum2(mapping.GetBinNumber()); + std::vector count(mapping.GetBinNumber()); + + sum[0] = 2 + 3 + 4 + 5; + sum2[0] = 4 + 9 + 16 + 25; + count[0] = 4; + + sum[1] = 1 + 1; + sum2[1] = 1 + 1; + count[1] = 2; + + sum[2] = 1; + sum2[2] = 1; + count[2] = 1; + + sum[3] = 0; + sum2[3] = 0; + count[3] = 0; + + REQUIRE_NOTHROW(profile.Add(sum, sum2, count)); + + auto ret_mean = profile.GetResult(); + auto ret_stddev = profile.GetStd(); + auto ret_count = profile.GetPixelCount(); + REQUIRE(ret_mean.size() == mapping.GetBinNumber()); + REQUIRE(ret_stddev.size() == mapping.GetBinNumber()); + REQUIRE(ret_count.size() == mapping.GetBinNumber()); + + CHECK(ret_mean[0] == Catch::Approx(3.5)); + CHECK(ret_stddev[0] == Catch::Approx(std::sqrt(5.0/ 3.0))); + CHECK(ret_count[0] == 4); + + CHECK(ret_mean[1] == Catch::Approx(1.0)); + CHECK(ret_stddev[1] == 0.0f); + CHECK(ret_count[1] == 2); + + CHECK(ret_mean[2] == Catch::Approx(1.0)); + CHECK(std::isnan(ret_stddev[2])); + CHECK(ret_count[2] == 1); + + CHECK(std::isnan(ret_mean[3])); + CHECK(std::isnan(ret_stddev[3])); + CHECK(ret_count[3] == 0); +} + TEST_CASE("AzimuthalIntegrationProfile_operatorAdd","[AzimuthalIntegration]") { DiffractionExperiment x(DetJF4M()); x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000); @@ -208,14 +266,15 @@ TEST_CASE("AzimuthalIntegrationProfile_operatorAdd","[AzimuthalIntegration]") { AzimuthalIntegrationProfile profile0(mapping), profile1(mapping); - std::vector sum(mapping.GetBinNumber()); + std::vector sum(mapping.GetBinNumber()), sum2(mapping.GetBinNumber()); std::vector count(mapping.GetBinNumber()); for (int i = 0; i < mapping.GetBinNumber(); i++) { sum[i] = (i + 1) * i * 4; + sum2[i] = (i+ 1) * i * 5; count[i] = i + 1; } - REQUIRE_NOTHROW(profile0.Add(sum, count)); + REQUIRE_NOTHROW(profile0.Add(sum, sum2, count)); REQUIRE_NOTHROW(profile1 += profile0); auto plot = profile1.GetPlot(); @@ -240,13 +299,15 @@ TEST_CASE("AzimuthalIntegrationProfile_GetMeanValueOfBins","[AzimuthalIntegratio AzimuthalIntegrationProfile profile(mapping); std::vector sum(mapping.GetBinNumber()); + std::vector sum2(mapping.GetBinNumber()); std::vector count(mapping.GetBinNumber()); for (int i = 0; i < mapping.GetBinNumber(); i++) { sum[i] = i * i * 4; + sum2[i] = i * i * i * i * 4; count[i] = i; } - REQUIRE_NOTHROW(profile.Add(sum, count)); + REQUIRE_NOTHROW(profile.Add(sum, sum2, count)); REQUIRE(profile.GetMeanValueOfBins(0,2) == Catch::Approx((sum[0] + sum[1] + sum[2]) / double(count[0] + count[1] + count[2]))); REQUIRE(profile.GetMeanValueOfBins(5,7) == Catch::Approx((sum[5] + sum[6] + sum[7]) / double (count[5] + count[6] + count[7]))); @@ -276,6 +337,7 @@ TEST_CASE("AzimuthalIntegrationProfile_GetResult1D","[AzimuthalIntegration]") { REQUIRE(mapping.GetBinNumber() == 9); std::vector sum(mapping.GetBinNumber(), 0.0f); + std::vector sum2(mapping.GetBinNumber(), 20.0f); std::vector count(mapping.GetBinNumber(), 0); // Layout is [azimuth][q], flattened: @@ -299,7 +361,7 @@ TEST_CASE("AzimuthalIntegrationProfile_GetResult1D","[AzimuthalIntegration]") { sum[7] = 31; count[7] = 1; // az2 q1 sum[8] = 32; count[8] = 1; // az2 q2 - REQUIRE_NOTHROW(profile.Add(sum, count)); + REQUIRE_NOTHROW(profile.Add(sum, sum2, count)); auto result_1d = profile.GetResult1D(); diff --git a/tests/CBORTest.cpp b/tests/CBORTest.cpp index 0c871a4c..8690f51a 100644 --- a/tests/CBORTest.cpp +++ b/tests/CBORTest.cpp @@ -797,7 +797,9 @@ TEST_CASE("CBORSerialize_Image_Rad_Int_Profile", "[CBOR]") { DataMessage message{ .number = 789, .image = image, - .az_int_profile = {4.0, 5.0, 7.0, 12.0, 13.25, 0.125} + .az_int_profile = {4.0, 5.0, 7.0, 12.0, 13.25, 0.125}, + .az_int_profile_std = {1.0, 0.0, NAN, NAN, 8.0, 12.0}, + .az_int_profile_count = {3,2,0,1,5,7} }; REQUIRE_NOTHROW(serializer.SerializeImage(message)); @@ -813,6 +815,15 @@ TEST_CASE("CBORSerialize_Image_Rad_Int_Profile", "[CBOR]") { REQUIRE(image_array.image.GetCompressedSize() == test.size() * sizeof(uint16_t)); REQUIRE(memcmp(image_array.image.GetCompressed(), test.data(), test.size() * sizeof(uint16_t)) == 0); REQUIRE(image_array.az_int_profile == message.az_int_profile); + REQUIRE(image_array.az_int_profile_std.size() == message.az_int_profile_std.size()); + for (int i = 0; i < image_array.az_int_profile_std.size(); i++) { + if (std::isnan(message.az_int_profile_std[i])) + CHECK(std::isnan(image_array.az_int_profile_std[i])); + else + CHECK(message.az_int_profile_std[i] == Catch::Approx(image_array.az_int_profile_std[i])); + } + REQUIRE(image_array.az_int_profile_std == message.az_int_profile_std); + REQUIRE(image_array.az_int_profile_count == message.az_int_profile_count); } TEST_CASE("CBORSerialize_Image_Spots", "[CBOR]") { diff --git a/writer/HDF5DataFilePluginAzInt.cpp b/writer/HDF5DataFilePluginAzInt.cpp index 1f032e04..5274aae9 100644 --- a/writer/HDF5DataFilePluginAzInt.cpp +++ b/writer/HDF5DataFilePluginAzInt.cpp @@ -25,6 +25,14 @@ void HDF5DataFilePluginAzInt::OpenFile(HDF5File &data_file, const DataMessage &m data_file.SaveVector("/entry/azint/bin_to_phi", az_int_bin_to_phi, dim); az_int_image.reserve(images_per_file * azimuthal_bins * q_bins); + if (!msg.az_int_profile_count.empty()) { + save_profile_count = true; + az_int_image_count.reserve(images_per_file * azimuthal_bins * q_bins); + } + if (!msg.az_int_profile_std.empty()) { + save_profile_std = true; + az_int_image_std.reserve(images_per_file * azimuthal_bins * q_bins); + } } void HDF5DataFilePluginAzInt::Write(const DataMessage &msg, uint64_t image_number) { @@ -33,17 +41,41 @@ void HDF5DataFilePluginAzInt::Write(const DataMessage &msg, uint64_t image_numbe if (image_number >= max_image_number || (max_image_number == 0)) { max_image_number = image_number; - az_int_image.resize((max_image_number + 1) * azimuthal_bins * q_bins); + az_int_image.resize((max_image_number + 1) * azimuthal_bins * q_bins, NAN); + if (save_profile_count) + az_int_image_count.resize((max_image_number + 1) * azimuthal_bins * q_bins, 0); + if (save_profile_std) + az_int_image_std.resize((max_image_number + 1) * azimuthal_bins * q_bins, NAN); } if (!msg.az_int_profile.empty() && (msg.az_int_profile.size() == azimuthal_bins * q_bins)) { for (int i = 0; i < azimuthal_bins * q_bins; i++) az_int_image[image_number * azimuthal_bins * q_bins + i] = msg.az_int_profile[i]; } + + if (save_profile_count + && !msg.az_int_profile_count.empty() + && (msg.az_int_profile_count.size() == azimuthal_bins * q_bins)) { + for (int i = 0; i < azimuthal_bins * q_bins; i++) + az_int_image_count[image_number * azimuthal_bins * q_bins + i] = msg.az_int_profile_count[i]; + } + + if (save_profile_std + && !msg.az_int_profile_std.empty() + && (msg.az_int_profile_std.size() == azimuthal_bins * q_bins)) { + for (int i = 0; i < azimuthal_bins * q_bins; i++) + az_int_image_std[image_number * azimuthal_bins * q_bins + i] = msg.az_int_profile_std[i]; + } } void HDF5DataFilePluginAzInt::WriteFinal(HDF5File &data_file) { if (!az_int_image.empty()) data_file.SaveVector("/entry/azint/image", az_int_image, {(hsize_t) (max_image_number + 1), azimuthal_bins, q_bins}); + if (!az_int_image_std.empty()) + data_file.SaveVector("/entry/azint/image_std", az_int_image_std, + {(hsize_t) (max_image_number + 1), azimuthal_bins, q_bins}); + if (!az_int_image_count.empty()) + data_file.SaveVector("/entry/azint/image_count", az_int_image_count, + {(hsize_t) (max_image_number + 1), azimuthal_bins, q_bins}); } diff --git a/writer/HDF5DataFilePluginAzInt.h b/writer/HDF5DataFilePluginAzInt.h index dc4e3c9a..07b82942 100644 --- a/writer/HDF5DataFilePluginAzInt.h +++ b/writer/HDF5DataFilePluginAzInt.h @@ -14,6 +14,11 @@ class HDF5DataFilePluginAzInt : public HDF5DataFilePlugin { std::vector az_int_bin_to_two_theta; std::vector az_int_bin_to_phi; std::vector az_int_image; + std::vector az_int_image_std; + std::vector az_int_image_count; + + bool save_profile_count = false; + bool save_profile_std = false; // std::unique_ptr dataset; public: