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

JFJochReceiver: Per-data-file radial integration profile

Browse Source
This commit is contained in:
leonarski_f
2023-05-06 10:28:45 +02:00
parent 59fbeab176
commit 0d2fa850a2
8 changed files with 190 additions and 104 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
grpc/jfjoch.proto
View File

@@ -270,6 +270,7 @@ enum PlotType {
message PlotRequest {
PlotType type = 1;
uint64 binning = 2;
optional uint64 file_number = 3;
}
message ReceiverDataProcessingPlots {

2
image_analysis/CMakeLists.txt
View File

@@ -4,7 +4,7 @@ ADD_LIBRARY(ImageAnalysis STATIC
GPUImageAnalysis.h
RadialIntegration.cpp RadialIntegration.h
RadialIntegrationMapping.cpp RadialIntegrationMapping.h
StrongPixelSet.cpp StrongPixelSet.h GPUImageAnalysis_Alt.cpp)
StrongPixelSet.cpp StrongPixelSet.h GPUImageAnalysis_Alt.cpp RadialIntegrationProfile.cpp RadialIntegrationProfile.h)
TARGET_LINK_LIBRARIES(ImageAnalysis CommonFunctions)

37
image_analysis/RadialIntegrationProfile.cpp Normal file
View File

@@ -0,0 +1,37 @@
// Copyright (2019-2023) Paul Scherrer Institute
// SPDX-License-Identifier: GPL-3.0-or-later
#include "RadialIntegrationProfile.h"
#include "../common/JFJochException.h"
RadialIntegrationProfile::RadialIntegrationProfile(RadialIntegrationMapping &mapping)
: bin_to_q(mapping.GetBinToQ()),
sum(mapping.GetBinNumber(), 0),
count(mapping.GetBinNumber(), 0) {
}
void RadialIntegrationProfile::Add(const std::vector<int32_t> &in_sum, const std::vector<int32_t> &in_count) {
std::unique_lock<std::mutex> 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];
}
} else if (!in_sum.empty() && !in_count.empty())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Mismatch in size of sum/count datasets");
}
void RadialIntegrationProfile::GetPlot(JFJochProtoBuf::Plot &plot) const {
std::unique_lock<std::mutex> ul(m);
std::vector<float> rad_int_profile(sum.size());
for (int i = 0; i < sum.size(); i++) {
if (count[i] > 0)
rad_int_profile[i] = static_cast<float>(sum[i]) / static_cast<float>(count[i]);
}
*plot.mutable_x() = {bin_to_q.begin(), bin_to_q.end()};
*plot.mutable_y() = {rad_int_profile.begin(), rad_int_profile.end()};
}

24
image_analysis/RadialIntegrationProfile.h Normal file
View File

@@ -0,0 +1,24 @@
// Copyright (2019-2023) Paul Scherrer Institute
// SPDX-License-Identifier: GPL-3.0-or-later
#ifndef JUNGFRAUJOCH_RADIALINTEGRATIONPROFILE_H
#define JUNGFRAUJOCH_RADIALINTEGRATIONPROFILE_H
#include <vector>
#include <mutex>
#include <jfjoch.pb.h>
#include "RadialIntegrationMapping.h"
class RadialIntegrationProfile {
mutable std::mutex m;
std::vector<int64_t> sum;
std::vector<int64_t> count;
std::vector<float> bin_to_q;
public:
explicit RadialIntegrationProfile(RadialIntegrationMapping &mapping);
void Add(const std::vector<int32_t> &sum, const std::vector<int32_t> &count);
void GetPlot(JFJochProtoBuf::Plot &plot) const;
};
#endif //JUNGFRAUJOCH_RADIALINTEGRATIONPROFILE_H

136
python/jfjoch_pb2.py
View File

File diff suppressed because one or more lines are too long

45
receiver/JFJochReceiver.cpp
View File

