From 545cf3b13eaaf376306b492c6f298baf6f8bac0c Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 12 Oct 2026 16:11:14 +0200 Subject: [PATCH] rugnux pre-scan: spot width on per-spot squares, engines built during the read, tiers looked ahead, measurements memoised The spot-width stage of the pre-scan copied every preprocessed frame back to the host (72 MB a frame on a 16M detector, page-locked host mirrors per worker) for measurements that read a few pixels about each spot. Now: - MeasureSpotsOnTiles: the flux curves and the footprint read (2*24+1)^2 squares about each spot, gathered on the card (ImagePreprocessorBufferGPU::GatherTiles) - no copy_image_to_host, no host mirror. The flux curves are bounded at R_BKG_OUT = 20 px. The footprint keeps its arithmetic; a per-iteration check that its window is inside the square, and a spot whose window walks off it is measured again from scratch on an 80 px square, then on the whole frame. Where the squares would be more pixels than the frame, the frame is read. Same pixels in the same order, so the same bits. Harness over 34 rotation sets (~1500 frames): 0 differences; misses at 24 on a few spots of most frames, at 54 on up to 80% of frames of a few sets, at 80 on none. - The azimuthal mapping and the per-worker spot engines are built while the frames are read. - The tiers no longer wait for each other: workers run on through the visit order, each tier's answer is drawn from exactly its own frames, and frames measured past the tier the answer stops at are thrown away - a frame is measured the same whichever tier it is read for. - PreScanMemo keeps the spot measurement, the defective pixels and the background gain under the pre-scan key plus the spot-finding and azimuthal settings, so a copy of the run that pre-scans again on the same inputs reads no frames. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- .../bragg_integration/SpotFootprint.cpp | 131 ++++---- .../bragg_integration/SpotFootprint.h | 14 + .../bragg_integration/SpotWidth.cpp | 99 +++++- image_analysis/bragg_integration/SpotWidth.h | 33 ++ .../ImagePreprocessorBuffer.cpp | 18 ++ .../ImagePreprocessorBuffer.h | 10 + .../ImagePreprocessorBufferGPU.cu | 57 ++++ .../ImagePreprocessorBufferGPU.h | 11 + rugnux/Rugnux.h | 35 ++- rugnux/RugnuxPreScan.cpp | 296 ++++++++++++------ tests/CMakeLists.txt | 1 + tests/SpotTilesTest.cpp | 145 +++++++++ 12 files changed, 691 insertions(+), 159 deletions(-) create mode 100644 tests/SpotTilesTest.cpp diff --git a/image_analysis/bragg_integration/SpotFootprint.cpp b/image_analysis/bragg_integration/SpotFootprint.cpp index 584d973da..88a327f24 100644 --- a/image_analysis/bragg_integration/SpotFootprint.cpp +++ b/image_analysis/bragg_integration/SpotFootprint.cpp @@ -30,69 +30,82 @@ float median_of(std::vector &v) { } // namespace +FootprintOutcome MeasureFootprintSpot(const FootprintTile &tile, int width, int height, 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; + 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 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; + + // 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. + float cx = x, cy = y; + float s2r = 1.0f, s2t = 1.0f; + for (int it = 0; it < ITERATIONS; ++it) { + const float wr = std::clamp(WINDOW_NSIGMA * std::sqrt(s2r), WINDOW_MIN_PX, WINDOW_MAX_PX); + const float wt = std::clamp(WINDOW_NSIGMA * std::sqrt(s2t), WINDOW_MIN_PX, WINDOW_MAX_PX); + 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; + // 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; + const auto pixel = [&](int px, int py) { + return tile.pixels[static_cast(py - tile.y0) * tile.w + (px - tile.x0)]; + }; + + ring.clear(); + for (int py = y0; py <= y1; ++py) + for (int px = x0; px <= x1; ++px) { + const float dx = px - cx, dy = py - cy; + const float rad = dx * ux + dy * uy, tn = -dx * uy + dy * ux; + const float e = rad * rad / (wr * wr) + tn * tn / (wt * wt); + const int32_t v = pixel(px, py); + if (e >= BKG_INNER * BKG_INNER && e < BKG_OUTER * BKG_OUTER && valid(v)) + ring.push_back(static_cast(v)); + } + if (ring.size() < 10) return FootprintOutcome::LeftOut; + const double bkg = median_of(ring); + + double w = 0.0, mr = 0.0, mt = 0.0, m2r = 0.0, m2t = 0.0; + for (int py = y0; py <= y1; ++py) + for (int px = x0; px <= x1; ++px) { + const float dx = px - cx, dy = py - cy; + const float rad = dx * ux + dy * uy, tn = -dx * uy + dy * ux; + if (rad * rad / (wr * wr) + tn * tn / (wt * wt) >= 1.0f) continue; + const int32_t v = pixel(px, py); + if (!valid(v)) return FootprintOutcome::LeftOut; + const double net = static_cast(v) - bkg; + w += net; + mr += net * rad; + mt += net * tn; + m2r += net * rad * rad; + m2t += net * tn * tn; + } + if (!(w > 0.0)) return FootprintOutcome::LeftOut; + const double cr = mr / w, ct = mt / w; + s2r = static_cast(std::max(0.25, m2r / w - cr * cr)); + s2t = static_cast(std::max(0.25, m2t / w - ct * ct)); + cx += static_cast(cr * ux - ct * uy); + cy += static_cast(cr * uy + ct * ux); + } + spot = {r, std::sqrt(s2r), std::sqrt(s2t)}; + return FootprintOutcome::Measured; +} + void MeasureFootprintSpots(const int32_t *img, int width, int height, float beam_x, float beam_y, const std::vector &x, const std::vector &y, std::vector &out) { - const int half = static_cast(std::ceil(BKG_OUTER * WINDOW_MAX_PX)) + 1; - std::vector ring; + const FootprintTile frame{img, 0, 0, width, height}; for (size_t s = 0; s < x.size(); ++s) { - const float rx = x[s] - beam_x, ry = y[s] - beam_y; - const float r = std::sqrt(rx * rx + ry * ry); - if (!(r > 1.0f)) continue; - const float ux = rx / r, uy = ry / r; - const int ix = static_cast(std::lround(x[s])), iy = static_cast(std::lround(y[s])); - if (ix - half < 0 || iy - half < 0 || ix + half >= width || iy + half >= height) continue; - - // 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. - float cx = x[s], cy = y[s]; - float s2r = 1.0f, s2t = 1.0f; - bool ok = true; - for (int it = 0; it < ITERATIONS && ok; ++it) { - const float wr = std::clamp(WINDOW_NSIGMA * std::sqrt(s2r), WINDOW_MIN_PX, WINDOW_MAX_PX); - const float wt = std::clamp(WINDOW_NSIGMA * std::sqrt(s2t), WINDOW_MIN_PX, WINDOW_MAX_PX); - 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) { ok = false; break; } - - ring.clear(); - for (int py = y0; py <= y1; ++py) - for (int px = x0; px <= x1; ++px) { - const float dx = px - cx, dy = py - cy; - const float rad = dx * ux + dy * uy, tn = -dx * uy + dy * ux; - const float e = rad * rad / (wr * wr) + tn * tn / (wt * wt); - const int32_t v = img[static_cast(py) * width + px]; - if (e >= BKG_INNER * BKG_INNER && e < BKG_OUTER * BKG_OUTER && valid(v)) - ring.push_back(static_cast(v)); - } - if (ring.size() < 10) { ok = false; break; } - const double bkg = median_of(ring); - - double w = 0.0, mr = 0.0, mt = 0.0, m2r = 0.0, m2t = 0.0; - for (int py = y0; py <= y1 && ok; ++py) - for (int px = x0; px <= x1; ++px) { - const float dx = px - cx, dy = py - cy; - const float rad = dx * ux + dy * uy, tn = -dx * uy + dy * ux; - if (rad * rad / (wr * wr) + tn * tn / (wt * wt) >= 1.0f) continue; - const int32_t v = img[static_cast(py) * width + px]; - if (!valid(v)) { ok = false; break; } - const double net = static_cast(v) - bkg; - w += net; - mr += net * rad; - mt += net * tn; - m2r += net * rad * rad; - m2t += net * tn * tn; - } - if (!ok || !(w > 0.0)) { ok = false; break; } - const double cr = mr / w, ct = mt / w; - s2r = static_cast(std::max(0.25, m2r / w - cr * cr)); - s2t = static_cast(std::max(0.25, m2t / w - ct * ct)); - cx += static_cast(cr * ux - ct * uy); - cy += static_cast(cr * uy + ct * ux); - } - if (ok) - out.push_back({r, std::sqrt(s2r), std::sqrt(s2t)}); + FootprintSpot spot; + if (MeasureFootprintSpot(frame, width, height, 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 9b54c8186..d549abbba 100644 --- a/image_analysis/bragg_integration/SpotFootprint.h +++ b/image_analysis/bragg_integration/SpotFootprint.h @@ -57,6 +57,20 @@ void MeasureFootprintSpots(const int32_t *img, int width, int height, float beam const std::vector &x, const std::vector &y, std::vector &out); +// The pixels [x0, x0 + w) x [y0, y0 + h) of a frame, rows w apart: the whole frame, or a square cut +// out of it about one spot. +struct FootprintTile { + const int32_t *pixels = nullptr; + int x0 = 0, y0 = 0, w = 0, h = 0; +}; +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, + 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 // holding too few spots takes the nearest bin that has enough; with no such bin the table is empty. SpotFootprint FootprintFromSpots(const std::vector &spots, float r_max); diff --git a/image_analysis/bragg_integration/SpotWidth.cpp b/image_analysis/bragg_integration/SpotWidth.cpp index b9600e010..327f1896c 100644 --- a/image_analysis/bragg_integration/SpotWidth.cpp +++ b/image_analysis/bragg_integration/SpotWidth.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -107,7 +108,7 @@ constexpr int R_MOMENT = 9; // clipped: a median is biased low on sparse Poisson counts, and a bias here adds the same flux at // 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 width, int cx, int cy, double mx, double my, +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) { SpotShape shape; double sum = 0.0; @@ -127,7 +128,7 @@ SpotShape shape_of(const int32_t *centre_px, int width, int cx, int cy, double m for (int iter = 0; iter < 4; iter++) { double f = 0.0, sx = 0.0, sy = 0.0, qxx = 0.0, qyy = 0.0, qxy = 0.0; for (int dy = -R_MOMENT - 2; dy <= R_MOMENT + 2; dy++) { - const int32_t *row = centre_px + static_cast(dy) * width; + const int32_t *row = centre_px + static_cast(dy) * stride; for (int dx = -R_MOMENT - 2; dx <= R_MOMENT + 2; dx++) { const double ddx = dx - mx, ddy = dy - my; if (ddx * ddx + ddy * ddy >= R2) continue; @@ -221,10 +222,21 @@ void MeasureSpotFluxCurves(const ImagePreprocessorBuffer &image, int width, int const DiffractionGeometry &geometry, const std::vector &spots, std::vector &out) { - if (spots.empty()) return; - const int32_t *pixels = image.data(); if (pixels == nullptr) return; + const auto centre_px = [&](size_t i) { + 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); +} + +void MeasureSpotFluxCurves(const std::function ¢re_px, int stride, int width, int height, + const DiffractionGeometry &geometry, + const std::vector &spots, + std::vector &out) { + if (spots.empty()) return; + const float beam_x = geometry.GetBeamX_pxl(), beam_y = geometry.GetBeamY_pxl(); // Where every spot of this image sits, so isolation can be tested against all of them and not @@ -303,7 +315,7 @@ void MeasureSpotFluxCurves(const ImagePreprocessorBuffer &image, int width, int for (size_t k = 0; k < take; k++) { const Coord &c = centre[cand[k].index]; const int cx = static_cast(std::lround(c.x)), cy = static_cast(std::lround(c.y)); - const int32_t *centre_px = pixels + static_cast(cy) * width + cx; + const int32_t *centre = centre_px(cand[k].index); // The background under the spot, and a check that the whole aperture is readable: a hole // in it removes flux from one radius and not another, which is exactly the shape this @@ -312,7 +324,7 @@ void MeasureSpotFluxCurves(const ImagePreprocessorBuffer &image, int width, int bool readable = true; for (int dy = -R_BKG_OUT; dy <= R_BKG_OUT && readable; dy++) { const int half = HALF_BKG[std::abs(dy)], inner = INNER_BKG[std::abs(dy)]; - const int32_t *row = centre_px + static_cast(dy) * width; + const int32_t *row = centre + static_cast(dy) * stride; for (int dx = -half; dx <= half; dx++) { const int32_t px = row[dx]; if (!valid(px)) { readable = false; break; } @@ -329,7 +341,7 @@ void MeasureSpotFluxCurves(const ImagePreprocessorBuffer &image, int width, int int n_core = 0; for (int dy = -R_CENTROID; dy <= R_CENTROID; dy++) { const int half = HALF_CORE[std::abs(dy)]; - const int32_t *row = centre_px + static_cast(dy) * width; + const int32_t *row = centre + static_cast(dy) * stride; for (int dx = -half; dx <= half; dx++) { const double v = row[dx] - bkg; core += v; @@ -359,7 +371,7 @@ void MeasureSpotFluxCurves(const ImagePreprocessorBuffer &image, int width, int const double span = std::sqrt(R2_MAX - ddy2); const int lo = std::max(-R_MAX, static_cast(std::floor(mx - span))); const int hi = std::min(R_MAX, static_cast(std::ceil(mx + span))); - const int32_t *row = centre_px + static_cast(dy) * width; + const int32_t *row = centre + static_cast(dy) * stride; for (int dx = lo; dx <= hi; dx++) { const double ddx = dx - mx, rc2 = ddx * ddx + ddy2; if (rc2 > R2_MAX) continue; @@ -378,12 +390,81 @@ void MeasureSpotFluxCurves(const ImagePreprocessorBuffer &image, int width, int 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_px, width, cx, cy, mx, my, ring, geometry); + curve.shape = shape_of(centre, stride, cx, cy, mx, my, ring, geometry); out.push_back(curve); } } } +SpotTileMisses MeasureSpotsOnTiles(ImagePreprocessorBuffer &image, int width, int height, + const DiffractionGeometry &geometry, const std::vector &spots, + std::vector &curves, std::vector &footprint, + int half, int footprint_half) { + SpotTileMisses misses; + // The pixel each spot rounds to, which both measurements centre on. + std::vector tx(spots.size()), ty(spots.size()); + for (size_t s = 0; s < spots.size(); s++) { + const Coord c = spots[s].RawCoord(); + tx[s] = static_cast(std::lround(c.x)); + ty[s] = static_cast(std::lround(c.y)); + } + const size_t frame_pixels = static_cast(width) * height; + const auto square = [](int h) { return static_cast(2 * h + 1) * (2 * h + 1); }; + const bool tiled = spots.size() * square(half) <= frame_pixels; + const int32_t *whole = tiled ? nullptr : image.HostImage(); + const int32_t *tiles = tiled ? image.GatherTiles(tx, ty, width, half) : nullptr; + const int side = 2 * half + 1; + 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); + + 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) { + return FootprintTile{t + k * square(h), tx[s] - h, ty[s] - h, 2 * h + 1, 2 * h + 1}; + }; + std::vector measured(spots.size()); + std::vector outcome(spots.size()); + std::vector off_tile; + for (size_t s = 0; s < spots.size(); s++) { + const Coord c = spots[s].RawCoord(); + 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]); + if (outcome[s] == FootprintOutcome::OffTile) + off_tile.push_back(s); + } + // Again from the start: on the larger squares, then on the frame. 0 stands for the frame. + for (const int h : {footprint_half, 0}) { + if (off_tile.empty()) break; + const bool larger = h > 0 && off_tile.size() * square(h) <= frame_pixels; + (larger ? misses.larger : misses.whole_frame) += off_tile.size(); + const int32_t *again = nullptr; + if (larger) { + std::vector mx, my; + for (const size_t s : off_tile) { + mx.push_back(tx[s]); + my.push_back(ty[s]); + } + again = image.GatherTiles(mx, my, width, h); + } else if (!whole) + whole = image.HostImage(); + std::vector still; + for (size_t k = 0; k < off_tile.size(); k++) { + 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]); + if (outcome[s] == FootprintOutcome::OffTile) + still.push_back(s); + } + off_tile = std::move(still); + } + for (size_t s = 0; s < spots.size(); s++) + if (outcome[s] == FootprintOutcome::Measured) + footprint.push_back(measured[s]); + return misses; +} + float spot_width::R80Fit::At(double d_A) const { return static_cast(std::clamp(c0 + c1 / d_A, lo, hi)); } diff --git a/image_analysis/bragg_integration/SpotWidth.h b/image_analysis/bragg_integration/SpotWidth.h index 6186a747f..085e06527 100644 --- a/image_analysis/bragg_integration/SpotWidth.h +++ b/image_analysis/bragg_integration/SpotWidth.h @@ -26,12 +26,14 @@ // ============================================================================= #include +#include #include #include #include "../../common/DiffractionGeometry.h" #include "../../common/DiffractionSpot.h" #include "../image_preprocessing/ImagePreprocessorBuffer.h" +#include "SpotFootprint.h" namespace spot_width { // The encircled-flux curve runs to R_MAX px and is normalised at R_NORM, clear of the noisy tail. @@ -189,3 +191,34 @@ void MeasureSpotFluxCurves(const ImagePreprocessorBuffer &image, int width, int const DiffractionGeometry &geometry, const std::vector &spots, std::vector &out); + +// The same, with spot i's pixels read about centre_px(i) - the pixel its position rounds to - with +// 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 std::vector &spots, + std::vector &out); + +// The squares MeasureSpotsOnTiles reads about each spot. The flux curves stay within +// spot_width::R_BKG_OUT of their spot; the footprint's window re-centres without a bound, so a spot +// whose window walks off its square is measured again from the start on a square of +// FOOTPRINT_TILE_HALF, and on the whole frame where even that is left. Over the rotation battery's +// pre-scans a window leaves the first square on a few spots of most frames, and the second on none. +constexpr int SPOT_TILE_HALF = 24; +constexpr int FOOTPRINT_TILE_HALF = 80; + +// Spots whose footprint had to be measured again on a larger square, and on the whole frame. +struct SpotTileMisses { + size_t larger = 0; + size_t whole_frame = 0; +}; + +// MeasureSpotFluxCurves and MeasureFootprintSpots of one frame, reading the (2*half+1)^2 pixels about +// each spot (ImagePreprocessorBuffer::GatherTiles) rather than the frame - on the GPU only those cross +// 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, + std::vector &curves, std::vector &footprint, + int half = SPOT_TILE_HALF, int footprint_half = FOOTPRINT_TILE_HALF); diff --git a/image_analysis/image_preprocessing/ImagePreprocessorBuffer.cpp b/image_analysis/image_preprocessing/ImagePreprocessorBuffer.cpp index 46cae1d23..01dd73c43 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorBuffer.cpp +++ b/image_analysis/image_preprocessing/ImagePreprocessorBuffer.cpp @@ -22,3 +22,21 @@ int32_t *ImagePreprocessorBuffer::getGPUBuffer() { const int32_t *ImagePreprocessorBuffer::getGPUBuffer() const { return nullptr; } + +const int32_t *ImagePreprocessorBuffer::GatherTiles(const std::vector &x, const std::vector &y, + int width, int half) { + const int side = 2 * half + 1; + const int height = static_cast(npixels / width); + tiles.resize(x.size() * side * side); + int32_t *out = tiles.data(); + for (size_t i = 0; i < x.size(); i++) + for (int py = y[i] - half; py <= y[i] + half; py++) + for (int px = x[i] - half; px <= x[i] + half; px++) + *out++ = (px < 0 || py < 0 || px >= width || py >= height) + ? INT32_MIN : buffer[static_cast(py) * width + px]; + return tiles.data(); +} + +const int32_t *ImagePreprocessorBuffer::HostImage() { + return buffer.data(); +} diff --git a/image_analysis/image_preprocessing/ImagePreprocessorBuffer.h b/image_analysis/image_preprocessing/ImagePreprocessorBuffer.h index c57305152..40fb715be 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorBuffer.h +++ b/image_analysis/image_preprocessing/ImagePreprocessorBuffer.h @@ -6,6 +6,7 @@ #include #include #include +#include class ImagePreprocessorBuffer { protected: @@ -43,4 +44,13 @@ public: // Values of the pixels with the given flat indices. Not just operator[] in a loop: on the GPU the // preprocessed image can live in device memory only, and then this gathers them there. virtual void Gather(const std::vector &npixel, std::vector &values) const; + + // The (2*half+1)^2 pixels about each (x[i], y[i]) of the image, `width` pixels wide: square after + // square, each row by row, with a pixel off the image read as INT32_MIN. Valid until the next call. + virtual const int32_t *GatherTiles(const std::vector &x, const std::vector &y, + int width, int half); + // The whole image on the host - on the GPU copied there first, where it was not already. + virtual const int32_t *HostImage(); +private: + std::vector tiles; }; diff --git a/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu b/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu index 8f8fe5381..6c5981ada 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu +++ b/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu @@ -3,6 +3,9 @@ #include "ImagePreprocessorBufferGPU.h" +#include +#include + static inline void cuda_err(cudaError_t val) { if (val != cudaSuccess) throw JFJochException(JFJochExceptionCategory::GPUCUDAError, cudaGetErrorString(val)); @@ -16,6 +19,21 @@ __global__ void gather_kernel(const int32_t *__restrict__ image, values[i] = image[npixel[i]]; } +__global__ void gather_tiles_kernel(const int32_t *__restrict__ image, int width, int height, + const int32_t *__restrict__ centre, int half, + int32_t *__restrict__ tiles, size_t count) { + const int side = 2 * half + 1; + const size_t area = static_cast(side) * side; + for (size_t i = blockIdx.x * blockDim.x + threadIdx.x; i < count; i += static_cast(blockDim.x) * gridDim.x) { + const size_t tile = i / area; + const int p = static_cast(i - tile * area); + const int px = centre[2 * tile] - half + p % side; + const int py = centre[2 * tile + 1] - half + p / side; + tiles[i] = (px < 0 || py < 0 || px >= width || py >= height) + ? INT32_MIN : image[static_cast(py) * width + px]; + } +} + ImagePreprocessorBufferGPU::ImagePreprocessorBufferGPU(size_t npixel, bool host_mirror) : ImagePreprocessorBuffer(npixel, host_mirror), gpu_image(npixel), @@ -49,3 +67,42 @@ void ImagePreprocessorBufferGPU::Gather(const std::vector &npixel, std cudaMemcpyDeviceToHost, gather_stream); cudaStreamSynchronize(gather_stream); } + +const int32_t *ImagePreprocessorBufferGPU::GatherTiles(const std::vector &x, const std::vector &y, + int width, int half) { + const size_t side = 2 * half + 1; + const size_t count = x.size() * side * side; + if (count == 0) + return nullptr; + if (count > tile_capacity) { + tile_capacity = std::max(count, 2 * tile_capacity); + gpu_tiles = CudaDevicePtr(tile_capacity); + host_tiles.resize(tile_capacity); + } + if (2 * x.size() > tile_centre.size()) { + tile_centre.resize(std::max(2 * x.size(), 2 * tile_centre.size())); + gpu_tile_centre = CudaDevicePtr(tile_centre.size()); + } + for (size_t i = 0; i < x.size(); i++) { + tile_centre[2 * i] = x[i]; + tile_centre[2 * i + 1] = y[i]; + } + cuda_err(cudaMemcpyAsync(gpu_tile_centre.get(), tile_centre.data(), 2 * x.size() * sizeof(int32_t), + cudaMemcpyHostToDevice, gather_stream)); + const int blocks = static_cast(std::min((count + 255) / 256, 4096)); + gather_tiles_kernel<<>>(gpu_image.get(), width, static_cast(npixels / width), + gpu_tile_centre.get(), half, gpu_tiles.get(), count); + cuda_err(cudaGetLastError()); + cuda_err(cudaMemcpyAsync(host_tiles.data(), gpu_tiles.get(), count * sizeof(int32_t), + cudaMemcpyDeviceToHost, gather_stream)); + cuda_err(cudaStreamSynchronize(gather_stream)); + return host_tiles.data(); +} + +const int32_t *ImagePreprocessorBufferGPU::HostImage() { + buffer.resize(npixels); + cuda_err(cudaMemcpyAsync(buffer.data(), gpu_image.get(), npixels * sizeof(int32_t), + cudaMemcpyDeviceToHost, gather_stream)); + cuda_err(cudaStreamSynchronize(gather_stream)); + return buffer.data(); +} diff --git a/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h b/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h index 45fbb7bc6..757eccb04 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h +++ b/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h @@ -25,6 +25,14 @@ class ImagePreprocessorBufferGPU : public ImagePreprocessorBuffer { // image is final by the time a gather is asked for. The NULL stream would serialise all workers. CudaStream gather_stream; + // GatherTiles: the squares' centres, the squares on the card, and the page-locked landing for them. + // Grown to the largest frame's spots. + size_t tile_capacity = 0; // pixels of gpu_tiles and host_tiles + CudaDevicePtr gpu_tile_centre; + CudaDevicePtr gpu_tiles; + PinnedVector tile_centre; + PinnedVector host_tiles; + public: // host_mirror = false skips the host copy of the preprocessed image entirely (and with it the // page-locking): pass it when every engine reading this buffer runs on the device. @@ -33,4 +41,7 @@ public: int32_t *getGPUBuffer() override; const int32_t *getGPUBuffer() const override; void Gather(const std::vector &npixel, std::vector &values) const override; + const int32_t *GatherTiles(const std::vector &x, const std::vector &y, + int width, int half) override; + const int32_t *HostImage() override; }; diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index df0d1fa21..73ac65de4 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -42,6 +42,8 @@ #include "../image_analysis/rotation_indexer/RotationIndexer.h" // RotationIndexerResult #include "../image_analysis/geom_refinement/PowderCalibration.h" // CalibrationMethod, CalibrationResult #include "../image_analysis/WriteReflections.h" // OutputProvenance +#include "../image_analysis/bragg_integration/SpotWidth.h" // FluxCurve +#include "../image_analysis/bragg_integration/SpotFootprint.h" // FootprintSpot class JFJochReader; class JFJochHDF5Reader; @@ -1122,9 +1124,40 @@ class Rugnux { // The pre-scan's defective-pixel mask (HotPixelFinder), read on `sample` with the rings drawn // about `ring_centre` where the background measured one. frame(i, scratch) is sample frame i, as // the pre-scan kept it or read into scratch; null where it cannot be read. + // What the pre-scan measures on its frames besides the beam stop - the spots, the defective pixels + // and the background gain - kept under everything they read (measure_key in PreScan), for the same + // reason as the beam stop below: a copy of the run that pre-scans again on the same inputs would + // measure the same again. Shared with the copies the run makes of itself. + struct PreScanMemo { + static constexpr size_t MAX_ENTRIES = 4; + struct Spots { + std::vector width_curves, shape_curves; + std::vector footprint; + std::vector spot_q; + size_t width_images = 0; + }; + struct Defects { + size_t pixels = 0; + std::vector> mask; // the mask's non-zero pixels + size_t hot = 0, error = 0, dead = 0; + uint32_t frames = 0; + }; + struct Entry { + std::vector key; + size_t frames_read = 0; + std::optional spots; + std::optional defects; + std::optional> gain; + }; + std::mutex m; + std::vector entries; + }; + std::shared_ptr prescan_memo_ = std::make_shared(); + // `defects` is the answer kept from an earlier pre-scan, or empty and then filled in here. void MaskDefectivePixels(int start_image, const std::vector &sample, const std::optional> &ring_centre, - const std::function &frame); + const std::function &frame, + std::optional &defects); // The beam-stop detection of the pre-scan, kept under everything it reads (the key in PreScan): a // detection that finds no shadow leaves the mask as it was, and every copy of the run that // pre-scans after it would detect again on the same frames, geometry and mask. Shared with the diff --git a/rugnux/RugnuxPreScan.cpp b/rugnux/RugnuxPreScan.cpp index 095e7c1d3..300d2f2c1 100644 --- a/rugnux/RugnuxPreScan.cpp +++ b/rugnux/RugnuxPreScan.cpp @@ -224,7 +224,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // answer kept under them instead (BeamStopMemo). std::vector shadow_key; std::optional shadow_memo; - if (want_shadow) { + { const auto &g = experiment_.GetDiffractionGeometry(); const uint64_t mask_checksum = pixel_mask_.GetBinaryMaskChecksum(); shadow_key = {static_cast(start_image), static_cast(config_.stride), @@ -239,6 +239,8 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru for (const float v : g.GetDetectorMatrix().arr()) shadow_key.push_back(v); shadow_key.insert(shadow_key.end(), shadow_sample.begin(), shadow_sample.end()); + } + if (want_shadow) { std::lock_guard lock(beam_stop_memo_->m); for (const auto &e : beam_stop_memo_->entries) if (e.key == shadow_key) @@ -246,6 +248,42 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru } const bool find_shadow = want_shadow && !shadow_memo; + // The spot measurement, the defective pixels and the background gain are kept the same way, under + // the same key and everything else they read: the settings the spots are found with, the + // azimuthal binning the adaptive finder thresholds on, and whether the shadow is detected first + // (the defective pixels are compared with it masked). A copy of the run that pre-scans again on + // the same inputs takes them from here and, where nothing else wants the frames, reads none. + std::vector measure_key = shadow_key; + { + const auto &f = config_.spot_finding; + const auto &a = experiment_.GetAzimuthalIntegrationSettings(); + measure_key.insert(measure_key.end(), + {static_cast(want_shadow), static_cast(experiment_.GetSaturationLimit()), + goniometer ? goniometer->GetIncrement_deg() : 0.0, + a.GetLowQ_recipA(), a.GetHighQ_recipA(), a.GetQSpacing_recipA(), + static_cast(a.GetAzimuthalBinCount()), static_cast(a.IsSolidAngleCorrection()), + static_cast(a.IsPolarizationCorrection()), + static_cast(f.enable), f.signal_to_noise_threshold, + static_cast(f.photon_count_threshold), + static_cast(f.min_pix_per_spot.value_or(-1)), static_cast(f.max_pix_per_spot), + f.high_resolution_limit.value_or(-1.0f), f.low_resolution_limit.value_or(-1.0f), + f.cutoff_spot_count_low_res, f.high_res_gap_Q_recipA.value_or(-1.0f), + f.ice_ring_width_Q_recipA, static_cast(f.adaptive_threshold), + f.false_pixels_per_frame, static_cast(experiment_.IsDetectIceRings()), + static_cast(f.measured_ring_q_recipA.size())}); + measure_key.insert(measure_key.end(), f.measured_ring_q_recipA.begin(), f.measured_ring_q_recipA.end()); + } + std::optional measured; + { + std::lock_guard lock(prescan_memo_->m); + for (const auto &e : prescan_memo_->entries) + if (e.key == measure_key) + measured = e; + } + const bool measure_width = want_width && !(measured && measured->spots); + const bool measure_hot = want_hot && !(measured && measured->defects); + const bool measure_gain = want_gain && !(measured && measured->gain); + std::optional finder; if (find_shadow) finder.emplace(experiment_, pixel_mask_); @@ -253,7 +291,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // The spot width is read in the same pass as the projection rather than in a pre-pass of its own. // The image is preprocessed a second time for it: the shadow is accumulated from raw counts and // the spot finder works on the converted image. - const bool want_spots = want_width; + const bool want_spots = measure_width; // The spots are measured in the background while the shadow, the defective pixels and the beam // centre are found (see below), and those write the mask. So the spot measurement reads its own // copies of the experiment and of the mask as they are now - which is what it always read, since @@ -279,14 +317,14 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // On the card where there is one, as in the image loops: the frame is decoded and preprocessed // there and the spots are found and extracted there. The device finders give the host's spot // list to the bit - integer ring sums, the same connected components in the same order - // (AdaptiveSpotFinderGPU, SpotExtractorGPU) - so this is a choice of where, not of what. The - // preprocessed image comes back only for the spot width, which reads pixels around the spots. + // (AdaptiveSpotFinderGPU, SpotExtractorGPU) - so this is a choice of where, not of what. Of the + // preprocessed image only the squares about the spots come back (SPOT_TILE_HALF). if (want_spots && get_gpu_count() > 0) { auto stream = std::make_shared(); w.preprocessor = std::make_unique(prescan_x, prescan_mask, stream, - /*copy_image_to_host=*/want_width); + /*copy_image_to_host=*/false); w.preprocessed = std::make_unique(prescan_x.GetPixelsNum(), - /*host_mirror=*/want_width); + /*host_mirror=*/false); if (config_.spot_finding.adaptive_threshold) w.spot_finder = std::make_unique(*prescan_mapping, stream); else @@ -345,23 +383,12 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru } } // The width goes first: it wants the whole spot list, to know which spots are isolated. - if (for_width) { - MeasureSpotFluxCurves(*w.preprocessed, static_cast(prescan_x.GetXPixelsNumConv()), - static_cast(prescan_x.GetYPixelsNumConv()), - prescan_x.GetDiffractionGeometry(), spots, curves); - // The footprint is measured on every spot, not only the isolated ones the width takes: - // it is the far, wide spots it is for (SpotFootprint.h). - std::vector fx, fy; - for (const auto &spot : spots) { - const Coord c = spot.RawCoord(); - fx.push_back(c.x); - fy.push_back(c.y); - } - const auto &g = prescan_x.GetDiffractionGeometry(); - MeasureFootprintSpots(w.preprocessed->data(), static_cast(prescan_x.GetXPixelsNumConv()), - static_cast(prescan_x.GetYPixelsNumConv()), - g.GetBeamX_pxl(), g.GetBeamY_pxl(), fx, fy, footprint); - } + // The footprint is measured on every spot, not only the isolated ones the width takes: it is + // the far, wide spots it is for (SpotFootprint.h). + if (for_width) + MeasureSpotsOnTiles(*w.preprocessed, static_cast(prescan_x.GetXPixelsNumConv()), + static_cast(prescan_x.GetYPixelsNumConv()), prescan_x.GetDiffractionGeometry(), + spots, curves, footprint); }; std::vector prescan_spot_q; @@ -394,12 +421,12 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru const size_t nworkers = std::min(std::max(config_.nthreads, 1), std::min(PRESCAN_MAX_WORKERS, ordinals.size())); // Per frame, in sample order, so the gain does not depend on how the workers interleaved. - std::vector gain_levels(want_gain ? ordinals.size() : 0); + std::vector gain_levels(measure_gain ? ordinals.size() : 0); const size_t gain_frame_stride = (ordinals.size() + background_gain::GAIN_MAX_FRAMES - 1) / background_gain::GAIN_MAX_FRAMES; // The sample as read here, by position, for the spot measurement and the defective pixels below // (PRESCAN_KEPT_BYTES). - const bool keep_frames = want_spots || want_hot; + const bool keep_frames = want_spots || measure_hot; std::vector kept(keep_frames ? ordinals.size() : 0); std::vector is_kept(ordinals.size(), 0); std::atomic kept_bytes{0}; @@ -418,7 +445,21 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru return nullptr; }; - if (find_shadow || want_gain || keep_frames) { + // The spot measurement's engines - the azimuthal mapping its adaptive finder thresholds on, and a + // preprocessor and a spot finder per worker - read only the experiment and the mask as they are + // now, so they are built while the frames are read rather than after. + std::vector workers; + std::future engines_built; + if (want_spots) { + workers.resize(nworkers); + engines_built = std::async(std::launch::async, [&] { + prescan_mapping = std::make_unique(prescan_x, prescan_mask, + *azint_geometry_); + OnWorkers(nworkers, [&](size_t t) { workers[t] = make_worker(); }); + }); + } + + if (find_shadow || measure_gain || keep_frames) { TimingMark mark(logger, "pre-scan read and projection"); std::atomic next{0}; OnWorkers(nworkers, [&](size_t) { @@ -458,7 +499,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru msg.original_number = image_idx; finder->AddImage(msg, shadow_buffer); } - if (want_gain && i % gain_frame_stride == 0) { + if (measure_gain && i % gain_frame_stride == 0) { const auto &img = raw_image.image; const size_t n = img.GetWidth() * img.GetHeight(); const uint8_t *ptr = img.GetUncompressedPtr(shadow_buffer); @@ -491,7 +532,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru logger.Debug("Pre-scan: {} of {} frames kept for the steps after the projection, {:.0f} MB", std::count(is_kept.begin(), is_kept.end(), 1), ordinals.size(), kept_bytes.load() / 1e6); } else - frames_read = shadow_memo->frames_read; + frames_read = measured ? measured->frames_read : shadow_memo->frames_read; // Applied before anything below copies the integration settings (bragg_before_adaptive_), so every // pass and every fallback integrates with it. At or below 1 - an ADSC or a marCCD, whose pixel @@ -499,7 +540,12 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // exactly as it is. if (want_gain) { background_gain_measured_ = true; - const auto g = background_gain::Gain(gain_levels); + if (measure_gain) { + if (!measured) + measured.emplace(); + measured->gain = background_gain::Gain(gain_levels); + } + const std::optional g = *measured->gain; if (g) { BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings(); bis.BackgroundPixelGain(std::max(1.0f, *g)); @@ -518,8 +564,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru const auto spot_measurement = [&] { TimingMark mark(logger, "pre-scan spot measurement"); - prescan_mapping = std::make_unique(prescan_x, prescan_mask, - *azint_geometry_); + engines_built.get(); std::vector> curves_of(ordinals.size()); std::vector> spot_q_of(ordinals.size()); std::vector> footprint_of(ordinals.size()); @@ -527,16 +572,13 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // out of the powder measurement, which is taken on the frames it always was. std::vector shape_only(ordinals.size(), 0); - // The order the frames are visited in, and where each width tier ends in it. Without a width - // to measure there is one tier and the sample is read in its own order. + // The order the frames are visited in, and where each width tier ends in it. std::vector visit; std::vector tier_end; visit.reserve(ordinals.size()); { std::vector taken(ordinals.size(), 0); - const size_t n_tier = want_width ? WIDTH_TIER_STRIDE.size() : 1; - for (size_t tier = 0; tier < n_tier; tier++) { - const size_t stride = want_width ? WIDTH_TIER_STRIDE[tier] : 1; + for (const int stride : WIDTH_TIER_STRIDE) { for (size_t j = 0; j < ordinals.size(); j += stride) if (!taken[j]) { taken[j] = 1; @@ -545,68 +587,97 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru tier_end.push_back(visit.size()); } } + // Where each frame is in that order, and the tier it is in. + std::vector position(ordinals.size()); + for (size_t v = 0; v < visit.size(); v++) + position[visit[v]] = v; + const auto tier_of = [&](size_t v) { + return static_cast(std::upper_bound(tier_end.begin(), tier_end.end(), v) - tier_end.begin()); + }; - // Kept across the tiers rather than rebuilt for each. A worker's preprocessed image is a - // buffer the size of the detector, so each is still built on its own worker the first time - // that worker runs and its pages are touched by the thread that reads them. - std::vector workers(nworkers); - std::vector worker_built(nworkers, 0); - - // Each tier runs on fresh threads (OnWorkers), worker t on the card its engine was built on. - - bool width_settled = !want_width; - // The bandwidth reads the same spots, but its slope needs more of them than the width does: - // the tiers carry on past the width's for it alone, until the pool is large enough. - bool shapes_settled = !want_width; - std::optional previous_r80; - size_t phase_begin = 0; - for (const size_t phase_end : tier_end) { - std::atomic next{phase_begin}; + // The workers go through the frames in that order without waiting at the end of a tier: a frame + // is measured the same whichever tier it is read for, so the next tier's frames are started + // while the last of this one's finish and the answer below is drawn. What that answer reads is + // the frames of the tiers so far and nothing else, and the frames measured past the tier it + // stops at are thrown away - so how far ahead the workers got changes nothing. + std::mutex m; + std::condition_variable cv; + std::vector done_in_tier(tier_end.size(), 0); + bool failed = false; + std::atomic next{0}; + std::atomic stop{false}; + auto measuring = std::async(std::launch::async, [&] { OnWorkers(nworkers, [&](size_t t) { - if (!worker_built[t]) { - workers[t] = make_worker(); - worker_built[t] = 1; - } PreScanWorker &w = workers[t]; - for (size_t v = next.fetch_add(1); v < phase_end; v = next.fetch_add(1)) { - const size_t i = visit[v]; - const int ordinal = ordinals[i]; - const int image_idx = start_image + ordinal * config_.stride; - const bool for_width = !width_settled; - const bool for_shape = width_settled && !shapes_settled; - if (!(for_width || for_shape)) - continue; // nothing to measure on this frame - const CompressedImage *image = frame(i, w.raw_image); - if (!image) continue; - - if (for_shape) shape_only[i] = 1; - find_spots(w, *image, image_idx, for_width || for_shape, curves_of[i], - spot_q_of[i], footprint_of[i]); + try { + for (size_t v = next.fetch_add(1); v < visit.size() && !stop.load(); v = next.fetch_add(1)) { + const size_t i = visit[v]; + const int image_idx = start_image + ordinals[i] * config_.stride; + if (const CompressedImage *image = frame(i, w.raw_image)) + find_spots(w, *image, image_idx, true, curves_of[i], spot_q_of[i], footprint_of[i]); + { + std::lock_guard lock(m); + done_in_tier[tier_of(v)]++; + } + cv.notify_all(); + } + } catch (...) { + { + std::lock_guard lock(m); + failed = true; + } + cv.notify_all(); + throw; } }); - phase_begin = phase_end; + }); + + bool width_settled = false; + // The bandwidth reads the same spots, but its slope needs more of them than the width does: + // the tiers carry on past the width's for it alone, until the pool is large enough. + bool shapes_settled = false; + std::optional previous_r80; + size_t measured_end = 0; // the frames, in visit order, the answer is taken on + for (size_t tier = 0; tier < tier_end.size() && !shapes_settled; tier++) { + const size_t tier_begin = tier ? tier_end[tier - 1] : 0; + { + std::unique_lock lock(m); + cv.wait(lock, [&] { return failed || done_in_tier[tier] == tier_end[tier] - tier_begin; }); + if (failed) break; + } + measured_end = tier_end[tier]; + // A tier begun once the width had settled is read for the bandwidth alone. + if (width_settled) + for (size_t v = tier_begin; v < measured_end; v++) + shape_only[visit[v]] = 1; // What this tier makes of the width, pooled in sample order so the answer does not // depend on how the workers interleaved. A tier too sparse to measure settles nothing: // the test is against what the smaller sample actually said, not against the radius it // fell back to. - if (shapes_settled) continue; + std::vector pooled; + for (size_t i = 0; i < ordinals.size(); i++) + if (position[i] < measured_end) + pooled.insert(pooled.end(), curves_of[i].begin(), curves_of[i].end()); if (!width_settled) { - width_images = phase_end; - width_curves.clear(); - for (const auto &c : curves_of) - width_curves.insert(width_curves.end(), c.begin(), c.end()); + width_images = measured_end; + width_curves = pooled; const auto r80 = spot_width::R80AtReference(width_curves); width_settled = r80 && previous_r80 && spot_width::WidthSettled(*r80, *previous_r80); previous_r80 = r80; } - shape_curves.clear(); - for (const auto &c : curves_of) - shape_curves.insert(shape_curves.end(), c.begin(), c.end()); + shape_curves = std::move(pooled); shapes_settled = width_settled && shape_curves.size() >= spot_width::BANDWIDTH_POOL_SPOTS; } + stop = true; + measuring.get(); OnWorkers(nworkers, [&](size_t t) { PreScanWorker released = std::move(workers[t]); }); + for (size_t i = 0; i < ordinals.size(); i++) + if (position[i] >= measured_end) { + spot_q_of[i].clear(); + footprint_of[i].clear(); + } for (size_t i = 0; i < ordinals.size(); i++) if (!shape_only[i]) prescan_spot_q.insert(prescan_spot_q.end(), spot_q_of[i].begin(), spot_q_of[i].end()); @@ -623,6 +694,21 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru const auto finish_spot_measurement = [&] { if (spots_measured.valid()) spots_measured.get(); + if (want_width && measure_width) { + if (!measured) + measured.emplace(); + measured->spots = PreScanMemo::Spots{width_curves, shape_curves, footprint_spots, prescan_spot_q, + width_images}; + } else if (want_width) { + logger.Info("Spot width: the spots kept from an earlier pre-scan of this run on the same frames, " + "geometry, mask and spot finding"); + const auto &k = *measured->spots; + width_curves = k.width_curves; + shape_curves = k.shape_curves; + footprint_spots = k.footprint; + prescan_spot_q = k.spot_q; + width_images = k.width_images; + } // The integration radius from the recorded spot width. r1 is the aperture the integrator learns // the profile width over and a second moment over a disk of radius a saturates at a^2/4, so at the @@ -885,8 +971,11 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru // projecting every frame for this alone is the largest per-frame cost of the pre-scan. On the // second pass of a two-pass run the shadow is already in the mask, so this does not run again - // by then the centre has been post-refined and the file's value is no longer what is in it. - if (want_hot) - MaskDefectivePixels(start_image, shadow_sample, ring_centre, frame); + if (want_hot) { + if (!measured) + measured.emplace(); + MaskDefectivePixels(start_image, shadow_sample, ring_centre, frame, measured->defects); + } // The centre is fitted with the shadow masked but NOT the transmitting pieces the shadow finder // adds after it. Their dimming is multiplicative in azimuth, which the fit's per-sector amplitude @@ -986,15 +1075,45 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru }).share(); finish_spot_measurement(); + + if (measured) { + measured->key = measure_key; + measured->frames_read = frames_read.load(); + std::lock_guard lock(prescan_memo_->m); + auto &entries = prescan_memo_->entries; + std::erase_if(entries, [&](const auto &e) { return e.key == measure_key; }); + if (entries.size() == PreScanMemo::MAX_ENTRIES) + entries.erase(entries.begin()); + entries.push_back(std::move(*measured)); + } } void Rugnux::MaskDefectivePixels(int start_image, const std::vector &sample, const std::optional> &ring_centre, - const std::function &frame) { + const std::function &frame, + std::optional &defects) { Logger logger("Rugnux"); TimingMark mark(logger, "pre-scan defective pixels"); hot_pixels_measured_ = true; + const auto apply = [&] { + if (defects->frames == 0) + return; + std::vector mask(defects->pixels, 0); + for (const auto &[i, v] : defects->mask) + mask[i] = v; + pixel_mask_.LoadHotPixelMask(experiment_, mask); + const size_t masked = defects->hot + defects->error + defects->dead; + logger.Info("Defective pixels: {} masked from {} pre-scan frames - {} lit above their ring on more " + "frames than one reflection or chance explains, {} holding the detector's error value, {} dead", + masked, defects->frames, defects->hot, defects->error, defects->dead); + Note(fmt::format("{} defective pixel{} masked", masked, masked == 1 ? "" : "s")); + }; + if (defects) { + apply(); + return; + } + // The rings each pixel is compared against have to be drawn about the beam, and the centre the // file states can be far enough out that a ring crosses the background's own radial fall-off - the // pixels on its bright side then read high on every frame. The centre the beam-stop projection @@ -1076,12 +1195,9 @@ void Rugnux::MaskDefectivePixels(int start_image, const std::vector &sample ? oscillation * (sample.back() - sample.front()) / static_cast(sample.size() - 1) : 0.0; const auto hot = finder.GetMask(oscillation, spacing, nthreads); - if (hot.frames == 0) - return; - pixel_mask_.LoadHotPixelMask(experiment_, hot.mask); - const size_t masked = hot.hot + hot.error + hot.dead; - logger.Info("Defective pixels: {} masked from {} pre-scan frames - {} lit above their ring on more " - "frames than one reflection or chance explains, {} holding the detector's error value, {} dead", - masked, hot.frames, hot.hot, hot.error, hot.dead); - Note(fmt::format("{} defective pixel{} masked", masked, masked == 1 ? "" : "s")); + defects = PreScanMemo::Defects{hot.mask.size(), {}, hot.hot, hot.error, hot.dead, hot.frames}; + for (size_t i = 0; i < hot.mask.size(); i++) + if (hot.mask[i] != 0) + defects->mask.emplace_back(static_cast(i), hot.mask[i]); + apply(); } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 864c19593..1ff1356f4 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -82,6 +82,7 @@ ADD_EXECUTABLE(jfjoch_test BackgroundGainTest.cpp BraggStencilTest.cpp SpotFootprintTest.cpp + SpotTilesTest.cpp LossyFilterTest.cpp ImageBufferTest.cpp PixelMaskTest.cpp diff --git a/tests/SpotTilesTest.cpp b/tests/SpotTilesTest.cpp new file mode 100644 index 000000000..7da0047a0 --- /dev/null +++ b/tests/SpotTilesTest.cpp @@ -0,0 +1,145 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include + +#include "../common/CUDAWrapper.h" +#include "../common/DiffractionExperiment.h" +#include "../common/PixelMask.h" +#include "../writer/HDF5Objects.h" +#include "../image_analysis/bragg_integration/SpotWidth.h" +#include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h" +#include "../image_analysis/spot_finding/ImageSpotFinderCPU.h" +#ifdef JFJOCH_USE_CUDA +#include "../image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h" +#endif + +// The pre-scan measures the spot width and the footprint on squares cut out about each spot instead +// of the whole frame (MeasureSpotsOnTiles). That is only allowed because it gives the whole frame's +// answer to the bit, so this checks exactly that on real diffraction frames - including the spots +// whose footprint window walks off its square and is measured again, on a larger square or on the +// frame, and the frame whose squares would be more than the frame itself. +namespace { + +struct Frame { + DiffractionExperiment x; + ImagePreprocessorBuffer image; + std::vector spots; +}; + +std::vector ReadFrames() { + HDF5ReadOnlyFile data("../../tests/test_data/compression_benchmark.h5"); + HDF5DataSet dataset(data, "/entry/data/data"); + HDF5DataSpace file_space(dataset); + const auto dim = file_space.GetDimensions(); + + DiffractionExperiment x(DetDECTRIS(dim[2], dim[1], "Test", {})); + x.PixelSigned(true).DetectorDistance_mm(75).BeamY_pxl(1136).BeamX_pxl(1090).IncidentEnergy_keV(12.4); + const PixelMask mask(x); + ImagePreprocessorCPU preprocessor(x, mask); + ImageSpotFinderCPU finder(static_cast(x.GetXPixelsNumConv()), static_cast(x.GetYPixelsNumConv())); + + std::vector frames; + std::vector raw(dim[1] * dim[2]); + for (hsize_t i = 0; i < dim[0]; i += 4) { + dataset.ReadVector(raw, {i, 0, 0}, {1, dim[1], dim[2]}); + Frame f{x, ImagePreprocessorBuffer(x.GetPixelsNum()), {}}; + preprocessor.Analyze(f.image, reinterpret_cast(raw.data()), CompressedImageMode::Int16); + f.spots = finder.Run(f.image, SpotFindingSettings{}); + frames.push_back(std::move(f)); + } + return frames; +} + +bool Same(const std::vector &a, const std::vector &b) { + if (a.size() != b.size()) return false; + for (size_t k = 0; k < a.size(); k++) { + const auto &p = a[k], &q = b[k]; + if (p.d_A != q.d_A || p.c != q.c || p.shape.valid != q.shape.valid || p.shape.two_theta != q.shape.two_theta + || p.shape.jr != q.shape.jr || p.shape.jt != q.shape.jt || p.shape.m_rad != q.shape.m_rad + || p.shape.m_tan != q.shape.m_tan || p.shape.cos_psi != q.shape.cos_psi + || p.shape.par_u != q.shape.par_u || p.shape.par_v != q.shape.par_v) + return false; + } + return true; +} + +bool Same(const std::vector &a, const std::vector &b) { + if (a.size() != b.size()) return false; + for (size_t k = 0; k < a.size(); k++) + if (a[k].r_px != b[k].r_px || a[k].sigma_rad != b[k].sigma_rad || a[k].sigma_tan != b[k].sigma_tan) + return false; + return true; +} + +// Every frame through MeasureSpotsOnTiles on `image(f)`, against the whole frame; the misses summed. +template +SpotTileMisses CheckFrames(std::vector &frames, Image &&image, int half, int footprint_half) { + SpotTileMisses total; + size_t curves = 0, footprints = 0; + for (auto &f : frames) { + const int width = static_cast(f.x.GetXPixelsNumConv()), height = static_cast(f.x.GetYPixelsNumConv()); + 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); + std::vector sx, sy; + for (const auto &s : f.spots) { + sx.push_back(s.RawCoord().x); + sy.push_back(s.RawCoord().y); + } + 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, + half, footprint_half); + CHECK(Same(c_frame, c_tiles)); + CHECK(Same(f_frame, f_tiles)); + total.larger += misses.larger; + total.whole_frame += misses.whole_frame; + curves += c_frame.size(); + footprints += f_frame.size(); + } + // Something was measured, or the comparison says nothing. + CHECK(curves > 0); + CHECK(footprints > 100); + return total; +} + +} // namespace + +TEST_CASE("SpotTiles_MatchWholeFrame", "[SpotWidth]") { + auto frames = ReadFrames(); + REQUIRE(!frames.empty()); + const auto host = [](Frame &f) -> ImagePreprocessorBuffer & { return f.image; }; + + // As the pre-scan reads them. + CheckFrames(frames, host, SPOT_TILE_HALF, FOOTPRINT_TILE_HALF); + // Squares too small for the footprint: its spots are measured again on the larger squares... + CHECK(CheckFrames(frames, host, spot_width::R_BKG_OUT, spot_width::R_BKG_OUT + 1).larger > 0); + // ... and on the frame, where the larger squares would be more than it. + CHECK(CheckFrames(frames, host, spot_width::R_BKG_OUT, 4000).whole_frame > 0); + // Squares so large that the frame is read instead from the start. + CHECK(CheckFrames(frames, host, 1000, FOOTPRINT_TILE_HALF).whole_frame == 0); +} + +#ifdef JFJOCH_USE_CUDA +// The same with the squares cut out on the card, from a device image the host never holds. +TEST_CASE("SpotTiles_MatchWholeFrame_GPU", "[SpotWidth][GPU]") { + if (get_gpu_count() == 0) + SKIP("No GPU"); + auto frames = ReadFrames(); + REQUIRE(!frames.empty()); + ImagePreprocessorBufferGPU device(frames.front().x.GetPixelsNum(), /*host_mirror=*/false); + const auto upload = [&](Frame &f) -> ImagePreprocessorBuffer & { + REQUIRE(cudaMemcpy(device.getGPUBuffer(), f.image.data(), f.image.size() * sizeof(int32_t), + cudaMemcpyHostToDevice) == cudaSuccess); + return device; + }; + CheckFrames(frames, upload, SPOT_TILE_HALF, FOOTPRINT_TILE_HALF); + CHECK(CheckFrames(frames, upload, spot_width::R_BKG_OUT, spot_width::R_BKG_OUT + 1).larger > 0); + CHECK(CheckFrames(frames, upload, spot_width::R_BKG_OUT, 4000).whole_frame > 0); + CheckFrames(frames, upload, 1000, FOOTPRINT_TILE_HALF); +} +#endif