From 0b5bbec1fc36bce9ff2e5bcf6bd416875234bea5 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 20 Oct 2023 18:00:29 +0200 Subject: [PATCH] AcquisitionDevice: Setup rad. int. mapping automatically --- image_analysis/CMakeLists.txt | 1 + .../RadialIntegrationMappingFPGA.cpp | 78 +++++++++++++++++++ image_analysis/RadialIntegrationMappingFPGA.h | 25 ++++++ receiver/AcquisitionDevice.cpp | 6 ++ receiver/AcquisitionDevice.h | 1 + tests/RadialIntegrationTest.cpp | 55 ++++++++++++- 6 files changed, 165 insertions(+), 1 deletion(-) create mode 100644 image_analysis/RadialIntegrationMappingFPGA.cpp create mode 100644 image_analysis/RadialIntegrationMappingFPGA.h diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index 0977a096..862094d5 100644 --- a/image_analysis/CMakeLists.txt +++ b/image_analysis/CMakeLists.txt @@ -3,6 +3,7 @@ ADD_LIBRARY(ImageAnalysis STATIC IndexerWrapper.cpp IndexerWrapper.h GPUImageAnalysis.h RadialIntegration.cpp RadialIntegration.h + RadialIntegrationMappingFPGA.cpp RadialIntegrationMappingFPGA.h RadialIntegrationMapping.cpp RadialIntegrationMapping.h StrongPixelSet.cpp StrongPixelSet.h GPUImageAnalysis.cpp RadialIntegrationProfile.cpp RadialIntegrationProfile.h PredictSpotsOnDetector.h) diff --git a/image_analysis/RadialIntegrationMappingFPGA.cpp b/image_analysis/RadialIntegrationMappingFPGA.cpp new file mode 100644 index 00000000..3a341677 --- /dev/null +++ b/image_analysis/RadialIntegrationMappingFPGA.cpp @@ -0,0 +1,78 @@ +// Copyright (2019-2023) Paul Scherrer Institute + +#include +#include "RadialIntegrationMappingFPGA.h" +#include "../common/RawToConvertedGeometry.h" +#include "../common/DiffractionGeometry.h" + +RadialIntegrationMappingFPGA::RadialIntegrationMappingFPGA(const DiffractionExperiment &experiment, uint16_t data_stream) : + low_q(experiment.GetLowQForRadialInt_recipA()), + high_q(experiment.GetHighQForRadialInt_recipA()), + q_spacing(experiment.GetQSpacingForRadialInt_recipA()), + pixel_to_bin(experiment.GetModulesNum(data_stream) * RAW_MODULE_SIZE), + max_bin_number(0) +{ + + if (q_spacing < 0.0) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Q-spacing has to be >= 0.0"); + + if (low_q >= high_q) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Low-Q must be smaller than high-Q"); + + if (low_q <= 0) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Low-Q must be > 0"); + + if (std::ceil((high_q-low_q) / q_spacing ) >= FPGA_INTEGRATION_BIN_COUNT) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Cannot accommodate more than " + std::to_string(FPGA_INTEGRATION_BIN_COUNT) + " rings"); + + for (int m = 0; m < experiment.GetModulesNum(data_stream); m++) { + for (int pxl = 0; pxl < RAW_MODULE_SIZE; pxl++) { + auto [x,y] = RawToConvertedCoordinate(experiment, m + experiment.GetFirstModuleOfDataStream(data_stream), pxl); + double pixel_q = 2 * M_PI / PxlToRes(experiment, x, y); + if ((pixel_q < low_q) || (pixel_q >= high_q)) + pixel_to_bin[m * RAW_MODULE_SIZE + pxl] = UINT16_MAX; + else + pixel_to_bin[m * RAW_MODULE_SIZE + pxl] = std::floor((pixel_q - low_q) / q_spacing); + } + } + + // In principle max bin number is equal to std::floor((high_q - low_q) / q_spacing), but it might be lower than this + // depending on detector distance and beam center position + for (const auto &i: pixel_to_bin) { + if ((i != UINT16_MAX) && (i > max_bin_number)) + max_bin_number = i; + } + + bin_to_q.resize(max_bin_number + 1); + for (int i = 0; i < max_bin_number + 1; i++) + bin_to_q[i] = static_cast(q_spacing * (i + 0.5) + low_q); + + solid_angle_corr.resize(max_bin_number + 1); + for (int i = 0; i < max_bin_number + 1; i++) + solid_angle_corr[i] = 1.0f / CalcRadIntSolidAngleCorr(experiment, bin_to_q[i]); + +} + +uint16_t RadialIntegrationMappingFPGA::GetBinNumber() const { + return max_bin_number + 1; +} + +const std::vector &RadialIntegrationMappingFPGA::GetPixelToBinMapping() const { + return pixel_to_bin; +} + +const std::vector &RadialIntegrationMappingFPGA::GetBinToQ() const { + return bin_to_q; +} + +const std::vector &RadialIntegrationMappingFPGA::GetSolidAngleCorr() const { + return solid_angle_corr; +} + +double RadialIntegrationMappingFPGA::QToBin(double q) const { + return std::min(static_cast(max_bin_number), std::max(0.0, (q - low_q) / q_spacing)); +} diff --git a/image_analysis/RadialIntegrationMappingFPGA.h b/image_analysis/RadialIntegrationMappingFPGA.h new file mode 100644 index 00000000..4957db8e --- /dev/null +++ b/image_analysis/RadialIntegrationMappingFPGA.h @@ -0,0 +1,25 @@ +// Copyright (2019-2023) Paul Scherrer Institute + +#ifndef JUNGFRAUJOCH_RADIALINTEGRATIONMAPPINGFPGA_H +#define JUNGFRAUJOCH_RADIALINTEGRATIONMAPPINGFPGA_H + +#include +#include +#include "../common/DiffractionExperiment.h" + +class RadialIntegrationMappingFPGA { + const double low_q, high_q, q_spacing; + std::vector bin_to_q; + std::vector solid_angle_corr; + std::vector pixel_to_bin; + uint16_t max_bin_number; +public: + explicit RadialIntegrationMappingFPGA(const DiffractionExperiment& experiment, uint16_t data_stream); + [[nodiscard]] uint16_t GetBinNumber() const; + [[nodiscard]] const std::vector &GetPixelToBinMapping() const; + [[nodiscard]] const std::vector &GetBinToQ() const; + [[nodiscard]] const std::vector &GetSolidAngleCorr() const; + [[nodiscard]] double QToBin(double q) const; +}; + +#endif //JUNGFRAUJOCH_RADIALINTEGRATIONMAPPINGFPGA_H diff --git a/receiver/AcquisitionDevice.cpp b/receiver/AcquisitionDevice.cpp index 067a7da5..d85a1baa 100644 --- a/receiver/AcquisitionDevice.cpp +++ b/receiver/AcquisitionDevice.cpp @@ -13,6 +13,7 @@ #include "../common/JFJochException.h" #include "AcquisitionDevice.h" #include "../common/NetworkAddressConvert.h" +#include "../image_analysis/RadialIntegrationMappingFPGA.h" void *mmap_acquisition_buffer(size_t size, int16_t numa_node) { void *ret = mmap(nullptr, size, PROT_READ | PROT_WRITE, MAP_PRIVATE | MAP_ANONYMOUS, -1, 0); @@ -43,6 +44,11 @@ void AcquisitionDevice::PrepareAction(const DiffractionExperiment &experiment) { throw(JFJochException(JFJochExceptionCategory::InputParameterAboveMax, "Number of modules exceeds max possible for FPGA")); + if (experiment.GetDetectorMode() == DetectorMode::Conversion) { + RadialIntegrationMappingFPGA map(experiment, data_stream); + InitializeIntegrationMap(experiment, map.GetPixelToBinMapping()); + } + counters.Reset(experiment, data_stream); } diff --git a/receiver/AcquisitionDevice.h b/receiver/AcquisitionDevice.h index 7643efa9..897bcf95 100644 --- a/receiver/AcquisitionDevice.h +++ b/receiver/AcquisitionDevice.h @@ -105,6 +105,7 @@ public: // Calibration virtual void InitializeCalibration(const DiffractionExperiment &experiment, const JFCalibration &calib); virtual void InitializeIntegrationMap(const DiffractionExperiment &experiment, const std::vector &v); + void InitializeIntegrationMap(const DiffractionExperiment& experiment); const AcquisitionCounters& Counters() const; virtual void SetSpotFinderParameters(int16_t count_threshold, double snr_threshold) {}; diff --git a/tests/RadialIntegrationTest.cpp b/tests/RadialIntegrationTest.cpp index 963ab7fa..e7b861a5 100644 --- a/tests/RadialIntegrationTest.cpp +++ b/tests/RadialIntegrationTest.cpp @@ -3,7 +3,7 @@ #include #include "../image_analysis/RadialIntegrationProfile.h" -#include "../image_analysis/RadialIntegrationMapping.h" +#include "../image_analysis/RadialIntegrationMappingFPGA.h" #include "../image_analysis/RadialIntegration.h" TEST_CASE("RadialIntegrationMapping_Constructor","[RadialIntegration]") { @@ -15,6 +15,15 @@ TEST_CASE("RadialIntegrationMapping_Constructor","[RadialIntegration]") { REQUIRE_NOTHROW(radial = std::make_unique(x)); } +TEST_CASE("RadialIntegrationMappingFPGA_Constructor","[RadialIntegration]") { + DiffractionExperiment x; + + std::unique_ptr radial; + + x.QSpacingForRadialInt_recipA(0.1).LowQForRadialInt_recipA(0.1).HighQForRadialInt_recipA(5); + REQUIRE_NOTHROW(radial = std::make_unique(x, 0)); +} + TEST_CASE("RadialIntegrationMapping_Constructor_Mask","[RadialIntegration]") { DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36)); @@ -41,6 +50,14 @@ TEST_CASE("RadialIntegrationMapping_GetBinNumber","[RadialIntegration]") { REQUIRE(mapping.GetBinNumber() == 39); } +TEST_CASE("RadialIntegrationMappingFPGA_GetBinNumber","[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); + RadialIntegrationMappingFPGA mapping(x, 0); + REQUIRE(mapping.GetBinNumber() == 39); +} + TEST_CASE("RadialIntegrationMapping_GetBinNumber_DetectorLimit","[RadialIntegration]") { DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36)); x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000); @@ -50,6 +67,15 @@ TEST_CASE("RadialIntegrationMapping_GetBinNumber_DetectorLimit","[RadialIntegrat REQUIRE(mapping.GetBinNumber() < 80); } +TEST_CASE("RadialIntegrationMappingFPGA_GetBinNumber_DetectorLimit","[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(9.9); + RadialIntegrationMappingFPGA mapping(x, 0); + REQUIRE(mapping.GetBinNumber() < 80); +} + TEST_CASE("RadialIntegrationMapping_GetBinToQ","[RadialIntegration]") { DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36)); x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000); @@ -64,6 +90,20 @@ TEST_CASE("RadialIntegrationMapping_GetBinToQ","[RadialIntegration]") { REQUIRE(bin_to_q[38] == Approx(3.95)); } +TEST_CASE("RadialIntegrationMappingFPGA_GetBinToQ","[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); + + RadialIntegrationMappingFPGA mapping(x, 0); + auto bin_to_q = mapping.GetBinToQ(); + + REQUIRE(bin_to_q[0] == Approx(0.15)); + REQUIRE(bin_to_q[1] == Approx(0.25)); + REQUIRE(bin_to_q[15] == Approx(1.65)); + REQUIRE(bin_to_q[38] == Approx(3.95)); +} + TEST_CASE("RadialIntegrationMapping_QToBin","[RadialIntegration]") { DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36)); x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000); @@ -77,6 +117,19 @@ TEST_CASE("RadialIntegrationMapping_QToBin","[RadialIntegration]") { REQUIRE(mapping.QToBin(50.0) == Approx(38)); } +TEST_CASE("RadialIntegrationMappingFPGA_QToBin","[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); + + RadialIntegrationMappingFPGA mapping(x, 0); + + REQUIRE(mapping.QToBin(0.0) == 0); + REQUIRE(std::floor(mapping.QToBin(0.200001)) == 1); + REQUIRE(mapping.QToBin(0.6) == Approx(5)); + 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);