@@ -35,7 +35,6 @@ JFJochReceiver::JFJochReceiver(const JFJochProtoBuf::ReceiverInput &settings,
frame_transformation_nthreads((experiment.GetSummation() >= threaded_summation_threshold) ?
2 : in_forward_and_sum_nthreads),
preview_publisher(in_preview_publisher),
rad_int_profile_window(experiment.GetSpotFindingBin()),
ndatastreams(experiment.GetDataStreamsNum()),
data_acquisition_ready(ndatastreams),
frame_transformation_ready((experiment.GetImageNum() > 0) ? frame_transformation_nthreads : 0)
@@ -76,6 +75,11 @@ JFJochReceiver::JFJochReceiver(const JFJochProtoBuf::ReceiverInput &settings,
logger.Info("GPU support missing");
rad_int_mapping = std::make_unique<RadialIntegrationMapping>(experiment, one_byte_mask.data());
rad_int_profile = std::make_unique<RadialIntegrationProfile>(*rad_int_mapping);
for (int i = 0; i < experiment.GetDataFileCount(); i++)
rad_int_profile_per_file.emplace_back(std::make_unique<RadialIntegrationProfile>(*rad_int_mapping));
spot_finder_mask = calib->CalculateOneByteMask(experiment);
}
@@ -446,7 +450,13 @@ int64_t JFJochReceiver::FrameTransformationThread() {
bkg_estimate.AddElement(image_number, bkg_estimate_val);
spot_finder->GetRadialIntegrationProfile(message.rad_int_profile);
AddRadialIntegrationProfile(message.rad_int_profile);
rad_int_profile->Add(spot_finder->GetRadialIntegrationSum(),
spot_finder->GetRadialIntegrationCount());
if (image_number % experiment.GetDataFileCount() < rad_int_profile_per_file.size())
rad_int_profile_per_file[image_number % experiment.GetDataFileCount()]
->Add(spot_finder->GetRadialIntegrationSum(),
spot_finder->GetRadialIntegrationCount());
}
if (push_images_to_writer) {
@@ -636,29 +646,6 @@ JFJochProtoBuf::DataProcessingSettings JFJochReceiver::GetDataProcessingSettings
return data_processing_settings;
}
void JFJochReceiver::AddRadialIntegrationProfile(const std::vector<float> &result) {
std::unique_lock<std::mutex> ul(rad_int_profile_mutex);
if (rad_int_profile.empty())
rad_int_profile = result;
else if (rad_int_profile.size() == result.size()) {
for (int i = 0; i < rad_int_profile.size(); i++)
rad_int_profile[i] += (result[i] - rad_int_profile[i]) / rad_int_profile_window;
} else {
// Throw exception?
Abort();
}
}
void JFJochReceiver::GetRadialIntegrationProfile(JFJochProtoBuf::Plot &plot) {
std::unique_lock<std::mutex> ul(rad_int_profile_mutex);
const auto &bin_to_q = rad_int_mapping->GetBinToQ();
if (!rad_int_profile.empty()) {
*plot.mutable_x() = {bin_to_q.begin(), bin_to_q.end()};
*plot.mutable_y() = {rad_int_profile.begin(), rad_int_profile.end()};
}
}
JFJochProtoBuf::Plot JFJochReceiver::GetPlots(const JFJochProtoBuf::PlotRequest &request) {
JFJochProtoBuf::Plot ret;
auto nbins = experiment.GetSpotFindingBin();
@@ -667,7 +654,13 @@ JFJochProtoBuf::Plot JFJochReceiver::GetPlots(const JFJochProtoBuf::PlotRequest
switch (request.type()) {
case JFJochProtoBuf::RAD_INT:
GetRadialIntegrationProfile(ret);
if (request.has_file_number()) {
if (request.file_number() < rad_int_profile_per_file.size())
rad_int_profile_per_file[request.file_number()]->GetPlot(ret);
} else {
if (rad_int_profile)
rad_int_profile->GetPlot(ret);
}
break;
case JFJochProtoBuf::SPOT_COUNT:
spot_count.GetPlot(ret, nbins);

14
receiver/JFJochReceiver.h
View File

@@ -23,7 +23,8 @@
#include "../common/ThreadSafeFIFO.h"
#include "../common/ZMQPreviewPublisher.h"
#include "../image_analysis/RadialIntegration.h"
#include "../image_analysis/RadialIntegrationMapping.h"
#include "../image_analysis/RadialIntegrationProfile.h"
#include "../common/StatusVector.h"
#include "../jungfrau/JFConversionFixedPoint.h"
@@ -41,10 +42,8 @@ class JFJochReceiver {
std::vector<std::future<int64_t>> data_acquisition_futures;
std::unique_ptr<RadialIntegrationMapping> rad_int_mapping;
std::vector<float> rad_int_profile;
const size_t rad_int_profile_window;
std::mutex rad_int_profile_mutex;
std::unique_ptr<RadialIntegrationProfile> rad_int_profile;
std::vector<std::unique_ptr<RadialIntegrationProfile>> rad_int_profile_per_file;
std::vector<uint8_t> spot_finder_mask;
@@ -99,10 +98,7 @@ class JFJochReceiver {
FrameTransformation &transformation, DataMessage &message);
void Abort(const JFJochException &e);
void FinalizeMeasurement();
JFJochProtoBuf::DataProcessingSettings GetDataProcessingSettings();
void AddRadialIntegrationProfile(const std::vector<float> &result);
void GetRadialIntegrationProfile(JFJochProtoBuf::Plot &plots);
JFJochProtoBuf::DataProcessingSettings GetDataProcessingSettings();
void UpdateMaxImage(int64_t image_number);
public:

35
tests/RadialIntegrationTest.cpp
View File

@@ -2,6 +2,9 @@
// SPDX-License-Identifier: GPL-3.0-or-later
#include <catch2/catch.hpp>
#include "../image_analysis/RadialIntegrationProfile.h"
#include "../image_analysis/RadialIntegrationMapping.h"
#include "../image_analysis/RadialIntegration.h"
TEST_CASE("RadialIntegrationMapping_Constructor","[RadialIntegration]") {
@@ -75,6 +78,38 @@ TEST_CASE("RadialIntegrationMapping_QToBin","[RadialIntegration]") {
REQUIRE(mapping.QToBin(50.0) == Approx(38));
}
TEST_CASE("RadialIntegrationProfile","[RadialIntegration]") {
DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36));
x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
x.QSpacingForRadialInt_recipA(0.1).LowQForRadialInt_recipA(0.1).HighQForRadialInt_recipA(4);
RadialIntegrationMapping mapping(x);
RadialIntegrationProfile profile(mapping);
std::vector<int32_t> sum(mapping.GetBinNumber());
std::vector<int32_t> count(mapping.GetBinNumber());
for (int i = 0; i < mapping.GetBinNumber(); i++) {
sum[i] = i * i * 4;
count[i] = i;
}
REQUIRE_NOTHROW(profile.Add(sum, count));
std::vector<int32_t> sum_wr(mapping.GetBinNumber() - 1);
REQUIRE_THROWS(profile.Add(sum_wr, count));
JFJochProtoBuf::Plot plot;
profile.GetPlot(plot);
REQUIRE(plot.x_size() == mapping.GetBinNumber());
REQUIRE(plot.y_size() == mapping.GetBinNumber());
for (int i = 0; i < mapping.GetBinNumber(); i++) {
REQUIRE(plot.x(i) == Approx(mapping.GetBinToQ()[i]));
REQUIRE(plot.y(i) == Approx(i * 4));
}
}
TEST_CASE("RadialIntegration_Process","[RadialIntegration]") {
std::vector<uint16_t> pixel_to_bin = {0,1,2,4,3,1,2,3};
std::vector<int16_t> test_image = {7,6,5,4,3,2,1,0};