From d6166f6fa5413fa86f549a27ba269c24cb6e4636 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 21 Oct 2026 10:40:20 +0200 Subject: [PATCH] Panel layout stage 2: prediction and Bragg integration through the panel table Everything sparse that took a radial quantity from an image coordinate now takes it from the reflection's (or spot's) position on the detector plane, through its panel (DetectorModules). With an empty module table every changed expression reduces to the old one in the same floating-point order, so flat detectors - every rotation set, EIGER, CCDs, assembled JUNGFRAU - integrate bit for bit as before. - Prediction (CPU, GPU, rotation CPU and GPU): bounds test against PlaneBounds instead of the image rectangle. - IntegrateOnImage: a prediction in a gap or off every panel is dropped (PlaneToPanel), no longer passed through onto an unrelated stored pixel. The panel id stays local: the engine and the way back find it from the image position (the panel boxes tile the image), and a prediction that rounding to float put on another panel's edge is dropped so the two agree. Nothing is added to Reflection. - Bragg engines (CPU, CUDA, Metal) take a BraggPanel per reflection (BraggStencil.h): its plane position, the panel's pixel axes on the plane and its box. The stencil's radial direction is taken on the plane and expressed in the panel's pixel axes; r0, the radial background curve, the Gaussian's elongation, the learned moments, tan(2theta) and Finalize's polarisation use the plane position; r_max / n_rad come from PlaneBounds. The neighbour mask and the owner map are clamped to the reflection's own panel; inside the walked box a pixel off the reflection's panel is read as a gap (unreadable, statically masked - not an overload), which is what the assembled image has there, so a disk cut by a panel edge is partial rather than silently "full". - Risk R4: the ProfileEmpirical grid is learned and applied in the plane's axes (the panel axes taken to the nearest quarter turn, BraggProfileCell), so turned panels pool one shape. Without it the stacked layout read up to 0.3 % off the assembled one in the harness. - SpotWidth / SpotFootprint: aperture bounds against the spot's panel box, beam distance, resolution and the shape's 2theta Jacobian on the plane, moments read in the panel's axes. - GetDetectorMaxQ_recipA, the far-corner radii (first-pass tilt report, metric-symmetry arm, footprint table range) from PlaneBounds; pre-scan spot resolution on the plane; LogHeaderBeamCenter: geometric centre of the panels' plane box, and the masked-pixel check on the panel the beam lands on (none when it is in a gap). Harness: PanelLayout_BraggIntegration is now a real [PanelLayout] check, over six settings (default, bandwidth-elongated ring, radial background, overlap exclude, box sum, empirical profile, all with polarisation), stacked vs assembled to 1e-5; PanelLayout_BraggIntegrationGPU runs the same on the CUDA engine. Its spots avoid exact half pixels: the grid centres on the rounded position, and a half rounds to the other neighbour on a turned panel. The Metal engine (.metal/.mm) mirrors the CPU change exactly but is UNCOMPILED here (no macOS): it needs a build and BraggIntegrationEngineMetal_BitIdenticalToCPU on a Mac. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu --- common/DiffractionExperiment.cpp | 8 +- image_analysis/MXAnalysisWithoutFPGA.cpp | 29 ++-- .../BraggIntegrationEngine.cpp | 22 ++- .../BraggIntegrationEngine.h | 12 +- .../BraggIntegrationEngineCPU.cpp | 127 +++++++++++------- .../BraggIntegrationEngineCPU.h | 13 +- .../BraggIntegrationEngineGPU.cu | 96 ++++++++----- .../BraggIntegrationEngineGPU.h | 2 + .../BraggIntegrationEngineMetal.metal | 55 +++++--- .../BraggIntegrationEngineMetal.mm | 19 ++- .../bragg_integration/BraggPanel.cpp | 18 +++ image_analysis/bragg_integration/BraggPanel.h | 16 +++ .../bragg_integration/BraggStencil.h | 49 ++++++- .../bragg_integration/CMakeLists.txt | 2 + .../bragg_integration/SpotFootprint.cpp | 15 ++- .../bragg_integration/SpotFootprint.h | 12 +- .../bragg_integration/SpotWidth.cpp | 51 ++++--- image_analysis/bragg_integration/SpotWidth.h | 11 +- .../bragg_prediction/BraggPrediction.cpp | 9 +- .../bragg_prediction/BraggPredictionGPU.cu | 9 +- .../bragg_prediction/BraggPredictionGPU.h | 3 +- .../bragg_prediction/BraggPredictionRot.cpp | 9 +- .../bragg_prediction/BraggPredictionRotGPU.cu | 9 +- .../bragg_prediction/BraggPredictionRotGPU.h | 3 +- rugnux/Rugnux.cpp | 15 ++- rugnux/RugnuxFirstPassAccept.cpp | 7 +- rugnux/RugnuxPasses.cpp | 11 +- rugnux/RugnuxPreScan.cpp | 17 ++- tests/BandwidthEstimateTest.cpp | 2 +- tests/PanelLayoutTest.cpp | 90 +++++++++++-- tests/RugnuxTest.cpp | 2 +- tests/SpotTilesTest.cpp | 4 +- 32 files changed, 531 insertions(+), 216 deletions(-) create mode 100644 image_analysis/bragg_integration/BraggPanel.cpp create mode 100644 image_analysis/bragg_integration/BraggPanel.h diff --git a/common/DiffractionExperiment.cpp b/common/DiffractionExperiment.cpp index db1ce112b..f81b52253 100644 --- a/common/DiffractionExperiment.cpp +++ b/common/DiffractionExperiment.cpp @@ -1419,13 +1419,13 @@ float DiffractionExperiment::GetDetectorMaxResolution_A() const { float DiffractionExperiment::GetDetectorMaxQ_recipA() const { const DiffractionGeometry geom = GetDiffractionGeometry(); - const auto width = static_cast(GetXPixelsNumConv()); - const auto height = static_cast(GetYPixelsNumConv()); + // The plane the panels span; without panels the image. + const PlaneBox box = GetModules().PlaneBounds(GetXPixelsNumConv(), GetYPixelsNumConv()); // The largest scattering angle sits at one of the detector corners, wherever the beam centre is. float q = 0.0f; - for (const float x: {0.0f, width}) - for (const float y: {0.0f, height}) + for (const float x: {static_cast(box.x0), static_cast(box.x1)}) + for (const float y: {static_cast(box.y0), static_cast(box.y1)}) q = std::max(q, geom.PxlToQ(x, y)); return q; } diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 0e50fbc00..b969a0092 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -513,20 +513,31 @@ std::vector MXAnalysisWithoutFPGA::IntegrateOnImage(const std::vecto return bragg_engine->Run(*preprocessor_buffer, predicted, npredicted, image_number); // Predictions are where the reflections fall on the detector plane; the engine reads the image, so - // each goes to where its module put those pixels, and comes back to the plane once integrated. - std::vector on_image(predicted.begin(), predicted.begin() + static_cast(npredicted)); - for (auto &r: on_image) { - const auto [x, y] = modules.PlaneToImage(r.predicted_x, r.predicted_y); - r.predicted_x = static_cast(x); - r.predicted_y = static_cast(y); + // each goes to where its panel put those pixels, and comes back to the plane once integrated. One in + // a gap, or off every panel, has no pixels and is not integrated. The engine finds a reflection's + // panel from its image position, as the way back below does; a position that rounding to float put + // on another panel's edge is dropped, so that the two agree. + std::vector on_image; + on_image.reserve(npredicted); + for (size_t i = 0; i < npredicted; i++) { + const auto p = modules.PlaneToPanel(predicted[i].predicted_x, predicted[i].predicted_y); + if (!p) + continue; + Reflection r = predicted[i]; + r.predicted_x = static_cast(p->x); + r.predicted_y = static_cast(p->y); + if (modules.ModuleOfImagePixel(r.predicted_x, r.predicted_y) != p->panel) + continue; + on_image.push_back(r); } - auto ret = bragg_engine->Run(*preprocessor_buffer, on_image, npredicted, image_number); + auto ret = bragg_engine->Run(*preprocessor_buffer, on_image, on_image.size(), image_number); for (auto &r: ret) { - const auto [px, py] = modules.ImageToPlane(r.predicted_x, r.predicted_y); + const int64_t panel = modules.ModuleOfImagePixel(r.predicted_x, r.predicted_y); + 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); if (std::isfinite(r.observed_x) && std::isfinite(r.observed_y)) { - const auto [ox, oy] = modules.ImageToPlane(r.observed_x, r.observed_y); + 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); } diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp index 28a40d186..bfcb4401f 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp @@ -51,7 +51,8 @@ std::map> radial_kernel_cache; } // namespace BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &experiment) - : geom(experiment.GetDiffractionGeometry()) { + : geom(experiment.GetDiffractionGeometry()), modules(experiment.GetModules()), + plane_bounds(modules.PlaneBounds(experiment.GetXPixelsNum(), experiment.GetYPixelsNum())) { const auto settings = experiment.GetBraggIntegrationSettings(); const auto &det = experiment.GetDetectorSetup(); @@ -148,8 +149,8 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe // bkg_radial is raised after construction only by the auto mode (MXAnalysisWithoutFPGA), which // requires bkg_radial_auto, and the GPU allocates its curve buffers under the same condition. // n_kern is the largest row BraggStencilKernelIndex can select, plus one. - r_max = std::hypot(std::max(beam_x, static_cast(xpixel) - beam_x), - std::max(beam_y, static_cast(ypixel) - beam_y)); + r_max = std::hypot(std::max(beam_x - plane_bounds.x0, plane_bounds.x1 - beam_x), + std::max(beam_y - plane_bounds.y0, plane_bounds.y1 - beam_y)); bkg_radial_built = bkg_radial || bkg_radial_auto; const float grow_max = bkg_radial_built ? BraggStencilMaxGrow_px(static_cast(r_max), stencil) : 0.0f; @@ -217,8 +218,20 @@ void BraggIntegrationEngine::BuildRadialKernel(float grow) { k_diff.push_back(static_cast(hist_ann[k] / sa - hist_disk[k] / sd)); } +std::vector BraggIntegrationEngine::Panels(const std::vector &predicted, + size_t npredicted) const { + std::vector out(npredicted); + for (size_t i = 0; i < npredicted; ++i) { + const float x = predicted[i].predicted_x, y = predicted[i].predicted_y; + out[i] = MakeBraggPanel(modules, modules.ModuleOfImagePixel(x, y), static_cast(xpixel), + static_cast(ypixel), x, y); + } + return out; +} + std::vector BraggIntegrationEngine::Finalize(const std::vector &predicted, size_t npredicted, + const std::vector &panels, const std::vector &results, int64_t image_number) const { // Sized to what is kept, not to what was predicted: this vector is the one every image holds to the @@ -247,7 +260,8 @@ std::vector BraggIntegrationEngine::Finalize(const std::vector(image_number); out.push_back(refl); diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.h b/image_analysis/bragg_integration/BraggIntegrationEngine.h index a5edef530..bec850887 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.h @@ -42,6 +42,7 @@ #include "../../common/DiffractionGeometry.h" #include "../../common/Reflection.h" #include "../image_preprocessing/ImagePreprocessorBuffer.h" +#include "BraggPanel.h" #include "BraggStencil.h" namespace bragg_engine { @@ -239,12 +240,20 @@ protected: DiffractionGeometry geom; // kept for the per-reflection polarization correction std::optional polarization; + // The panels of the image (DetectorModules), empty for a flat detector, and the box on the plane + // they span - where the far corner is, for anything sized by the largest radius on the detector. + DetectorModules modules; + PlaneBox plane_bounds; + // The panel of each prediction, found from its position in the image (BraggPanel). + std::vector Panels(const std::vector &predicted, size_t npredicted) const; + // Accumulated over every image this engine has run. An engine belongs to one worker thread, so // these are plain counters and the caller sums over the engines it made. BraggIntegrationCounts counts; // Assemble output reflections from the per-reflection fit results (polarization + scale corr). std::vector Finalize(const std::vector &predicted, size_t npredicted, + const std::vector &panels, const std::vector &results, int64_t image_number) const; @@ -253,7 +262,8 @@ public: virtual ~BraggIntegrationEngine() = default; // predicted[0..npredicted) are the reflections to extract; image is the preprocessed int32 - // frame (image.size() == npixel). Returns only the observed reflections. + // frame (image.size() == npixel). Returns only the observed reflections. Positions in and out are + // the image's: a reflection is on the panel whose box holds its predicted position. virtual std::vector Run(const ImagePreprocessorBuffer &image, const std::vector &predicted, size_t npredicted, int64_t image_number) = 0; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index bdf16f9ca..e99c37f81 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -7,6 +7,8 @@ #include #include #include +#include +#include #include #if defined(_MSC_VER) && (defined(_M_X64) || defined(_M_IX86)) @@ -23,6 +25,11 @@ namespace { // The engine reads pixels in the INT32_MIN(masked)/INT32_MAX(saturated) convention. inline bool valid(int32_t v) { return v != INT32_MIN && v != INT32_MAX; } +// A direction on the plane in a panel's pixel axes (MakeBraggStencil). Unchanged without panels. +inline std::pair to_panel_axes(const BraggPanel &pn, double ux, double uy) { + return {ux * pn.fx + uy * pn.fy, ux * pn.sx + uy * pn.sy}; +} + // Ask for the cache lines from first to last ahead of their use. Only a hint: it changes no result. inline void prefetch(const void *p) { #if defined(_MSC_VER) && (defined(_M_X64) || defined(_M_IX86)) @@ -85,24 +92,26 @@ BraggIntegrationEngineCPU::BraggIntegrationEngineCPU(const DiffractionExperiment } void BraggIntegrationEngineCPU::MarkNeighbours(const std::vector &predicted, + const std::vector &panels, const std::vector &order) { - const int W = static_cast(xpixel), H = static_cast(ypixel); // --- Reflection mask: mark the r2 signal region of every predicted reflection, so a neighbour's // signal is excluded from this reflection's r2..r3 background ring. MASK_FLUX where the // neighbour puts flux on this frame (NEIGHBOUR_MASK_MIN_PARTIALITY) - only those are excluded // - and MASK_TAIL where only the tail of a rocking curve does, which is counted for the density // guard but read as background. The region is the INNER stencil ellipse, taken in each // neighbour's OWN frame - an elongated neighbour whose streak is still marked as a disk would - // leak its tails into this reflection's ring. --- + // leak its tails into this reflection's ring. Only on the neighbour's own panel: beyond it the + // image holds pixels that are elsewhere on the detector. --- refl_mask.Clear(); for (const uint32_t i : order) { const auto &r = predicted[i]; + const BraggPanel &pn = panels[i]; const uint8_t level = r.partiality >= NEIGHBOUR_MASK_MIN_PARTIALITY ? MASK_FLUX : MASK_TAIL; - const BraggStencil st = MakeBraggStencil(r.predicted_x, r.predicted_y, stencil); - const int x0 = std::max(0, static_cast(std::floor(r.predicted_x - st.ex_in - 1.0f))); - const int x1 = std::min(W - 1, static_cast(std::ceil(r.predicted_x + st.ex_in + 1.0f))); - const int y0 = std::max(0, static_cast(std::floor(r.predicted_y - st.ey_in - 1.0f))); - const int y1 = std::min(H - 1, static_cast(std::ceil(r.predicted_y + st.ey_in + 1.0f))); + const BraggStencil st = MakeBraggStencil(pn, stencil); + const int x0 = std::max(pn.x0, static_cast(std::floor(r.predicted_x - st.ex_in - 1.0f))); + const int x1 = std::min(pn.x1, static_cast(std::ceil(r.predicted_x + st.ex_in + 1.0f))); + const int y0 = std::max(pn.y0, static_cast(std::floor(r.predicted_y - st.ey_in - 1.0f))); + const int y1 = std::min(pn.y1, static_cast(std::ceil(r.predicted_y + st.ey_in + 1.0f))); for (int y = y0; y <= y1; ++y) for (int x = x0; x <= x1; ++x) { const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y); @@ -119,10 +128,12 @@ void BraggIntegrationEngineCPU::CountNeighbourCrowding(const std::vector(xpixel), H = static_cast(ypixel); counts.predicted += npredicted; const std::vector order = FrameOrder(predicted, npredicted); - MarkNeighbours(predicted, order); + const std::vector panels = Panels(predicted, npredicted); + MarkNeighbours(predicted, panels, order); for (const uint32_t i : order) { const auto &r = predicted[i]; - const BraggStencil st = MakeBraggStencil(r.predicted_x, r.predicted_y, stencil); + const BraggPanel &pn = panels[i]; + const BraggStencil st = MakeBraggStencil(pn, stencil); const int x0 = std::max(0, static_cast(std::floor(r.predicted_x - st.ex_out - 1.0f))); const int x1 = std::min(W - 1, static_cast(std::ceil(r.predicted_x + st.ex_out + 1.0f))); const int y0 = std::max(0, static_cast(std::floor(r.predicted_y - st.ey_out - 1.0f))); @@ -135,7 +146,8 @@ void BraggIntegrationEngineCPU::CountNeighbourCrowding(const std::vector= r3_sq) continue; const size_t idx = static_cast(y) * W + x; - if (idx / 32 < static_mask.size() && ((static_mask[idx / 32] >> (idx % 32)) & 1U)) + if (!BraggPanelHolds(pn, x, y) + || (idx / 32 < static_mask.size() && ((static_mask[idx / 32] >> (idx % 32)) & 1U))) continue; ++n_ring; if (refl_mask.Get(x, y) == 0) @@ -151,7 +163,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, const std::vector &predicted, size_t npredicted, int64_t image_number) { if (npredicted == 0) - return Finalize(predicted, npredicted, std::vector(npredicted), image_number); + return Finalize(predicted, npredicted, {}, std::vector(npredicted), image_number); counts.predicted += npredicted; const int W = static_cast(xpixel), H = static_cast(ypixel); @@ -159,7 +171,8 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, // Every loop below whose results do not depend on the order walks the reflections in this one. const std::vector order = FrameOrder(predicted, npredicted); - MarkNeighbours(predicted, order); + const std::vector panels = Panels(predicted, npredicted); + MarkNeighbours(predicted, panels, order); owner.Clear(); // --- Signal-region ownership: a pixel inside two reflections' signal regions belongs to the @@ -173,10 +186,11 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, const float claim_sq = claim * claim; for (const uint32_t i : order) { const auto &r = predicted[i]; - const int x0 = std::max(0, static_cast(std::floor(r.predicted_x - claim))); - const int x1 = std::min(W - 1, static_cast(std::ceil(r.predicted_x + claim))); - const int y0 = std::max(0, static_cast(std::floor(r.predicted_y - claim))); - const int y1 = std::min(H - 1, static_cast(std::ceil(r.predicted_y + claim))); + const BraggPanel &pn = panels[i]; + const int x0 = std::max(pn.x0, static_cast(std::floor(r.predicted_x - claim))); + const int x1 = std::min(pn.x1, static_cast(std::ceil(r.predicted_x + claim))); + const int y0 = std::max(pn.y0, static_cast(std::floor(r.predicted_y - claim))); + const int y1 = std::min(pn.y1, static_cast(std::ceil(r.predicted_y + claim))); for (int y = y0; y <= y1; ++y) for (int x = x0; x <= x1; ++x) { const float dx = x - r.predicted_x, dy = y - r.predicted_y; @@ -224,15 +238,16 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, // prefetch.) if (k + 1 < npredicted) { const Reflection &n = predicted[order[k + 1]]; - const Box b = box_of(n, MakeBraggStencil(n.predicted_x, n.predicted_y, stencil)); + const Box b = box_of(n, MakeBraggStencil(panels[order[k + 1]], stencil)); for (int y = b.y0; y <= b.y1 && b.x0 <= b.x1; ++y) prefetch_row(img.address(static_cast(y) * W + b.x0), img.address(static_cast(y) * W + b.x1)); } const auto &r = predicted[i]; + const BraggPanel &pn = panels[i]; // This reflection's stencil: the r1 signal disk, and the r2..r3 ring elongated radially by // the analytic smear. - const BraggStencil st = MakeBraggStencil(r.predicted_x, r.predicted_y, stencil); + const BraggStencil st = MakeBraggStencil(pn, stencil); const Box box = box_of(r, st); const int x0 = box.x0, x1 = box.x1, y0 = box.y0, y1 = box.y1; @@ -247,7 +262,9 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, for (int y = y0; y <= y1; ++y) for (int x = x0; x <= x1; ++x) { const auto d = BraggStencilDistances(st, x - r.predicted_x, y - r.predicted_y); - const int32_t px = img[y * W + x]; + // Off the reflection's panel is a gap: unreadable, and not lost to the flux. + const bool on = BraggPanelHolds(pn, x, y); + const int32_t px = on ? img[y * W + x] : INT32_MIN; if (do_clip && d.inner >= r2_sq && d.outer < r3_sq && valid(px) && refl_mask.Get(x, y) != MASK_FLUX) clip_ring.push_back({px, d.rad}); if (BraggStencilInSignal(st, x - r.predicted_x, y - r.predicted_y, d.signal, r1_sq)) { @@ -262,7 +279,7 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, } ++s.n_inner; if (px == INT32_MAX - || (px == INT32_MIN && !static_masked(static_cast(y) * W + x))) + || (px == INT32_MIN && on && !static_masked(static_cast(y) * W + x))) s.overloaded = 1; if (!valid(px)) continue; s.I_sum += px; @@ -306,17 +323,18 @@ std::vector BraggIntegrationEngineCPU::RunImpl(const Sampler &img, FinishPassA(a, i, r, st, s, clip, bkg_vals); } - return ProfilePasses(img, owner, predicted, npredicted, a, image_number); + return ProfilePasses(img, owner, predicted, npredicted, panels, a, image_number); } BraggIntegrationEngineCPU::PassA BraggIntegrationEngineCPU::StartPassA(size_t npredicted) const { - const int W = static_cast(xpixel), H = static_cast(ypixel); PassA a; a.rough.resize(npredicted); // A pixel's radius is the reflection's radius plus the pixel's projection on the beam->reflection - // direction, so no per-pixel sqrt is needed. - a.n_rad = static_cast(std::ceil(std::hypot(std::max(beam_x, W - beam_x), - std::max(beam_y, H - beam_y)))) + 2; + // direction, so no per-pixel sqrt is needed. Out to the far corner of the plane the panels span. + const float bx0 = static_cast(plane_bounds.x0), bx1 = static_cast(plane_bounds.x1); + const float by0 = static_cast(plane_bounds.y0), by1 = static_cast(plane_bounds.y1); + a.n_rad = static_cast(std::ceil(std::hypot(std::max(beam_x - bx0, bx1 - beam_x), + std::max(beam_y - by0, by1 - beam_y)))) + 2; if (bkg_radial) { a.rad_sum.assign(a.n_rad, 0); a.rad_cnt.assign(a.n_rad, 0); @@ -469,7 +487,7 @@ bool BraggIntegrationEngineCPU::FitKept(const BraggFitOut &o) const { template BraggFitOut BraggIntegrationEngineCPU::FitProfile(const Sampler &img, const OwnerMap &owner_map, - const BraggFitSetup &s, size_t i, + const BraggFitSetup &s, const BraggPanel &pn, size_t i, const std::vector &shell_profiles, std::vector &cells, std::vector &fit_pixels) const { @@ -480,9 +498,15 @@ BraggFitOut BraggIntegrationEngineCPU::FitProfile(const Sampler &img, const Owne BraggFitOut o{}; const float *P; - if (s.shell >= 0) - P = shell_profiles.data() + static_cast(s.shell) * GG; - else { + if (s.shell >= 0) { + // The shell's learned profile is in the plane's axes; laid out here in the panel's. + const float *learned = shell_profiles.data() + static_cast(s.shell) * GG; + cells.resize(GG); + for (int dy = -R; dy <= R; ++dy) + for (int dx = -R; dx <= R; ++dx) + cells[(dy + R) * G + (dx + R)] = learned[BraggProfileCell(pn, dx, dy, R)]; + P = cells.data(); + } else { // The grid's image rows, asked for now (as in pass A): building the profile below takes long // enough for them to arrive before the fit reads them. const int gx0 = std::max(0, s.cx - Rf), gx1 = std::min(W - 1, s.cx + Rf); @@ -510,7 +534,7 @@ BraggFitOut BraggIntegrationEngineCPU::FitProfile(const Sampler &img, const Owne // pixels, what Reject cuts on. m_read / m_all is the fraction of the r1 disk the // summation seed the runaway guard compares against actually summed; with nothing missing // and nothing excluded it is 1 and the guard is untouched. --- - const BraggStencil st = MakeBraggStencil(s.px, s.py, stencil); + const BraggStencil st = MakeBraggStencil(pn, stencil); fit_pixels.clear(); for (int dy = -Rf; dy <= Rf; ++dy) for (int dx = -Rf; dx <= Rf; ++dx) { @@ -527,7 +551,7 @@ BraggFitOut BraggIntegrationEngineCPU::FitProfile(const Sampler &img, const Owne o.p_grid += Pp; o.p_peak = std::max(o.p_peak, Pp); if (in_disk) o.m_all += Pp; - const int32_t v = img[y * W + x]; + const int32_t v = BraggPanelHolds(pn, x, y) ? img[y * W + x] : INT32_MIN; if (!valid(v)) { o.p_lost_peak = std::max(o.p_lost_peak, Pp); continue; @@ -561,7 +585,8 @@ BraggFitOut BraggIntegrationEngineCPU::FitProfile(const Sampler &img, const Owne template std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler &img, const OwnerMap &owner_map, const std::vector &predicted, - size_t npredicted, PassA &a, + size_t npredicted, + const std::vector &panels, PassA &a, int64_t image_number) { const int W = static_cast(xpixel), H = static_cast(ypixel); std::vector results(npredicted); @@ -572,7 +597,6 @@ std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler & return !use_owner || BraggOwnedBy(owner_map.Get(x, y), static_cast(i)); }; - auto grid_idx = [this](int dx, int dy) { return (dy + R) * G + (dx + R); }; // --- Radial background curvature correction. The annulus mean is blind to the curvature of the // radial background (a linear background cancels between the concentric disk and annulus), so @@ -613,7 +637,7 @@ std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler & static_cast(rh.obs_x), static_cast(rh.obs_y), static_cast(rh.var_bkg), true, rh.has_obs, rh.overloaded}; } - return Finalize(predicted, npredicted, results, image_number); + return Finalize(predicted, npredicted, panels, results, image_number); } auto shell_of = [&](float d) { @@ -626,7 +650,7 @@ std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler & // --- Learn the profile per shell (+ global) from the strong spots. --- // Two things are learned. The empirical profile is the average grid in the DETECTOR frame, which is - // where it is applied. The width is a pair of second moments taken in each spot's OWN radial / + // where it is applied - in the plane's axes, which a turned panel's pixels are mapped onto. The width is a pair of second moments taken in each spot's OWN radial / // tangential frame: a grid stacked in the detector frame is azimuthally averaged, so its is // sigma_r^2 + sigma_t^2 with no way back, and a radially smeared spot reads as a wide TANGENTIAL // one. Rotating each contribution into the spot's frame keeps the two apart. @@ -641,19 +665,22 @@ std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler & for (size_t i = 0; i < npredicted; ++i) { const auto &rh = rough[i]; if (!rh.ok || !rh.strong || rh.I <= 0.0) continue; - const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y; + const BraggPanel &pn = panels[i]; + const double rx = pn.plane_x - beam_x, ry = pn.plane_y - beam_y; const double Rpx = std::hypot(rx, ry); - const double ux = Rpx > 1e-6 ? rx / Rpx : 1.0, uy = Rpx > 1e-6 ? ry / Rpx : 0.0; - for (int dy = -R; dy <= R; ++dy) + const auto [ux, uy] = to_panel_axes(pn, Rpx > 1e-6 ? rx / Rpx : 1.0, Rpx > 1e-6 ? ry / Rpx : 0.0); + for (int dy = -R; dy <= R; ++dy) for (int dx = -R; dx <= R; ++dx) { const int x = rh.cx + dx, y = rh.cy + dy; - if (x < 0 || y < 0 || x >= W || y >= H) continue; + if (x < 0 || y < 0 || x >= W || y >= H || !BraggPanelHolds(pn, x, y)) continue; const int32_t px = img[y * W + x]; if (!valid(px)) continue; if (exclude && !clean(x, y, i)) continue; const double v = (static_cast(px) - rh.bkg) / rh.I; - shell_grid[rh.shell][grid_idx(dx, dy)] += v; - global_grid[grid_idx(dx, dy)] += v; + // The grid in the plane's axes (BraggProfileCell); the moments below are axis-free. + const int cell = BraggProfileCell(pn, dx, dy, R); + shell_grid[rh.shell][cell] += v; + global_grid[cell] += v; if (dx * dx + dy * dy >= r1_sq) continue; const double rad = dx * ux + dy * uy, tn = -dx * uy + dy * ux; shell_mom[rh.shell].rad += v * rad * rad; @@ -734,7 +761,8 @@ std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler & continue; } const int sh = rh.shell < 0 ? 0 : rh.shell; - const BraggStencil st = MakeBraggStencil(predicted[i].predicted_x, predicted[i].predicted_y, stencil); + const BraggPanel &pn = panels[i]; + const BraggStencil st = MakeBraggStencil(pn, stencil); fs.fit = 1; fs.px = predicted[i].predicted_x; fs.py = predicted[i].predicted_y; @@ -744,7 +772,7 @@ std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler & fs.shell = sh; fs.footprint = st.fp_s2t > 0.0f; if (!empirical) { - const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y; + const double rx = pn.plane_x - beam_x, ry = pn.plane_y - beam_y; const double Rpx = std::hypot(rx, ry); const double tan2t = Rpx / F_px; double s2t = shell_sigma2[sh].tan; @@ -758,7 +786,7 @@ std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler & const double radial_extra = std::max(shell_sigma2[sh].rad - s2t, sbw * sbw + c_radial * tan2t * tan2t); if (Rpx > 1e-6 && radial_extra > 0.25) { - ux = rx / Rpx; uy = ry / Rpx; + std::tie(ux, uy) = to_panel_axes(pn, rx / Rpx, ry / Rpx); s2r = s2t + radial_extra; elong = true; } @@ -796,7 +824,7 @@ std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler & std::vector fit_pixels; for (const uint32_t i : FrameOrder(predicted, npredicted)) if (setups[i].fit) - fits[i] = FitProfile(img, owner_map, setups[i], i, shell_profiles, cells_buf, fit_pixels); + fits[i] = FitProfile(img, owner_map, setups[i], panels[i], i, shell_profiles, cells_buf, fit_pixels); } for (size_t i = 0; i < npredicted; ++i) { @@ -841,18 +869,19 @@ std::vector BraggIntegrationEngineCPU::ProfilePasses(const Sampler & static_cast(var_bkg), true, rh.has_obs, rh.overloaded}; } - return Finalize(predicted, npredicted, results, image_number); + return Finalize(predicted, npredicted, panels, results, image_number); } template std::vector BraggIntegrationEngineCPU::ProfilePasses( const BraggIntegrationEngineCPU::BufferSampler &, const BraggIntegrationEngineCPU::FlatOwnerMap &, - const std::vector &, size_t, BraggIntegrationEngineCPU::PassA &, int64_t); + const std::vector &, size_t, const std::vector &, + BraggIntegrationEngineCPU::PassA &, int64_t); std::vector BraggIntegrationEngineCPU::Run(const ImagePreprocessorBuffer &image, const std::vector &predicted, size_t npredicted, int64_t image_number) { if (image.size() != npixel) - return Finalize(predicted, npredicted, std::vector(npredicted), image_number); + return Finalize(predicted, npredicted, {}, std::vector(npredicted), image_number); return RunImpl(BufferSampler{image.data()}, predicted, npredicted, image_number); } @@ -860,7 +889,7 @@ std::vector BraggIntegrationEngineCPU::Run(const CompressedImage &im const std::vector &predicted, size_t npredicted, int64_t image_number) { if (image.GetWidth() * image.GetHeight() != npixel) - return Finalize(predicted, npredicted, std::vector(npredicted), image_number); + return Finalize(predicted, npredicted, {}, std::vector(npredicted), image_number); std::vector scratch; const auto *ptr = image.GetUncompressedPtr(scratch); diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h index 7ad908cf5..237fa4ad3 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.h @@ -56,9 +56,10 @@ class BraggIntegrationEngineCPU : public BraggIntegrationEngine { TiledFrame refl_mask; TiledFrame owner; - // Marks every predicted reflection's r2 region in refl_mask (MASK_FLUX / MASK_TAIL), visiting them in - // `order` (any order gives the same mask). - void MarkNeighbours(const std::vector &predicted, const std::vector &order); + // Marks every predicted reflection's r2 region in refl_mask (MASK_FLUX / MASK_TAIL), on its own panel, + // visiting them in `order` (any order gives the same mask). + void MarkNeighbours(const std::vector &predicted, const std::vector &panels, + const std::vector &order); template std::vector RunImpl(const Sampler &img, const std::vector &predicted, @@ -69,7 +70,8 @@ class BraggIntegrationEngineCPU : public BraggIntegrationEngine { // fit. fit_pixels is scratch: the pixels the fit reads, gathered once. struct FitPixel { float Pp, signal; }; template - BraggFitOut FitProfile(const Sampler &img, const OwnerMap &owner_map, const BraggFitSetup &s, size_t i, + BraggFitOut FitProfile(const Sampler &img, const OwnerMap &owner_map, const BraggFitSetup &s, + const BraggPanel &pn, size_t i, const std::vector &shell_profiles, std::vector &cells, std::vector &fit_pixels) const; // Whether enough of a reflection's profile is readable, and its own, for the fit to be kept. @@ -140,7 +142,8 @@ protected: template std::vector ProfilePasses(const Sampler &img, const OwnerMap &owner_map, const std::vector &predicted, size_t npredicted, - PassA &a, int64_t image_number); + const std::vector &panels, PassA &a, + int64_t image_number); public: BraggIntegrationEngineCPU(const DiffractionExperiment &experiment, const PixelMask &mask); diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index c024db81c..2c1428552 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -67,18 +67,20 @@ constexpr int COUNT_SLOTS = 5; // `value` is what the launch writes: MASK_TAIL from every reflection, which also claims the owner map; // then MASK_FLUX from the ones that put flux on this frame (`marks`), a later launch so that it wins // over the first where two regions meet, as the CPU engine's maximum does; and 0 to clear. -__global__ void mark_mask(const float *px_x, const float *px_y, const uint8_t *marks, uint8_t *mask, - uint32_t *owner, BraggGpuParams p, int n, int value) { +__global__ void mark_mask(const float *px_x, const float *px_y, const BraggPanel *panel, const uint8_t *marks, + uint8_t *mask, uint32_t *owner, BraggGpuParams p, int n, int value) { const int i = blockIdx.x; if (i >= n) return; if (value == MASK_FLUX && !marks[i]) return; const bool clear = value == 0; const float cx = px_x[i], cy = px_y[i]; - const BraggStencil st = MakeBraggStencil(cx, cy, p.stencil); - const int x0 = max(0, (int) floorf(cx - st.ex_in - 1.0f)); - const int x1 = min(p.W - 1, (int) ceilf(cx + st.ex_in + 1.0f)); - const int y0 = max(0, (int) floorf(cy - st.ey_in - 1.0f)); - const int y1 = min(p.H - 1, (int) ceilf(cy + st.ey_in + 1.0f)); + // Only on the reflection's own panel (see the CPU engine). + const BraggPanel pn = panel[i]; + const BraggStencil st = MakeBraggStencil(pn, p.stencil); + const int x0 = max(pn.x0, (int) floorf(cx - st.ex_in - 1.0f)); + const int x1 = min(pn.x1, (int) ceilf(cx + st.ex_in + 1.0f)); + const int y0 = max(pn.y0, (int) floorf(cy - st.ey_in - 1.0f)); + const int y1 = min(pn.y1, (int) ceilf(cy + st.ey_in + 1.0f)); const int bw = x1 - x0 + 1, bh = y1 - y0 + 1; if (bw <= 0 || bh <= 0) return; for (int t = threadIdx.x; t < bw * bh; t += blockDim.x) { @@ -159,7 +161,7 @@ __device__ __forceinline__ float profile_float(unsigned long long v) { } // --- Pass A box-sum: rough I / background / centroid / strong flag, one block per reflection. --- -__global__ void boxsum(const float *px_x, const float *px_y, const float *dd, +__global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *panel, const float *dd, const int32_t *img, const uint8_t *mask, const uint32_t *owner, const uint8_t *static_mask, BraggGpuParams p, int n, int *cx_o, int *cy_o, float *I_o, float *sigma_o, float *bkg_o, @@ -194,8 +196,9 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, unsigned long long *s_radv = s_rad; // signed value carried as two's complement int *s_radn = (int *) (s_rad + p.rad_w); __shared__ int s_radbase; + const BraggPanel pn = panel[i]; if (threadIdx.x == 0) { - const float rx = px_x[i] - p.beam_x, ry = px_y[i] - p.beam_y; + const float rx = pn.plane_x - p.beam_x, ry = pn.plane_y - p.beam_y; s_r0 = sqrtf(rx * rx + ry * ry); s_Isum = 0; s_Ix = 0; s_Iy = 0; s_x = 0; s_y = 0; s_ninner = 0; s_ninner_valid = 0; s_nbkg = 0; s_bkgsum = 0; @@ -206,7 +209,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, const float cx = px_x[i], cy = px_y[i]; // Every thread builds the same stencil from the same inputs; it is a few flops against a whole // bounding box of pixel reads, so it costs less than staging it through shared memory. - const BraggStencil st = MakeBraggStencil(cx, cy, p.stencil); + const BraggStencil st = MakeBraggStencil(pn, p.stencil); const int x0 = max(0, (int) floorf(cx - st.ex_out - 1.0f)); const int x1 = min(p.W - 1, (int) ceilf(cx + st.ex_out + 1.0f)); const int y0 = max(0, (int) floorf(cy - st.ey_out - 1.0f)); @@ -232,12 +235,15 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, else if (p.exclude_px) continue; } ++l_ni; - const int32_t px = img[y * p.W + x]; + // Off the reflection's panel is a gap: unreadable, and not lost to the flux (see the CPU engine). + const bool on = BraggPanelHolds(pn, x, y); + const int32_t px = on ? img[y * p.W + x] : INT32_MIN; // Saturated, or unreadable on this frame beyond the run's mask: an overload (see the CPU engine). - if (px == INT32_MAX || (px == INT32_MIN && !static_mask[y * p.W + x])) ++l_nsat; + if (px == INT32_MAX || (px == INT32_MIN && on && !static_mask[y * p.W + x])) ++l_nsat; if (valid(px)) { l_Isum += px; l_Ix += (long long) x * px; l_Iy += (long long) y * px; l_x += x; l_y += y; ++l_niv; } } else if (d.inner >= p.r2_sq && d.outer < p.r3_sq) { + if (!BraggPanelHolds(pn, x, y)) continue; const int32_t px = img[y * p.W + x]; if (!valid(px)) continue; ++l_nring; @@ -307,7 +313,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, const int x = x0 + t % bw, y = y0 + t / bw; const BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy); if (!(d.inner >= p.r2_sq && d.outer < p.r3_sq)) continue; - if (mask[y * p.W + x] == MASK_FLUX) continue; + if (!BraggPanelHolds(pn, x, y) || mask[y * p.W + x] == MASK_FLUX) continue; const int32_t px = img[y * p.W + x]; if (!valid(px)) continue; const int slot = atomicAdd(&s_bn, 1); @@ -349,7 +355,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const float *dd, const int x = x0 + t % bw, y = y0 + t / bw; const BraggStencilDist d = BraggStencilDistances(st, (float) x - cx, (float) y - cy); if (!(d.inner >= p.r2_sq && d.outer < p.r3_sq)) continue; - if (mask[y * p.W + x] == MASK_FLUX) continue; + if (!BraggPanelHolds(pn, x, y) || mask[y * p.W + x] == MASK_FLUX) continue; const int32_t px = img[y * p.W + x]; if (!valid(px)) continue; if ((long long) px <= s_thr_i) { @@ -482,8 +488,7 @@ __global__ void reset(unsigned long long *shell_grid, unsigned long long *global // the spot's own radial/tangential frame, because a detector-frame stack is azimuthally averaged // and cannot tell a radially smeared spot from a tangentially wide one (see the CPU engine). // One block per reflection. --- -__global__ void learn_profile(const int32_t *img, const uint32_t *owner, - const float *px_x, const float *px_y, +__global__ void learn_profile(const int32_t *img, const uint32_t *owner, const BraggPanel *panel, const int *cx_a, const int *cy_a, const float *dd, const unsigned long long *invd2mm, const float *I_a, const float *bkg_a, const uint8_t *ok_a, const uint8_t *strong_a, @@ -502,9 +507,12 @@ __global__ void learn_profile(const int32_t *img, const uint32_t *owner, __syncthreads(); const int sh = s_sh; const float bkg = bkg_a[i]; - const float rx = px_x[i] - p.beam_x, ry = px_y[i] - p.beam_y; + const BraggPanel pn = panel[i]; + const float rx = pn.plane_x - p.beam_x, ry = pn.plane_y - p.beam_y; const float Rpx = sqrtf(rx * rx + ry * ry); - const float ux = Rpx > 1e-6f ? rx / Rpx : 1.0f, uy = Rpx > 1e-6f ? ry / Rpx : 0.0f; + // The radial direction on the plane, in the panel's pixel axes (see the CPU engine). + const float pux = Rpx > 1e-6f ? rx / Rpx : 1.0f, puy = Rpx > 1e-6f ? ry / Rpx : 0.0f; + const float ux = pux * pn.fx + puy * pn.fy, uy = pux * pn.sx + puy * pn.sy; unsigned long long *sg = shell_grid + (size_t) sh * p.GG; __shared__ unsigned long long s_rad, s_tan, s_w; if (threadIdx.x == 0) { s_rad = 0; s_tan = 0; s_w = 0; } @@ -513,14 +521,16 @@ __global__ void learn_profile(const int32_t *img, const uint32_t *owner, for (int k = threadIdx.x; k < p.GG; k += blockDim.x) { const int dx = k % p.G - p.R, dy = k / p.G - p.R; const int x = cx + dx, y = cy + dy; - if (x < 0 || y < 0 || x >= p.W || y >= p.H) continue; + if (x < 0 || y < 0 || x >= p.W || y >= p.H || !BraggPanelHolds(pn, x, y)) continue; const int32_t px = img[y * p.W + x]; if (!valid(px)) continue; if (p.overlap == 2 && !BraggOwnedBy(owner[y * p.W + x], i)) continue; const float v = ((float) px - bkg) / I; const unsigned long long vq = profile_fixed(v); - atomicAdd(&sg[k], vq); - atomicAdd(&global_grid[k], vq); + // The grid in the plane's axes (BraggProfileCell, see the CPU engine). + const int cell = BraggProfileCell(pn, dx, dy, p.R); + atomicAdd(&sg[cell], vq); + atomicAdd(&global_grid[cell], vq); if ((float) (dx * dx + dy * dy) >= p.r1_sq) continue; const float rad = dx * ux + dy * uy, tn = -dx * uy + dy * ux; l_rad += v * rad * rad; l_tan += v * tn * tn; l_w += v; @@ -618,6 +628,7 @@ __global__ void radial_correct(const unsigned long long *rad_sum, const int *rad // Otwinowski & Minor, Methods Enzymol 276, 307-326 (1997). // One block per reflection; the (possibly elongated) profile is built in shared memory. --- __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x, const float *px_y, + const BraggPanel *panel, const int *cx_a, const int *cy_a, const float *dd, const unsigned long long *invd2mm, const float *I_seed, const float *sigma_seed, const float *bkg_a, const float *bkgvar_a, const float *varbkg_seed, const uint8_t *ok_a, const uint8_t *overloaded_a, @@ -640,7 +651,8 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x if (overloaded_a[i]) return; const int cx = cx_a[i], cy = cy_a[i]; - const BraggStencil st = MakeBraggStencil(px_x[i], px_y[i], p.stencil); + const BraggPanel pn = panel[i]; + const BraggStencil st = MakeBraggStencil(pn, p.stencil); // As in learn_profile: one thread per block, not one per thread. Two double-precision divisions // on a card whose double throughput is a sixty-fourth of its single is worth doing once. __shared__ int s_sh; @@ -654,12 +666,14 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x const float *Psrc = use_shell ? (shell_P + (size_t) sh * p.GG) : global_P; if (threadIdx.x == 0) { s_Rf = p.R; s_Gf = p.G; } __syncthreads(); - for (int k = threadIdx.x; k < p.GG; k += blockDim.x) Pbuf[k] = Psrc[k]; + // The learned profile is in the plane's axes; laid out here in the panel's. + for (int k = threadIdx.x; k < p.GG; k += blockDim.x) + Pbuf[k] = Psrc[BraggProfileCell(pn, k % p.G - p.R, k / p.G - p.R, p.R)]; __syncthreads(); } else { const int si = use_shell ? sh : N_SHELL; // N_SHELL is the global slot float s2t = sigma2_t[si]; - const float rx = px_x[i] - p.beam_x, ry = px_y[i] - p.beam_y; + const float rx = pn.plane_x - p.beam_x, ry = pn.plane_y - p.beam_y; const float Rpx = sqrtf(rx * rx + ry * ry); const float tan2t = Rpx / p.F_px; float s2r = s2t, ux = 1.0f, uy = 0.0f; @@ -669,7 +683,14 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x // parallax/capture term (see the CPU engine). const float sbw = p.bw_sigma * Rpx; const float radial_extra = fmaxf(sigma2_r[si] - s2t, sbw * sbw + p.c_radial * tan2t * tan2t); - if (Rpx > 1e-6f && radial_extra > 0.25f) { ux = rx / Rpx; uy = ry / Rpx; s2r = s2t + radial_extra; elong = true; } + if (Rpx > 1e-6f && radial_extra > 0.25f) { + // On the plane, in the panel's pixel axes (see the CPU engine). + const float pux = rx / Rpx, puy = ry / Rpx; + ux = pux * pn.fx + puy * pn.fy; + uy = pux * pn.sx + puy * pn.sy; + s2r = s2t + radial_extra; + elong = true; + } } // A spot that outgrows the r1 disk takes the measured footprint's widths (see the CPU engine). if (st.fp_s2t > 0.0f) { @@ -728,7 +749,7 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x l_pgrid += Pp; l_ppeak = fmaxf(l_ppeak, Pp); if (in_disk) l_mall += Pp; - if (!valid(img[y * p.W + x])) { + if (!BraggPanelHolds(pn, x, y) || !valid(img[y * p.W + x])) { l_plost = fmaxf(l_plost, Pp); continue; } @@ -772,7 +793,7 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x const float Pp = Pbuf[k]; if (Pp <= 0.0f) continue; const int x = cx + (k % Gf - Rf), y = cy + (k / Gf - Rf); - if (x < 0 || y < 0 || x >= p.W || y >= p.H) continue; + if (x < 0 || y < 0 || x >= p.W || y >= p.H || !BraggPanelHolds(pn, x, y)) continue; const int32_t px = img[y * p.W + x]; if (!valid(px)) continue; const float v = fmaxf(B + Ihere * Pp, (float) WEIGHT_VARIANCE_MIN_FRACTION * B); @@ -905,6 +926,7 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) { d_px_x = CudaDevicePtr(new_capacity); d_px_y = CudaDevicePtr(new_capacity); d_d = CudaDevicePtr(new_capacity); + d_panel = CudaDevicePtr(new_capacity); d_mark = CudaDevicePtr(new_capacity); d_cx = CudaDevicePtr(new_capacity); d_cy = CudaDevicePtr(new_capacity); @@ -927,6 +949,7 @@ void BraggIntegrationEngineGPU::EnsureCapacity(size_t n) { h_px_x = CudaHostPtr(new_capacity); h_px_y = CudaHostPtr(new_capacity); h_d = CudaHostPtr(new_capacity); + h_panel = CudaHostPtr(new_capacity); h_mark = CudaHostPtr(new_capacity); h_I = CudaHostPtr(new_capacity); h_sigma = CudaHostPtr(new_capacity); @@ -945,7 +968,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu size_t npredicted, int64_t image_number) { std::vector results(npredicted); if (image.size() != npixel || npredicted == 0) - return Finalize(predicted, npredicted, results, image_number); + return Finalize(predicted, npredicted, {}, results, image_number); counts.predicted += npredicted; const int32_t *img = image.getGPUBuffer(); @@ -955,7 +978,9 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu EnsureCapacity(npredicted); const int n = static_cast(npredicted); + const std::vector panels = Panels(predicted, npredicted); for (size_t i = 0; i < npredicted; ++i) { + h_panel[i] = panels[i]; h_px_x[i] = predicted[i].predicted_x; h_px_y[i] = predicted[i].predicted_y; h_d[i] = predicted[i].d; @@ -964,6 +989,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu cuda_err(cudaMemcpyAsync(d_px_x, h_px_x.get(), sizeof(float) * npredicted, cudaMemcpyHostToDevice, *stream)); cuda_err(cudaMemcpyAsync(d_px_y, h_px_y.get(), sizeof(float) * npredicted, cudaMemcpyHostToDevice, *stream)); cuda_err(cudaMemcpyAsync(d_d, h_d.get(), sizeof(float) * npredicted, cudaMemcpyHostToDevice, *stream)); + cuda_err(cudaMemcpyAsync(d_panel, h_panel.get(), sizeof(BraggPanel) * npredicted, cudaMemcpyHostToDevice, *stream)); cuda_err(cudaMemcpyAsync(d_mark, h_mark.get(), sizeof(uint8_t) * npredicted, cudaMemcpyHostToDevice, *stream)); BraggGpuParams p{ @@ -1015,11 +1041,11 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu reset<<<32, 256, 0, *stream>>>(d_shell_grid, d_global_grid, d_mom, d_shell_n, d_global_n, d_invd2, GG); cuda_err(cudaGetLastError()); dirty = true; - mark_mask<<>>(d_px_x, d_px_y, d_mark, d_mask, d_owner, p, n, MASK_TAIL); + mark_mask<<>>(d_px_x, d_px_y, d_panel, d_mark, d_mask, d_owner, p, n, MASK_TAIL); cuda_err(cudaGetLastError()); - mark_mask<<>>(d_px_x, d_px_y, d_mark, d_mask, d_owner, p, n, MASK_FLUX); + mark_mask<<>>(d_px_x, d_px_y, d_panel, d_mark, d_mask, d_owner, p, n, MASK_FLUX); cuda_err(cudaGetLastError()); - boxsum<<>>(d_px_x, d_px_y, d_d, img, d_mask, d_owner, + boxsum<<>>(d_px_x, d_px_y, d_panel, d_d, img, d_mask, d_owner, d_static_mask->get(), p, n, d_cx, d_cy, d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_obs_x, d_obs_y, d_ok, d_strong, d_has_obs, d_overloaded, d_invd2, @@ -1038,7 +1064,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu if (mode != IntegratorMode::BoxSum) { // Pass B: learn (shell computed inline) -> build -> fit. - learn_profile<<>>(img, d_owner, d_px_x, d_px_y, d_cx, d_cy, d_d, d_invd2, + learn_profile<<>>(img, d_owner, d_panel, d_cx, d_cy, d_d, d_invd2, d_I, d_bkg, d_ok, d_strong, d_shell_grid, d_global_grid, d_mom, d_shell_n, d_global_n, p, n); cuda_err(cudaGetLastError()); @@ -1046,7 +1072,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu d_shell_grid, d_global_grid, d_mom, d_shell_n, d_global_n, d_shell_P, d_global_P, d_sigma2_r, d_sigma2_t, p); cuda_err(cudaGetLastError()); - fit<<>>(img, d_owner, d_px_x, d_px_y, d_cx, d_cy, d_d, d_invd2, + fit<<>>(img, d_owner, d_px_x, d_px_y, d_panel, d_cx, d_cy, d_d, d_invd2, d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_ok, d_overloaded, d_shell_P, d_global_P, d_sigma2_r, d_sigma2_t, d_shell_n, @@ -1057,7 +1083,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu // Nothing below reads the mask or the owner map, so put them back now, over the boxes that were // marked rather than over the frame - where that is the cheaper of the two. if (clear_by_box) { - mark_mask<<>>(d_px_x, d_px_y, d_mark, d_mask, d_owner, p, n, 0); + mark_mask<<>>(d_px_x, d_px_y, d_panel, d_mark, d_mask, d_owner, p, n, 0); cuda_err(cudaGetLastError()); } @@ -1091,7 +1117,7 @@ std::vector BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu results[i].has_observed = true; } } - return Finalize(predicted, npredicted, results, image_number); + return Finalize(predicted, npredicted, panels, results, image_number); } BraggIntegrationCounts BraggIntegrationEngineGPU::Counts() const { diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h index 3e0cdd2f3..8a4de8253 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.h @@ -38,6 +38,7 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { // --- per-reflection device arrays (grown by EnsureCapacity) --- CudaDevicePtr d_px_x, d_px_y, d_d; + CudaDevicePtr d_panel; // the panel of each reflection (BraggStencil.h) CudaDevicePtr d_mark; // the reflection marks its signal region in d_mask CudaDevicePtr d_cx, d_cy; CudaDevicePtr d_I, d_sigma, d_bkg, d_bkg_var, d_var_bkg, d_obs_x, d_obs_y; @@ -78,6 +79,7 @@ class BraggIntegrationEngineGPU : public BraggIntegrationEngine { // driver has staged it through a bounce buffer, so eight of them in a row are eight serialised // round-trips rather than eight queued transfers. CudaHostPtr h_px_x, h_px_y, h_d; + CudaHostPtr h_panel; CudaHostPtr h_mark; CudaHostPtr h_I, h_sigma, h_bkg, h_var_bkg, h_obs_x, h_obs_y; CudaHostPtr h_ok, h_has_obs, h_overloaded; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineMetal.metal b/image_analysis/bragg_integration/BraggIntegrationEngineMetal.metal index 4d2632fa3..46d12341f 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineMetal.metal +++ b/image_analysis/bragg_integration/BraggIntegrationEngineMetal.metal @@ -31,16 +31,19 @@ kernel void mark_neighbours(device const float2 *xy [[buffer(0)]], constant char *stencil_raw [[buffer(3)]], constant BraggPassAFrame &f [[buffer(4)]], constant uint8_t &want [[buffer(5)]], + device const BraggPanel *panels [[buffer(6)]], uint i [[thread_position_in_grid]]) { if (i >= f.n || (want != 0 && level[i] != want)) return; const BraggStencilParams sp = stencil_params(stencil_raw); const float px = xy[i].x, py = xy[i].y; - const BraggStencil st = MakeBraggStencil(px, py, sp); - const int x0 = max(0, (int) floor(px - st.ex_in - 1.0f)); - const int x1 = min(f.width - 1, (int) ceil(px + st.ex_in + 1.0f)); - const int y0 = max(0, (int) floor(py - st.ey_in - 1.0f)); - const int y1 = min(f.height - 1, (int) ceil(py + st.ey_in + 1.0f)); + // Only on the reflection's own panel, as the CPU engine marks. + const BraggPanel pn = panels[i]; + const BraggStencil st = MakeBraggStencil(pn, sp); + const int x0 = max(pn.x0, (int) floor(px - st.ex_in - 1.0f)); + const int x1 = min(pn.x1, (int) ceil(px + st.ex_in + 1.0f)); + const int y0 = max(pn.y0, (int) floor(py - st.ey_in - 1.0f)); + const int y1 = min(pn.y1, (int) ceil(py + st.ey_in + 1.0f)); for (int y = y0; y <= y1; ++y) for (int x = x0; x <= x1; ++x) { const BraggStencilDist d = BraggStencilDistances(st, x - px, y - py); @@ -56,15 +59,17 @@ kernel void claim_owner(device const float2 *xy [[buffer(0)]], device atomic_uint *owner [[buffer(1)]], constant BraggPassAFrame &f [[buffer(2)]], constant uint &release [[buffer(3)]], + device const BraggPanel *panels [[buffer(4)]], uint i [[thread_position_in_grid]]) { if (i >= f.n) return; const float px = xy[i].x, py = xy[i].y; const float claim_sq = f.claim * f.claim; - const int x0 = max(0, (int) floor(px - f.claim)); - const int x1 = min(f.width - 1, (int) ceil(px + f.claim)); - const int y0 = max(0, (int) floor(py - f.claim)); - const int y1 = min(f.height - 1, (int) ceil(py + f.claim)); + const BraggPanel pn = panels[i]; + const int x0 = max(pn.x0, (int) floor(px - f.claim)); + const int x1 = min(pn.x1, (int) ceil(px + f.claim)); + const int y0 = max(pn.y0, (int) floor(py - f.claim)); + const int y1 = min(pn.y1, (int) ceil(py + f.claim)); for (int y = y0; y <= y1; ++y) for (int x = x0; x <= x1; ++x) { const float dx = x - px, dy = y - py; @@ -88,12 +93,14 @@ kernel void box_sums(device const float2 *xy [[buffer(0)]], constant char *stencil_raw [[buffer(5)]], constant BraggPassAFrame &f [[buffer(6)]], device BraggBoxSums *sums [[buffer(7)]], + device const BraggPanel *panels [[buffer(8)]], uint i [[thread_position_in_grid]]) { if (i >= f.n) return; const BraggStencilParams sp = stencil_params(stencil_raw); const float px = xy[i].x, py = xy[i].y; - const BraggStencil st = MakeBraggStencil(px, py, sp); + const BraggPanel pn = panels[i]; + const BraggStencil st = MakeBraggStencil(pn, sp); const int x0 = max(0, (int) floor(px - st.ex_out - 1.0f)); const int x1 = min(f.width - 1, (int) ceil(px + st.ex_out + 1.0f)); const int y0 = max(0, (int) floor(py - st.ey_out - 1.0f)); @@ -103,7 +110,9 @@ kernel void box_sums(device const float2 *xy [[buffer(0)]], for (int x = x0; x <= x1; ++x) { const BraggStencilDist d = BraggStencilDistances(st, x - px, y - py); const int idx = y * f.width + x; - const int32_t v = img[idx]; + // Off the reflection's panel is a gap: unreadable, and not lost to the flux. + const bool on = BraggPanelHolds(pn, x, y); + const int32_t v = on ? img[idx] : PX_MASKED; if (BraggStencilInSignal(st, x - px, y - py, d.signal, f.r1_sq)) { ++s.n_disk; if (f.use_owner) { @@ -114,7 +123,7 @@ kernel void box_sums(device const float2 *xy [[buffer(0)]], // Unreadable here but not in the run's mask: lost to the flux on this frame. const bool masked = (uint) idx / 32 < f.static_mask_words && ((static_mask[(uint) idx / 32] >> ((uint) idx % 32)) & 1u); - if (v == PX_SATURATED || (v == PX_MASKED && !masked)) + if (v == PX_SATURATED || (v == PX_MASKED && on && !masked)) s.overloaded = 1; if (!valid_px(v)) continue; s.I_sum += v; @@ -157,12 +166,14 @@ kernel void clip_sums(device const float2 *xy [[buffer(0)]], device BraggClipSums *out [[buffer(7)]], device atomic_uint *rad_sum [[buffer(8)]], device atomic_int *rad_cnt [[buffer(9)]], + device const BraggPanel *panels [[buffer(10)]], uint i [[thread_position_in_grid]]) { if (i >= f.n || !clips[i]) return; const BraggStencilParams sp = stencil_params(stencil_raw); const float px = xy[i].x, py = xy[i].y; - const BraggStencil st = MakeBraggStencil(px, py, sp); + const BraggPanel pn = panels[i]; + const BraggStencil st = MakeBraggStencil(pn, sp); const int x0 = max(0, (int) floor(px - st.ex_out - 1.0f)); const int x1 = min(f.width - 1, (int) ceil(px + st.ex_out + 1.0f)); const int y0 = max(0, (int) floor(py - st.ey_out - 1.0f)); @@ -173,7 +184,7 @@ kernel void clip_sums(device const float2 *xy [[buffer(0)]], for (int x = x0; x <= x1; ++x) { const BraggStencilDist d = BraggStencilDistances(st, x - px, y - py); const int idx = y * f.width + x; - const int32_t v = img[idx]; + const int32_t v = BraggPanelHolds(pn, x, y) ? img[idx] : PX_MASKED; if (!(d.inner >= f.r2_sq && d.outer < f.r3_sq && valid_px(v) && mask[idx] != MASK_FLUX && v <= t)) continue; c.sum += v; @@ -198,6 +209,7 @@ kernel void profile_fit(device const BraggFitSetup *setups [[buffer(0)]], constant char *stencil_raw [[buffer(5)]], constant BraggPassBFrame &f [[buffer(6)]], device BraggFitOut *out [[buffer(7)]], + device const BraggPanel *panels [[buffer(8)]], uint i [[thread_position_in_grid]]) { if (i >= f.n || !setups[i].fit) return; @@ -205,8 +217,11 @@ kernel void profile_fit(device const BraggFitSetup *setups [[buffer(0)]], const int Rf = s.Rf, Gf = 2 * Rf + 1; BraggFitOut o = {}; + // A learned profile is in the plane's axes (BraggProfileCell); a Gaussian is built in the panel's. + const BraggPanel pn = panels[i]; + const bool learned = s.shell >= 0; device const float *P; - if (s.shell >= 0) + if (learned) P = shell_profiles + s.shell * (2 * f.R + 1) * (2 * f.R + 1); else { device float *g = cells + s.cells; @@ -222,10 +237,10 @@ kernel void profile_fit(device const BraggFitSetup *setups [[buffer(0)]], P = g; } - const BraggStencil st = MakeBraggStencil(s.px, s.py, stencil_params(stencil_raw)); + const BraggStencil st = MakeBraggStencil(pn, stencil_params(stencil_raw)); for (int dy = -Rf; dy <= Rf; ++dy) for (int dx = -Rf; dx <= Rf; ++dx) { - const float Pp = P[(dy + Rf) * Gf + (dx + Rf)]; + const float Pp = P[learned ? BraggProfileCell(pn, dx, dy, Rf) : (dy + Rf) * Gf + (dx + Rf)]; if (Pp <= 0.0f) continue; const int x = s.cx + dx, y = s.cy + dy; if (x < 0 || y < 0 || x >= f.width || y >= f.height) continue; @@ -236,7 +251,7 @@ kernel void profile_fit(device const BraggFitSetup *setups [[buffer(0)]], o.p_grid += Pp; o.p_peak = o.p_peak < Pp ? Pp : o.p_peak; // std::max if (in_disk) o.m_all += Pp; - if (!valid_px(img[y * f.width + x])) { + if (!BraggPanelHolds(pn, x, y) || !valid_px(img[y * f.width + x])) { o.p_lost_peak = o.p_lost_peak < Pp ? Pp : o.p_lost_peak; continue; } @@ -254,10 +269,10 @@ kernel void profile_fit(device const BraggFitSetup *setups [[buffer(0)]], o.wsum = 0.0f; for (int dy = -Rf; dy <= Rf; ++dy) for (int dx = -Rf; dx <= Rf; ++dx) { - const float Pp = P[(dy + Rf) * Gf + (dx + Rf)]; + const float Pp = P[learned ? BraggProfileCell(pn, dx, dy, Rf) : (dy + Rf) * Gf + (dx + Rf)]; if (Pp <= 0.0f) continue; const int x = s.cx + dx, y = s.cy + dy; - if (x < 0 || y < 0 || x >= f.width || y >= f.height) continue; + if (x < 0 || y < 0 || x >= f.width || y >= f.height || !BraggPanelHolds(pn, x, y)) continue; const int32_t v = img[y * f.width + x]; if (!valid_px(v)) continue; if (f.exclude && f.use_owner && !BraggOwnedBy(owner[y * f.width + x], (int) i)) continue; diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineMetal.mm b/image_analysis/bragg_integration/BraggIntegrationEngineMetal.mm index 85ffafbe7..d086742e2 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineMetal.mm +++ b/image_analysis/bragg_integration/BraggIntegrationEngineMetal.mm @@ -109,7 +109,7 @@ struct BraggIntegrationEngineMetal::Device { id rad_sum = nil, rad_cnt = nil; // the radial background curve, as the clip kernel adds it up // Per reflection, grown with the frames. size_t capacity = 0; - id xy = nil, level = nil, sums = nil, clips = nil, thr = nil, clip_out = nil; + id xy = nil, panels = nil, level = nil, sums = nil, clips = nil, thr = nil, clip_out = nil; id fit_setups = nil, fit_out = nil; // Pass B's Gaussians, grown with the frames, and the learned profiles. size_t cells_capacity = 0, profiles_capacity = 0; @@ -168,7 +168,7 @@ std::vector BraggIntegrationEngineMetal::Run(const ImagePreprocessor counts.predicted += npredicted; const bool use_owner = overlap != OverlapMode::Off; - if (dev->unmark != nil) { // it reads xy, which is overwritten below + if (dev->unmark != nil) { // it reads xy and the panels, which are overwritten below id unmark = dev->unmark; dev->unmark = nil; dev->clear = false; // unless it completed @@ -188,6 +188,7 @@ std::vector BraggIntegrationEngineMetal::Run(const ImagePreprocessor if (dev->capacity < npredicted) { const size_t n = npredicted + npredicted / 2; dev->xy = buffer(n * 2 * sizeof(float)); + dev->panels = buffer(n * sizeof(BraggPanel)); dev->level = buffer(n); dev->sums = buffer(n * sizeof(BraggBoxSums)); dev->clips = buffer(n); @@ -200,6 +201,8 @@ std::vector BraggIntegrationEngineMetal::Run(const ImagePreprocessor auto *xy = static_cast(dev->xy.contents); auto *level = static_cast(dev->level.contents); + const std::vector panels = Panels(predicted, npredicted); + std::copy(panels.begin(), panels.end(), static_cast(dev->panels.contents)); for (size_t i = 0; i < npredicted; ++i) { xy[2 * i] = predicted[i].predicted_x; xy[2 * i + 1] = predicted[i].predicted_y; @@ -243,6 +246,7 @@ std::vector BraggIntegrationEngineMetal::Run(const ImagePreprocessor [enc setBytes:&stencil length:sizeof(stencil) atIndex:3]; [enc setBytes:&f length:sizeof(f) atIndex:4]; [enc setBytes:&want length:sizeof(want) atIndex:5]; + [enc setBuffer:dev->panels offset:0 atIndex:6]; [enc dispatchThreads:grid threadsPerThreadgroup:group]; } if (use_owner) { @@ -251,6 +255,7 @@ std::vector BraggIntegrationEngineMetal::Run(const ImagePreprocessor [enc setBuffer:dev->owner offset:0 atIndex:1]; [enc setBytes:&f length:sizeof(f) atIndex:2]; [enc setBytes:&claim_not_release length:sizeof(claim_not_release) atIndex:3]; + [enc setBuffer:dev->panels offset:0 atIndex:4]; [enc dispatchThreads:grid threadsPerThreadgroup:group]; } [enc setComputePipelineState:k.box]; @@ -262,6 +267,7 @@ std::vector BraggIntegrationEngineMetal::Run(const ImagePreprocessor [enc setBytes:&stencil length:sizeof(stencil) atIndex:5]; [enc setBytes:&f length:sizeof(f) atIndex:6]; [enc setBuffer:dev->sums offset:0 atIndex:7]; + [enc setBuffer:dev->panels offset:0 atIndex:8]; [enc dispatchThreads:grid threadsPerThreadgroup:group]; [enc endEncoding]; Wait(cb); @@ -302,6 +308,7 @@ std::vector BraggIntegrationEngineMetal::Run(const ImagePreprocessor [enc setBuffer:dev->clip_out offset:0 atIndex:7]; [enc setBuffer:(bkg_radial ? dev->rad_sum : dev->clips) offset:0 atIndex:8]; // not written without the curve [enc setBuffer:(bkg_radial ? dev->rad_cnt : dev->clips) offset:0 atIndex:9]; + [enc setBuffer:dev->panels offset:0 atIndex:10]; [enc dispatchThreads:grid threadsPerThreadgroup:group]; [enc endEncoding]; Wait(cb); @@ -322,12 +329,13 @@ std::vector BraggIntegrationEngineMetal::Run(const ImagePreprocessor std::vector no_trim; // the trimmed mean is the CPU engine's alone (Available) for (size_t i = 0; i < npredicted; ++i) { const Reflection &r = predicted[i]; - FinishPassA(a, i, r, MakeBraggStencil(r.predicted_x, r.predicted_y, stencil), sums[i], + FinishPassA(a, i, r, MakeBraggStencil(panels[i], stencil), sums[i], clips[i] ? clip_out[i] : BraggClipSums{}, no_trim); } const FlatOwnerMap owner_map{use_owner ? static_cast(dev->owner.contents) : nullptr, xpixel}; dev->img = img; - auto ret = ProfilePasses(BufferSampler{image.data()}, owner_map, predicted, npredicted, a, image_number); + auto ret = ProfilePasses(BufferSampler{image.data()}, owner_map, predicted, npredicted, panels, a, + image_number); dev->img = nil; // --- Unmark this frame's regions and give back its claims, while the CPU goes on --- @@ -341,6 +349,7 @@ std::vector BraggIntegrationEngineMetal::Run(const ImagePreprocessor [enc setBytes:&stencil length:sizeof(stencil) atIndex:3]; [enc setBytes:&f length:sizeof(f) atIndex:4]; [enc setBytes:&unmark length:sizeof(unmark) atIndex:5]; + [enc setBuffer:dev->panels offset:0 atIndex:6]; [enc dispatchThreads:grid threadsPerThreadgroup:group]; if (use_owner) { const uint32_t release = 1; @@ -349,6 +358,7 @@ std::vector BraggIntegrationEngineMetal::Run(const ImagePreprocessor [enc setBuffer:dev->owner offset:0 atIndex:1]; [enc setBytes:&f length:sizeof(f) atIndex:2]; [enc setBytes:&release length:sizeof(release) atIndex:3]; + [enc setBuffer:dev->panels offset:0 atIndex:4]; [enc dispatchThreads:grid threadsPerThreadgroup:group]; } [enc endEncoding]; @@ -395,6 +405,7 @@ bool BraggIntegrationEngineMetal::FitOnDevice(const std::vector & [enc setBytes:&stencil length:sizeof(stencil) atIndex:5]; [enc setBytes:&f length:sizeof(f) atIndex:6]; [enc setBuffer:dev->fit_out offset:0 atIndex:7]; + [enc setBuffer:dev->panels offset:0 atIndex:8]; [enc dispatchThreads:MTLSizeMake(n, 1, 1) threadsPerThreadgroup:MTLSizeMake(64, 1, 1)]; [enc endEncoding]; Wait(cb); diff --git a/image_analysis/bragg_integration/BraggPanel.cpp b/image_analysis/bragg_integration/BraggPanel.cpp new file mode 100644 index 000000000..876067b18 --- /dev/null +++ b/image_analysis/bragg_integration/BraggPanel.cpp @@ -0,0 +1,18 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "BraggPanel.h" + +BraggPanel MakeBraggPanel(const DetectorModules &modules, int64_t panel, int64_t width, int64_t height, + float x, float y) { + if (modules.empty()) + return {x, y, 1.0f, 0.0f, 0.0f, 1.0f, 0, 0, static_cast(width - 1), static_cast(height - 1)}; + if (panel < 0) + return {x, y, 1.0f, 0.0f, 0.0f, 1.0f, 0, 0, -1, -1}; + const DetectorModule &m = modules[panel]; + const auto [px, py] = m.ImageToPlane(x, y); + return {static_cast(px), static_cast(py), + static_cast(m.fx), static_cast(m.fy), static_cast(m.sx), static_cast(m.sy), + static_cast(m.x0), static_cast(m.y0), + static_cast(m.x0 + m.width - 1), static_cast(m.y0 + m.height - 1)}; +} diff --git a/image_analysis/bragg_integration/BraggPanel.h b/image_analysis/bragg_integration/BraggPanel.h new file mode 100644 index 000000000..9866a16a8 --- /dev/null +++ b/image_analysis/bragg_integration/BraggPanel.h @@ -0,0 +1,16 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#include "../../common/DetectorModules.h" +#include "BraggStencil.h" + +// Panel `panel` of the module table, for a point at (x, y) of the image (width x height): its plane +// position, the panel's pixel axes on the plane and its box in the image (BraggPanel). With an empty +// table the image is one flat panel, the plane the image itself. Panel -1 of a table that has panels - +// a point in no panel's box - is a panel with no pixels. +BraggPanel MakeBraggPanel(const DetectorModules &modules, int64_t panel, int64_t width, int64_t height, + float x, float y); diff --git a/image_analysis/bragg_integration/BraggStencil.h b/image_analysis/bragg_integration/BraggStencil.h index 00030799d..2dc44b8ca 100644 --- a/image_analysis/bragg_integration/BraggStencil.h +++ b/image_analysis/bragg_integration/BraggStencil.h @@ -171,15 +171,40 @@ BRAGG_STENCIL_HD float BraggStencilMaxGrow_px(float r_max, const BRAGG_STENCIL_T return g; } -BRAGG_STENCIL_HD BraggStencil MakeBraggStencil(float px_x, float px_y, const BRAGG_STENCIL_THREAD BraggStencilParams &p) { +// The panel a reflection is on (DetectorModules), as the stencil needs it. Everything radial about a +// reflection - its distance from the beam, the direction the ring is elongated in, the radial +// background - is a property of the detector PLANE, while the pixels it is walked over are the image's. +// A panel's image pixels lie on the plane at a fixed turn (and maybe mirror), so the radial direction is +// taken on the plane and expressed in the panel's own pixel axes; a panel is rigid, so a distance in its +// pixels is the same distance on the plane. A pixel outside the panel's box is not the reflection's +// neighbour on the plane, whatever is stored there: the engines read it as a gap. +// Without panels the image is the plane: the whole image, plane == image, and every expression below +// is the one it always was. +struct BraggPanel { + float plane_x, plane_y; // the reflection on the plane + float fx, fy, sx, sy; // the plane step of one image pixel along x, and of one along y + int32_t x0, y0, x1, y1; // the panel's pixels in the image, inclusive +}; + +BRAGG_STENCIL_HD bool BraggPanelHolds(const BRAGG_STENCIL_THREAD BraggPanel &pn, int x, int y) { + return x >= pn.x0 && x <= pn.x1 && y >= pn.y0 && y <= pn.y1; +} + +BRAGG_STENCIL_HD BraggStencil MakeBraggStencil(const BRAGG_STENCIL_THREAD BraggPanel &pn, + const BRAGG_STENCIL_THREAD BraggStencilParams &p) { BraggStencil s; - const float rx = px_x - p.beam_x, ry = px_y - p.beam_y; + const float rx = pn.plane_x - p.beam_x, ry = pn.plane_y - p.beam_y; const float r0 = sqrtf(rx * rx + ry * ry); s.r0 = r0; if (r0 > 1e-6f) { s.ux = rx / r0; s.uy = ry / r0; } + // The radial direction in the panel's pixel axes. A pixel step (dx, dy) is f dx + s dy on the plane, + // so its radial part is dx (u.f) + dy (u.s). The identity of a flat detector leaves u as it is. + const float ux = s.ux, uy = s.uy; + s.ux = ux * pn.fx + uy * pn.fy; + s.uy = ux * pn.sx + uy * pn.sy; s.grow = BraggStencilGrow_px(r0, p); // A spot wider than the r1 disk can resolve keeps its background ring clear of itself: the ring @@ -223,6 +248,26 @@ BRAGG_STENCIL_HD BraggStencil MakeBraggStencil(float px_x, float px_y, const BRA return s; } +// The cell of a learned profile (ProfileEmpirical) for grid offset (dx, dy) in a panel's pixels, on a grid +// of half-size R. The learned profile is kept in the plane's axes, so that reflections on panels turned +// against each other pool one shape rather than four orientations of it. A grid of whole pixels can only +// follow a quarter turn (and a mirror), so each panel axis is taken as the plane axis nearest to it; the +// small tilt of a placed panel is left out. Without panels, the cell of (dx, dy) itself. +BRAGG_STENCIL_HD int BraggProfileCell(const BRAGG_STENCIL_THREAD BraggPanel &pn, int dx, int dy, int R) { + const bool f_on_x = pn.fx * pn.fx >= pn.fy * pn.fy, s_on_x = pn.sx * pn.sx > pn.sy * pn.sy; + const int f = f_on_x ? (pn.fx > 0.0f ? 1 : -1) : (pn.fy > 0.0f ? 1 : -1); + const int s = s_on_x ? (pn.sx > 0.0f ? 1 : -1) : (pn.sy > 0.0f ? 1 : -1); + const int px = (f_on_x ? f * dx : 0) + (s_on_x ? s * dy : 0); + const int py = (f_on_x ? 0 : f * dx) + (s_on_x ? 0 : s * dy); + return (py + R) * (2 * R + 1) + (px + R); +} + +// A reflection at (x, y) of a flat detector, where the image is the plane. +BRAGG_STENCIL_HD BraggStencil MakeBraggStencil(float x, float y, const BRAGG_STENCIL_THREAD BraggStencilParams &p) { + const BraggPanel pn = {x, y, 1.0f, 0.0f, 0.0f, 1.0f, 0, 0, 0, 0}; + return MakeBraggStencil(pn, p); +} + // Which of the n_kern radial-background kernels this reflection uses: one per whole pixel of growth. // Taken from the stencil, so the kernel and the ring it corrects are built from the same radius - // deriving it from a separately computed r0 let the two engines round differently and pick diff --git a/image_analysis/bragg_integration/CMakeLists.txt b/image_analysis/bragg_integration/CMakeLists.txt index 2a74ecd82..a52ca23ea 100644 --- a/image_analysis/bragg_integration/CMakeLists.txt +++ b/image_analysis/bragg_integration/CMakeLists.txt @@ -3,6 +3,8 @@ ADD_LIBRARY(JFJochBraggIntegration STATIC BraggIntegrationEngine.h BraggIntegrationEngineCPU.cpp BraggIntegrationEngineCPU.h + BraggPanel.cpp + BraggPanel.h BraggPassA.h BraggPassB.h BraggExp.h diff --git a/image_analysis/bragg_integration/SpotFootprint.cpp b/image_analysis/bragg_integration/SpotFootprint.cpp index 88a327f24..67a5de7be 100644 --- a/image_analysis/bragg_integration/SpotFootprint.cpp +++ b/image_analysis/bragg_integration/SpotFootprint.cpp @@ -30,16 +30,18 @@ float median_of(std::vector &v) { } // namespace -FootprintOutcome MeasureFootprintSpot(const FootprintTile &tile, int width, int height, float beam_x, float beam_y, +FootprintOutcome MeasureFootprintSpot(const FootprintTile &tile, const BraggPanel &pn, float beam_x, float beam_y, float x, float y, FootprintSpot &spot) { const int half = static_cast(std::ceil(BKG_OUTER * WINDOW_MAX_PX)) + 1; std::vector ring; - const float rx = x - beam_x, ry = y - beam_y; + // The distance from the beam on the plane, and the radial direction in the panel's pixel axes. + const float rx = pn.plane_x - beam_x, ry = pn.plane_y - beam_y; const float r = std::sqrt(rx * rx + ry * ry); if (!(r > 1.0f)) return FootprintOutcome::LeftOut; - const float ux = rx / r, uy = ry / r; + const float ux = rx / r * pn.fx + ry / r * pn.fy, uy = rx / r * pn.sx + ry / r * pn.sy; const int ix = static_cast(std::lround(x)), iy = static_cast(std::lround(y)); - if (ix - half < 0 || iy - half < 0 || ix + half >= width || iy + half >= height) return FootprintOutcome::LeftOut; + if (ix - half < pn.x0 || iy - half < pn.y0 || ix + half > pn.x1 || iy + half > pn.y1) + return FootprintOutcome::LeftOut; // Start from a compact spot at the prediction; each round re-centres on the signal and takes // the window to three of the widths just measured. @@ -51,7 +53,7 @@ FootprintOutcome MeasureFootprintSpot(const FootprintTile &tile, int width, int const float reach = BKG_OUTER * std::max(wr, wt); const int x0 = static_cast(std::floor(cx - reach)), x1 = static_cast(std::ceil(cx + reach)); const int y0 = static_cast(std::floor(cy - reach)), y1 = static_cast(std::ceil(cy + reach)); - if (x0 < 0 || y0 < 0 || x1 >= width || y1 >= height) return FootprintOutcome::LeftOut; + if (x0 < pn.x0 || y0 < pn.y0 || x1 > pn.x1 || y1 > pn.y1) return FootprintOutcome::LeftOut; // The re-centring is not bounded, so the window can walk off a tile the frame still holds. if (x0 < tile.x0 || y0 < tile.y0 || x1 >= tile.x0 + tile.w || y1 >= tile.y0 + tile.h) return FootprintOutcome::OffTile; @@ -104,7 +106,8 @@ void MeasureFootprintSpots(const int32_t *img, int width, int height, float beam const FootprintTile frame{img, 0, 0, width, height}; for (size_t s = 0; s < x.size(); ++s) { FootprintSpot spot; - if (MeasureFootprintSpot(frame, width, height, beam_x, beam_y, x[s], y[s], spot) == FootprintOutcome::Measured) + const BraggPanel pn = {x[s], y[s], 1.0f, 0.0f, 0.0f, 1.0f, 0, 0, width - 1, height - 1}; + if (MeasureFootprintSpot(frame, pn, beam_x, beam_y, x[s], y[s], spot) == FootprintOutcome::Measured) out.push_back(spot); } } diff --git a/image_analysis/bragg_integration/SpotFootprint.h b/image_analysis/bragg_integration/SpotFootprint.h index d549abbba..f84f03935 100644 --- a/image_analysis/bragg_integration/SpotFootprint.h +++ b/image_analysis/bragg_integration/SpotFootprint.h @@ -42,6 +42,7 @@ #include "../../common/BraggIntegrationSettings.h" #include "../../common/Reflection.h" #include "../../common/SpotToSave.h" +#include "BraggStencil.h" // One spot's measured widths, and its distance from the beam centre [px]. struct FootprintSpot { @@ -51,7 +52,7 @@ struct FootprintSpot { }; // Measure the spots at the given (sub-pixel) positions on one preprocessed frame (INT32_MIN masked, -// INT32_MAX saturated). A spot whose window holds an unreadable pixel, or no signal above its +// INT32_MAX saturated) of a flat detector. A spot whose window holds an unreadable pixel, or no signal above its // background, is left out. Results are appended in the order of the positions. void MeasureFootprintSpots(const int32_t *img, int width, int height, float beam_x, float beam_y, const std::vector &x, const std::vector &y, @@ -65,10 +66,11 @@ struct FootprintTile { }; enum class FootprintOutcome { Measured, LeftOut, OffTile }; -// One spot of MeasureFootprintSpots, reading the frame (width x height) through `tile`. OffTile: the -// window walked off the tile while still on the frame, so the spot has to be measured again, from -// the start, on a tile that holds more of the frame - the answer is then the one the whole frame gives. -FootprintOutcome MeasureFootprintSpot(const FootprintTile &tile, int width, int height, float beam_x, float beam_y, +// One spot at (x, y) of the image, on panel `pn` (BraggPanel: its plane position, the panel's pixel axes +// and its box), reading the frame through `tile`. The window stays on the panel; OffTile: it walked off +// the tile while still on the panel, so the spot has to be measured again, from the start, on a tile +// that holds more of the frame - the answer is then the one the whole frame gives. +FootprintOutcome MeasureFootprintSpot(const FootprintTile &tile, const BraggPanel &pn, float beam_x, float beam_y, float x, float y, FootprintSpot &spot); // Tabulate the spots: the median widths in bins of distance from the beam centre out to r_max. A bin diff --git a/image_analysis/bragg_integration/SpotWidth.cpp b/image_analysis/bragg_integration/SpotWidth.cpp index 327f1896c..7c83fc243 100644 --- a/image_analysis/bragg_integration/SpotWidth.cpp +++ b/image_analysis/bragg_integration/SpotWidth.cpp @@ -14,6 +14,7 @@ #include "../../common/JFJochMath.h" #include "../SensorAbsorption.h" +#include "BraggPanel.h" using namespace spot_width; @@ -109,7 +110,8 @@ constexpr int R_MOMENT = 9; // every radius of the disk - an extra second moment that does not cancel between the two directions // once they are weighted by the obliquity. SpotShape shape_of(const int32_t *centre_px, int stride, int cx, int cy, double mx, double my, - const std::vector &ring, const DiffractionGeometry &geometry) { + const std::vector &ring, const DiffractionGeometry &geometry, + const DetectorModules &modules, int64_t panel, const BraggPanel &pn) { SpotShape shape; double sum = 0.0; for (const int32_t v : ring) sum += v; @@ -157,8 +159,10 @@ SpotShape shape_of(const int32_t *centre_px, int stride, int cx, int cy, double } if (!(flux > 0.0)) return shape; - // Where one radian of 2theta, and of the angle across the scattering plane, moves the spot. - const float x = static_cast(cx + mx), y = static_cast(cy + my); + // Where one radian of 2theta, and of the angle across the scattering plane, moves the spot - on + // the plane, where the geometry is. + const auto [plane_x, plane_y] = modules.ImageToPlane(panel, cx + mx, cy + my); + const float x = static_cast(plane_x), y = static_cast(plane_y); const Coord s0(0, 0, 1); const Coord s1 = geometry.LabCoord(x, y).Normalize(); const float two_theta = std::acos(std::clamp(s1 * s0, -1.0f, 1.0f)); @@ -198,8 +202,11 @@ SpotShape shape_of(const int32_t *centre_px, int stride, int cx, int cy, double shape.two_theta = two_theta; shape.jr = static_cast(jr); shape.jt = static_cast(jt); - shape.m_rad = static_cast(ux * ux * sxx + 2.0 * ux * uy * sxy + uy * uy * syy); - shape.m_tan = static_cast(vx * vx * sxx + 2.0 * vx * vy * sxy + vy * vy * syy); + // The moments were taken in the panel's pixel axes, so the two directions are read there. + const double iux = ux * pn.fx + uy * pn.fy, iuy = ux * pn.sx + uy * pn.sy; + const double ivx = vx * pn.fx + vy * pn.fy, ivy = vx * pn.sx + vy * pn.sy; + shape.m_rad = static_cast(iux * iux * sxx + 2.0 * iux * iuy * sxy + iuy * iuy * syy); + shape.m_tan = static_cast(ivx * ivx * sxx + 2.0 * ivx * ivy * sxy + ivy * ivy * syy); shape.cos_psi = cos_psi; shape.par_u = static_cast(par_u); shape.par_v = static_cast(par_v); @@ -219,7 +226,7 @@ double median_of(std::vector &v) { } // namespace void MeasureSpotFluxCurves(const ImagePreprocessorBuffer &image, int width, int height, - const DiffractionGeometry &geometry, + const DiffractionGeometry &geometry, const DetectorModules &modules, const std::vector &spots, std::vector &out) { const int32_t *pixels = image.data(); @@ -228,11 +235,11 @@ void MeasureSpotFluxCurves(const ImagePreprocessorBuffer &image, int width, int const Coord c = spots[i].RawCoord(); return pixels + static_cast(std::lround(c.y)) * width + std::lround(c.x); }; - MeasureSpotFluxCurves(centre_px, width, width, height, geometry, spots, out); + MeasureSpotFluxCurves(centre_px, width, width, height, geometry, modules, spots, out); } void MeasureSpotFluxCurves(const std::function ¢re_px, int stride, int width, int height, - const DiffractionGeometry &geometry, + const DiffractionGeometry &geometry, const DetectorModules &modules, const std::vector &spots, std::vector &out) { if (spots.empty()) return; @@ -284,23 +291,26 @@ void MeasureSpotFluxCurves(const std::function ¢re_ return true; }; - // Candidates that pass the geometric gates, by band, strongest first. - struct Candidate { size_t index; int64_t count; float d_A; }; + // Candidates that pass the geometric gates, by band, strongest first. The aperture has to lie on the + // spot's own panel; distances from the beam and the resolution are taken on the plane. + struct Candidate { size_t index; int64_t count; float d_A; int64_t panel; BraggPanel pn; }; std::array, N_BAND> candidates; constexpr double MIN_BEAM_DISTANCE_PX2 = static_cast(MIN_BEAM_DISTANCE_PX) * MIN_BEAM_DISTANCE_PX; for (size_t i = 0; i < spots.size(); i++) { const Coord &c = centre[i]; const int cx = static_cast(std::lround(c.x)), cy = static_cast(std::lround(c.y)); - if (cx < R_BKG_OUT || cy < R_BKG_OUT || cx >= width - R_BKG_OUT || cy >= height - R_BKG_OUT) + const int64_t panel = modules.ModuleOfImagePixel(c.x, c.y); + const BraggPanel pn = MakeBraggPanel(modules, panel, width, height, c.x, c.y); + if (cx - R_BKG_OUT < pn.x0 || cy - R_BKG_OUT < pn.y0 || cx + R_BKG_OUT > pn.x1 || cy + R_BKG_OUT > pn.y1) continue; - const double bx = c.x - beam_x, by = c.y - beam_y; + const double bx = pn.plane_x - beam_x, by = pn.plane_y - beam_y; if (bx * bx + by * by < MIN_BEAM_DISTANCE_PX2) continue; - const float d_A = geometry.PxlToRes(c.x, c.y); + const float d_A = geometry.PxlToRes(pn.plane_x, pn.plane_y); const int band = band_of(d_A); if (band < 0) continue; if (!isolated(i)) continue; - candidates[band].push_back({i, spots[i].Count(), d_A}); + candidates[band].push_back({i, spots[i].Count(), d_A, panel, pn}); } // The ring is gathered as the counts it is - the median of an int list is the same number, and @@ -390,14 +400,15 @@ void MeasureSpotFluxCurves(const std::function ¢re_ if (!(curve.c[R_NORM - 1] > 0.0f) || !(curve.c[R_MAX - 1] > 0.0f)) continue; const float norm = curve.c[R_NORM - 1]; for (float &v : curve.c) v /= norm; - curve.shape = shape_of(centre, stride, cx, cy, mx, my, ring, geometry); + curve.shape = shape_of(centre, stride, cx, cy, mx, my, ring, geometry, modules, cand[k].panel, cand[k].pn); out.push_back(curve); } } } SpotTileMisses MeasureSpotsOnTiles(ImagePreprocessorBuffer &image, int width, int height, - const DiffractionGeometry &geometry, const std::vector &spots, + const DiffractionGeometry &geometry, const DetectorModules &modules, + const std::vector &spots, std::vector &curves, std::vector &footprint, int half, int footprint_half) { SpotTileMisses misses; @@ -417,7 +428,7 @@ SpotTileMisses MeasureSpotsOnTiles(ImagePreprocessorBuffer &image, int width, in MeasureSpotFluxCurves([&](size_t s) { return tiled ? tiles + s * square(half) + static_cast(half) * (side + 1) : whole + static_cast(ty[s]) * width + tx[s]; - }, tiled ? side : width, width, height, geometry, spots, curves); + }, tiled ? side : width, width, height, geometry, modules, spots, curves); const float beam_x = geometry.GetBeamX_pxl(), beam_y = geometry.GetBeamY_pxl(); const auto tile_of = [&](const int32_t *t, size_t k, size_t s, int h) { @@ -425,11 +436,13 @@ SpotTileMisses MeasureSpotsOnTiles(ImagePreprocessorBuffer &image, int width, in }; std::vector measured(spots.size()); std::vector outcome(spots.size()); + std::vector panel(spots.size()); std::vector off_tile; for (size_t s = 0; s < spots.size(); s++) { const Coord c = spots[s].RawCoord(); + panel[s] = MakeBraggPanel(modules, modules.ModuleOfImagePixel(c.x, c.y), width, height, c.x, c.y); outcome[s] = MeasureFootprintSpot(tiled ? tile_of(tiles, s, s, half) : FootprintTile{whole, 0, 0, width, height}, - width, height, beam_x, beam_y, c.x, c.y, measured[s]); + panel[s], beam_x, beam_y, c.x, c.y, measured[s]); if (outcome[s] == FootprintOutcome::OffTile) off_tile.push_back(s); } @@ -453,7 +466,7 @@ SpotTileMisses MeasureSpotsOnTiles(ImagePreprocessorBuffer &image, int width, in const size_t s = off_tile[k]; const Coord c = spots[s].RawCoord(); outcome[s] = MeasureFootprintSpot(larger ? tile_of(again, k, s, h) : FootprintTile{whole, 0, 0, width, height}, - width, height, beam_x, beam_y, c.x, c.y, measured[s]); + panel[s], beam_x, beam_y, c.x, c.y, measured[s]); if (outcome[s] == FootprintOutcome::OffTile) still.push_back(s); } diff --git a/image_analysis/bragg_integration/SpotWidth.h b/image_analysis/bragg_integration/SpotWidth.h index 085e06527..1635968a3 100644 --- a/image_analysis/bragg_integration/SpotWidth.h +++ b/image_analysis/bragg_integration/SpotWidth.h @@ -30,6 +30,7 @@ #include #include +#include "../../common/DetectorModules.h" #include "../../common/DiffractionGeometry.h" #include "../../common/DiffractionSpot.h" #include "../image_preprocessing/ImagePreprocessorBuffer.h" @@ -186,9 +187,10 @@ struct BandwidthEstimate { // Every isolated, strong, fully readable spot of one image, as a normalised encircled-flux curve // with its shape, appended to `out`. Called on one worker's own vector, so the pool is independent of the order the -// workers finish in. +// workers finish in. The spots are in the image; `modules` places its panels on the plane the geometry +// describes (empty: the image is the plane), and a spot's aperture has to lie on its own panel. void MeasureSpotFluxCurves(const ImagePreprocessorBuffer &image, int width, int height, - const DiffractionGeometry &geometry, + const DiffractionGeometry &geometry, const DetectorModules &modules, const std::vector &spots, std::vector &out); @@ -196,7 +198,7 @@ void MeasureSpotFluxCurves(const ImagePreprocessorBuffer &image, int width, int // rows `stride` apart: R_BKG_OUT pixels of them every way that the image itself has. The frame is // width x height. void MeasureSpotFluxCurves(const std::function ¢re_px, int stride, int width, int height, - const DiffractionGeometry &geometry, + const DiffractionGeometry &geometry, const DetectorModules &modules, const std::vector &spots, std::vector &out); @@ -219,6 +221,7 @@ struct SpotTileMisses { // back to the host. The answers are the whole frame's to the bit: each measurement reads the same // pixels in the same order. The whole frame is read where the squares would be more pixels than it. SpotTileMisses MeasureSpotsOnTiles(ImagePreprocessorBuffer &image, int width, int height, - const DiffractionGeometry &geometry, const std::vector &spots, + const DiffractionGeometry &geometry, const DetectorModules &modules, + const std::vector &spots, std::vector &curves, std::vector &footprint, int half = SPOT_TILE_HALF, int footprint_half = FOOTPRINT_TILE_HALF); diff --git a/image_analysis/bragg_prediction/BraggPrediction.cpp b/image_analysis/bragg_prediction/BraggPrediction.cpp index 53fa376c2..a5aa7b7c5 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.cpp +++ b/image_analysis/bragg_prediction/BraggPrediction.cpp @@ -71,8 +71,11 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal const BraggPredictionSettings &settings) { const auto geom = experiment.GetDiffractionGeometry(); - const auto det_width_pxl = static_cast(experiment.GetXPixelsNum()); - const auto det_height_pxl = static_cast(experiment.GetYPixelsNum()); + // The plane the panels span (without panels, the image): a quick reject, a reflection in a gap + // between panels is dropped when it is integrated. + const PlaneBox bounds = experiment.GetModules().PlaneBounds(experiment.GetXPixelsNum(), experiment.GetYPixelsNum()); + const auto det_x0 = static_cast(bounds.x0), det_x1 = static_cast(bounds.x1); + const auto det_y0 = static_cast(bounds.y0), det_y1 = static_cast(bounds.y1); const float one_over_dmax = 1.0f / settings.high_res_A; const float one_over_dmax_sq = one_over_dmax * one_over_dmax; @@ -192,7 +195,7 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal float x = beam_x + F * S_rot_x / S_rot_z; float y = beam_y + F * S_rot_y / S_rot_z; - if ((x < 0) || (x >= det_width_pxl) || (y < 0) || (y >= det_height_pxl)) + if ((x < det_x0) || (x >= det_x1) || (y < det_y0) || (y >= det_y1)) continue; // Sensor quantum efficiency at this reflection's angle of incidence on the diff --git a/image_analysis/bragg_prediction/BraggPredictionGPU.cu b/image_analysis/bragg_prediction/BraggPredictionGPU.cu index 9d1689f95..09477cdbf 100644 --- a/image_analysis/bragg_prediction/BraggPredictionGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionGPU.cu @@ -122,7 +122,7 @@ namespace { float coeff = C.coeff_const / Srz; float x = C.beam_x + Srx * coeff; float y = C.beam_y + Sry * coeff; - if (x < 0.0f || x >= C.det_width_pxl || y < 0.0f || y >= C.det_height_pxl) return false; + if (x < C.det_x0 || x >= C.det_x1 || y < C.det_y0 || y >= C.det_y1) return false; // Sensor quantum efficiency and the air flight path at this reflection's angle of incidence // on the detector. Sr is the diffracted direction in the detector's own frame, so Srz over @@ -189,8 +189,11 @@ namespace { float bandwidth_sigma) { KernelConsts kc{}; auto geom = experiment.GetDiffractionGeometry(); - kc.det_width_pxl = static_cast(experiment.GetXPixelsNum()); - kc.det_height_pxl = static_cast(experiment.GetYPixelsNum()); + const PlaneBox bounds = experiment.GetModules().PlaneBounds(experiment.GetXPixelsNum(), experiment.GetYPixelsNum()); + kc.det_x0 = static_cast(bounds.x0); + kc.det_y0 = static_cast(bounds.y0); + kc.det_x1 = static_cast(bounds.x1); + kc.det_y1 = static_cast(bounds.y1); kc.beam_x = geom.GetBeamX_pxl(); kc.beam_y = geom.GetBeamY_pxl(); kc.coeff_const = geom.GetDetectorDistance_mm() / geom.GetPixelSize_mm(); diff --git a/image_analysis/bragg_prediction/BraggPredictionGPU.h b/image_analysis/bragg_prediction/BraggPredictionGPU.h index 7f38ba8f6..5d75cf4e9 100644 --- a/image_analysis/bragg_prediction/BraggPredictionGPU.h +++ b/image_analysis/bragg_prediction/BraggPredictionGPU.h @@ -10,8 +10,7 @@ #include "../indexing/CUDAMemHelpers.h" struct KernelConsts { - float det_width_pxl; - float det_height_pxl; + float det_x0, det_y0, det_x1, det_y1; // the plane the panels span (DetectorModules::PlaneBounds) float beam_x; float beam_y; float coeff_const; diff --git a/image_analysis/bragg_prediction/BraggPredictionRot.cpp b/image_analysis/bragg_prediction/BraggPredictionRot.cpp index 929b90a66..e4cab889b 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRot.cpp +++ b/image_analysis/bragg_prediction/BraggPredictionRot.cpp @@ -16,8 +16,11 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys const BraggPredictionSettings &settings) { const auto geom = experiment.GetDiffractionGeometry(); - const auto det_width_pxl = static_cast(experiment.GetXPixelsNum()); - const auto det_height_pxl = static_cast(experiment.GetYPixelsNum()); + // The plane the panels span (without panels, the image): a quick reject, a reflection in a gap + // between panels is dropped when it is integrated. + const PlaneBox bounds = experiment.GetModules().PlaneBounds(experiment.GetXPixelsNum(), experiment.GetYPixelsNum()); + const auto det_x0 = static_cast(bounds.x0), det_x1 = static_cast(bounds.x1); + const auto det_y0 = static_cast(bounds.y0), det_y1 = static_cast(bounds.y1); const float one_over_dmax = 1.0f / settings.high_res_A; const float one_over_dmax_sq = one_over_dmax * one_over_dmax; @@ -238,7 +241,7 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys float x = beam_x + F * S_rot_x / S_rot_z; float y = beam_y + F * S_rot_y / S_rot_z; - if ((x < 0) || (x >= det_width_pxl) || (y < 0) || (y >= det_height_pxl)) + if ((x < det_x0) || (x >= det_x1) || (y < det_y0) || (y >= det_y1)) continue; float dist_ewald_sphere = std::fabs(S.Length() - one_over_wavelength); diff --git a/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu b/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu index 0fb5e32a5..1e6f932f8 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu @@ -186,7 +186,7 @@ namespace { float x = C.beam_x + Srx * coeff; float y = C.beam_y + Sry * coeff; - if (x < 0.0f || x >= C.det_width_pxl || y < 0.0f || y >= C.det_height_pxl) + if (x < C.det_x0 || x >= C.det_x1 || y < C.det_y0 || y >= C.det_y1) continue; float dist_ewald = fabsf(sqrtf(Sx * Sx + Sy * Sy + Sz * Sz) - C.one_over_wavelength); @@ -350,8 +350,11 @@ namespace { KernelConstsRot kc{}; auto geom = experiment.GetDiffractionGeometry(); - kc.det_width_pxl = static_cast(experiment.GetXPixelsNum()); - kc.det_height_pxl = static_cast(experiment.GetYPixelsNum()); + const PlaneBox bounds = experiment.GetModules().PlaneBounds(experiment.GetXPixelsNum(), experiment.GetYPixelsNum()); + kc.det_x0 = static_cast(bounds.x0); + kc.det_y0 = static_cast(bounds.y0); + kc.det_x1 = static_cast(bounds.x1); + kc.det_y1 = static_cast(bounds.y1); kc.beam_x = geom.GetBeamX_pxl(); kc.beam_y = geom.GetBeamY_pxl(); kc.coeff_const = geom.GetDetectorDistance_mm() / geom.GetPixelSize_mm(); diff --git a/image_analysis/bragg_prediction/BraggPredictionRotGPU.h b/image_analysis/bragg_prediction/BraggPredictionRotGPU.h index 63ea227cc..17c17a95b 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRotGPU.h +++ b/image_analysis/bragg_prediction/BraggPredictionRotGPU.h @@ -10,8 +10,7 @@ #include "../indexing/CUDAMemHelpers.h" struct KernelConstsRot { - float det_width_pxl; - float det_height_pxl; + float det_x0, det_y0, det_x1, det_y1; // the plane the panels span (DetectorModules::PlaneBounds) float beam_x; float beam_y; float coeff_const; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 037f225a1..eaa72cb87 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -286,17 +286,24 @@ static void LogHeaderBeamCenter(const DiffractionExperiment &experiment, const P const float beam_y = experiment.GetBeamY_pxl(); const auto width = static_cast(experiment.GetXPixelsNumConv()); const auto height = static_cast(experiment.GetYPixelsNumConv()); + // The detector is the plane its panels span; without panels, the image. + const auto &modules = experiment.GetModules(); + const PlaneBox box = modules.PlaneBounds(width, height); - if (std::hypot(beam_x - static_cast(width) / 2.0f, - beam_y - static_cast(height) / 2.0f) < 1.0f) + if (std::hypot(beam_x - (static_cast(box.x0) + static_cast(box.x1)) / 2.0f, + beam_y - (static_cast(box.y0) + static_cast(box.y1)) / 2.0f) < 1.0f) logger.Info("Beam centre ({:.2f},{:.2f}) is the geometric centre of the detector - a value " "written rather than measured", beam_x, beam_y); else if (beam_x == std::round(beam_x) && beam_y == std::round(beam_y)) logger.Info("Beam centre ({:.0f},{:.0f}) is a whole number of pixels in both coordinates - a " "value written rather than measured", beam_x, beam_y); - const int x = static_cast(std::lround(beam_x)); - const int y = static_cast(std::lround(beam_y)); + // The pixel the beam lands on, if it lands on a panel at all: in the gap between panels it does not. + const auto on = modules.PlaneToPanel(beam_x, beam_y); + if (!on) + return; + const int x = static_cast(std::lround(on->x)); + const int y = static_cast(std::lround(on->y)); if (x >= 0 && y >= 0 && x < width && y < height && pixel_mask.GetMask()[static_cast(y) * width + x] != 0) logger.Warning("Beam centre ({:.2f},{:.2f}) lands on a masked pixel. A real beam does not, so " diff --git a/rugnux/RugnuxFirstPassAccept.cpp b/rugnux/RugnuxFirstPassAccept.cpp index 943e4e448..6c1bef3bf 100644 --- a/rugnux/RugnuxFirstPassAccept.cpp +++ b/rugnux/RugnuxFirstPassAccept.cpp @@ -406,9 +406,12 @@ void Rugnux::FirstPassRun::ReportTiltWalk(const FirstPass &best) { const double start_rot1 = tw.refused ? g.GetPoniRot1_rad() : rx.experiment_.GetPoniRot1_rad(); const double start_rot2 = tw.refused ? g.GetPoniRot2_rad() : rx.experiment_.GetPoniRot2_rad(); const double walk_rad = std::hypot(tw.tilt_rad[0] - start_rot1, tw.tilt_rad[1] - start_rot2); + // The far corner of the plane the panels span (without panels, the image). + const PlaneBox box = rx.experiment_.GetModules().PlaneBounds(rx.experiment_.GetXPixelsNumConv(), + rx.experiment_.GetYPixelsNumConv()); const double corner_px = std::hypot( - std::max(g.GetBeamX_pxl(), rx.experiment_.GetXPixelsNumConv() - g.GetBeamX_pxl()), - std::max(g.GetBeamY_pxl(), rx.experiment_.GetYPixelsNumConv() - g.GetBeamY_pxl())); + std::max(g.GetBeamX_pxl() - static_cast(box.x0), static_cast(box.x1) - g.GetBeamX_pxl()), + std::max(g.GetBeamY_pxl() - static_cast(box.y0), static_cast(box.y1) - g.GetBeamY_pxl())); const double lever_px = g.GetDetectorDistance_mm() / g.GetPixelSize_mm(); const std::string what = fmt::format( "Detector tilt: the free fit walked to {:.4f},{:.4f} deg, {:.2f} deg from the " diff --git a/rugnux/RugnuxPasses.cpp b/rugnux/RugnuxPasses.cpp index 7256ab640..2dc4d1d36 100644 --- a/rugnux/RugnuxPasses.cpp +++ b/rugnux/RugnuxPasses.cpp @@ -1000,10 +1000,13 @@ std::optional Rugnux::RefinedPasses(RugnuxObserver *observer, Pro std::shared_ptr demoted_recheck; // the demoted arm again on all the frames, for a near tie if (!cancelled_ && !config_.fast_live && (length_arm || angle_arm) && pass1.rotation_lattice_type) { split = length_arm ? *pass1.lattice_length_split : *pass1.lattice_angle_departure_deg * PI / 180.0; - const double dx = std::max(experiment_.GetBeamX_pxl(), - static_cast(experiment_.GetXPixelsNum()) - experiment_.GetBeamX_pxl()); - const double dy = std::max(experiment_.GetBeamY_pxl(), - static_cast(experiment_.GetYPixelsNum()) - experiment_.GetBeamY_pxl()); + // The far corner of the plane the panels span (without panels, the image). + const PlaneBox box = experiment_.GetModules().PlaneBounds(experiment_.GetXPixelsNum(), + experiment_.GetYPixelsNum()); + const double dx = std::max(experiment_.GetBeamX_pxl() - static_cast(box.x0), + static_cast(box.x1) - experiment_.GetBeamX_pxl()); + const double dy = std::max(experiment_.GetBeamY_pxl() - static_cast(box.y0), + static_cast(box.y1) - experiment_.GetBeamY_pxl()); const double r_corner_px = std::hypot(dx, dy); const double miss_px = 0.5 * split * r_corner_px; const double disc_px = experiment_.GetBraggIntegrationSettings().GetR1(); diff --git a/rugnux/RugnuxPreScan.cpp b/rugnux/RugnuxPreScan.cpp index 0ae59a893..42f489d29 100644 --- a/rugnux/RugnuxPreScan.cpp +++ b/rugnux/RugnuxPreScan.cpp @@ -383,9 +383,12 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // sharpens them against the crystal's own rows. { const auto &g = prescan_x.GetDiffractionGeometry(); + const auto &modules = prescan_x.GetModules(); for (const auto &spot : spots) { + // The spot is in the image; its resolution is where its panel puts it on the plane. const Coord c = spot.RawCoord(); - const float d = g.PxlToRes(c.x, c.y); + const auto [px, py] = modules.ImageToPlane(modules.ModuleOfImagePixel(c.x, c.y), c.x, c.y); + const float d = g.PxlToRes(static_cast(px), static_cast(py)); if (d > 0.0f && std::isfinite(d)) spot_q.push_back(2 * PI / d); } @@ -396,7 +399,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru if (for_width) MeasureSpotsOnTiles(*w.preprocessed, static_cast(prescan_x.GetXPixelsNumConv()), static_cast(prescan_x.GetYPixelsNumConv()), prescan_x.GetDiffractionGeometry(), - spots, curves, footprint); + prescan_x.GetModules(), spots, curves, footprint); }; std::vector prescan_spot_q; @@ -765,11 +768,15 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // longer holds them (SpotFootprint.h). It changes the integration only where a spot // outgrows the r1 disk. const auto &g = experiment_.GetDiffractionGeometry(); + // Out to the far corner of the plane the panels span (without panels, the image). const float bx = g.GetBeamX_pxl(), by = g.GetBeamY_pxl(); - const float W = static_cast(experiment_.GetXPixelsNum()); - const float H = static_cast(experiment_.GetYPixelsNum()); + const PlaneBox box = experiment_.GetModules().PlaneBounds(experiment_.GetXPixelsNum(), + experiment_.GetYPixelsNum()); + const float x0 = static_cast(box.x0), x1 = static_cast(box.x1); + const float y0 = static_cast(box.y0), y1 = static_cast(box.y1); const SpotFootprint fp = FootprintFromSpots(footprint_spots, - std::hypot(std::max(bx, W - bx), std::max(by, H - by))); + std::hypot(std::max(bx - x0, x1 - bx), + std::max(by - y0, y1 - by))); if (!fp.empty()) prescan_footprint_ = fp; // It acts only where a spot outgrows the r1 disk (BraggStencil.h), so a pattern of compact diff --git a/tests/BandwidthEstimateTest.cpp b/tests/BandwidthEstimateTest.cpp index 56d68ebf9..c8979b63c 100644 --- a/tests/BandwidthEstimateTest.cpp +++ b/tests/BandwidthEstimateTest.cpp @@ -71,7 +71,7 @@ namespace { ImagePreprocessorBuffer image(static_cast(W) * H); std::vector spots; PaintSpots(image.getBuffer(), spots, bw_fwhm, offset, rng); - MeasureSpotFluxCurves(image, W, H, geometry, spots, curves); + MeasureSpotFluxCurves(image, W, H, geometry, DetectorModules(), spots, curves); } return spot_width::EstimateBandwidth(curves, sensor_absorption::AttenuationLength_um("Si", LAMBDA_A), THICKNESS_UM, PIXEL_MM * 1000.0); diff --git a/tests/PanelLayoutTest.cpp b/tests/PanelLayoutTest.cpp index 4c8f5a484..084e4a384 100644 --- a/tests/PanelLayoutTest.cpp +++ b/tests/PanelLayoutTest.cpp @@ -24,6 +24,7 @@ #include #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/spot_finding/ImageSpotFinderCPU.h" #include "../image_analysis/spot_finding/SpotUtils.h" @@ -148,7 +153,8 @@ void CheckSameSpots(const std::vector &a, const std::vector Integrate(const DiffractionExperiment &x, const PixelMask &mask, - const ImagePreprocessorBuffer &image, const std::vector &spots) { + const ImagePreprocessorBuffer &image, const std::vector &spots, + bool gpu = false) { const auto &modules = x.GetModules(); const auto geom = x.GetDiffractionGeometry(); std::vector predicted; @@ -166,9 +172,25 @@ std::map Integrate(const DiffractionExperiment &x, const PixelM predicted.push_back(r); panel_of[r.h] = p->panel; } - BraggIntegrationEngineCPU engine(x, mask); + std::vector integrated; +#ifdef JFJOCH_USE_CUDA + if (gpu) { + auto stream = std::make_shared(); + 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 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(px); @@ -293,20 +315,60 @@ TEST_CASE("PanelLayout_ResolutionMask", "[.][PanelLayoutTodo]") { 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 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. diff --git a/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index c85c7b86e..0aa0f08bf 100644 --- a/tests/RugnuxTest.cpp +++ b/tests/RugnuxTest.cpp @@ -466,7 +466,7 @@ TEST_CASE("SpotWidth_Gaussian", "[process]") { PaintGaussianSpots(image, W, sigma, 20000.0, spots); std::vector curves; - MeasureSpotFluxCurves(image, W, H, geometry, spots, curves); + MeasureSpotFluxCurves(image, W, H, geometry, DetectorModules(), spots, curves); REQUIRE(curves.size() >= 45); const auto r80 = spot_width::R80AtReference(curves); diff --git a/tests/SpotTilesTest.cpp b/tests/SpotTilesTest.cpp index d5c2f5d12..0cb1d0eb2 100644 --- a/tests/SpotTilesTest.cpp +++ b/tests/SpotTilesTest.cpp @@ -84,7 +84,7 @@ SpotTileMisses CheckFrames(std::vector &frames, Image &&image, int half, const auto &g = f.x.GetDiffractionGeometry(); std::vector c_frame, c_tiles; std::vector f_frame, f_tiles; - MeasureSpotFluxCurves(f.image, width, height, g, f.spots, c_frame); + MeasureSpotFluxCurves(f.image, width, height, g, DetectorModules(), f.spots, c_frame); std::vector sx, sy; for (const auto &s : f.spots) { sx.push_back(s.RawCoord().x); @@ -92,7 +92,7 @@ SpotTileMisses CheckFrames(std::vector &frames, Image &&image, int half, } MeasureFootprintSpots(f.image.data(), width, height, g.GetBeamX_pxl(), g.GetBeamY_pxl(), sx, sy, f_frame); - const auto misses = MeasureSpotsOnTiles(image(f), width, height, g, f.spots, c_tiles, f_tiles, + const auto misses = MeasureSpotsOnTiles(image(f), width, height, g, DetectorModules(), f.spots, c_tiles, f_tiles, half, footprint_half); CHECK(Same(c_frame, c_tiles)); CHECK(Same(f_frame, f_tiles));