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