mx/Jungfraujoch
3
2 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 92 Wiki Activity
Files
rc184
Jungfraujoch/tests/BraggIntegrationEngineMetalTest.cpp
T
leonarski_fandClaude Opus 5.5 5c28330de6 Azimuthal integration: Q bins equally spaced in log10(Q), or linear by bin count 
The Q bins can now be laid out in equal steps of log10(Q) between low and high Q,
dense at low Q and coarse at high Q, as small-angle work needs; such binning is
given by its number of bins. Linear binning takes either its spacing (as before)
or its number of bins, not both. Linear binning by spacing bins every pixel and
reports every bin centre exactly as before.

The Q axis (bin_to_q, written to CBOR and HDF5) holds each bin's centre on its
binning scale: the geometric mean of the edges for log10 bins. The ring background
of the Bragg radial correction and the ice-ring score read the bin centres instead
of assuming a constant spacing. A q rounding onto the far edge of the last bin now
stays in that bin instead of spilling into the next azimuthal row.

rugnux: --azim-q-binning linear|log10 and --azim-q-bins; --mode calibration
refuses log10 bins, as its ring search works in equal Q steps.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK
2026-11-10 16:55:48 +01:00

514 lines
25 KiB
C++
Raw Permalink Blame History

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include <catch2/catch_all.hpp>
#include "LinearQCentres.h"
#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 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 {
struct Frame {
std::vector<int32_t> image;
std::vector<Reflection> predicted;
};
// A crowded random frame: predictions at random sub-pixel positions (so the pixel classification meets
// every kind of rounding near the stencil boundaries), a quarter of them on whole and half pixels (so it
// meets exact ties with the radii too), with random partialities (rocking-curve tails and flux, which
// mask the rings differently), about a third of them with a spot, overlapping where they fall close. On
// a background that rises with the radius (the radial correction has a curve to fit) with noise that
// takes it below zero near the beam, as a pedestal-subtracted detector does (a negative pixel is what
// carries into the high word of the GPU's 64-bit radial sums), hot pixels in the rings (so the
// background clip rejects something), and unreadable pixels: some in the run's mask, some saturated,
// some unreadable on this frame alone (both read as an overload), inside spots as well as between them.
Frame RandomFrame(const DiffractionExperiment &x, uint32_t seed, const std::vector<size_t> &masked) {
const int W = static_cast<int>(x.GetXPixelsNum()), H = static_cast<int>(x.GetYPixelsNum());
const float beam_x = x.GetBeamX_pxl(), beam_y = x.GetBeamY_pxl();
std::mt19937 rng(seed);
std::uniform_real_distribution<float> U(0.0f, 1.0f);
Frame f;
f.image.resize(static_cast<size_t>(W) * H);
for (int y = 0; y < H; ++y)
for (int x_ = 0; x_ < W; ++x_) {
const float r = std::hypot(x_ - beam_x, y - beam_y);
const float mean = 0.5f + 0.01f * r;
f.image[static_cast<size_t>(y) * W + x_] = static_cast<int32_t>(std::lround(mean + (U(rng) - 0.5f) * 8.0f));
}
const int nrefl = W * H / 400;
int hkl = 1;
for (int n = 0; n < nrefl; ++n) {
Reflection r{};
r.h = hkl; r.k = -hkl; r.l = hkl % 7; ++hkl;
r.predicted_x = 2.0f + U(rng) * (W - 4);
r.predicted_y = 2.0f + U(rng) * (H - 4);
if (n % 4 == 1) {
r.predicted_x = std::round(r.predicted_x * 2.0f) / 2.0f;
r.predicted_y = std::round(r.predicted_y);
}
const float u = U(rng);
r.partiality = u < 0.3f ? 0.01f * u : u; // a third of them rocking-curve tails
r.d = 1.2f + 3.0f * U(rng);
r.prescaling_corr = 1.0f;
f.predicted.push_back(r);
if (U(rng) < 0.35f) {
const float amp = 20.0f + 3000.0f * U(rng) * U(rng);
const float s = 0.8f + 1.4f * U(rng);
for (int dy = -8; dy <= 8; ++dy)
for (int dx = -8; dx <= 8; ++dx) {
const int px = static_cast<int>(std::lround(r.predicted_x)) + dx;
const int py = static_cast<int>(std::lround(r.predicted_y)) + dy;
if (px < 0 || py < 0 || px >= W || py >= H) continue;
const float ex = px - r.predicted_x, ey = py - r.predicted_y;
f.image[static_cast<size_t>(py) * W + px] += static_cast<int32_t>(
std::lround(amp * std::exp(-(ex * ex + ey * ey) / (2.0f * s * s))));
}
}
}
std::uniform_int_distribution<size_t> P(0, f.image.size() - 1);
for (int k = 0; k < W * H / 500; ++k)
f.image[P(rng)] += 200 + static_cast<int32_t>(5000 * U(rng)); // hot pixels
for (int k = 0; k < W * H / 3000; ++k)
f.image[P(rng)] = U(rng) < 0.5f ? INT32_MIN : INT32_MAX; // unreadable on this frame
for (size_t i = 0; i < f.predicted.size(); i += 13) { // ... inside spots too
const auto &r = f.predicted[i];
const size_t idx = static_cast<size_t>(std::lround(r.predicted_y)) * W + std::lround(r.predicted_x) + i % 3;
if (idx < f.image.size()) f.image[idx] = i % 2 ? INT32_MIN : INT32_MAX;
}
for (size_t idx : masked)
f.image[idx] = INT32_MIN; // what the run's mask holds
return f;
}
struct Case {
IntegratorMode mode = IntegratorMode::ProfileGaussian;
OverlapMode overlap = OverlapMode::Exclude;
bool radial = false;
std::optional<float> bandwidth;
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, LinearQCentres(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).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;
fp.bin_px = 100.0f;
for (int b = 0; b < 8; ++b) {
fp.sigma_rad.push_back(0.9f + 0.35f * b);
fp.sigma_tan.push_back(0.8f + 0.3f * b);
}
s.Footprint(fp);
}
x.ImportBraggIntegrationSettings(s);
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, LinearQCentres(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) {
REQUIRE(metal.size() == cpu.size());
for (size_t i = 0; i < cpu.size(); ++i) {
const Reflection &a = cpu[i], &b = metal[i];
INFO("reflection " << i << " hkl " << a.h << " " << a.k << " " << a.l
<< " I " << a.I << " / " << b.I << " bkg " << a.bkg << " / " << b.bkg);
REQUIRE(a.h == b.h);
REQUIRE(a.k == b.k);
REQUIRE(a.l == b.l);
REQUIRE(Bits(a.image_number) == Bits(b.image_number));
REQUIRE(Bits(a.I) == Bits(b.I));
REQUIRE(Bits(a.sigma) == Bits(b.sigma));
REQUIRE(Bits(a.bkg) == Bits(b.bkg));
REQUIRE(Bits(a.var_bkg) == Bits(b.var_bkg));
REQUIRE(Bits(a.observed_x) == Bits(b.observed_x));
REQUIRE(Bits(a.observed_y) == Bits(b.observed_y));
REQUIRE(Bits(a.image_scale_corr) == Bits(b.image_scale_corr));
REQUIRE(Bits(a.prescaling_corr) == Bits(b.prescaling_corr));
REQUIRE(a.observed == b.observed);
REQUIRE(a.overloaded == b.overloaded);
REQUIRE(a.clipped == b.clipped);
}
}
void RequireSame(const BraggIntegrationCounts &a, const BraggIntegrationCounts &b) {
CHECK(a.predicted == b.predicted);
CHECK(a.bkg_starved == b.bkg_starved);
CHECK(a.bkg_starved_by_neighbour == b.bkg_starved_by_neighbour);
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
// masks, its buffers) has to be left as the CPU engine leaves its own - each compared exactly.
// Returns how many reflections came out, so a case can show it was not vacuous.
size_t CompareFrames(const DetectorSetup &det, const Case &c) {
const DiffractionExperiment x = MakeExperiment(det, c);
REQUIRE(BraggIntegrationEngineMetal::Available(x));
const size_t W = x.GetXPixelsNum();
// A few pixels of the run's mask, under spots and between them: unreadable there, but not overloads.
PixelMask mask(x);
std::vector<uint32_t> user(x.GetPixelsNum(), 0);
std::vector<size_t> masked;
for (size_t k = 0; k < x.GetPixelsNum(); k += 7919) {
user[k] = 1;
masked.push_back(k);
}
mask.LoadUserMask(x, user);
BraggIntegrationEngineCPU cpu(x, mask);
BraggIntegrationEngineMetal metal(x, mask);
if (c.curve) {
HandOverCurve(cpu);
HandOverCurve(metal);
}
ImagePreprocessorBuffer image(x.GetPixelsNum());
size_t n_out = 0;
for (uint32_t frame = 0; frame < 3; ++frame) {
INFO("frame " << frame);
const Frame f = RandomFrame(x, 1000 * frame + static_cast<uint32_t>(W), masked);
std::copy(f.image.begin(), f.image.end(), image.data());
// The auto radial mode switches per image; here it is on for the radial cases, every image.
cpu.BackgroundRadial(c.radial);
metal.BackgroundRadial(c.radial);
const auto out_cpu = cpu.Run(image, f.predicted, f.predicted.size(), frame);
const auto out_metal = metal.Run(image, f.predicted, f.predicted.size(), frame);
RequireSame(out_cpu, out_metal);
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;
}
} // namespace
TEST_CASE("BraggIntegrationEngineMetal_BitIdenticalToCPU", "[portable][Integration][Metal]") {
// A JUNGFRAU 1M, and a small odd-sized detector, where most boxes meet an edge.
const DetectorSetup large = DetJF(2);
const DetectorSetup small = DetDECTRIS(211, 197, "small", {});
std::vector<std::pair<std::string, Case>> cases;
for (IntegratorMode mode : {IntegratorMode::BoxSum, IntegratorMode::ProfileGaussian, IntegratorMode::ProfileEmpirical})
for (OverlapMode overlap : {OverlapMode::Off, OverlapMode::Reject, OverlapMode::Exclude})
for (bool radial : {false, true}) {
Case c;
c.mode = mode;
c.overlap = overlap;
c.radial = radial;
cases.emplace_back("mode " + std::to_string(static_cast<int>(mode)) + " overlap "
+ std::to_string(static_cast<int>(overlap)) + (radial ? " radial" : ""), c);
}
// Rings elongated by the bandwidth, and the measured footprint: per-reflection stencils.
Case broadband;
broadband.bandwidth = 0.01f;
broadband.stencil_k = 3.0f;
broadband.radial = true;
cases.emplace_back("broadband stencil radial", broadband);
Case footprint;
footprint.footprint = true;
footprint.radial = true;
cases.emplace_back("footprint radial", footprint);
Case footprint_box = footprint;
footprint_box.mode = IntegratorMode::BoxSum;
cases.emplace_back("footprint box sum", footprint_box);
// No clip: the plain ring mean, and no clip kernel at all.
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) {
CHECK(CompareFrames(large, c) > 1000);
CHECK(CompareFrames(small, c) > 50);
}
}
}
// 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, LinearQCentres(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]") {
Case trim;
trim.clip_nsigma = 0.0f;
DiffractionExperiment x = MakeExperiment(DetJF(2), trim);
BraggIntegrationSettings s = x.GetBraggIntegrationSettings();
s.BackgroundTrimFraction(0.1f);
x.ImportBraggIntegrationSettings(s);
CHECK(!BraggIntegrationEngineMetal::Available(x));
}
#endif
Reference in new issue View Git Blame Copy Permalink