mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests 1 Actions 2 Packages Projects Releases 55 Wiki Activity

AzimuthalIntegrationProfile: Add standard deviation + save pixel count + improve handling of special cases (missing image)

Browse Source
This commit is contained in:
leonarski_f
2026-06-02 15:50:11 +02:00
parent 75de40f52b
commit 953089009f
13 changed files with 175 additions and 16 deletions
Download Patch File Download Diff File Expand all files Collapse all files
common/AzimuthalIntegrationProfile.cpp
+32 -3
View File
@@ -7,8 +7,15 @@
inline float sum_to_count(float sum, uint64_t count) {
if (count == 0)
return NAN;
else
return sum / (static_cast<float>(count));
return sum / (static_cast<float>(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<float>(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<float> &in_sum, const std::vector<uint32_t> &in_count) {
void AzimuthalIntegrationProfile::Add(const std::vector<float> &in_sum,
const std::vector<float> &in_sum2,
const std::vector<uint32_t> &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<float> AzimuthalIntegrationProfile::GetResult() const {
return rad_int_profile;
}
std::vector<float> AzimuthalIntegrationProfile::GetStd() const {
std::unique_lock ul(m);
std::vector<float> 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<uint64_t> AzimuthalIntegrationProfile::GetPixelCount() const {
std::unique_lock ul(m);
return count;
}
std::vector<float> 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];
}
common/AzimuthalIntegrationProfile.h
+6 -1
View File
@@ -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<float> &sum, const std::vector<uint32_t> &count);
void Add(const std::vector<float> &sum,
const std::vector<float> &sum2,
const std::vector<uint32_t> &count);
void Add(int64_t bin, int64_t value);
std::vector<float> GetResult() const;
std::vector<float> GetStd() const;
std::vector<uint64_t> GetPixelCount() const;
std::vector<float> GetResult1D() const;
float GetMeanValueOfBins(uint16_t min_bin, uint16_t max_bin) const;
common/JFJochMessages.h
+3
View File
@@ -107,6 +107,9 @@ struct DataMessage {
std::vector<float> spot_plot_one_over_d_square;
std::vector<float> az_int_profile;
std::vector<float> az_int_profile_std;
std::vector<uint64_t> az_int_profile_count;
std::optional<float> bkg_estimate;
std::optional<bool> indexing_result;
frame_serialize/CBORStream2Deserializer.cpp
+4
View File
@@ -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") {
frame_serialize/CBORStream2Serializer.cpp
+4
View File
@@ -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());
image_analysis/MXAnalysisWithoutFPGA.cpp
+3
View File
@@ -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());
}
image_analysis/azint/AzIntEngineCPU.cpp
+1 -1
View File
@@ -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);
}
image_analysis/azint/AzIntEngineGPU.cu
+1 -1
View File
@@ -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);
}
receiver/JFJochReceiverFPGA.cpp
+1
View File
@@ -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);
tests/AzimuthalIntegrationTest.cpp
+70 -8
View File
@@ -171,18 +171,19 @@ TEST_CASE("AzimuthalIntegrationProfile","[AzimuthalIntegration]") {
AzimuthalIntegrationProfile profile(mapping);
std::vector<float> sum(mapping.GetBinNumber());
std::vector<float> sum(mapping.GetBinNumber()), sum2(mapping.GetBinNumber());
std::vector<uint32_t> 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<float> 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<float> sum(mapping.GetBinNumber()), sum2(mapping.GetBinNumber());
std::vector<uint32_t> 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<float> sum(mapping.GetBinNumber());
std::vector<float> sum(mapping.GetBinNumber()), sum2(mapping.GetBinNumber());
std::vector<uint32_t> 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<float> sum(mapping.GetBinNumber());
std::vector<float> sum2(mapping.GetBinNumber());
std::vector<uint32_t> 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<float> sum(mapping.GetBinNumber(), 0.0f);
std::vector<float> sum2(mapping.GetBinNumber(), 20.0f);
std::vector<uint32_t> 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();
tests/CBORTest.cpp
+12 -1
View File
@@ -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]") {
writer/HDF5DataFilePluginAzInt.cpp
+33 -1
View File
@@ -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});
}
writer/HDF5DataFilePluginAzInt.h
+5
View File
@@ -14,6 +14,11 @@ class HDF5DataFilePluginAzInt : public HDF5DataFilePlugin {
std::vector<float> az_int_bin_to_two_theta;
std::vector<float> az_int_bin_to_phi;
std::vector<float> az_int_image;
std::vector<float> az_int_image_std;
std::vector<uint64_t> az_int_image_count;
bool save_profile_count = false;
bool save_profile_std = false;
// std::unique_ptr<HDF5DataSet> dataset;
public: