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

Bragg integration: option to use azimuthal integration profile

Browse Source
This commit is contained in:
leonarski_f
2026-06-10 14:41:36 +02:00
parent e051eed033
commit e4230bc14e
15 changed files with 218 additions and 115 deletions
Download Patch File Download Diff File Expand all files Collapse all files
common/BraggIntegrationSettings.cpp
+9
View File
@@ -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;
}
common/BraggIntegrationSettings.h
+3 -2
View File
@@ -12,14 +12,14 @@ class BraggIntegrationSettings {
float d_min_limit_A = 1.0;
std::optional<float> 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<float> 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;
};
image_analysis/IndexAndRefine.cpp
+7 -7
View File
@@ -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<float>(integration_end_time - integration_start_time).count();
@@ -314,7 +314,7 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
msg.bragg_prediction_time_s = std::chrono::duration<float>(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<float>(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;
image_analysis/IndexAndRefine.h
+4 -4
View File
@@ -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<ScaleOnTheFly> 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<MergedReflection> &reference);
ScalingResult ScaleAllImages(const std::vector<MergedReflection> &reference, size_t nthreads = 0);
image_analysis/MXAnalysisAfterFPGA.cpp
+7 -4
View File
@@ -25,10 +25,12 @@ double stddev(const std::vector<float> &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) {
image_analysis/MXAnalysisAfterFPGA.h
+3 -2
View File
@@ -28,9 +28,10 @@ class MXAnalysisAfterFPGA {
std::chrono::duration<double, std::micro> 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);
};
image_analysis/MXAnalysisWithoutFPGA.cpp
+1 -1
View File
@@ -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;
image_analysis/bragg_integration/BraggIntegrate2D.cpp
+157 -62
View File
@@ -8,65 +8,62 @@
#include <vector>
namespace {
template<class T>
float Median(std::vector<T> &values) {
if (values.empty())
return 0.0f;
template<class T>
float Median(std::vector<T> &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<float>(values[middle]);
if (values.size() % 2 == 1)
return static_cast<float>(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<float>(lower + upper);
}
return 0.5f * static_cast<float>(lower + upper);
}
void MarkReflectionMask(std::vector<uint8_t> &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<uint8_t> &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<int64_t>(xpixel))
x1 = static_cast<int64_t>(xpixel) - 1;
if (y1 >= static_cast<int64_t>(ypixel))
y1 = static_cast<int64_t>(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<int64_t>(xpixel))
x1 = static_cast<int64_t>(xpixel) - 1;
if (y1 >= static_cast<int64_t>(ypixel))
y1 = static_cast<int64_t>(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<uint8_t> BuildReflectionMask(const std::vector<Reflection> &predicted,
size_t npredicted,
size_t xpixel, size_t ypixel,
float r_2, float r_2_sq) {
std::vector<uint8_t> mask(xpixel * ypixel, 0);
std::vector<uint8_t> BuildReflectionMask(const std::vector<Reflection> &predicted,
size_t npredicted,
size_t xpixel, size_t ypixel,
float r_2, float r_2_sq) {
std::vector<uint8_t> 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<class T>
@@ -75,7 +72,6 @@ void IntegrateReflection(Reflection &r, const T *image, const std::vector<uint8_
int64_t special_value, int64_t saturation,
float r_3, float r_1_sq, float r_2_sq, float r_3_sq,
float minimum_sigma_in_regards_to_i) {
int64_t x0 = std::floor(r.predicted_x - r_3 - 1.0);
int64_t x1 = std::ceil(r.predicted_x + r_3 + 1.0);
int64_t y0 = std::floor(r.predicted_y - r_3 - 1.0);
@@ -146,15 +142,94 @@ void IntegrateReflection(Reflection &r, const T *image, const std::vector<uint8_
}
}
template<class T>
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<uint16_t> &pixel_to_bin,
const std::vector<float> &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<int64_t>(xpixel - 1), x1);
y0 = std::max(0L, y0);
y1 = std::min(static_cast<int64_t>(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<T> bkg_values;
bkg_values.reserve(static_cast<size_t>((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<float>(I_npixel_integrated);
r.I = static_cast<float>(I_sum) - bkg;
if (I_sum > 0) {
r.observed_x = static_cast<float>(I_sum_x) / static_cast<float>(I_sum);
r.observed_y = static_cast<float>(I_sum_y) / static_cast<float>(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<float>(I_sum)));
r.observed = true;
} else {
r.I = 0;
r.bkg = 0;
r.sigma = NAN;
r.observed = false;
}
}
template<class T>
std::vector<Reflection> IntegrateInternal(const DiffractionExperiment &experiment,
const CompressedImage &image,
const AzimuthalIntegrationMapping &mapping,
const AzimuthalIntegrationProfile &profile,
const std::vector<Reflection> &predicted,
size_t npredicted,
int64_t special_value,
int64_t saturation,
int64_t image_number) {
std::vector<Reflection> ret;
ret.reserve(npredicted);
@@ -164,24 +239,36 @@ std::vector<Reflection> IntegrateInternal(const DiffractionExperiment &experimen
std::vector<uint8_t> buffer;
auto ptr = reinterpret_cast<const T *>(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<float>(image_number);
ret.emplace_back(r);
@@ -193,6 +280,8 @@ std::vector<Reflection> IntegrateInternal(const DiffractionExperiment &experimen
std::vector<Reflection> BraggIntegrate2D(const DiffractionExperiment &experiment,
const CompressedImage &image,
const std::vector<Reflection> &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<Reflection> BraggIntegrate2D(const DiffractionExperiment &experiment
switch (image.GetMode()) {
case CompressedImageMode::Int8:
return IntegrateInternal<int8_t>(experiment, image, predicted, npredicted, INT8_MIN, INT8_MAX, image_number);
return IntegrateInternal<int8_t>(experiment, image, mapping, profile, predicted, npredicted,
INT8_MIN, INT8_MAX, image_number);
case CompressedImageMode::Int16:
return IntegrateInternal<int16_t>(experiment, image, predicted, npredicted, INT16_MIN, INT16_MAX, image_number);
return IntegrateInternal<int16_t>(experiment, image, mapping, profile, predicted, npredicted,
INT16_MIN, INT16_MAX, image_number);
case CompressedImageMode::Int32:
return IntegrateInternal<int32_t>(experiment, image, predicted, npredicted, INT32_MIN, INT32_MAX, image_number);
return IntegrateInternal<int32_t>(experiment, image, mapping, profile, predicted, npredicted,
INT32_MIN, INT32_MAX, image_number);
case CompressedImageMode::Uint8:
return IntegrateInternal<uint8_t>(experiment, image, predicted, npredicted, UINT8_MAX, UINT8_MAX, image_number);
return IntegrateInternal<uint8_t>(experiment, image, mapping, profile, predicted, npredicted,
UINT8_MAX, UINT8_MAX, image_number);
case CompressedImageMode::Uint16:
return IntegrateInternal<uint16_t>(experiment, image, predicted, npredicted, UINT16_MAX, UINT16_MAX, image_number);
return IntegrateInternal<uint16_t>(experiment, image, mapping, profile, predicted, npredicted,
UINT16_MAX, UINT16_MAX, image_number);
case CompressedImageMode::Uint32:
return IntegrateInternal<uint32_t>(experiment, image, predicted, npredicted, UINT32_MAX, UINT32_MAX, image_number);
return IntegrateInternal<uint32_t>(experiment, image, mapping, profile, predicted, npredicted,
UINT32_MAX, UINT32_MAX, image_number);
default:
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Image mode not supported");
image_analysis/bragg_integration/BraggIntegrate2D.h
+5 -5
View File
@@ -1,17 +1,17 @@
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#ifndef JFJOCH_BRAGGINTEGRATE2D_H
#define JFJOCH_BRAGGINTEGRATE2D_H
#pragma once
#include <vector>
#include "../../common/DiffractionExperiment.h"
#include "../../common/Reflection.h"
#include "../../common/AzimuthalIntegrationProfile.h"
std::vector<Reflection> BraggIntegrate2D(const DiffractionExperiment &experiment,
const CompressedImage &image,
const std::vector<Reflection> &predicted,
const AzimuthalIntegrationMapping &mapping,
const AzimuthalIntegrationProfile &profile,
size_t npredicted,
int64_t image_number);
#endif //JFJOCH_BRAGGINTEGRATE2D_H
int64_t image_number);
receiver/JFJochReceiver.cpp
+10 -15
View File
@@ -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<AzimuthalIntegrationMapping>(experiment, pixel_mask);
auto end_time_point = std::chrono::steady_clock::now();
auto duration = std::chrono::duration<float>(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();
receiver/JFJochReceiver.h
+1 -1
View File
@@ -80,7 +80,7 @@ protected:
std::vector<std::unique_ptr<ADUHistogram>> adu_histogram_module;
PixelMask pixel_mask;
std::unique_ptr<AzimuthalIntegrationMapping> az_int_mapping;
AzimuthalIntegrationMapping az_int_mapping;
std::optional<uint64_t> max_delay;
std::mutex max_delay_mutex;
receiver/JFJochReceiverFPGA.cpp
+4 -7
View File
@@ -297,7 +297,7 @@ void JFJochReceiverFPGA::FrameTransformationThread(uint32_t threadid) {
try {
numa_policy.Bind(threadid);
analyzer = std::make_unique<MXAnalysisAfterFPGA>(experiment, indexer);
analyzer = std::make_unique<MXAnalysisAfterFPGA>(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);
}
receiver/JFJochReceiverLite.cpp
+2 -2
View File
@@ -242,7 +242,7 @@ void JFJochReceiverLite::DataAnalysisThread(uint32_t id) {
measurement_started.wait();
try {
analysis = std::make_unique<MXAnalysisWithoutFPGA>(experiment, *az_int_mapping, pixel_mask, indexer);
analysis = std::make_unique<MXAnalysisWithoutFPGA>(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();
tests/BraggIntegrate2DTest.cpp
+5 -1
View File
@@ -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<int>(height); y++) {
for (int x = 0; x < static_cast<int>(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<Reflection> 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);
tools/jfjoch_scale.cpp
-2
View File
@@ -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();