diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index da48248f1..3eac9e404 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -14,6 +14,7 @@ ADD_EXECUTABLE(jfjoch_portable_test EXCLUDE_FROM_ALL XFELReaderTest.cpp SaturationCeilingTest.cpp DetectorModulesTest.cpp + PanelLayoutTest.cpp CellPopulationTest.cpp MultiReaderTest.cpp ModuleRefinementTest.cpp @@ -111,6 +112,7 @@ ADD_EXECUTABLE(jfjoch_test XFELReaderTest.cpp SaturationCeilingTest.cpp DetectorModulesTest.cpp + PanelLayoutTest.cpp CellPopulationTest.cpp ModuleRefinementTest.cpp LiveCollectionTest.cpp diff --git a/tests/PanelLayoutTest.cpp b/tests/PanelLayoutTest.cpp new file mode 100644 index 000000000..4c8f5a484 --- /dev/null +++ b/tests/PanelLayoutTest.cpp @@ -0,0 +1,323 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +// Equivalence of the stored (stacked) and the assembled layout of a panel detector. +// +// One synthetic detector, recorded twice: as an assembled image (panels where they are on the plane, gaps +// masked) and as the stacked array a panel detector stores (panels one under the other, in four +// orientations, no gaps), with the DetectorModules that place the stacked panels on the plane. Both carry +// the same value in every physical pixel, so anything computed on the plane must come out the same. +// +// [PanelLayout] cases check what is already layout-correct. The [PanelLayoutTodo] cases are hidden: they +// state the equivalence each later stage has to deliver, and fail until it does. A stage that converts a +// consumer retags its case [PanelLayout]: +// stage 1 - dense maps: PanelLayout_AzintMapping, PanelLayout_ResolutionMask +// stage 2 - sparse/Bragg: PanelLayout_BraggIntegration +// stage 3 - per-panel spot finding: PanelLayout_SpotsNeverCrossPanels; the raster algorithms +// (ShadowFinder, BeamCenterFFT, BeamCenterFromBackground, DiagnosticOutput) add their case here, +// comparing the result on the stacked layout with the result on the assembled image. + +#include + +#include +#include +#include + +#include "../common/AzimuthalIntegrationMapping.h" +#include "../common/BraggIntegrationSettings.h" +#include "../common/DetectorModules.h" +#include "../common/DetectorSetup.h" +#include "../common/DiffractionExperiment.h" +#include "../common/PixelMask.h" +#include "../common/Reflection.h" +#include "../image_analysis/bragg_integration/BraggIntegrationEngineCPU.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h" +#include "../image_analysis/spot_finding/ImageSpotFinderCPU.h" +#include "../image_analysis/spot_finding/SpotUtils.h" + +namespace { + +// One panel as stored, and the plane the four of them are placed on. +constexpr int64_t PW = 96, PH = 64, NPANELS = 4; +constexpr int64_t PLANE_W = 220, PLANE_H = 200; + +// Panel p is stored at rows [p * PH, (p + 1) * PH) and placed 2 x 2 on the plane with gaps of 8 and 16 px. +DetectorModules StackedPanels() { + auto panel = [](int64_t p, double ox, double oy, double fx, double fy, double sx, double sy) { + DetectorModule m; + m.x0 = 0; + m.y0 = p * PH; + m.width = PW; + m.height = PH; + m.ox = ox; m.oy = oy; + m.fx = fx; m.fy = fy; + m.sx = sx; m.sy = sy; + return m; + }; + return DetectorModules({ + panel(0, 10, 10, 1, 0, 0, 1), // as stored: x 10..105, y 10..73 + panel(1, 209, 73, -1, 0, 0, -1), // turned 180: x 114..209, y 10..73 + panel(2, 73, 90, 0, 1, -1, 0), // turned 90: x 10..73, y 90..185 + panel(3, 120, 185, 0, -1, 1, 0) // turned -90: x 120..183, y 90..185 + }); +} + +struct Spot { + double x, y; // on the plane +}; + +// Spots well inside their panels: more than the spot finder's 15 px window from every panel edge. +const std::vector INTERIOR_SPOTS = {{50.3, 40.7}, {80.5, 45.3}, {160.6, 30.2}, {40.4, 140.8}, {150.2, 130.5}}; + +DiffractionExperiment Experiment(int64_t width, int64_t height, const DetectorModules &modules) { + DiffractionExperiment x(DetDECTRIS(width, height, "PanelLayout", "")); + // Beam in the gap cross; at 30 mm the plane reaches about 2 A in its corners. + x.DetectorDistance_mm(30.0f).IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(110.0f).BeamY_pxl(82.0f); + x.ImportBraggIntegrationSettings(BraggIntegrationSettings()); + x.Modules(modules); + return x; +} + +// The same detector, assembled and stacked. +struct Scene { + DetectorModules modules = StackedPanels(); + DiffractionExperiment assembled = Experiment(PLANE_W, PLANE_H, {}); + DiffractionExperiment stacked = Experiment(PW, NPANELS * PH, modules); + ImagePreprocessorBuffer assembled_image{PLANE_W * PLANE_H}; + ImagePreprocessorBuffer stacked_image{PW * NPANELS * PH}; + PixelMask assembled_mask; + PixelMask stacked_mask{PW, NPANELS * PH}; + + // Assembled pixel at the plane position of a stacked pixel (the placements are whole-pixel). + static size_t AssembledIndex(double px, double py) { + return static_cast(std::lround(py)) * PLANE_W + static_cast(std::lround(px)); + } + + explicit Scene(const std::vector &spots) { + auto value = [&](double px, double py) { + double v = 2.0; + for (const auto &s: spots) { + const double dx = px - s.x, dy = py - s.y; + v += 800.0 * std::exp(-(dx * dx + dy * dy) / (2.0 * 1.3 * 1.3)); + } + return static_cast(std::lround(v)); + }; + // Gaps are what the preprocessor makes of a masked pixel. + for (size_t i = 0; i < assembled_image.size(); i++) + assembled_image[i] = INT32_MIN; + std::vector gaps(PLANE_W * PLANE_H, 1u << PixelMask::ModuleGapPixelBit); + modules.ForEachPixel(PW, NPANELS * PH, [&](int64_t, size_t i, double px, double py) { + stacked_image[i] = value(px, py); + assembled_image[AssembledIndex(px, py)] = value(px, py); + gaps[AssembledIndex(px, py)] = 0; + }); + assembled_mask = PixelMask(gaps); + } +}; + +std::vector FindSpots(const DiffractionExperiment &x, const ImagePreprocessorBuffer &image, + float snr = 0.0f) { + SpotFindingSettings settings; + settings.signal_to_noise_threshold = snr; + settings.photon_count_threshold = 20; + settings.min_pix_per_spot = 1; + settings.high_res_gap_Q_recipA = std::nullopt; + ImageSpotFinderCPU finder(static_cast(x.GetXPixelsNum()), static_cast(x.GetYPixelsNum())); + finder.SetResolutionMask(std::vector(x.GetPixelsNum(), false)); + DataMessage msg{}; + msg.number = 0; + SpotAnalyze(x, settings, finder.Run(image, settings), msg); + return msg.spots; +} + +// Every spot of `a` has one in `b` at the same plane position, and the counts agree. +void CheckSameSpots(const std::vector &a, const std::vector &b) { + REQUIRE(a.size() == b.size()); + for (const auto &s: a) { + size_t matched = 0; + for (const auto &t: b) + if (std::hypot(s.x - t.x, s.y - t.y) < 1e-3) { + CHECK(s.intensity == t.intensity); + CHECK(s.d_A == Catch::Approx(t.d_A)); + matched++; + } + CHECK(matched == 1); + } +} + +// 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 Integrate(const DiffractionExperiment &x, const PixelMask &mask, + const ImagePreprocessorBuffer &image, const std::vector &spots) { + const auto &modules = x.GetModules(); + const auto geom = x.GetDiffractionGeometry(); + std::vector predicted; + std::map panel_of; + for (size_t i = 0; i < spots.size(); i++) { + const auto p = modules.PlaneToPanel(spots[i].x, spots[i].y); + REQUIRE(p.has_value()); + Reflection r{}; + r.h = static_cast(i) + 1; r.k = 0; r.l = 1; + r.predicted_x = static_cast(p->x); + r.predicted_y = static_cast(p->y); + r.d = geom.PxlToRes(static_cast(spots[i].x), static_cast(spots[i].y)); + r.prescaling_corr = 1.0f; + r.partiality = 1.0f; + predicted.push_back(r); + panel_of[r.h] = p->panel; + } + BraggIntegrationEngineCPU engine(x, mask); + std::map out; + for (auto r: engine.Run(image, predicted, predicted.size(), 0)) { + 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(px); + r.predicted_y = static_cast(py); + const auto [ox, oy] = modules.ImageToPlane(panel, r.observed_x, r.observed_y); + r.observed_x = static_cast(ox); + r.observed_y = static_cast(oy); + out[r.h] = r; + } + return out; +} + +} // namespace + +// The panel table itself: every stored pixel is visited once and lands on a free plane pixel, plane -> +// panel inverts it, a gap is on no panel, and the plane bounds hold every panel. +TEST_CASE("PanelLayout_PanelTable", "[PanelLayout]") { + const DetectorModules modules = StackedPanels(); + const PlaneBox box = modules.PlaneBounds(PW, NPANELS * PH); + CHECK(box.x0 == 9.5); + CHECK(box.y0 == 9.5); + CHECK(box.x1 == 209.5); + CHECK(box.y1 == 185.5); + + std::vector visits(PW * NPANELS * PH, 0), on_plane(PLANE_W * PLANE_H, 0); + modules.ForEachPixel(PW, NPANELS * PH, [&](int64_t panel, size_t i, double px, double py) { + visits[i]++; + on_plane[Scene::AssembledIndex(px, py)]++; + CHECK(panel == static_cast(i / (PW * PH))); + const auto p = modules.PlaneToPanel(px, py); + REQUIRE(p.has_value()); + CHECK(p->panel == panel); + CHECK(static_cast(std::lround(p->y)) * PW + static_cast(std::lround(p->x)) == i); + const auto [qx, qy] = modules.PlaneToImage(panel, px, py); + CHECK(qx == Catch::Approx(p->x)); + CHECK(qy == Catch::Approx(p->y)); + }); + for (int v: visits) + CHECK(v == 1); + for (int v: on_plane) + CHECK(v <= 1); + + CHECK_FALSE(modules.PlaneToPanel(110.0, 40.0).has_value()); // gap between panels 0 and 1 + CHECK_FALSE(modules.PlaneToPanel(110.0, 82.0).has_value()); // gap cross, where the beam is + + // Empty table: the identity, the whole image one panel -1. + const DetectorModules none; + const auto p = none.PlaneToPanel(12.5, -3.0); + REQUIRE(p.has_value()); + CHECK(p->panel == -1); + CHECK(p->x == 12.5); + CHECK(p->y == -3.0); + CHECK(none.ImageToPlane(-1, 4.0, 5.0) == std::pair{4.0, 5.0}); + const PlaneBox image = none.PlaneBounds(PW, PH); + CHECK(image.x0 == 0.0); + CHECK(image.x1 == static_cast(PW)); + size_t n = 0; + none.ForEachPixel(PW, PH, [&](int64_t panel, size_t i, double px, double py) { + CHECK(panel == -1); + CHECK(i == n++); + CHECK(px == static_cast(i % PW)); + CHECK(py == static_cast(i / PW)); + }); + CHECK(n == static_cast(PW * PH)); +} + +// Spot centroids: found in the image, mapped to the plane by SpotAnalyze. A count threshold on interior +// spots, so the finder's neighbourhood never reaches another panel. +TEST_CASE("PanelLayout_SpotPlaneCoordinates", "[PanelLayout]") { + const Scene scene(INTERIOR_SPOTS); + const auto a = FindSpots(scene.assembled, scene.assembled_image); + const auto s = FindSpots(scene.stacked, scene.stacked_image); + CHECK(a.size() == INTERIOR_SPOTS.size()); + CheckSameSpots(a, s); +} + +// Predictions go plane -> panel -> image and come back to where they were on the plane; the integrated +// spot is observed where it is on the plane in both layouts. +TEST_CASE("PanelLayout_IntegrationRoundTrip", "[PanelLayout]") { + 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(a.size() == INTERIOR_SPOTS.size()); + REQUIRE(s.size() == INTERIOR_SPOTS.size()); + for (const auto &[h, r]: s) { + const Spot &spot = INTERIOR_SPOTS[h - 1]; + CHECK(r.predicted_x == Catch::Approx(spot.x).margin(1e-3)); + CHECK(r.predicted_y == Catch::Approx(spot.y).margin(1e-3)); + CHECK(r.observed_x == Catch::Approx(a.at(h).observed_x).margin(1e-3)); + CHECK(r.observed_y == Catch::Approx(a.at(h).observed_y).margin(1e-3)); + } +} + +// TODO stage 1: the azimuthal-integration map built per stored pixel through the panel table - same bin, +// resolution and correction as the assembled pixel it lands on. +TEST_CASE("PanelLayout_AzintMapping", "[.][PanelLayoutTodo]") { + const Scene scene(INTERIOR_SPOTS); + const AzimuthalIntegrationMapping a(scene.assembled, scene.assembled_mask); + const AzimuthalIntegrationMapping s(scene.stacked, scene.stacked_mask); + size_t bin = 0, resolution = 0, correction = 0; // pixels that differ + scene.modules.ForEachPixel(PW, NPANELS * PH, [&](int64_t, size_t i, double px, double py) { + const size_t j = Scene::AssembledIndex(px, py); + bin += s.GetPixelToBin()[i] != a.GetPixelToBin()[j]; + resolution += std::fabs(s.Resolution()[i] - a.Resolution()[j]) > 1e-5f * a.Resolution()[j]; + correction += std::fabs(s.Corrections()[i] - a.Corrections()[j]) > 1e-5f * a.Corrections()[j]; + }); + CHECK(bin == 0); + CHECK(resolution == 0); + CHECK(correction == 0); +} + +// TODO stage 1: the spot finders' resolution mask, per stored pixel through the panel table. +TEST_CASE("PanelLayout_ResolutionMask", "[.][PanelLayoutTodo]") { + const Scene scene(INTERIOR_SPOTS); + const auto a = AzimuthalIntegrationMapping(scene.assembled, scene.assembled_mask).ResolutionMaskBits(3.0f, 8.0f); + const auto s = AzimuthalIntegrationMapping(scene.stacked, scene.stacked_mask).ResolutionMaskBits(3.0f, 8.0f); + auto bit = [](const std::vector &v, size_t i) { return (v[i / 32] >> (i % 32)) & 1u; }; + size_t differ = 0; + scene.modules.ForEachPixel(PW, NPANELS * PH, [&](int64_t, size_t i, double px, double py) { + differ += bit(*s, i) != bit(*a, Scene::AssembledIndex(px, py)); + }); + 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)); + } +} + +// 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. +TEST_CASE("PanelLayout_SpotsNeverCrossPanels", "[.][PanelLayoutTodo]") { + // Panel 0 stored row 63 and panel 1 stored row 0, both at stored column 30. + const std::vector edge = {{40.0, 73.0}, {179.0, 73.0}}; + const Scene scene(edge); + CheckSameSpots(FindSpots(scene.assembled, scene.assembled_image), + FindSpots(scene.stacked, scene.stacked_image)); + + const Scene interior(INTERIOR_SPOTS); + CheckSameSpots(FindSpots(interior.assembled, interior.assembled_image, 3.0f), + FindSpots(interior.stacked, interior.stacked_image, 3.0f)); +}