diff --git a/common/DiffractionGeometry.h b/common/DiffractionGeometry.h index b89b82d1..1005e991 100644 --- a/common/DiffractionGeometry.h +++ b/common/DiffractionGeometry.h @@ -5,6 +5,7 @@ #include "DiffractionExperiment.h" #include +#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 CalcRadIntCorr(const DiffractionExperiment& experiment return corr; } +inline std::vector CalcRadIntCorrRawCoord(const DiffractionExperiment& experiment) { + std::vector 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(x), + static_cast(y)); + if (experiment.GetApplyPolarizationCorr()) + corr[m * RAW_MODULE_SIZE + i] *= CalcRadIntPolarizationCorr(experiment, + static_cast(x), + static_cast(y)); + } + } + + return corr; +} #endif //JUNGFRAUJOCH_DIFFRACTIONGEOMETRY_H diff --git a/image_analysis/GPUImageAnalysis.cpp b/image_analysis/GPUImageAnalysis.cpp index 7f7fd441..e336b79f 100644 --- a/image_analysis/GPUImageAnalysis.cpp +++ b/image_analysis/GPUImageAnalysis.cpp @@ -55,10 +55,6 @@ std::vector GPUImageAnalysis::GetRadialIntegrationCount() const { return {}; } -float GPUImageAnalysis::GetRadialIntegrationRangeValue(uint16_t min_bin, uint16_t max_bin) { - return 0; -} - void GPUImageAnalysis::LoadRadialIntegrationCorr(const std::vector& v) {} #endif diff --git a/image_analysis/GPUImageAnalysis.cu b/image_analysis/GPUImageAnalysis.cu index 45226352..d2dabedf 100644 --- a/image_analysis/GPUImageAnalysis.cu +++ b/image_analysis/GPUImageAnalysis.cu @@ -486,25 +486,3 @@ std::vector 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; -} diff --git a/image_analysis/GPUImageAnalysis.h b/image_analysis/GPUImageAnalysis.h index 683132a4..aea1e19d 100644 --- a/image_analysis/GPUImageAnalysis.h +++ b/image_analysis/GPUImageAnalysis.h @@ -72,7 +72,6 @@ public: void GetRadialIntegrationProfile(std::vector &result); [[nodiscard]] std::vector GetRadialIntegrationSum() const; [[nodiscard]] std::vector GetRadialIntegrationCount() const; - [[nodiscard]] float GetRadialIntegrationRangeValue(uint16_t min_bin, uint16_t max_bin); static bool GPUPresent(); void RegisterBuffer(); diff --git a/image_analysis/RadialIntegrationProfile.cpp b/image_analysis/RadialIntegrationProfile.cpp index 599ead54..282a931a 100644 --- a/image_analysis/RadialIntegrationProfile.cpp +++ b/image_analysis/RadialIntegrationProfile.cpp @@ -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 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 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 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(static_cast(result.integration_result[i].sum) / (1LU << 24)); + count[i] += result.integration_result[i].count; + } +} \ No newline at end of file diff --git a/image_analysis/RadialIntegrationProfile.h b/image_analysis/RadialIntegrationProfile.h index 75dcec40..fc9d4326 100644 --- a/image_analysis/RadialIntegrationProfile.h +++ b/image_analysis/RadialIntegrationProfile.h @@ -8,6 +8,7 @@ #include #include "RadialIntegrationMapping.h" +#include "../receiver/AcquisitionDevice.h" class RadialIntegrationProfile { mutable std::mutex m; @@ -16,9 +17,12 @@ class RadialIntegrationProfile { std::vector bin_to_q; public: explicit RadialIntegrationProfile(RadialIntegrationMapping &mapping, const DiffractionExperiment &experiment); + void Add(const DeviceOutput &result); void Add(const std::vector &sum, const std::vector &count); std::vector 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 diff --git a/receiver/AcquisitionDevice.cpp b/receiver/AcquisitionDevice.cpp index 56e44b35..dcc4a284 100644 --- a/receiver/AcquisitionDevice.cpp +++ b/receiver/AcquisitionDevice.cpp @@ -180,7 +180,7 @@ void AcquisitionDevice::InitializeIntegrationMap(const DiffractionExperiment &ex void AcquisitionDevice::InitializeIntegrationMap(const DiffractionExperiment &experiment, const std::vector &v, - const std::vector &weights) {} + const std::vector &weights) {} void AcquisitionDevice::MapBuffersStandard(size_t c2h_buffer_count, size_t h2c_buffer_count, int16_t numa_node) { try { diff --git a/receiver/AcquisitionDevice.h b/receiver/AcquisitionDevice.h index ff24925a..27d68275 100644 --- a/receiver/AcquisitionDevice.h +++ b/receiver/AcquisitionDevice.h @@ -106,7 +106,7 @@ public: virtual void InitializeCalibration(const DiffractionExperiment &experiment, const JFCalibration &calib); virtual void InitializeIntegrationMap(const DiffractionExperiment &experiment, const std::vector &v); virtual void InitializeIntegrationMap(const DiffractionExperiment &experiment, const std::vector &v, - const std::vector &weights); + const std::vector &weights); const AcquisitionCounters& Counters() const; virtual void SetSpotFinderParameters(int16_t count_threshold, double snr_threshold) {}; diff --git a/receiver/FPGAAcquisitionDevice.cpp b/receiver/FPGAAcquisitionDevice.cpp index 8631deea..2e97c673 100644 --- a/receiver/FPGAAcquisitionDevice.cpp +++ b/receiver/FPGAAcquisitionDevice.cpp @@ -63,32 +63,13 @@ void FPGAAcquisitionDevice::SendWorkRequestThread() { void FPGAAcquisitionDevice::InitializeIntegrationMap(const DiffractionExperiment &experiment, const std::vector &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 weights(experiment.GetModulesNum() * RAW_MODULE_SIZE, 1.0); + InitializeIntegrationMap(experiment, v, weights); } void FPGAAcquisitionDevice::InitializeIntegrationMap(const DiffractionExperiment &experiment, const std::vector &v, - const std::vector &weights) { + const std::vector &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); diff --git a/receiver/FPGAAcquisitionDevice.h b/receiver/FPGAAcquisitionDevice.h index ec2dd06e..21af4bd0 100644 --- a/receiver/FPGAAcquisitionDevice.h +++ b/receiver/FPGAAcquisitionDevice.h @@ -47,7 +47,7 @@ public: void InitializeCalibration(const DiffractionExperiment &experiment, const JFCalibration &calib) override; void InitializeIntegrationMap(const DiffractionExperiment &experiment, const std::vector &v) override; void InitializeIntegrationMap(const DiffractionExperiment &experiment, const std::vector &v, - const std::vector &weights) override; + const std::vector &weights) override; void SetInternalGeneratorFrame(const std::vector &v); void SetInternalGeneratorFrame(); diff --git a/receiver/JFJochReceiver.cpp b/receiver/JFJochReceiver.cpp index 2cd4a811..d8179130 100644 --- a/receiver/JFJochReceiver.cpp +++ b/receiver/JFJochReceiver.cpp @@ -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) diff --git a/tests/HDF5WritingTest.cpp b/tests/HDF5WritingTest.cpp index 7198119b..9b8111f1 100644 --- a/tests/HDF5WritingTest.cpp +++ b/tests/HDF5WritingTest.cpp @@ -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 rad_int_profile(mapping.GetBinNumber(), 4.0); std::vector rad_int_avg(mapping.GetBinNumber(), 0.33); diff --git a/tests/RadialIntegrationTest.cpp b/tests/RadialIntegrationTest.cpp index e8453deb..42f69480 100644 --- a/tests/RadialIntegrationTest.cpp +++ b/tests/RadialIntegrationTest.cpp @@ -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 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 sum(mapping.GetBinNumber()); + std::vector 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 sum(mapping.GetBinNumber()); + std::vector 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 pixel_to_bin = {0,1,2,4,3,1,2,3}; std::vector 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 pixel_to_bin = {0,1,2,4,3,1,2,3}; - std::vector test_image = {7,6,5,4,3,2,1,0}; - std::vector 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 -}