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

RadialIntegrationProfile: Extra routines to handle GPU/CPU/FPGA workflows in more versatile way

Browse Source
This commit is contained in:
leonarski_f
2023-10-21 17:14:17 +02:00
parent a7706546b7
commit b4ab3087f1
13 changed files with 144 additions and 77 deletions
Download Patch File Download Diff File Expand all files Collapse all files
common/DiffractionGeometry.h
+20
View File
@@ -5,6 +5,7 @@
#include "DiffractionExperiment.h"
#include <cmath>
#include "RawToConvertedGeometry.h"
inline Coord DetectorToRecip(const DiffractionExperiment &experiment, float x, float y) {
return experiment.LabCoord(x, y).Normalize() / experiment.GetWavelength_A() - experiment.GetScatteringVector();
@@ -96,4 +97,23 @@ inline std::vector<float> CalcRadIntCorr(const DiffractionExperiment& experiment
return corr;
}
inline std::vector<float> CalcRadIntCorrRawCoord(const DiffractionExperiment& experiment) {
std::vector<float> corr(experiment.GetModulesNum() * RAW_MODULE_SIZE, 1.0);
for (int m = 0; m < experiment.GetModulesNum(); m++) {
for (int i = 0; i < RAW_MODULE_SIZE; i++) {
auto [x,y] = RawToConvertedCoordinate(experiment, m, i);
if (experiment.GetApplySolidAngleCorr())
corr[m * RAW_MODULE_SIZE + i] *= CalcRadIntSolidAngleCorr(experiment,
static_cast<float>(x),
static_cast<float>(y));
if (experiment.GetApplyPolarizationCorr())
corr[m * RAW_MODULE_SIZE + i] *= CalcRadIntPolarizationCorr(experiment,
static_cast<float>(x),
static_cast<float>(y));
}
}
return corr;
}
#endif //JUNGFRAUJOCH_DIFFRACTIONGEOMETRY_H
image_analysis/GPUImageAnalysis.cpp
-4
View File
@@ -55,10 +55,6 @@ std::vector<float> GPUImageAnalysis::GetRadialIntegrationCount() const {
return {};
}
float GPUImageAnalysis::GetRadialIntegrationRangeValue(uint16_t min_bin, uint16_t max_bin) {
return 0;
}
void GPUImageAnalysis::LoadRadialIntegrationCorr(const std::vector<float>& v) {}
#endif
image_analysis/GPUImageAnalysis.cu
-22
View File
@@ -486,25 +486,3 @@ std::vector<float> GPUImageAnalysis::GetRadialIntegrationCount() const {
memcpy(out.data(), host_rad_integration_count, rad_integration_nbins * sizeof(int32_t));
return out;
}
float GPUImageAnalysis::GetRadialIntegrationRangeValue(uint16_t min_bin, uint16_t max_bin) {
if (rad_integration_nbins == 0)
throw JFJochException(JFJochExceptionCategory::SpotFinderError, "Radial integration not initialized");
cuda_err(cudaStreamSynchronize(cudastream->v));
float ret_sum = 0;
float ret_count = 0;
for (int i = std::min(rad_integration_nbins,min_bin);
i <= std::min((uint16_t)(rad_integration_nbins-1),max_bin);
i++) {
ret_sum += host_rad_integration_sum[i];
ret_count += host_rad_integration_count[i];
}
if (ret_count == 0)
return 0;
else
return ret_sum / ret_count;
}
image_analysis/GPUImageAnalysis.h
-1
View File
@@ -72,7 +72,6 @@ public:
void GetRadialIntegrationProfile(std::vector<float> &result);
[[nodiscard]] std::vector<float> GetRadialIntegrationSum() const;
[[nodiscard]] std::vector<float> GetRadialIntegrationCount() const;
[[nodiscard]] float GetRadialIntegrationRangeValue(uint16_t min_bin, uint16_t max_bin);
static bool GPUPresent();
void RegisterBuffer();
image_analysis/RadialIntegrationProfile.cpp
+49
View File
@@ -39,3 +39,52 @@ void RadialIntegrationProfile::GetPlot(JFJochProtoBuf::Plot &plot) const {
*plot.mutable_x() = {bin_to_q.begin(), bin_to_q.end()};
*plot.mutable_y() = {rad_int_profile.begin(), rad_int_profile.end()};
}
float RadialIntegrationProfile::GetMeanValueOfBins(uint16_t min_bin, uint16_t max_bin) const {
std::unique_lock<std::mutex> ul(m);
float ret_sum = 0;
float ret_count = 0;
for (int i = std::min((uint16_t)sum.size(),min_bin);
i <= std::min((uint16_t)(sum.size()-1),max_bin);
i++) {
ret_sum += sum[i];
ret_count += count[i];
}
if (ret_count == 0)
return 0;
else
return ret_sum / ret_count;
}
RadialIntegrationProfile &RadialIntegrationProfile::operator+=(const RadialIntegrationProfile &other) {
std::unique_lock<std::mutex> ul(m);
// It is expected that "other" has exclusive access
if ((other.bin_to_q != bin_to_q) || (sum.size() != other.sum.size())) {
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Error combining two radial integration profiles");
}
for (int i = 0; i < sum.size(); i++) {
sum[i] += other.sum[i];
count[i] += other.count[i];
}
return *this;
}
void RadialIntegrationProfile::Add(const DeviceOutput &result) {
std::unique_lock<std::mutex> ul(m);
if (sum.size() > FPGA_INTEGRATION_BIN_COUNT )
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Error in getting result from FPGA");
for (int i = 0; i < sum.size(); i++) {
sum[i] += static_cast<float>(static_cast<double>(result.integration_result[i].sum) / (1LU << 24));
count[i] += result.integration_result[i].count;
}
}
image_analysis/RadialIntegrationProfile.h
+4
View File
@@ -8,6 +8,7 @@
#include <jfjoch.pb.h>
#include "RadialIntegrationMapping.h"
#include "../receiver/AcquisitionDevice.h"
class RadialIntegrationProfile {
mutable std::mutex m;
@@ -16,9 +17,12 @@ class RadialIntegrationProfile {
std::vector<float> bin_to_q;
public:
explicit RadialIntegrationProfile(RadialIntegrationMapping &mapping, const DiffractionExperiment &experiment);
void Add(const DeviceOutput &result);
void Add(const std::vector<float> &sum, const std::vector<float> &count);
std::vector<float> GetResult() const;
float GetMeanValueOfBins(uint16_t min_bin, uint16_t max_bin) const;
void GetPlot(JFJochProtoBuf::Plot &plot) const;
RadialIntegrationProfile& operator+=(const RadialIntegrationProfile& profile);
};
#endif //JUNGFRAUJOCH_RADIALINTEGRATIONPROFILE_H
receiver/AcquisitionDevice.cpp
+1 -1
View File
@@ -180,7 +180,7 @@ void AcquisitionDevice::InitializeIntegrationMap(const DiffractionExperiment &ex
void AcquisitionDevice::InitializeIntegrationMap(const DiffractionExperiment &experiment,
const std::vector<uint16_t> &v,
const std::vector<double> &weights) {}
const std::vector<float> &weights) {}
void AcquisitionDevice::MapBuffersStandard(size_t c2h_buffer_count, size_t h2c_buffer_count, int16_t numa_node) {
try {
receiver/AcquisitionDevice.h
+1 -1
View File
@@ -106,7 +106,7 @@ public:
virtual void InitializeCalibration(const DiffractionExperiment &experiment, const JFCalibration &calib);
virtual void InitializeIntegrationMap(const DiffractionExperiment &experiment, const std::vector<uint16_t> &v);
virtual void InitializeIntegrationMap(const DiffractionExperiment &experiment, const std::vector<uint16_t> &v,
const std::vector<double> &weights);
const std::vector<float> &weights);
const AcquisitionCounters& Counters() const;
virtual void SetSpotFinderParameters(int16_t count_threshold, double snr_threshold) {};
receiver/FPGAAcquisitionDevice.cpp
+4 -23
View File
@@ -63,32 +63,13 @@ void FPGAAcquisitionDevice::SendWorkRequestThread() {
void FPGAAcquisitionDevice::InitializeIntegrationMap(const DiffractionExperiment &experiment,
const std::vector<uint16_t> &v) {
auto offset = experiment.GetFirstModuleOfDataStream(data_stream);
if (v.size() != experiment.GetModulesNum() * RAW_MODULE_SIZE)
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Mismatch regarding integration map array");
size_t modules = experiment.GetModulesNum(data_stream);
if (modules > 2 * buffer_device.size())
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
"Not enough host/FPGA buffers to load all integration map values");
for (int m = 0; m < modules; m++)
memcpy(buffer_device[m], v.data() + (offset + m) * RAW_MODULE_SIZE, RAW_MODULE_SIZE * sizeof(uint16_t));
for (int m = 0; m < modules; m++) {
for (int i = 0; i < RAW_MODULE_SIZE; i++) {
buffer_device[modules + m][i] = to_fixed(1.0, 15);
}
}
HW_LoadIntegrationMap(modules);
std::vector<float> weights(experiment.GetModulesNum() * RAW_MODULE_SIZE, 1.0);
InitializeIntegrationMap(experiment, v, weights);
}
void FPGAAcquisitionDevice::InitializeIntegrationMap(const DiffractionExperiment &experiment,
const std::vector<uint16_t> &v,
const std::vector<double> &weights) {
const std::vector<float> &weights) {
auto offset = experiment.GetFirstModuleOfDataStream(data_stream);
if (v.size() != experiment.GetModulesNum() * RAW_MODULE_SIZE)
@@ -111,7 +92,7 @@ void FPGAAcquisitionDevice::InitializeIntegrationMap(const DiffractionExperiment
for (int m = 0; m < modules; m++) {
for (int i = 0; i < RAW_MODULE_SIZE; i++) {
buffer_device[modules + m][i] = to_fixed(weights[m * RAW_MODULE_SIZE + i], 15);
buffer_device[modules + m][i] = to_fixed(weights[(offset + m) * RAW_MODULE_SIZE + i], 15);
}
}
HW_LoadIntegrationMap(modules);
receiver/FPGAAcquisitionDevice.h
+1 -1
View File
@@ -47,7 +47,7 @@ public:
void InitializeCalibration(const DiffractionExperiment &experiment, const JFCalibration &calib) override;
void InitializeIntegrationMap(const DiffractionExperiment &experiment, const std::vector<uint16_t> &v) override;
void InitializeIntegrationMap(const DiffractionExperiment &experiment, const std::vector<uint16_t> &v,
const std::vector<double> &weights) override;
const std::vector<float> &weights) override;
void SetInternalGeneratorFrame(const std::vector<uint16_t> &v);
void SetInternalGeneratorFrame();
receiver/JFJochReceiver.cpp
+11 -8
View File
@@ -548,21 +548,24 @@ void JFJochReceiver::FrameTransformationThread() {
}
if (send_bkg_estimate) {
RadialIntegrationProfile rad_int_profile_image(*rad_int_mapping, experiment);
rad_int_profile_image.Add(spot_finder->GetRadialIntegrationSum(),
spot_finder->GetRadialIntegrationCount());
uint16_t rad_int_min_bin = std::floor(
rad_int_mapping->QToBin(local_data_processing_settings.bkg_estimate_low_q()));
uint16_t rad_int_max_bin = std::ceil(
rad_int_mapping->QToBin(local_data_processing_settings.bkg_estimate_high_q()));
float bkg_estimate_val = spot_finder->GetRadialIntegrationRangeValue(rad_int_min_bin, rad_int_max_bin);
bkg_estimate.AddElement(image_number, bkg_estimate_val);
spot_finder->GetRadialIntegrationProfile(message.rad_int_profile);
rad_int_profile->Add(spot_finder->GetRadialIntegrationSum(),
spot_finder->GetRadialIntegrationCount());
float bkg_estimate_val = rad_int_profile_image.GetMeanValueOfBins(rad_int_min_bin, rad_int_max_bin);
bkg_estimate.AddElement(image_number, bkg_estimate_val);
message.rad_int_profile = rad_int_profile_image.GetResult();
*rad_int_profile += rad_int_profile_image;
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());
*rad_int_profile_per_file[image_number % experiment.GetDataFileCount()] += rad_int_profile_image;
}
if (send_preview)
tests/HDF5WritingTest.cpp
-1
View File
@@ -230,7 +230,6 @@ TEST_CASE("HDF5Writer_Rad_Int_Profile", "[HDF5][Full]") {
x.QSpacingForRadialInt_recipA(0.1).LowQForRadialInt_recipA(0.1).HighQForRadialInt_recipA(4);
RadialIntegrationMapping mapping(x);
RadialIntegrationProfile profile(mapping, x);
std::vector<float> rad_int_profile(mapping.GetBinNumber(), 4.0);
std::vector<float> rad_int_avg(mapping.GetBinNumber(), 0.33);
tests/RadialIntegrationTest.cpp
+53 -15
View File
@@ -147,6 +147,7 @@ TEST_CASE("RadialIntegrationProfile","[RadialIntegration]") {
count[i] = i;
}
REQUIRE_NOTHROW(profile.Add(sum, count));
REQUIRE_NOTHROW(profile.Add(sum, count));
std::vector<float> sum_wr(mapping.GetBinNumber() - 1);
REQUIRE_THROWS(profile.Add(sum_wr, count));
@@ -162,6 +163,58 @@ TEST_CASE("RadialIntegrationProfile","[RadialIntegration]") {
}
}
TEST_CASE("RadialIntegrationProfile_operatorAdd","[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 profile0(mapping, x), profile1(mapping, x);
std::vector<float> sum(mapping.GetBinNumber());
std::vector<float> count(mapping.GetBinNumber());
for (int i = 0; i < mapping.GetBinNumber(); i++) {
sum[i] = i * i * 4;
count[i] = i;
}
REQUIRE_NOTHROW(profile0.Add(sum, count));
REQUIRE_NOTHROW(profile1 += profile0);
JFJochProtoBuf::Plot plot;
profile1.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("RadialIntegrationProfile_GetMeanValueOfBins","[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, x);
std::vector<float> sum(mapping.GetBinNumber());
std::vector<float> 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));
REQUIRE(profile.GetMeanValueOfBins(0,2) == Approx ((sum[0] + sum[1] + sum[2]) / (count[0] + count[1] + count[2])));
REQUIRE(profile.GetMeanValueOfBins(5,7) == Approx ((sum[5] + sum[6] + sum[7]) / (count[5] + count[6] + count[7])));
}
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};
@@ -311,19 +364,4 @@ TEST_CASE("RadialIntegrationGPU_GetResult","[RadialIntegration]") {
REQUIRE(result[3] == (3 + 0) / 2.0);
}
TEST_CASE("RadialIntegrationGPU_GetRangeValue","[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};
std::vector<uint8_t> one_byte_mask = {1,1,1,1,1,1,1,1};
GPUImageAnalysis image_analysis(8, 1, one_byte_mask, pixel_to_bin, 4);
image_analysis.SetInputBuffer(test_image.data());
image_analysis.LoadDataToGPU();
image_analysis.RunRadialIntegration();
REQUIRE(image_analysis.GetRadialIntegrationRangeValue(0, 7) == Approx((0+1+2+3+5+6+7) / 7.0));
REQUIRE(image_analysis.GetRadialIntegrationRangeValue(2, 2) == Approx((5+1) / 2.0));
REQUIRE(image_analysis.GetRadialIntegrationRangeValue(2, 3) == Approx((5+3+1+0) / 4.0));
REQUIRE(image_analysis.GetRadialIntegrationRangeValue(15, 15) == 0); // Empty set
}