Merge branch 'panel-s2' into panel-s5

tests/PanelLayoutTest.cpp

@@ -24,6 +24,7 @@
 #include <vector>
 
 #include "../common/AzimuthalIntegrationMapping.h"
+#include "../common/CUDAWrapper.h"
 #include "../common/BraggIntegrationSettings.h"
 #include "../common/DetectorModules.h"
 #include "../common/DetectorSetup.h"
@@ -31,6 +32,10 @@
 #include "../common/PixelMask.h"
 #include "../common/Reflection.h"
 #include "../image_analysis/bragg_integration/BraggIntegrationEngineCPU.h"
+#ifdef JFJOCH_USE_CUDA
+#include "../image_analysis/bragg_integration/BraggIntegrationEngineGPU.h"
+#include "../image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h"
+#endif
 #include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h"
 #include "../image_analysis/roi/ROIIntegrationCPU.h"
 #include "../image_analysis/spot_finding/ImageSpotFinderCPU.h"
@@ -149,7 +154,8 @@ void CheckSameSpots(const std::vector<SpotToSave> &a, const std::vector<SpotToSa
 // Integrate predictions placed on the plane: on the assembled image directly, on the stacked one through
 // the panel table - PlaneToPanel in, ImageToPlane of the same panel out. Keyed by h.
 std::map<int, Reflection> Integrate(const DiffractionExperiment &x, const PixelMask &mask,
-                                    const ImagePreprocessorBuffer &image, const std::vector<Spot> &spots) {
+                                    const ImagePreprocessorBuffer &image, const std::vector<Spot> &spots,
+                                    bool gpu = false) {
     const auto &modules = x.GetModules();
     const auto geom = x.GetDiffractionGeometry();
     std::vector<Reflection> predicted;
@@ -167,9 +173,25 @@ std::map<int, Reflection> Integrate(const DiffractionExperiment &x, const PixelM
         predicted.push_back(r);
         panel_of[r.h] = p->panel;
     }
-    BraggIntegrationEngineCPU engine(x, mask);
+    std::vector<Reflection> integrated;
+#ifdef JFJOCH_USE_CUDA
+    if (gpu) {
+        auto stream = std::make_shared<CudaStream>();
+        ImagePreprocessorBufferGPU gpu_image(image.size());
+        for (size_t i = 0; i < image.size(); i++)
+            gpu_image[i] = image[i];
+        REQUIRE(cudaMemcpyAsync(gpu_image.getGPUBuffer(), gpu_image.getBuffer().data(), image.size() * sizeof(int32_t),
+                                cudaMemcpyHostToDevice, *stream) == cudaSuccess);
+        BraggIntegrationEngineGPU engine(x, stream, mask);
+        integrated = engine.Run(gpu_image, predicted, predicted.size(), 0);
+    }
+#endif
+    if (!gpu) {
+        BraggIntegrationEngineCPU engine(x, mask);
+        integrated = engine.Run(image, predicted, predicted.size(), 0);
+    }
     std::map<int, Reflection> out;
-    for (auto r: engine.Run(image, predicted, predicted.size(), 0)) {
+    for (auto r: integrated) {
         const int64_t panel = panel_of.at(r.h);
         const auto [px, py] = modules.ImageToPlane(panel, r.predicted_x, r.predicted_y);
         r.predicted_x = static_cast<float>(px);
@@ -315,20 +337,60 @@ TEST_CASE("PanelLayout_ResolutionMask", "[PanelLayout]") {
     CHECK(differ == 0);
 }
 
-// TODO stage 2: Bragg integration on the stacked layout measures what it measures on the assembled image -
-// stencil direction, radial background and polarisation taken from the plane position, not the stored one.
-TEST_CASE("PanelLayout_BraggIntegration", "[.][PanelLayoutTodo]") {
-    const Scene scene(INTERIOR_SPOTS);
-    const auto a = Integrate(scene.assembled, scene.assembled_mask, scene.assembled_image, INTERIOR_SPOTS);
-    const auto s = Integrate(scene.stacked, scene.stacked_mask, scene.stacked_image, INTERIOR_SPOTS);
-    REQUIRE(s.size() == a.size());
-    for (const auto &[h, r]: s) {
-        CHECK(r.I == Catch::Approx(a.at(h).I).epsilon(1e-5));
-        CHECK(r.sigma == Catch::Approx(a.at(h).sigma).epsilon(1e-5));
-        CHECK(r.bkg == Catch::Approx(a.at(h).bkg).epsilon(1e-5));
+// Bragg integration on the stacked layout measures what it measures on the assembled image - stencil
+// direction, radial background and polarisation taken from the plane position, not the stored one.
+// Also with what reads the radial direction hardest: polarisation, a ring elongated by a bandwidth, the
+// radial background correction, the owner map, and the box sum. The interior spots moved off the exact
+// half pixel: there the profile grid centres on the rounded position, and rounding a half picks the other
+// neighbour once a panel is turned.
+void CheckBraggIntegration(bool gpu) {
+    const std::vector<Spot> spots = {{50.3, 40.7}, {80.4, 45.3}, {160.6, 30.2}, {40.4, 140.8}, {150.2, 130.6}};
+    const Scene scene(spots);
+    auto setup = [](DiffractionExperiment x, int variant) {
+        BraggIntegrationSettings bis = x.GetBraggIntegrationSettings();
+        x.PolarizationFactor(0.99f);
+        if (variant == 1) {
+            x.BandwidthFWHM(0.02f);
+            bis.StencilKSigma(2.0f);
+        }
+        if (variant == 2)
+            bis.BackgroundRadialCorrection(true);
+        if (variant == 3)
+            bis.Overlap(OverlapMode::Exclude);
+        if (variant == 4)
+            bis.Integrator(IntegratorMode::BoxSum);
+        if (variant == 5)
+            bis.Integrator(IntegratorMode::ProfileEmpirical);
+        x.ImportBraggIntegrationSettings(bis);
+        return x;
+    };
+    for (int variant = 0; variant < 6; variant++) {
+        CAPTURE(variant);
+        const auto a = Integrate(setup(scene.assembled, variant), scene.assembled_mask, scene.assembled_image, spots,
+                                 gpu);
+        const auto s = Integrate(setup(scene.stacked, variant), scene.stacked_mask, scene.stacked_image, spots, gpu);
+        REQUIRE(s.size() == a.size());
+        for (const auto &[h, r]: s) {
+            CAPTURE(h);
+            CHECK(r.I == Catch::Approx(a.at(h).I).epsilon(1e-5));
+            CHECK(r.sigma == Catch::Approx(a.at(h).sigma).epsilon(1e-5));
+            CHECK(r.bkg == Catch::Approx(a.at(h).bkg).epsilon(1e-5));
+            CHECK(r.image_scale_corr == Catch::Approx(a.at(h).image_scale_corr).epsilon(1e-5));
+        }
     }
 }
 
+TEST_CASE("PanelLayout_BraggIntegration", "[PanelLayout]") {
+    CheckBraggIntegration(false);
+}
+
+// The same on the CUDA engine, stacked against assembled on the GPU.
+TEST_CASE("PanelLayout_BraggIntegrationGPU", "[PanelLayout][GPU]") {
+    if (get_gpu_count() == 0)
+        SKIP("No GPU");
+    CheckBraggIntegration(true);
+}
+
 // TODO stage 3: spot finding per panel. Two half spots on the facing edges of panels 0 and 1 as stored -
 // adjacent rows of the stacked array, 139 px apart on the plane - stay two spots, and the local-background
 // (SNR) finder agrees with the assembled image away from panel edges.