diff --git a/common/DiffractionExperiment.cpp b/common/DiffractionExperiment.cpp index 3c0e059d7..1c2ec942b 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 e4e1efbc6..c54bfe0b6 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/roi/ROIIntegrationCPU.h" #include "../image_analysis/spot_finding/ImageSpotFinderCPU.h" @@ -149,7 +154,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; @@ -167,9 +173,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); @@ -315,20 +337,60 @@ TEST_CASE("PanelLayout_ResolutionMask", "[PanelLayout]") { CHECK(differ == 0); } -// TODO stage 2: Bragg integration on the stacked layout measures what it measures on the assembled image - -// stencil direction, radial background and polarisation taken from the plane position, not the stored one. -TEST_CASE("PanelLayout_BraggIntegration", "[.][PanelLayoutTodo]") { - const Scene scene(INTERIOR_SPOTS); - const auto a = Integrate(scene.assembled, scene.assembled_mask, scene.assembled_image, INTERIOR_SPOTS); - const auto s = Integrate(scene.stacked, scene.stacked_mask, scene.stacked_image, INTERIOR_SPOTS); - REQUIRE(s.size() == a.size()); - for (const auto &[h, r]: s) { - CHECK(r.I == Catch::Approx(a.at(h).I).epsilon(1e-5)); - CHECK(r.sigma == Catch::Approx(a.at(h).sigma).epsilon(1e-5)); - CHECK(r.bkg == Catch::Approx(a.at(h).bkg).epsilon(1e-5)); +// Bragg integration on the stacked layout measures what it measures on the assembled image - stencil +// direction, radial background and polarisation taken from the plane position, not the stored one. +// Also with what reads the radial direction hardest: polarisation, a ring elongated by a bandwidth, the +// radial background correction, the owner map, and the box sum. The interior spots moved off the exact +// half pixel: there the profile grid centres on the rounded position, and rounding a half picks the other +// neighbour once a panel is turned. +void CheckBraggIntegration(bool gpu) { + const std::vector 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));