v1.0.0-rc.179 (#90)

* Rugnux: Stills processing needs only the unit cell and space group beyond what the file records - detector distance, beam centre, tilt, panel placement, saturation level, gain and polarisation plane are found from the images - and is up to twice as fast.
* Rugnux: Stills can be indexed without a GPU, with a CPU implementation of the fast-feedback indexer (macOS, computers without a GPU).
* Rugnux: Rotation and stills data are integrated in the detector's own units with one gain model, improving data from CCD detectors.
* Web frontend: The Rugnux live page shows the placed model with its electron density and anomalous difference map in a Mol* viewer.

Reviewed-on: #90
Co-authored-by: Filip Leonarski <filip.leonarski@psi.ch>

tests/BraggIntegrationEngineMetalTest.cpp

@@ -3,26 +3,32 @@
 
 #include <catch2/catch_all.hpp>
 
-#ifdef __APPLE__
-
+#include <algorithm>
 #include <bit>
 #include <cmath>
 #include <cstdint>
+#include <map>
 #include <random>
+#include <tuple>
 #include <vector>
 
 #include "../common/BraggIntegrationSettings.h"
 #include "../common/DetectorSetup.h"
 #include "../common/DiffractionExperiment.h"
+#include "../common/JFJochMath.h"
 #include "../common/PixelMask.h"
 #include "../common/Reflection.h"
 #include "../image_analysis/bragg_integration/BraggIntegrationEngineCPU.h"
+#ifdef __APPLE__
 #include "../image_analysis/bragg_integration/BraggIntegrationEngineMetal.h"
+#endif
 
 // The Metal engine's whole claim is that it IS the CPU engine, bit for bit (BraggIntegrationEngineMetal.h).
 // So nothing here has a tolerance: every float of every reflection is compared as a bit pattern, and
 // every count as an integer. A Mac without a Metal device fails rather than skips - the macOS CI job is
-// the only place this can be checked, and a skip there would check nothing.
+// the only place this can be checked, and a skip there would check nothing. The frames and the cases
+// are built everywhere, and one test case checks on the CPU engine alone that they reach what they are
+// there for.
 
 namespace {
 
@@ -106,17 +112,34 @@ struct Case {
     float stencil_k = 0.0f;
     bool footprint = false;
     float clip_nsigma = 4.0f;
+    bool still = false;       // a still keeps its overloads as lower bounds (saturated pixels at the limit)
+    float gain = 1.0f;        // ADU per photon
+    bool curve = false;       // a radial curve handed over (RadialBackgroundFromRings), which the clip follows
 };
 
+// The radial curve a case hands over: rising steeply enough with the radius that the clip it makes the
+// ring follow (RadialClipOffset) moves by tens of counts across a ring.
+void HandOverCurve(BraggIntegrationEngine &engine) {
+    std::vector<float> rings(400);
+    for (size_t i = 0; i < rings.size(); ++i)
+        rings[i] = 1.0f + 0.01f * static_cast<float>(i * i);
+    engine.RadialBackgroundFromRings(rings, static_cast<int32_t>(rings.size()), 0.05f, 0.01f);
+}
+
 DiffractionExperiment MakeExperiment(const DetectorSetup &det, const Case &c) {
     DiffractionExperiment x(det);
     x.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV)
      .BeamX_pxl(0.4f * static_cast<float>(x.GetXPixelsNum()))
      .BeamY_pxl(0.55f * static_cast<float>(x.GetYPixelsNum()));
+    // A rotation drops an overloaded observation, a still keeps it as a lower bound (Finalize).
+    if (!c.still) {
+        x.Goniometer(GoniometerAxis("omega", 0.0f, 0.1f, Coord(-1, 0, 0), {}));
+        x.ImportIndexingSettings(IndexingSettings().RotationIndexing(true));
+    }
     x.BandwidthFWHM(c.bandwidth);
     BraggIntegrationSettings s;
     s.Integrator(c.mode).Overlap(c.overlap).BackgroundRadialCorrection(c.radial).StencilKSigma(c.stencil_k)
-     .BackgroundClipNSigma(c.clip_nsigma);
+     .BackgroundClipNSigma(c.clip_nsigma).PixelGain(c.gain);
     if (c.footprint) {
         // Wider than the r1 disk can hold from the middle out, so both stencils occur on one frame.
         SpotFootprint fp;
@@ -131,6 +154,162 @@ DiffractionExperiment MakeExperiment(const DetectorSetup &det, const Case &c) {
     return x;
 }
 
+// The beam-stop shadow of BraggIntegrationEngineCPU_RadialCurveOverThePixelsRead, on the same tilted
+// detector: the background rises towards the beam until the shadow, whose pixels are masked, cuts it
+// off, and the frame's curve holds that edge. Here with a few dozen reflections around the beam whose
+// rings reach into the shadow, some of them all the way, with noise and hot pixels (so the clip rejects
+// something). Where the curve is in use, the correction is summed over the ring pixels left
+// (BraggCurveCorrection) and the curve read at each reflection's 2 theta (CurveRadius), which the tilt
+// puts 17 px from its distance to the PONI.
+struct ShadowEdge {
+    DiffractionExperiment x{DetJF(2)};
+    Frame f;
+    std::vector<float> rings;
+    float low_q = 0.05f, dq = 0.0002f;
+    float beam_x = 0.0f, beam_y = 0.0f;   // where the beam lands
+};
+
+ShadowEdge MakeShadowEdge(IntegratorMode mode) {
+    ShadowEdge e;
+    e.x.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(400.0f).BeamY_pxl(400.0f)
+            .PoniRot1_rad(0.01f);
+    e.x.ImportBraggIntegrationSettings(BraggIntegrationSettings().Integrator(mode).BackgroundRadialCorrection(true));
+    const DiffractionGeometry geom = e.x.GetDiffractionGeometry();
+    const double F_px = 100.0 / 0.075, lambda = 1.0;
+    const auto radius = [&](double x, double y) {
+        return F_px * std::tan(2.0 * std::asin(lambda / (2.0 * geom.PxlToRes(static_cast<float>(x), static_cast<float>(y)))));
+    };
+    const double edge = 296.0;
+    const auto background = [](double r) { return 50.0 + 4.0 * std::max(0.0, 302.0 - r); };
+    const size_t W = e.x.GetXPixelsNum();
+    std::mt19937 rng(11);
+    std::uniform_real_distribution<float> U(0.0f, 1.0f);
+    e.f.image.resize(e.x.GetPixelsNum());
+    for (size_t i = 0; i < e.f.image.size(); ++i) {
+        const double r = radius(static_cast<double>(i % W), static_cast<double>(i / W));
+        e.f.image[i] = r < edge ? INT32_MIN
+                                : static_cast<int32_t>(std::lround(background(r) + (U(rng) - 0.5f) * 8.0f));
+    }
+    std::uniform_int_distribution<size_t> P(0, e.f.image.size() - 1);
+    for (size_t k = 0; k < e.f.image.size() / 3000; ++k) {
+        const size_t idx = P(rng);
+        if (e.f.image[idx] != INT32_MIN) e.f.image[idx] += 300;
+    }
+    e.rings.resize(10000);
+    for (size_t i = 0; i < e.rings.size(); ++i) {
+        const double q = e.low_q + (i + 0.5) * e.dq;
+        const double r = F_px * std::tan(2.0 * std::asin(q * lambda / (4.0 * PI)));
+        e.rings[i] = static_cast<float>(r < edge ? 10.0 : background(r));
+    }
+    std::tie(e.beam_x, e.beam_y) = geom.GetDirectBeam_pxl();
+    // Around the beam, 303 to 315 px out, at sub-pixel positions.
+    for (int k = 0; k < 48; ++k) {
+        const float rr = 303.37f + 4.0f * static_cast<float>(k % 4);
+        const float a = 2.0f * static_cast<float>(PI) * static_cast<float>(k) / 48.0f + 0.01f;
+        Reflection r{};
+        r.h = k + 1; r.k = 0; r.l = 1;
+        r.predicted_x = e.beam_x + rr * std::cos(a);
+        r.predicted_y = e.beam_y + rr * std::sin(a);
+        r.d = geom.PxlToRes(r.predicted_x, r.predicted_y);
+        r.prescaling_corr = 1.0f;
+        r.partiality = 1.0f;
+        e.f.predicted.push_back(r);
+    }
+    return e;
+}
+
+}   // namespace
+
+// Runs everywhere, on the CPU engine alone: the cases the Metal comparison below adds for the stills, the
+// gain and the handed-over curve do reach what they are there for, so that comparison cannot pass by
+// comparing nothing.
+TEST_CASE("BraggIntegrationEngineMetal_CasesReachTheirModes", "[portable][Integration]") {
+    Case c;
+    c.still = true;
+    c.gain = 4.0f;
+    c.radial = true;
+    c.curve = true;
+    const DiffractionExperiment x = MakeExperiment(DetJF(2), c);
+    const PixelMask mask(x);
+    BraggIntegrationEngineCPU cpu(x, mask);
+    HandOverCurve(cpu);
+    cpu.BackgroundRadial(true);
+    // The clip follows the curve: it moves by whole counts across a ring.
+    CHECK(cpu.RadialClipOffset(500.0f, 10.0f) > 10);
+    CHECK(cpu.RadialClipOffset(500.0f, -10.0f) < -10);
+    ImagePreprocessorBuffer image(x.GetPixelsNum());
+    const Frame f = RandomFrame(x, 7, {});
+    std::copy(f.image.begin(), f.image.end(), image.data());
+    const auto out = cpu.Run(image, f.predicted, f.predicted.size(), 0);
+    CHECK(out.size() > 1000);
+    // Saturated pixels counted at the limit: overloads kept as lower bounds.
+    CHECK(cpu.Counts().overloads_kept > 0);
+}
+
+// The same, for the beam-stop shadow the Metal comparison runs (MakeShadowEdge): the tilt moves the beam
+// well off the PONI, so a kernel reading the curve at the distance from the PONI would differ; the rings
+// next to the shadow have lost pixels to it, which the correction is taken over; and the correction
+// takes the background read off them back towards the one under the disk, in the profile modes (the
+// clip's kept pixels) and in the box sum (the whole ring) alike.
+TEST_CASE("BraggIntegrationEngineMetal_ShadowEdgeReachesItsModes", "[portable][Integration]") {
+    for (IntegratorMode mode : {IntegratorMode::BoxSum, IntegratorMode::ProfileGaussian}) {
+        INFO("mode " << static_cast<int>(mode));
+        ShadowEdge e = MakeShadowEdge(mode);
+        REQUIRE(std::hypot(e.beam_x - 400.0f, e.beam_y - 400.0f) > 10.0f);
+        const size_t W = e.x.GetXPixelsNum();
+        const PixelMask mask(e.x);
+        ImagePreprocessorBuffer image(e.x.GetPixelsNum());
+        std::copy(e.f.image.begin(), e.f.image.end(), image.data());
+        auto integrate = [&](bool correct) {
+            BraggIntegrationEngineCPU engine(e.x, mask);
+            engine.RadialBackgroundFromRings(e.rings, static_cast<int32_t>(e.rings.size()), e.low_q, e.dq);
+            engine.BackgroundRadial(correct);
+            return engine.Run(image, e.f.predicted, e.f.predicted.size(), 0);
+        };
+        const auto plain = integrate(false), corrected = integrate(true);
+        CHECK(plain.size() > 40);
+        std::map<int64_t, float> bkg_corrected;
+        for (const Reflection &r : corrected)
+            bkg_corrected[r.h] = r.bkg;
+        // The rings 303 px out reach into the shadow; over those reflections, how far each background is
+        // from the mean of its own disk (hot pixels left out).
+        const BraggIntegrationSettings bs = e.x.GetBraggIntegrationSettings();
+        const float r1_sq = bs.GetR1() * bs.GetR1(), r2_sq = bs.GetR2() * bs.GetR2(), r3_sq = bs.GetR3() * bs.GetR3();
+        double err_plain = 0.0, err_corrected = 0.0;
+        int n = 0, n_shadowed = 0;
+        for (const Reflection &r : plain) {
+            if (std::hypot(r.predicted_x - e.beam_x, r.predicted_y - e.beam_y) > 304.0f)
+                continue;
+            REQUIRE(bkg_corrected.count(r.h) == 1);
+            double truth = 0.0;
+            int n_disk = 0, n_masked = 0;
+            const int cx = static_cast<int>(std::lround(r.predicted_x)), cy = static_cast<int>(std::lround(r.predicted_y));
+            for (int dy = -12; dy <= 12; ++dy)
+                for (int dx = -12; dx <= 12; ++dx) {
+                    const int32_t v = e.f.image[static_cast<size_t>(cy + dy) * W + static_cast<size_t>(cx + dx)];
+                    const float ex = static_cast<float>(cx + dx) - r.predicted_x, ey = static_cast<float>(cy + dy) - r.predicted_y;
+                    const float d2 = ex * ex + ey * ey;
+                    if (d2 < r1_sq && v != INT32_MIN && v < 200) { truth += v; ++n_disk; }
+                    if (d2 >= r2_sq && d2 < r3_sq && v == INT32_MIN) ++n_masked;
+                }
+            if (n_masked > 0) ++n_shadowed;
+            REQUIRE(n_disk > 0);
+            truth /= n_disk;
+            err_plain += std::fabs(r.bkg - truth);
+            err_corrected += std::fabs(bkg_corrected[r.h] - truth);
+            ++n;
+        }
+        REQUIRE(n > 5);
+        CHECK(n_shadowed == n);
+        CHECK(err_plain / n > 2.0);
+        CHECK(err_corrected < 0.5 * err_plain);
+    }
+}
+
+#ifdef __APPLE__
+
+namespace {
+
 uint32_t Bits(float f) { return std::bit_cast<uint32_t>(f); }
 
 void RequireSame(const std::vector<Reflection> &cpu, const std::vector<Reflection> &metal) {
@@ -153,6 +332,7 @@ void RequireSame(const std::vector<Reflection> &cpu, const std::vector<Reflectio
         REQUIRE(Bits(a.prescaling_corr) == Bits(b.prescaling_corr));
         REQUIRE(a.observed == b.observed);
         REQUIRE(a.overloaded == b.overloaded);
+        REQUIRE(a.clipped == b.clipped);
     }
 }
 
@@ -163,6 +343,7 @@ void RequireSame(const BraggIntegrationCounts &a, const BraggIntegrationCounts &
     CHECK(a.profile_fallback == b.profile_fallback);
     CHECK(a.overlapped == b.overlapped);
     CHECK(a.overlapped_major == b.overlapped_major);
+    CHECK(a.overloads_kept == b.overloads_kept);
 }
 
 // Three frames through one pair of engines - so whatever the Metal engine keeps between frames (its
@@ -185,6 +366,10 @@ size_t CompareFrames(const DetectorSetup &det, const Case &c) {
 
     BraggIntegrationEngineCPU cpu(x, mask);
     BraggIntegrationEngineMetal metal(x, mask);
+    if (c.curve) {
+        HandOverCurve(cpu);
+        HandOverCurve(metal);
+    }
     ImagePreprocessorBuffer image(x.GetPixelsNum());
 
     size_t n_out = 0;
@@ -201,6 +386,9 @@ size_t CompareFrames(const DetectorSetup &det, const Case &c) {
         n_out += out_cpu.size();
     }
     RequireSame(cpu.Counts(), metal.Counts());
+    // Overloads kept as lower bounds count saturated pixels: the case is not vacuous there either.
+    if (c.still)
+        CHECK(cpu.Counts().overloads_kept > 0);
     // Not the CPU engine compared with itself: every frame went through the kernels.
     REQUIRE(metal.FramesOnGPU() == 3);
     return n_out;
@@ -241,6 +429,39 @@ TEST_CASE("BraggIntegrationEngineMetal_BitIdenticalToCPU", "[portable][Integrati
     Case no_clip;
     no_clip.clip_nsigma = 0.0f;
     cases.emplace_back("no clip", no_clip);
+    // A radial curve handed over, which the ring's clip follows (RadialClipOffset), under every overlap
+    // mode and every integrator.
+    for (IntegratorMode mode : {IntegratorMode::BoxSum, IntegratorMode::ProfileGaussian, IntegratorMode::ProfileEmpirical})
+        for (OverlapMode overlap : {OverlapMode::Off, OverlapMode::Exclude}) {
+            Case curve;
+            curve.mode = mode;
+            curve.overlap = overlap;
+            curve.radial = true;
+            curve.curve = true;
+            cases.emplace_back("curve mode " + std::to_string(static_cast<int>(mode)) + " overlap "
+                               + std::to_string(static_cast<int>(overlap)), curve);
+        }
+    // Stills: overloads kept as lower bounds, every saturated pixel counted. At a gain of one, and in
+    // detector units at two gains (the clip, the floors and the output's units all move with it), with
+    // and without the handed-over curve.
+    for (IntegratorMode mode : {IntegratorMode::BoxSum, IntegratorMode::ProfileGaussian, IntegratorMode::ProfileEmpirical})
+        for (float gain : {1.0f, 4.0f, 1.7f}) {
+            Case still;
+            still.mode = mode;
+            still.still = true;
+            still.gain = gain;
+            still.radial = gain != 1.0f;
+            still.curve = gain == 4.0f;
+            cases.emplace_back("still mode " + std::to_string(static_cast<int>(mode)) + " gain "
+                               + std::to_string(gain) + (still.radial ? " radial" : "")
+                               + (still.curve ? " curve" : ""), still);
+        }
+    // Footprint stencils on a still with a curve: the clip follows it on a per-reflection stencil.
+    Case footprint_still = footprint;
+    footprint_still.still = true;
+    footprint_still.curve = true;
+    footprint_still.gain = 4.0f;
+    cases.emplace_back("footprint still gain 4 curve", footprint_still);
 
     for (const auto &[name, c] : cases) {
         DYNAMIC_SECTION(name) {
@@ -250,6 +471,32 @@ TEST_CASE("BraggIntegrationEngineMetal_BitIdenticalToCPU", "[portable][Integrati
     }
 }
 
+// The beam-stop shadow on a tilted detector (MakeShadowEdge): the curve summed over rings that lost pixels
+// to the shadow, read at each reflection's 2 theta, with the clip (profile modes) and without (box sum).
+TEST_CASE("BraggIntegrationEngineMetal_ShadowEdgeBitIdenticalToCPU", "[portable][Integration][Metal]") {
+    for (IntegratorMode mode : {IntegratorMode::BoxSum, IntegratorMode::ProfileGaussian}) {
+        INFO("mode " << static_cast<int>(mode));
+        const ShadowEdge e = MakeShadowEdge(mode);
+        REQUIRE(BraggIntegrationEngineMetal::Available(e.x));
+        const PixelMask mask(e.x);
+        BraggIntegrationEngineCPU cpu(e.x, mask);
+        BraggIntegrationEngineMetal metal(e.x, mask);
+        for (BraggIntegrationEngine *engine : {static_cast<BraggIntegrationEngine *>(&cpu),
+                                               static_cast<BraggIntegrationEngine *>(&metal)}) {
+            engine->RadialBackgroundFromRings(e.rings, static_cast<int32_t>(e.rings.size()), e.low_q, e.dq);
+            engine->BackgroundRadial(true);
+        }
+        ImagePreprocessorBuffer image(e.x.GetPixelsNum());
+        std::copy(e.f.image.begin(), e.f.image.end(), image.data());
+        const auto out_cpu = cpu.Run(image, e.f.predicted, e.f.predicted.size(), 0);
+        const auto out_metal = metal.Run(image, e.f.predicted, e.f.predicted.size(), 0);
+        CHECK(out_cpu.size() > 40);
+        RequireSame(out_cpu, out_metal);
+        RequireSame(cpu.Counts(), metal.Counts());
+        REQUIRE(metal.FramesOnGPU() == 1);
+    }
+}
+
 // The trimmed-mean background sorts each ring, which only the CPU engine does: such a run keeps the CPU
 // engine (MXAnalysisWithoutFPGA asks Available).
 TEST_CASE("BraggIntegrationEngineMetal_TrimmedMeanStaysOnCPU", "[portable][Integration][Metal]") {